low pass filter and FFT for beginners with Python

I am new to signal processing and especially to FFT, hence I am not sure if I am doing the correct thing here and I am a bit confused with the result.

I have a discrete real function (measurement data) and want to set up a low pass filter on that. The tool of choice is Python with the numpy package. I follow this procedure:

• compute the fft of my function


• cut off high frequencies


• perform the inverse fft

Here is the code that I am using:

import numpy as np
sampling_length = 15.0*60.0 # measured every 15 minutes
Fs = 1.0/sampling_length
ls = range(len(data)) # data contains the function

freq = np.fft.fftfreq(len(data), d = sampling_length)
fft = np.fft.fft(data)
x = freq[:len(data)/2] 
for i in range(len(x)):
if x[i] > 0.005: # cut off all frequencies higher than 0.005
    fft[i] = 0.0
    fft[len(data)/2 + i] = 0.0
inverse = np.fft.ifft(fft)

Is this the correct procedure? The result inverse contains complex values, that confuses me.

Answer

The fact that the result is complex is to be expected. I want to point out a couple things:

You are applying a brick-wall frequency-domain filter to the data, attempting to zero out all FFT outputs that correspond to a frequency greater than 0.005 Hz, then inverse-transforming to get a time-domain signal again. In order for the result to be real, then the input to the inverse FFT must be conjugate symmetric. This means that for a length-$N$ FFT,

$$ X[k] = X^*[N-k], k = 1, 2, \ldots , \frac{N}{2} - 1\;\;\;\;\;\;\;(N\;\; even) $$

$$ X[k] = X^*[N-k], k = 1, 2, \ldots , \lfloor \frac{N}{2} \rfloor\;\;\;\;\;\;\;\;(N\;\; odd) $$

• Note that for $N$ even, $X[0]$ and $X[\frac{N}{2}]$ are not equal in general, but they are both real. For odd $N$, $X[0]$ must be real.

I see that you attempted to do something like this in your code above, but it is not quite correct. If you enforce the above condition on the signal that you pass to the inverse FFT, then you should get a real signal out.

My second point is more of a philosophical one: what you're doing will work, in that it will suppress the frequency-domain content that you don't want. However, this is not typically the way a lowpass filter would be implemented in practice. As I mentioned before, what you're doing is essentially applying a filter that has a brick-wall (i.e. perfectly rectangular) magnitude response. The impulse response of such a filter has a $sinc(x)$ shape. Since multiplication in the frequency domain is equivalent to (in the case of using the DFT, circular) convolution in the time domain, this operation is equivalent to convolving the time domain signal with a $sinc$ function.

Why is this a problem? Recall what the $sinc$ function looks like in the time domain (below image shamelessly borrowed from Wikipedia):

The $sinc$ function has very broad support in the time domain; it decays very slowly as you move in time away from its main lobe. For many applications, this is not a desirable property; when you convolve a signal with a $sinc$, the effects of the slowly-decaying sidelobes will often be apparent in the time-domain form of the filtered output signal. This sort of effect is often referred to as ringing. If you know what you're doing, there are some instances where this type of filtering might be appropriate, but in the general case, it's not what you want.

There are more practical means of applying lowpass filters, both in the time and frequency domains. Finite impulse response and infinite impulse response filters can be applied directly using their difference equation representation. Or, if your filter has a sufficiently-long impulse response, you can often obtain performance benefits using fast convolution techniques based on the FFT (applying the filter by multiplying in the frequency domain instead of convolution in the time domain), like the overlap-save and overlap-add methods.