```
def bandpass_ifft(X, Low_cutoff, High_cutoff, F_sample, M=None):
"""Bandpass filtering on a real signal using inverse FFT
Inputs
=======
X: 1-D numpy array of floats, the real time domain signal (time series) to be filtered
Low_cutoff: float, frequency components below this frequency will not pass the filter (physical frequency in unit of Hz)
High_cutoff: float, frequency components above this frequency will not pass the filter (physical frequency in unit of Hz)
F_sample: float, the sampling frequency of the signal (physical frequency in unit of Hz)
Notes
=====
1. The input signal must be real, not imaginary nor complex
2. The Filtered_signal will have only half of original amplitude. Use abs() to restore.
3. In Numpy/Scipy, the frequencies goes from 0 to F_sample/2 and then from negative F_sample to 0.
"""
import scipy, numpy
if M == None: # if the number of points for FFT is not specified
M = X.size # let M be the length of the time series
Spectrum = scipy.fft(X, n=M)
[Low_cutoff, High_cutoff, F_sample] = map(float, [Low_cutoff, High_cutoff, F_sample])
#Convert cutoff frequencies into points on spectrum
[Low_point, High_point] = map(lambda F: F/F_sample * M /2, [Low_cutoff, High_cutoff])# the division by 2 is because the spectrum is symmetric
Filtered_spectrum = [Spectrum[i] if i >= Low_point and i <= High_point else 0.0 for i in xrange(M)] # Filtering
Filtered_signal = scipy.ifft(Filtered_spectrum, n=M) # Construct filtered signal
return Spectrum, Filtered_spectrum, Filtered_signal, Low_point, High_point
```

Here is an example showing its usage.

First of all, let's generate a signal of different frequencies.

```
import numpy
N = 400 # signal length of number
x = numpy.arange(0, N, 1) # generate the time ticks
Sines = [numpy.sin(x*n)*(1-n) for n in [.9, .75, .5, .25, .12, .03, 0.025]]# different frequency components
y = numpy.sum(Sines, axis=0) # add them by column, low frequencies have higher amplitudes
import matplotlib.pyplot as plt
plt.plot(x, y) # visualize the data
```

The raw signal looks like this:

Then, let's assume the sampling frequency is 500Hz and we only want signals between 5 to 30 Hz to go thru the filter.

```
Low_cutoff, High_cutoff, F_sample = 5, 30, 500
Spectrum, Filtered_spectrum, Filtered_signal, Low_point, High_point = bandpass_ifft(y, Low_cutoff, High_cutoff, F_sample)
# Below is visualization
fig1 = plt.figure()
plt.stem(numpy.fft.fftfreq(N)*F_sample, abs(Spectrum))
plt.axis([0, F_sample/2, None, None]) # since the signal is real, the spectrum is symmetric|
plt.axvspan(Low_cutoff, High_cutoff,fc='g',alpha='.5') # The frequencies that we wanna get rid of
fig2 = plt.figure()
plt.plot(x, Filtered_signal)
```

In the first figure plotted, the green box shows frequencies that will be kept by the filter.

Here is the filtered signal. As you can see, high frequencies are gone. The signal looks very smooth.

Of course, this approach is very costly: FFT and iFFT are both $O(n\lg n)$ complexity algorithms. It is much more costly than convolving signal with filter terms (what people normally do in digital signal processing). Oh, and it cannot be done in real-time - you must wait until all data is sampled.

There is also a minor implementation issue here. Both Scipy and Numpy provides FFT/iFFT functions. But their speeds vary [1]. If you wanna go faster, try Python binding to some C/C++-based FFT libraries.

References:

1. http://stackoverflow.com/questions/6363154/what-is-the-difference-between-numpy-fft-and-scipy-fftpack

## 1 comment:

nice work, but it seems that this code set the all the power in [negative F_sample, 0] to zero, including the corresponding power btw [Low_cutoff, High_cutoff]. Doesn't this make some power loss in ifft?

Post a Comment