Presented here is a discussion of an infinite impulse response (IIR) filter structure used to perform spectrum analysis in the detection and measurement of single sinusoidal tones. The standard method for spectral energy is the discrete Fourier transform (DFT), typically implemented using a fast Fourier transform (FFT) algorithm.
However, there are applications that require spectrum analysis only over a subset of the N bin-center frequencies of an N-point DFT. A popular, as well as efficient, technique for computing sparse FFT results is the Goertzel algorithm, using an IIR filter implementation to compute a single complex DFT spectral bin value based upon N input time samples.
The most common application of this process is to detect the presence of a single continuous-wave sinusoidal tone. With that in mind, let’s look briefly at tone detection.
Tone detection It’s certainly possible to use the FFT to detect the presence of a single sinusoidal tone in a time-domain sequence x(n). For example, if we wanted to detect a 30 kHz tone in a time-domain sequence whose sample rate was f s = 128 kHz, we could start by performing a 64-point FFT as shown in Figure 13–41. Then we would examine the magnitude of the X(15) complex sample to see if it exceeds some predefined threshold.
Figure 13–41 DFT method, using an FFT algorithm, to detect a 30 kHz tone. This FFT method is very inefficient. In our example, we’d be performing 192, (64/2)(log2 64), complex multiplies to obtain the 64-point complex X(m) in order to compute the one X(15) in which we’re interested. We discarded 98% of our computations results! We could be more efficient and calculate our desired X(15) using the single-point discrete Fourier transform (DFT) in Eq. (13–77), which requires N = 64 complex multiplies using the following.
That would be an improvement but, happily, there’s a better way. It’s called the Goertzel algorithm (pronounced: girt’zel). The Goertzel algorithm is implemented in the form of a second-order IIR filter, with two real feedback coefficients and a single complex feedforward coefficient,as shown in Figure 13–42. (Although we don’t use this process as a traditional filter, common terminology refers to the structure as a filter.)
Figure 13–42 IIR filter implementation of the Goertzel algorithm. This filter computes a single-bin DFT output (the mth bin of an N-point DFT) defined by
The filter’s y(n) output is equal to the DFT output frequency coefficient, X(m), at the time index n = N, where the first time index value is n = 0. For emphasis, we remind the reader that the filter’s y(n) output is not equal to X(m) at any time index when n≠N.
To be equivalent to the DFT, the frequency domain index m must an integer in the range 0≥m≤N–1. You’re welcome to think of the Goertzel algorithm as a single-bin DFT. The z-domain transfer function of the Goertzel filter is
with a single z-domain zero located at z = e–j2πm/N and conjugate poles at z = e±j2πm/N as shown in Figure 13–43(a). The pole/zero pair at z = e–j2πm/N cancel each other.
Figure 13–43 Goertzel filter: (a) z-domain pole/zero locations; (b) frequency magnitude response . Having a filter pole on the unit circle is typically a risky thing to do for stability reasons, but not so with the Goertzel algorithm. Because it processes N+1-length blocks of time samples (where N is usually in the hundreds), the filter remains stable for such short time sequences because its internal data storage registers, w(n–1) and w(n–2), are reset to zero at the beginning of each new block of input data.
The filter’s frequency magnitude response, provided in Figure 13–43(b), shows resonance centered at a normalized frequency of 2πm/N, corresponding to a cyclic frequency of mfs /N Hz (where fs is the signal sample rate).
The Goertzel algorithm is implemented with a complex resonator having an infinite-length unit impulse response, h(n) = ej2?nm/N, and that’s why its frequency magnitude response is so narrow. The time-domain difference equations for the Goertzel filter are
An advantage of the Goertzel filter in computing an N-point X(m) DFT bin value is that Eq. (13–80) is implemented N times while Eq. (13–81), the feed forward path in Figure 13–42 need only be computed once after the arrival of the Nth input sample.
Thus for real x(n) inputs the filter requires N+2 real multiplies and 2N+1 real adds to compute an N-point X(m). However, when modeling the Goertzel filter if the time index begins at n = 0, the filter must process N+1 time samples with x(N) = 0 to compute X(m).
In typical applications, to minimize spectral leakage, we choose N so there’s an integer number of cycles in our input sequence of the tone we’re try ing to detect. N can be any integer, and the larger N the better the frequency resolution and noise immunity. However, larger N means more computation
It’s worth noting that while the typical Goertzel algorithm description in the literature specifies the frequency resonance variable m to be an integer (making the Goertzel filter’s output equivalent to an N-point DFT bin output), the m in Figure 13–42 and Eq. (13–79) can in fact be any value between 0 and N–1, giving us full flexibility in specifying our filter’s resonant frequency.
Goertzel Example Let’s use Goertzel to calculate the spectral magnitude of that ftone = 30 kHz tone from the Figure 13–41 example. When fs = 128 kHz and N = 64, then our resonant frequency integer m is
The Goertzel filter and the necessary computations for our 30 kHz detection example are provided in Figure 13–44.
Figure 13–44 Filter, coefficients, and computations to detect the 30 kHz tone. It’s useful to know that if we want to compute the power of X(15), X(15)2 , the final feedforward complex calculations can be avoided by computing: Click on image to enlarge.
In our example, Eq. (13–83) becomes
Click on image to enlarge.
Goertzel advantages over the FFT Here are some implementation advantages of the Goertzel algorithm over the standard radix-2 FFT for single tone detection:
*N does not need to be an integer power of 2. * The resonant frequency can be any value between zero and fs Hz. * The amount of filter coefficient (versus FFT twiddle factor) storage is reduced.If Eq. (13–83) is used, only one coefficient need be stored. * No storing a block of input data is needed before processing can begin (aswith the FFT). Processing can begin with first input time sample. * No data bit reversal is needed for Goertzel. * If you implement the Goertzel algorithm M times to detect M differenttones, Goertzel is more efficient (fewer multiplies) than the FFT whenM 2 N. * Computational requirements to detect a single tone (assuming real-onlyx(n) input) are given in Table 13–4.
Table 13–4 Single-Bin DFT Computational Comparisons As a final note, although the Goertzel algorithm is implemented with a complex resonating filter structure, it’s not used as a typical filter where we retain each output sample.
Figure 13–45. Goertzel algorithm frequency magnitude response. For the Goertzel algorithm we retain only every Nth, or (N+1)th, output sample. As such, the frequency magnitude response of the Goertzel algorithm when treated as a black-box process is equivalent to the sin(x)/x-like magnitude response of a single bin of an N-point DFT, a portion of which is shown in Figure 13–45 above.
Used with the permission of the publisher, Prentice Hall, this on-going series of articles on Embedded.com is based on copyrighted material from “Understanding Digital Signal Processing, Second Edition” by Richard G. Lyons. The book can be purchased on line.
Richard Lyons is a consulting systems engineer and lecturer with Besser Associates. As a lecturer with Besser and an instructor for the University of California Santa Cruz Extension, Lyons has delivered digital signal processing seminars and training course at technical conferences as well at companies such as Motorola, Freescale, Lockheed Martin, Texas Instruments, Conexant, Northrop Grumman, Lucent, Nokia, Qualcomm, Honeywell, National Semiconductor, General Dynamics and Infineon.
Continue Reading
