Search This Blog

Everyday DSP for Programmers: FIR Filters

Now that we have this shiny new tool called the DFT, it's time to look more closely at Finite Impulse Response (FIR) filters. We can use the DFT to create and analyze FIR filters so that we can better understand how they behave with different signals.

To review, an FIR filter is named as such because when it is applied to an impulse function, the response is complete (i.e. it has reached and will remain at zero) within a finite number of samples. An FIR filter consists of an array of values, called taps, and to apply the filter to a signal, the taps are overlapped with a block of samples with each sample multiplied by a tap. Then, the results are added together to get a single filtered sample. The filter is advanced by one sample and the process is repeated for the next filtered sample. This process is called convolution, as we covered way back in the post on averaging.

We would use a filter when we want to remove some part of the signal before analyzing what's left. Normally we want to remove high-frequency noise, and a filter that does this is called a low-pass filter because it only allows low frequencies to get through it. We can also design high-pass filters that remove low-frequency bias from a signal or band-pass filters that only allow certain frequency ranges of interest through. All of these types of filters have similar design processes, so we'll focus on the low-pass filter. Let's try to design a perfect low-pass FIR filter.

Brick Wall Filter

A filter with a steep cut-off at a certain frequency can be referred to as a brick wall filter because it looks just like that, a brick wall. We can specify such a filter in the frequency domain by simply setting the frequency magnitudes to the values that we want them to be. Setting a value of one means the filter will pass that frequency, and setting a value of zero means the filter will remove that frequency.

How does it work to specify the frequency response of a filter in the frequency domain? Well, it turns out that convolution in the time domain (the way we apply a filter to a signal) is equivalent to a point-by-point multiply in the frequency domain. If we converted our signal to the frequency domain with a DFT and multiplied each signal frequency with its corresponding filter frequency, we would end up with the frequency response of the filtered signal. Then, we can get the filtered signal back in the time domain by doing the inverse DFT, because the DFT is a completely reversible operation.

What if, instead of doing these DFTs and inverse DFTs on the signal to perform the filtering in the frequency domain, we wanted to do the filtering in the time domain? We can do this by simply converting the filter in the frequency domain back into the time domain with an inverse DFT. Then we have a set of filter taps to use with the time domain convolution operation. Let's see what kind of taps we get if we specify the filter in the frequency domain and then convert it back to the time domain.

I'll be using the same JavaScript DSP library I've been using, and the conversion of a filter in the frequency domain to the time domain can easily be done using the FFT.inverse(real, imag) function:
function GenFilter(size, fCutOff) {
  var firSpectrum = new Array().fill(1,0,fCutOff).fill(0,fCutOff,size);
  var fft = new FFT(size, sampleRate);
  var zeros = new Array(size).fill(0);
  var taps = Array.prototype.slice.call(fft.inverse(firSpectrum, zeros));
  return taps.slice(size/2,size).concat(taps.slice(0,size/2))
} 
The fCutOff variable is the cut-off frequency of the filter. The real values for the FFT inverse are the FIR spectrum we defined while the imaginary values can simply be set as an array of zeros.

The two halves of the array of taps are swapped because the inverse FFT calculates the taps in the time domain in such a way that they get split between the two ends of the time domain block. This is because a discrete frequency domain signal converts to a periodic time domain signal, and the phase is determined by the relation of the real and imaginary parts of the spectrum. We could have figured out what to set the real and imaginary parts to so that the time domain signal would come out the way we wanted, but it's easier to swap the two halves of the array after the fact.

We'll start with a filter that will pass everything, and you can click on the graph to have it sweep back through the frequencies, setting them to zero to see how the taps change as the filter passes only lower and lower frequencies:


How very interesting. A filter that passes all frequencies in the frequency domain is the same as the impulse function in the time domain. This should be obvious in retrospect, since an impulse function used as a filter would return each sample of the signal without modification as it was swept across the signal.

As the higher frequencies are set to zero, the impulse function stretches out into an ever-widening and shortening sinc function until it appears to finally spread out to a set of zero-valued taps when only the DC value is passed. The taps are not actually zero-valued, though. They're just very small.

We can continually scale the height of the sinc function to better see what happens to the taps at lower cut-off frequencies. In this next graph I do that as well as zoom in on the center of the taps to show more detail of the sinc function. The graph also pauses at intervals so you can see what's going on, so keep clicking to advance the animation:


Here you can really see how the sinc function widens and eventually stretches into a DC signal when all of the frequencies are removed except the zero frequency. So is it this easy to create a perfect brick wall filter simply by using an inverse DFT on a desired filter frequency response? I bet you're now thinking, no, and you'd be right. What's wrong with this assumption?

Brick Wall Filter in High-Res

Remember, the frequency spectrum of a DFT only shows the values of a discrete set of frequencies. We currently don't know what's happening at the frequencies between the frequency bins shown in the graphs. We can find out what's happening at those frequencies by calculating a forward DFT on the sinc function, but at a higher resolution than we did the inverse DFT. To get a higher resolution, we zero-pad the taps. After zero-padding, a specific brick wall filter looks like this:


Both the time and frequency domains are zoomed in to show extra detail, and the zero-padding of the taps is not shown. Imagine that there are an equal number of zeros as taps added to the right of the sinc function. The white curve shows the ideal frequency response. (Note the DC magnitude is half the value of the other low frequencies because the inverse DFT actually doubles the DC value, so it would otherwise add too much DC offset to the taps.) The blue curve shows the high resolution frequency response, which includes a lot of ripples near the cut-off frequency.

You can see how those ripples are generated in the frequency domain from the sinc function by clicking on the graph. The first click shows how the frequency domain would look if only the main peak of the sinc function were used as the filter with all other taps set to zero (or removed). The frequency response is one smooth transition from passing frequencies to removing frequencies with a single bounce in the higher frequency range.

As you keep clicking, more valleys and peaks of the sinc function are added back in, and more ripples get squeezed into the frequency response. To squeeze the ripples in so tightly that they disappear, you would need to have an infinitely long sinc function with taps going on forever in both directions. It's like the frequency compliment to the sine wave approximation of the square wave, which requires an infinite summation of sines to get the square wave. In fact, if you run a square wave through one of these brick wall filters, you will get a square wave with ripples, just like a sine wave approximation:


It's amazing how all of this stuff is connected together! So how can we tame those ripples?

A Gentle Filter

One of the reasons the brick wall filter creates ripples in sharp transitions is that the sinc function ends abruptly on both ends. If we can soften the ends a bit, we can reduce the ripples. The way to do this is to multiply the filter taps by a window function. A window function is an array of scale factors that adjusts the values of whatever it's multiplied by inside of its range, and forces zero values outside of its range. Our filter is already implicitly multiplied by a window function, albeit a degenerate one—the rectangular window. A rectangular window has a value of one over its range and zero everywhere else.

There are tons of window functions to choose from, but we'll use a common one that's fairly easy to calculate, called the Blackman window. The equation for the Blackman window is

w(k) = 0.42 - 0.5 cos(2πk/(N-1)) + 0.08 cos(4πk/(N-1))

Where k is the index of the window value and N is the number of window values. We can modify our filter generation code to include the window function:
function GenBlackmanFilter(N, fCutOff) {
  var firSpectrum = new Array().fill(1,0,fCutOff).fill(0,fCutOff,size);
  var fft = new FFT(N, sampleRate);
  var zeros = new Array(N).fill(0);
  var taps = Array.prototype.slice.call(fft.inverse(firSpectrum, zeros));
  taps = taps.slice(N/2,N).concat(taps.slice(0,N/2))
  return taps.map(function(x,k) {
    return x*(0.42 - 0.5*Math.cos(2*Math.PI*k/(N-1)) + 0.08*Math.cos(4*Math.PI*k/(N-1)))
  });
} 
When we apply the Blackman window to the brick wall filter, it smooths out the frequency response of the filter pretty well, as seen in the following graph (click to switch the window in and out):


This window alone won't be enough to remove all of the ripples, though. The ripples in the sinc function itself are also causing the ripples in the square wave. To really smooth out the filter response, we need to use only the main peak of the sinc function for the filter. If we select only those taps of the filter before applying the Blackman window, the resulting filter gives us a nice smoothing of sharp transitions, like so:


This gentle filter no longer has a perfectly sharp cutoff in the frequency domain, but it doesn't create ripples on sharp transitions in the time domain. This is an inherent trade-off to filter design, and depending on your application, you may want a sharper or more gentle filter.

Averaging as an FIR Filter

With all of this discussion of filtering, it would be interesting to see what a regular old moving average filter has for a frequency response. The moving average can be thought of as an FIR filter with every tap equal to 1/N. To calculate the DFT, we simply need to zero pad the taps and convert. It looks like this:


By now you should be familiar with this picture. Of course, the frequency spectrum turns out to be half of a sinc function (the other half is over at the high end of the frequency spectrum, not shown). We can even apply the Blackman window to smooth out all of the ripples in the frequency response, giving us another type of gentle filter. You can click the graph to see what this looks like.

So the moving average is really a special case of an FIR filter, and can be analyzed in a similar manner. Depending on the application, we can design a sharper brick wall filter that will create ripples on sharp transitions, or we can design a smoother gentle filter that doesn't differentiate between frequencies at the cutoff as well but doesn't create ripples. There's much more to filter design, but next week we'll move on to another DSP algorithm that detects when a specific frequency is present in a signal.


Other DSP posts:
Basic Signals
Transforms
Sampling Theory
Averaging
Step Response of Averaging
Signal Variance
Edge Detection
Signal Envelopes
Frequency Measurement
The Discrete Fourier Transform
The DFT in Use
Spectral Peak Detection
FIR Filters
Frequency Detection
DC and Impulsive Noise Removal 

No comments:

Post a Comment