Goertzel Algorithm
Table of Contents
What is The Goertzel Algorithm?
Although many engineers have never heard of it, the Goertzel algorithm can detect tones with much less CPU power than the Fast Fourier Transform. This article makes an effort to remedy that.
The Fast Fourier Transform (FFT) is a tool that is familiar to most engineers, and they wouldn’t have any trouble using a “canned” FFT routine to identify one or more tones in an audio signal. Many people are unaware that there is a much faster method available if you only need to find a few frequencies, though. The Goertzel algorithm is what it’s called.
Tone Detection
Many applications require tone detection, such as:
- DTMF (touch tone) decoding
- Call progress (dial tone, busy, and so on) decoding
- Frequency response measurements, which involve sending a tone and simultaneously reading the outcome, can provide useful information if you do so for a variety of frequencies. You can determine whether load coils (inductors) are present on a telephone line by looking at its frequency response curve, for instance.
Although there are ICs specifically designed for the aforementioned applications, performing these tasks in software is more affordable. Unfortunately, a lot of embedded systems lack the processing power necessary to carry out continuous real-time FFTs. The Goertzel algorithm is used in this situation.
I discuss a basic Goertzel and an optimized Goertzel in this article.
The fundamental Goertzel provides real and fictitious frequency components, just like a standard Discrete Fourier Transform (DFT) or FFT. And The real/imaginary pair can then be used to compute magnitude and phase, if necessary. The real and imaginary frequency components are not provided by the optimized Goertzel, despite it being faster (and simpler) than the basic Goertzel. Instead, it provides you with the squared relative magnitude. There is no way to determine the phase, but you can take the square root of this result to determine the relative magnitude (if necessary).
Goertzel Algorithm
First, a brief explanation of the algorithm: Each sample is processed in the middle. Every Nth sample, the actual tone detection takes place. (I’ll expand on N in a moment.)
You work with blocks of samples, just like the FFT. This does not require that you process the data in blocks, though. If you receive an interrupt for each sample, the numerical processing can be completed within the interrupt service routine (ISR) that collects the samples. Alternatively, if you receive buffers of samples, you can process them one batch at a time.
Before you can do the actual Goertzel, you must do some preliminary calculations:
- Decide on the sampling rate.
- Choose the block size, N .
- Precompute one cosine and one sine term.
- Precompute one coefficient.
To save RAM and ROM space, you can either precompute each of these once before having it hardcoded into your program, or you can compute them as you go.
Sampling Rate
The application might have already selected your sampling rate. For instance, an 8kHz (8,000 samples per second) sampling rate is frequently used in telecom applications. Alternately, your analog-to-digital converter (or CODEC) might be powere d by an uncontrollable external clock or crystal. The standard Nyquist guidelines still apply if you get to pick the sampling rate: it must be at least twice as high as the frequency you are most interested in. I say “at least” because it’s possible that a higher sampling frequency will produce even better results if you are detecting multiple frequencies. You want every frequency of interest to be an integer factor of the sampling rate, not the sampling rate itself.
Block Size
Size in Goertzel blocks N resembles the number of equivalent FFT points. It regulates the bin width, which is another name for frequency resolution. For instance, if N is 100 samples and your sampling rate is 8 kHz, your bin width will be 80 Hz. This would encourage you to increase N as much as you can to achieve the highest frequency resolution. The problem is that as N increases, it takes longer to detect each tone because you have to wait longer for all of the samples to arrive. For instance, 800 samples will accumulate in 100ms at 8kHz sampling. You must use compatible values of N if you want to be able to detect short-duration tones.
The correlation between the sampling rate and the target frequencies is the third factor affecting your choice of N. And The frequencies should ideally be in the middle of their respective bins. The target frequencies should therefore be integer multiples of sample_rate/N. The good news is that N need not be a power of two, unlike the FFT.
Precomputed Constants
It only takes five easy steps to calculate the constants you’ll need during processing once you’ve chosen your sampling rate and block size:
w = (2*π/N)*k
cosine = cos w
sine = sin w
coeff = 2 * cosine
For the per-sample processing you’re going to need three variables. Let’s call them Q0 , Q1 , and Q2 .
Q1 is just the value of Q0 last time. Q2 is just the value of Q0 two times ago (or Q1 last time).
Q1 and Q2 must be initialized to zero at the beginning of each block of samples. For every sample, you need to run the following three equations:
Q0 = coeff * Q1 – Q2 + sample
Q2 = Q1
Q1 = Q0
After running the per-sample equations N times, it’s time to see if the tone is present or not.
real = (Q1 – Q2 * cosine)
imag = (Q2 * sine)
magnitude2 = real2 + imag2
You can determine whether the tone was present or not using a straightforward threshold test of magnitude. Set Q2 and Q1 back to zero and begin the following block.
An Optimized Goertzel
The optimized Goertzel sacrifices phase information in order to require less computation than the basic one. Processing for each sample is identical, but processing at the end of a block is different. You calculate the following directly rather than first converting the real and imaginary components into relative magnitude squared:
magnitude2 = Q1 2 + Q2 2 -Q1 *Q2 *coeff