## DSP tutorial: Lowpass FIR filtering using FFT convolution

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.

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.

1. Initialize an overlap buffer with zeros, it’s length should be equal of the audio and filter buffer lengths.
2. 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.
3. Now that buffer has the resulting filtered audio.
4. Copy the resulting inverse FFT buffer’s second half to the overlap buffer.
5. 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

 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117 // see http://en.wikipedia.org/wiki/Sinc_function     private double sinc(double x) {         if (x == 0)             return 1;         return Math.sin(Math.PI*x)/(Math.PI*x);     }         LowPassFilter(TargetDataLine tdl) {         audio = new double[FFTRES/2];         filter = new double[FFTRES/2+1];         audioFFT = new double[FFTRES*2]; // *2 because FFT data has both real and imaginary parts         filterFFT = new double[FFTRES*2];         audioOverlap = new double[FFTRES/2];         audioFiltered = new double[FFTRES/2];                 this.tdl = tdl;         audioBytes = new byte[audio.length*2];         baos = new ByteArrayOutputStream();         fft = new DoubleFFT_1D(FFTRES);         System.out.println("FFT resolution: " + FFTRES);         System.out.println("Audio sample in buffer length: " + audio.length);         System.out.println("Filter length: " + filter.length);         // designing the windowed sinc filter         int cutFreq = 1000;         for (int i = 0; i < filter.length-1; i++) {             // see http://en.wikipedia.org/wiki/Sinc_filter             double sincFilter = (2*cutFreq/(double)SAMPLERATE)*sinc(2*cutFreq*(i-((filter.length-1)/2.0))/(double)SAMPLERATE);             // uncomment for a bandpass filter             //sincFilter = sincFilter-(2*(cutFreq+1000)/(double)SAMPLERATE)*sinc(2*(cutFreq+1000)*(i-((filter.length-1)/2.0))/(double)SAMPLERATE);             // uncomment for a highpass filter             /*if (i != (filter.length-1)/2)                 sincFilter *= -1;             else                 sincFilter = 1-sincFilter;*/             // applying a Hamming window, see http://stackoverflow.com/questions/5418951/what-is-the-hamming-window-for and http://en.wikipedia.org/wiki/Window_function             //filter[i] = (0.53836 - (0.46164 * Math.cos((Math.PI*2) * (double)i / (double)(filter.length-1)))) * idealFilter;             // applying a Hann window             //filter[i] = 0.5*(1-Math.cos( (2*Math.PI*i)/(double)(filter.length-1) )) * sincFilter;             // applying a Blackman window             filter[i] = (0.42-0.5*Math.cos((2*Math.PI*i)/(double)(filter.length-1))+0.08*Math.cos((4*Math.PI*i)/(double)(filter.length-1))) * sincFilter;         }         // clearing the audio overlap buffer         for (int i = 0; i < audioOverlap.length; i++)             audioOverlap[i] = 0;         // clearing the FFT buffer         for (int i = 0; i < filterFFT.length; i++)             filterFFT[i] = 0;         // copying time domain filter data to the FFT buffer         for (int i = 0; i < filter.length; i++)             filterFFT[2 * i] = filter[i];         fft.complexForward(filterFFT);     }     private void lowPassFilter(double[] audioDataIn, double[] audioDataOut) {         // zeroing out the audio FFT buffer         for (int i = 0; i < audioFFT.length; i++)             audioFFT[i] = 0;         // copying audio data to the FFT data buffer         for (int i = 0; i < audioDataIn.length; i++)             audioFFT[2 * i] = audioDataIn[i];         // calculating the fft of the data         fft.complexForward(audioFFT);         // pointwise multiplication of the filter and audio data in the frequency domain         for (int i = 0; i < filterFFT.length; i += 2) {             double temp = audioFFT[i] * filterFFT[i] - audioFFT[i+1]*filterFFT[i+1];             audioFFT[i+1] = audioFFT[i] * filterFFT[i+1] + audioFFT[i+1] * filterFFT[i]; // imaginary part             audioFFT[i] = temp; // real part         }         // built-in scaling hangs the thread, so we don't use it         fft.complexInverse(audioFFT, false);         // adding the first half of the audio FFT buffer to the overlap buffer         for (int i = 0; i < audioDataOut.length; i++)             audioDataOut[i] = (audioOverlap[i] + audioFFT[i * 2]) / 2000; // applying scaling         // copying the second half of the audio FFT buffer to the audio overlap buffer         for (int i = 0; i < audioOverlap.length; i++)             audioOverlap[i] = audioFFT[audioFFT.length/2 + i * 2];     }     @Override     public void run() {         tdl.start();         try {             while (!Thread.interrupted()) {                 // waiting for the buffer to get filled                 while (tdl.available() < audioBytes.length)                     Thread.sleep(0, 1); // without this, the audio will be choppy                 int bytesRead = tdl.read(audioBytes, 0, audioBytes.length);                 // converting frames stored as bytes to double values                 int samplesRead = bytesRead / 2;                 for (int i = 0; i < samplesRead; i++)                     audio[i] = ((audioBytes[i * 2] & 0xFF) | (audioBytes[i * 2 + 1] << 8)) / 32768.0;                 lowPassFilter(audio, audioFiltered);                 baos.write(getBytesFromDoubles(audioFiltered, audioFiltered.length), 0, audioFiltered.length * 2);             }         } catch (InterruptedException e) {         }         tdl.stop();         tdl.close();         writeWavFile(baos.toByteArray(), baos.size() / 2, "output.wav");     }  Michał 2012-03-23 23:01:26

