Fft bandpass filter in python

Question

I am trying to filter my data using fft. I have a noisy signal recorded at 500 Hz as a 1d array. My high frequency should be turned off at 20 Hz and my low frequency at 10 Hz. I tried:

fft=scipy.fft(signal) bp=fft[:] for i in range(len(bp)): if not 10<i<20: bp[i]=0 ibp=scipy.ifft(bp) 

Now I get complex numbers. So something must be wrong. What kind? How can I fix my code?

Answer 1

It is worth noting that the units of your bp will not necessarily be in Hz, but depending on the sampling rate of the signal, you should use scipy.fftpack.fftfreq to convert. Also, if your signal is real, you should use scipy.fftpack.rfft . Here is a minimal working example that filters out all frequencies less than a given amount:

 import numpy as np from scipy.fftpack import rfft, irfft, fftfreq time = np.linspace(0,10,2000) signal = np.cos(5*np.pi*time) + np.cos(7*np.pi*time) W = fftfreq(signal.size, d=time[1]-time[0]) f_signal = rfft(signal) # If our original signal time was in seconds, this is now in Hz cut_f_signal = f_signal.copy() cut_f_signal[(W<6)] = 0 cut_signal = irfft(cut_f_signal) 

We can build the evolution of a signal in real and rangefinding space:

 import pylab as plt plt.subplot(221) plt.plot(time,signal) plt.subplot(222) plt.plot(W,f_signal) plt.xlim(0,10) plt.subplot(223) plt.plot(W,cut_f_signal) plt.xlim(0,10) plt.subplot(224) plt.plot(time,cut_signal) plt.show() 

Answer 2

There is a fundamental flaw in what you are trying to do here: you are applying a rectangular window in the frequency domain that will result in a signal in the time domain that has been minimized using the sinc function. In other words, a large number of “rings” will occur in the time-domain signal due to step changes that you entered in the frequency domain. The correct way to filter such a frequency range is to apply a suitable window function in the frequency domain. Any good DSP introductory book should cover this.