Everyday DSP for Programmers: Frequency Detection

Up until now we've seen two ways to detect the frequency of a signal. If the signal has a dominant frequency, we can measure that frequency by counting zero crossings, and if there's more frequency content, we can use a DFT to measure the full frequency spectrum of the signal. If we're looking for frequency peaks, we can use spectral peak detection to calculate exact frequency values from the quantized frequency bins of the DFT.

Now we'll explore another way to detect a specific frequency in a signal using a type of filter called an Infinite Impulse Response (IIR) filter. There is a huge variety of IIR filters—the exponential average is an example of one that we've already seen—and in this case, the type of filter we'll use for detecting frequencies is called a complex resonator. We can use the DFT to help analyze the frequency response of the complex resonator, but first let's see why we would want to use it.

The Long Way (Using a DFT)

If we're looking for a specific frequency in a signal that's full of different frequencies, we can certainly use a DFT to convert the time-domain signal to the frequency domain and inspect the magnitude of the frequency bin of interest. That's a well-defined use of a DFT at this point, and it will give you the frequency content of your signal from the DC end up to half the sample rate, as shown in this graph (click the graph to sweep through frequencies):


But if we're only interested in a single frequency (or even a handful of frequencies), this is an awfully big hammer to use. For every block of samples of length N, we need to do O(N log N) operations in an FFT to calculate the frequency spectrum. We end up with N/2 frequency magnitudes when we really only need a couple of them, and we needed a lot of storage space for the N samples and the coefficients in the FFT algorithm. Can we do better? Absolutely.

The Short Way (Using a Complex Resonator)

Instead of calculating on a block of samples, like the moving average needs to do for averaging, we want a method of detecting a frequency that only uses a little bit of signal history, like the exponential average does. We are in luck. The method that we want to use actually looks fairly similar to an exponential average, but instead of multiplying the previous value by an exponential weight, like we do with the average, we multiply the previous value by a complex angle, like so

y(t) = x(t) + y(t-1)*(cos ω + i sin ω)

Where x(t) is the incoming sample, y(t-1) is the previous output of this filter, i is the imaginary number √(-1), and ω is the angle corresponding to the frequency we want to detect. To find this angle, we use the formula

ω = 2πk/N

Where k is the frequency bin corresponding to the desired frequency in an N-point DFT. But, we don't have to limit ourselves to the frequency bins of a DFT because we're not actually using a DFT. Another way to state this angle is

ω = 2πf0/fs

Where f0 is the desired frequency, and fs is the sampling frequency. This is another advantage of the complex resonator over the DFT method—we are not limited to quantized frequencies. We can detect any frequency in the range of 0 to fs. We can also implement this resonator as a simple JavaScript function that takes in samples, the previous output, and the resonator's pole, and returns the next filter output:
var iir_out = [0,0];
var angle = 2*Math.PI*16/512;
var atten = 0.995;
var pole = [Math.cos(angle)*atten, Math.sin(angle)*atten];

function complexResonator(x, y, pole) {
  return [y[0]*pole[0] - y[1]*pole[1] + x,
          y[0]*pole[1] + y[1]*pole[0]];
}

iir_out = complexResonator(sample, iir_out, pole);
var iir_mag = iir_out[0]*iir_out[0] + iir_out[1]*iir_out[1];
I could have used a complex number library for this code, but the calculation is simple enough that it should be clear what's going on. The pole is the complex angle that we're using to detect the frequency we want. In this example, the frequency to detect is 16 Hz, and the sample rate is 512 Hz. The output of the complex resonator (which is an IIR filter) is initialized to the complex number (0,0), and we can see that an attenuation factor is also being used in the definition of the pole. This is because the resonator has an infinite response, and if we didn't attenuate the output, as soon as it saw the desired frequency, it would resonate forever! We'll analyze the response more in a bit.

So this resonator works by multiplying the previous value of the filter by a complex angle and adding the new input sample to the result. This causes the output of the filter to naturally oscillate at the desired frequency, but it rejects all other frequencies. To plot the output of the filter, we take the magnitude of the complex output. The result when a signal is applied looks like this:


