I need to filter a signal. I want to keep frequencies between 0 and 51Hz. The following is the code that I am using, part of it is taken from these two questions (Python: Designing a time-series filter after Fourier analysis, Creating lowpass filter in SciPy - understanding methods and units):
def butter_lowpass(cutoff, fs, order=5):
nyq = 0.5 * fs
normal_cutoff = cutoff / nyq
b, a = butter(order, normal_cutoff, btype='low', analog=False)
return b, a
def butter_lowpass_filter(data, cutoff, fs, order=5):
b, a = butter_lowpass(cutoff, fs, order=order)
y = lfilter(b, a, data)
return y
# y is the original signal
# getting the unbiased signal
y = list(np.array(y)-sts.mean(y))
# saving the original signal
y_before = y
# x is the time vector
# for the spectrum of the original signal
yps_before = np.abs(np.fft.fft(y_before))**2
yfreqs_before = np.fft.fftfreq(6000, d = 0.001)
y_idx_before = np.argsort(yfreqs_before)
#Filtering
order = 8
fs = 1000.0
cutoff = 50.0
y = butter_lowpass_filter(y, cutoff, fs, order)
# for the spectrum of the filtered signal
yps = np.abs(np.fft.fft(y))**2
yfreqs = np.fft.fftfreq(6000, d = 0.001)
y_idx = np.argsort(yfreqs)
fig = plt.figure(figsize=(14,10))
fig.suptitle(file_name, fontsize=20)
plt.plot(yfreqs_before[y_idx_before], yps_before[y_idx_before], 'k-', label='original spectrum',linewidth=0.5)
plt.plot(yfreqs[y_idx], yps[y_idx], 'r-', linewidth=2, label='filtered spectrum')
plt.xlabel('Frequency [Hz]')
plt.yscale('log')
plt.grid()
plt.legend()
plt.show()
The result of this code is a filtered signal, however, that is the spectrum comparison:
As you can see the spectrum looks good after 100Hz, however, between 50Hz and about 100Hz there is still a component. Thus I tried to use a filter with higher order (20), but as output, I get a really strange spectrum:
So, I know that a filter can not and will never be perfect, but to me, this seems a little too much. In my experience, I was always able to obtain a pretty good filtered signal at my cutoff frequency. Any advice?