Goertzel and the chirp z-transform

Two loops, each computing one DFT bin, hear a 770 Hz tone. Watch the one tuned near 770 Hz climb while the other wanders below 20.

A loop that rings at one frequency

Your phone’s keypad sends two tones per key. The exchange must know which of 8 frequencies are present, 25 ms at a time. An FFT computes 256 bins to read 8.

This page computes one bin at a time. I start from the sum of The DFT (13.2), where $X[k]$ is the walking vector: each sample is multiplied by a probe that turns by $\Omega_k=2\pi k/N$ per step, and the products are added. Now run that sum as a loop. Start from $y_k[-1]=0$ and, for each new sample, write

This is a loop of the kind in Difference equations (6.1): each output comes from the last one. The factor $e^{j\Omega_k}$ is an arrow of length 1 (Complex exponentials, 3.4), so each step spins the running total by $\Omega_k$ and does not make it bigger or smaller. Then the new sample is added.

After $N$ samples the loop holds

Because $\Omega_kN=2\pi k$, the factor $e^{j\Omega_k(N-1)}$ equals $e^{-j\Omega_k}$, so $y_k[N-1]=e^{-j\Omega_k}X[k]$. The loop gives $X[k]$ turned by a fixed angle: $\lvert y_k[N-1]\rvert=\lvert X[k]\rvert$, and $X[k]=e^{j\Omega_k}y_k[N-1]$.

Why does it ring at one frequency? If the input matches the loop’s frequency, each new sample lands in step with the spinning total, and the total grows steadily. If not, the additions drift in and out of step, and the total stays small. It is pushing a child on a swing: in time with the swing, the swing climbs; out of time, your pushes cancel.

The picture at the top of the page runs two such loops on a 770 Hz tone, one tuned to bin 20 (780.5 Hz) and one to bin 18 (702.4 Hz). Both start with the first sample, at size 1. Step or pause on the halfway frame, $n=102$: the loop tuned near 770 Hz has grown to 49.93, and the other is at 7.54 and going nowhere. After 205 samples they are 91.18 and 13.65. These are $\lvert X[20]\rvert$ and $\lvert X[18]\rvert$, the same as the FFT gives.

A tone that sits exactly on a bin gives the largest possible answer. Unit-amplitude cosine on bin 20 gives $\lvert X[20]\rvert=N/2=102.5$. The 770 Hz tone is between bins 18 and 20, so the loop at 20 reaches 91.18, and the loop at 18 reaches only 13.65.

The real form and the cost

The complex loop needs a complex multiply per sample, which is 4 real multiplies (Complex numbers for signals, 3.3). A rearrangement avoids that. Use two delays and a real coefficient:

That is one real multiply per sample. This is Goertzel’s algorithm. After the last sample, two numbers are enough:

When you need only the size, no complex arithmetic remains:

For the 770 Hz tone and $k=20$ the coefficient is $2\cos\Omega_{20}=1.6359$, and the loop ends with $q[204]=-118.605$ and $q[203]=-157.493$. Then $X[20]=e^{j\Omega_{20}}q[204]-q[203]=60.483-68.236j$, which is what the FFT gives.

A loop like this, which keeps ringing at one frequency, is a resonator. Resonators, notches and combs (17.4) designs them properly, and Transfer functions, poles and zeros (16.3) explains why it rings. I leave the poles until then.

Now the cost. Eight bins on 205 samples take $8\times205=1640$ real multiplies, plus a few at the end of each bin. A 256-point FFT has $\tfrac N2\log_2N=1024$ complex multiplies (The FFT, 14.1), which is $4\times1024=4096$ real multiplies. Its bins also sit at multiples of 31.25 Hz, not where the keypad tones are.

Per bin, Goertzel costs $N$ real multiplies and the FFT costs $4\cdot\tfrac N2\log_2N=2N\log_2N$ in total. So Goertzel wins when you want fewer than about $2\log_2N$ bins: 16 of them for $N=256$.

Which key was pressed

Each key sends one row tone (697, 770, 852 or 941 Hz) and one column tone (1209, 1336, 1477 or 1633 Hz). This is the DTMF scheme of ITU-T Recommendations Q.23 and Q.24. With $f_s=8$ kHz and $N=205$, which is 25.6 ms of sound, the bins are $8000/205=39.0$ Hz apart. Rounding $k=Nf/f_s$ gives bins 18, 20, 22, 24, 31, 34, 38 and 42, at 702.4, 780.5, 858.5, 936.6, 1209.8, 1326.8, 1482.9 and 1639.0 Hz. Each is within 11 Hz of its tone.

Use eight loops, one per tone. The strongest row loop and the strongest column loop name the key. It is how a voicemail menu hears “press 1”. The next picture dials four keys by itself. Watch the bars: exactly two stand out each time, one in each group.

