14

What I try is to filter my data with fft. I have a noisy signal recorded with 500Hz as a 1d- array. My high-frequency should cut off with 20Hz and my low-frequency with 10Hz. What I have tried is: 

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

ibp=scipy.ifft(bp)

What I get now are complex numbers. So something must be wrong. What? How can I correct my code?

Paul R
  • 208,748
  • 37
  • 389
  • 560
men in black
  • 193
  • 1
  • 1
  • 9

2 Answers2

27

It's worth noting that the magnitude of the units of your bp are not necessarily going to be in Hz, but are dependent on the sampling frequency of signal, you should use scipy.fftpack.fftfreq for the conversion. Also if your signal is real you should be using scipy.fftpack.rfft. Here is a minimal working example that filters out all frequencies less than a specified 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 plot the evolution of the signal in real and fourier 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()

enter image description here

Community
  • 1
  • 1
Hooked
  • 84,485
  • 43
  • 192
  • 261
  • Thanks! But if I want a bandpass filter I have to do it twice? Once lowpass and once highpass? – men in black Oct 02 '13 at 09:19
  • 1
    @meninblack Sure. As you can see in the example I've only filtered the low freqs out with a rectangular cutoff (see the answer by Paul R for notes about this, and improvements). What I tried to do is give you a simple, minimal working answer that lets you see what went wrong (wrong units, etc...). – Hooked Oct 02 '13 at 12:56
9

There's a fundamental flaw in what you are trying to do here - you're applying a rectangular window in the frequency domain which will result in a time domain signal which has been convolved with a sinc function. In other words there will be a large amount of "ringing" in the time domain signal due to the step changes you have introduced in the frequency domain. The proper way to do this kind of frequency domain filtering is to apply a suitable window function in the frequency domain. Any good introductory DSP book should cover this.

Paul R
  • 208,748
  • 37
  • 389
  • 560
  • 2
    +1. Linked question: http://dsp.stackexchange.com/questions/6220/why-is-it-a-bad-idea-to-filter-by-zeroing-out-fft-bins – Basj Dec 29 '16 at 00:30