In the 7th example we already discussed how to low/band/high pass filter using FFT, but there’s a problem with the approach described there. As you can heard, the low pass filtered sound has a distorted ringing effect. This is okay for some applications (it doesn’t count when we are decoding an RTTY signal…), but in most cases you want to avoid it.

Why is the ringing?

(source) We are trying to suppress the frequency domain content that we don’t want. However, zeroing out bands in the frequency domain is not typically the way a filter would be implemented in practice. This way what we are doing is essentially applying a filter that has a brick-wall (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? The sinc function looks like this in the time domain:

(Image source: 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.

Convolution

To avoid this artifact but keep the filter FIR, you have to use convolution, which can be done more efficiently with FFT. The convolution theorem says: multiplication in the frequency domain corresponds to convolution in the time domain. So to apply your filter, you have to calculate the FFTs of the audio data and your filter, and you calculate their pointwise multiplication and then inverse FFT them.

Here’s the whole filtering procedure:

1. Design the filter using the sinc(x) function.
2. Apply a window function.
3. Calculate the FFT of the filter and store it.
4. Read a bunch of samples from the sound device to a buffer.
5. Calculate the FFT of the buffer.
6. Pointwise multiplicate the audio FFT buffer with the filter FFT buffer.
7. Calculate the inverse FFT of the resulting buffer.
8. Use the overlap-add method to get the filtered audio.
9. Go to step 3.

These things and the theory are described in this excellent article in detail. Check that page because there are lots of nice diagrams.

There’s a good article which implements the convolution and filter calculation in Matlab here.

Window function

(source) Whenever you do a finite Fourier transform, you’re implicitly applying it to an infinitely repeating signal. So, for instance, if the start and end of your finite sample don’t match then that will look just like a discontinuity in the signal, and show up as lots of high-frequency nonsense in the Fourier transform, which you don’t really want. And if your sample happens to be a beautiful sinusoid but an integer number of periods don’t happen to fit exactly into the finite sample, your FFT will show appreciable energy in all sorts of places nowhere near the real frequency. You don’t want any of that.

Windowing the data makes sure that the ends match up while keeping everything reasonably smooth; this greatly reduces the sort of “spectral leakage” described in the previous paragraph.

Overlap-add method

We are calculating convolution on real-time audio chunks at a time, not on a whole (pre-saved) .wav file. As the result of the convolution will be the sum of the length of the audio and filter data, we have to use some kind of overlapping.

0. Initialize an overlap buffer with zeros, it’s length should be equal of the audio and filter buffer lengths.
1. After you got the convolution after the inverse FFT, add each of the overlap buffer elements to the firsthalf of the resulting buffer’s elements one by one.
2. Now that buffer has the resulting filtered audio.
3. Copy the resulting inverse FFT buffer’s second half to the overlap buffer.
4. Go to step 1.

The overlap-add method is also described in the dspguide.com article in detail.

Now let’s see my implementation of FIR filtering with FFT convolution in Java:

Source code