Simple example of signal generation and application of a bandstop butterworth filter in python.

In [11]:

import numpy as np
import scipy
from scipy import signal
import matplotlib.pyplot as plt

Initial definitions:

In [12]:

freq1 = 60  # Hz
freq2 = 90  # Hz
sampling_rate = 1000  # samples per second
time = 1  # seconds
t = np.linspace(0, time, time * sampling_rate)

Added two sine signals of both frequencies:

In [13]:

y = np.sin(2*np.pi*freq1*t) + np.sin(2*np.pi*freq2*t)
# y + random noise
yn = y * 0.3 * np.random.rand(len(t))
dt = 250  # observation

In [14]:

# plot
plt.figure(tight_layout=True)

plt.subplot(211)
plt.plot(t[:dt], y[:dt])
plt.title('y(t)')
plt.grid()
plt.xlabel('t(s)')

plt.subplot(212)
plt.plot(t[:dt], yn[:dt])
plt.title('yn(t)')
plt.grid()
plt.xlabel('t(s)')

Out[14]:

Text(0.5, 0, 't(s)')

Check predominant frequencies in both (original and noise added) signals with Fast Fourier Transform (FFT):

In [15]:

# fft
Y = np.abs(np.fft.fft(y))
Yn = np.abs(np.fft.fft(yn))
f = np.arange(0, len(Y)) * sampling_rate/len(Y)
df = int(200 * len(y) / sampling_rate)

In [16]:

#fft plots
plt.figure(tight_layout=True)

plt.subplot(211)
plt.plot(f[:df], Y[:df])
plt.title('|Y(f)|')
plt.grid()
plt.xlabel('f(Hz)')

plt.subplot(212)
plt.plot(f[:df], Yn[:df])
plt.title('|Yn(f)|')
plt.grid()
plt.xlabel('f(Hz)')

Out[16]:

Text(0.5, 0, 'f(Hz)')

Define filter specifications, create filter and apply to both signals:

In [17]:

# filter bandstop - reject 60 Hztop)
lowcut = 57
highcut = 63
order = 4
nyq = 0.5 * sampling_rate
low = lowcut / nyq
high = highcut / nyq
b, a = signal.butter(order, [low, high], btype='bandstop')

y_filt = signal.filtfilt(b, a, y)
yn_filt = signal.filtfilt(b, a, yn)

In [18]:

# filtered signal
plt.figure(tight_layout=True)

plt.subplot(211)
plt.plot(t[:dt], y_filt[:dt])
plt.title('y_filt(t)')
plt.grid()
plt.xlabel('t(s)')

plt.subplot(212)
plt.plot(t[:dt], yn_filt[:dt])
plt.title('yn_filt(t)')
plt.grid()
plt.xlabel('t(s)')

Out[18]:

Text(0.5, 0, 't(s)')

Check again for predominant frequencies with FFT:

In [19]:

# fft
Y = np.abs(np.fft.fft(y_filt))
Yn = np.abs(np.fft.fft(yn_filt))
f = np.arange(0, len(Y)) * sampling_rate/len(Y)
df = int(200 * len(y_filt) / sampling_rate)

In [20]:

#fft plots
plt.figure(tight_layout=True)

plt.subplot(211)
plt.plot(f[:df], Y[:df])
plt.title('|Y(f) filered|')
plt.grid()
plt.xlabel('f(Hz)')

plt.subplot(212)
plt.plot(f[:df], Yn[:df])
plt.title('|Yn(f)| filered')
plt.grid()
plt.xlabel('f(Hz)')

Out[20]:

Text(0.5, 0, 'f(Hz)')