Hi,

Nice note about filtering signal in frequency domain. I’m only curious about one thing. Why scaling factor in:

audioDataOut[i] = (audioOverlap[i] + audioFFT[i * 2]) / 2000; // applying scaling

is exact 2000?

This is a value chosen empirically or I don’t understand something of this post and mentioned articles? :)

Greetings from Poland,
Michał  2012-03-24 07:30:22

Hello, yes, this scaling value was chosen empirically :)  dude 2012-04-05 14:18:25

I am abit curious, how would this be applied with the overlap add, if what you were processing was realtime? like a function that has input and out as parameters? e.X: void process(double *input, double *output); which processed each sample? I tried myself, but i just can’t seem to get my head around it  2012-04-05 23:41:30

Realtime processing of overlap add could be done on a buffer without problems.  Dude 2012-04-06 20:36:45

Yeah probably, i’ve been sitting with this one for ages.

If i have a fft buffer which is based on a struct with two entries.
Im and Re. It makes sense to creata a buffer that is 1024 large. No need to think about multiplying i * 2 to get the real part, because i can simply just access the re or im with [i].re or [i].im in an iteration.

So anyhow, having two fft buffers, one filter and one with the signal.

You multiply them both, and the perform the inverse FFT on the multiplied result. This multiplied result should also be 1024 right?

In that case i dont see why and where you should introduce this overlap add in a realtime scenario.

Since here, i will have to wait until i have saved 1024 successful samples before i can begin to compute.

I guess this is not your problem what so ever, but i just can’t figure this out. It bothers me, i hate not being able to solve stuff :(  saru 2012-08-14 06:13:19

iwant htis code in matlab could u please send it to me in my add ” saruhada@gmail.com  2012-08-14 07:50:39

Sorry, I don’t have this code in Matlab.  Brini 2012-09-12 22:05:27

I´m quite new in dsp world.

audio[i] = ((audioBytes[i * 2] & 0xFF) | (audioBytes[i * 2 + 1] << 8)) / 32768.0;

why do you divide by 32768?

Thanks for posting the code I was having big time trouble doing a low pass filter.  2012-09-12 23:02:40

That line converts the 16 bit input value (stored as two bytes) to double (between 0 and 1).  Brini 2012-09-14 21:01:30

Thanks.  Joe McDaniel 2015-04-27 17:04:44

Yes, but the max value of a 16-bit is 32767. The use of 32768.0 won’t cause any problems — you are simply decreasing the signal by 1/32768.  Thomas 2012-09-22 23:43:27

Hi

Thanks a lot for this!
I don’t know (don’t understand) how to set the bandwidth for the bandpass filter.

sincFilter-(2*(cutFreq+1000)/(double)SAMPLERATE)*sinc(2*(cutFreq+1000)*(i-((filter.length-1)/2.0))/(double)SAMPLERATE);

First I expected if going from cutFreq until cutFreq+1000, but it doesn’t seem to work like that…  Thomas 2012-09-23 00:03:11

Sorry, my fault: Seems exactly to work like that …  KernowPete 2013-03-01 13:30:13

Hi,

Thank you so much for this Tutorial, A really concise explanation. I have found simpler break downs of this process, for a DSP noob, very hard to find.

I have one question. I am trying to make a version of this to run in processing, due to my college assignment being implement a band pass filter in Processing.

Would you have any suggestions as to how I would need to change this into a processing version?

Thank you hugely as the step by step explanation has finaaly at least explained the full process to me in plain English!

Much appreciated  2013-03-01 13:38:35

Hello! I’m happy that you’ve found these tutorials useful. Sadly I don’t know Processing. :(  Cube 2013-04-02 23:08:48

Hi,

I found a site with dft filtering tutorial about a year ago, and
i cant seem to find it now, it had black background with green letters and some pictures of some robots drawing on paper, very interesting and fun to read. There was everything i needed for audio processing including all the math material complex numbers and stuff like that. So i was wandering if it was your website because your name (Nonoo) sounds very familiar to me?  2013-04-03 07:05:45

Hello, it wasn’t my website. Maybe it was dspguide.com :)  Cube 2013-04-10 01:38:23

It’s not dspguide.com, maybe the material is, but thanks anyway… :) although I really think it was your webpage, very good presentation;)  ali zaidi 2015-02-18 10:09:19

i am developing an Applicaton of speech Recognition using java. can someone help me to remove the background noise which comes  ASaf 2016-10-01 01:09:59

Hello,

Is there any way to calculate the scaling factor “2000” in order to get the original amplitude , based on cutoff frequency amplitude
is changing ?

Thanks  qaif 2016-10-03 13:01:02

Name (required)
E-mail (required - never shown publicly)
Webpage URL
Comment:
You may use <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong> in your comment.