The input signal is in blue, and it will ramp up and down in frequency between 8 Hz and 24 Hz. You'll know when the input hits 16 Hz because the red signal will start rising, and then it will fall again when the input signal is no longer 16 Hz. As long as a 16 Hz frequency is present, the complex resonator output will continue to increase, so care should be taken to not overflow the filter output. The output should be saturation limited at some value.

Frequency Response of an IIR Filter

Since introducing the DFT, we've seen the frequency response of a number of different signals and FIR filters, but how can we use the DFT to better understand IIR filters? An IIR filter essentially goes on forever because it has an infinite response to an input, while a DFT has to be done on a finite set of samples. The best we can do is take some number of samples of the response and perform the DFT on those samples to see what the initial frequency response of the IIR filter is.

The next issue that we have to deal with is how to run a DFT on an IIR filter. For an FIR filter it was easy. The taps of the filter provide the samples to use for the DFT, and it's frequency response is simply the spectrum that results from calculating the DFT of the filter taps. For the IIR filter, those taps don't exist. The key is in the name of the filters. They are named impulse response filters. For the FIR filter, the taps are equivalent to the impulse response. For the IIR filter, we can get a set of samples to run the DFT on by running an impulse function through the filter for some number of samples. Let's try this for the exponential average function and see what we get:


The time-domain signal on the top shows the impulse response of the exponential average. It's magnified in the y-direction to show more detail. The bottom plot shows how the frequency response evolves as more samples are generated from the filter with an impulse as input. The frequency spectrum is also magnified, both in the y-direction and the frequency direction, to show more detail. The DFT was also zero-padded so that it's shown in higher resolution.

It's easy to see why the exponential average works, since it only allows very low frequencies through, and the frequency response quickly falls off as the frequency increases. However, this response is for samples that are fairly far along in the filter. If you pause the animation near the beginning, you can see that for samples that have only been cycled through the exponential average a few times, the frequency response is much flatter, and more of the higher frequencies come through. This is why recent samples have a much greater effect on the exponential average than on other types of averaging.

Frequency Response of the Complex Resonator

Now that we know how to analyze the frequency response of an IIR filter, what does the complex resonator look like in the frequency spectrum? Here's a graph of the impulse response of both a pure and an attenuated (by 0.995) resonator:


The pure resonator is in blue, and the attenuated resonator is in red. We can clearly see that the pure resonator is going to go on generating that sine wave forever, and that's because it was hit with an impulse, which has infinite frequency content. An impulse will excite every single resonator, no matter what frequency it resonates at, and it will continue forever unless it is attenuated. This even applies to a DC resonator. If the resonant frequency is zero, an impulse will simply be converted to a constant DC value, generating a step function. The equation for the filter reduces to

y(t) = x(t) + y(t-1)

The frequency response in the graph shows that the pure resonator has a sharp peak at the resonant frequency. (Again, the DFT was zero-padded to show more resolution.) The attenuated resonator shows a much lower peak, and it's more smoothed out as well. This reveals a trade-off with attenuation. Like every other design decision, you don't get something for nothing. Less attenuation will give you a sharper filter that responds more strongly to a specific frequency, but it will also respond more slowly and hold onto that response for a much longer time. A higher attenuation filter will respond more quickly to the presence or absence of a frequency, but it will respond more generally to frequencies around the desired frequency. The only way to improve both the sharpness and the response time of the filter is to increase the sample rate, and that has its own trade-offs.

The complex resonator gives us a great way to detect specific frequencies much more efficiently than the DFT, and we can even use several resonators at once to detect multiple frequencies in a straightforward manner. The complex resonator is only one of a wide variety of IIR filters that perform all kinds of tasks. There is much more to them, and if you're interested in learning more about them or any of these other DSP topics, I would highly recommend picking up Richard Lyons' excellent book, Understanding Digital Signal Processing. I learned much of what I know about DSP from this book, and it's quite well-written.

We're almost done with this series. I'll wrap up next week with a couple efficient ways to remove the DC value or impulsive noise from a signal, and then it'll be time to write about something else.


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