Key 1 sends 697 Hz and 1209 Hz, and those two loops stand out. For key 5, 770 and 1336 Hz stand out at 45.9 and 46.9, and every other loop stays below 7. Each key needs 8 loops × 205 samples: 1640 multiplies, for 25.6 ms of sound.

Then it is yours. Press a key, or type it, and the bars regrow for that key. Try key 0: its bars are about 51.6 at 941 Hz and 46.4 at 1336 Hz, and the other six stay below 4.

The sixteen keys are told apart because every pair of tones sits on its own pair of bins. The bars of the other six loops are small but not zero, because a tone between bins leaks a little into the rest. I come back to leakage in Windowing and spectral leakage (15.1).

Zoom into a band

Now the opposite need. The DFT samples the spectrum at $N$ evenly spaced points (Sampling the spectrum, 13.1). Suppose you want to see the curve finely between 1000 and 1080 Hz. You could zero-pad to 8192 points (Zero-padding and resolution, 15.3), which computes 8192 bins to use 81.

Instead, compute only the points you want. Choose a lowest frequency $f_L$, a highest frequency $f_H$, a step $f_\text{step}$ and a number of points $N_\text{zoom}$. Evaluate the spectrum at

Each point is a DFT-like sum with $e^{-j2\pi(f_L+kf_\text{step})n/f_s}$ in it. More FFT algorithms (14.2) showed that Bluestein’s chirp trick turns sums of this form into one convolution, which an FFT computes. The same trick handles all the zoom points at once. The next picture zooms into 80 Hz around a 1037 Hz tone. Watch the peak readout move from the nearest FFT bin to the tone.

The FFT’s nearest bin is 33, at 1031.25 Hz: the tone is somewhere between the bins. The zoom’s 81 points, 1 Hz apart from 1000 to 1080 Hz, put the peak at 1037 Hz.

Compare the two peak readouts. The FFT’s nearest bin is 5.75 Hz from the tone at 1037 Hz, while the 1 Hz grid lands 0.06 Hz from the true top of the curve at 1037.06 Hz. The zoom did not find a new tone. Like a map app, zooming in shows the same coastline in more detail; it does not add islands. In the same way, two tones closer than about 31 Hz still blur into one hump here, because 256 samples limit what the curve can show.

Afterwards the band is yours. Drag it sideways, or use the arrow keys, and the zoom recomputes. When the band leaves out 1037 Hz, the highest point sits at an edge of the band.

The cost of one zoom with Bluestein’s trick is three FFTs of 512 points, the first power of two at least $256+81-1=336$. That is 6912 complex multiplies, against 53 248 for the 8192-point zero-padded FFT.

The general version evaluates the spectrum along any spiral, not only along the circle you used here. The spiral lives in the z-plane of The z-transform (16.1), and that is why the method is called the chirp z-transform.

Worked example

Every number below was recomputed in Python and compared with numpy.fft.fft.

1. DTMF bins ($f_s=8$ kHz, $N=205$, spacing 39.02 Hz): $k=\mathrm{round}(Nf/f_s)$ gives 18, 20, 22, 24, 31, 34, 38 and 42, at 702.4, 780.5, 858.5, 936.6, 1209.8, 1326.8, 1482.9 and 1639.0 Hz.
2. Goertzel on a 770 Hz tone (amplitude 1): $\lvert X[20]\rvert=91.183$, $\lvert X[18]\rvert=13.646$ and $\lvert X[24]\rvert=6.132$, each equal to the FFT bin. In real form for $k=20$: $2\cos\Omega_{20}=1.6359$, $q[204]=-118.605$, $q[203]=-157.493$, and $X[20]=60.483-68.236j$.
3. A tone exactly on bin 20 gives $\lvert X[20]\rvert=N/2=102.5$.
4. Cost. Eight bins by Goertzel: 1640 real multiplies, plus 3 per bin at the end. A 256-point FFT: 4096 real multiplies. A direct DFT of 8 bins on real input: $2\cdot205\cdot8=3280$ real multiplies.
5. Zoom. A 1037 Hz tone, $N=256$, $f_s=8$ kHz. FFT bins 32, 33 and 34 sit at 1000, 1031.25 and 1062.5 Hz, with $\lvert X\rvert$ of 19.16, 121.34 and 26.93. The zoom from 1000 to 1080 Hz, 81 points, peaks at 1037.0 Hz with $\lvert X\rvert=128.59$. SciPy’s zoom_fft (with endpoint=True) matches the direct sum to about $10^{-12}$.

Where you’ll meet this

Resonators and poles on the unit circle come in 16.3 and 17.4. Zero-padding and what resolution really means come in 15.3. Windowing, which tames the small leakage in the DTMF bars, is 15.1. Detection thresholds and false-alarm rates return in Chapter 25.

Reference card