21
$\begingroup$

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.

$\endgroup$
22
$\begingroup$

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):

plot of the sinc function

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.

$\endgroup$
  • $\begingroup$ The sinc function is ideal filtering, though, no? That's what all other filters are aiming for, but not achieving. It's bad for image processing because images aren't antialiased first, so it produces ringing which looks awful, but for audio or other signals which were antialias filtered before sampling, isn't it the best filter you can get? $\endgroup$ – endolith Nov 29 '11 at 15:08
  • 1
    $\begingroup$ Yes, my result was not conjugate symmetric. I corrected the code, now everything works fine. Thank you! $\endgroup$ – Till B Nov 29 '11 at 15:12
  • 3
    $\begingroup$ @endolith - a Sinc is an ideal interpolator for certains kinds of interpolation, but can be far from ideal as a filter for most kinds of common filter requirements, such as flatness of pass band response, stop band rejection, and etc. $\endgroup$ – hotpaw2 Nov 29 '11 at 15:27
  • $\begingroup$ +1 for the nice explanation on "why people don't implement the filter as the PO does" $\endgroup$ – Sibbs Gambling Aug 13 '13 at 6:05
  • $\begingroup$ You have to use a windowed sinc. If you are not time-constrained, this is the optimal filter, far better than Chebichev. $\endgroup$ – user13706 Jul 18 '15 at 21:09

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.