You probably have at least a nodding familiarity with the Fourier transform, a mathematical process for transforming a time-domain signal into a frequency domain signal. In particular, for computers, we don’t really have a nice equation so we use the discrete version of the transform which takes a series of measurements at regular intervals. If you need to understand the entire frequency spectrum of a signal or you want to filter portions of the signal, this is definitely the tool for the job. However, sometimes it is more than you need.
For example, consider tuning a guitar string. You only need to know if one frequency is present or if it isn’t. If you are decoding TouchTones, you only need to know if two of eight frequencies are present. You don’t care about anything else.
A Fourier transform can do either of those jobs. But if you go that route you are going to do a lot of math to compute things you don’t care about just so you can pick out the one or two pieces you do care about. That’s the idea behind the Goertzel. It is essentially a fast Fourier transform algorithm stripped down to compute just one frequency band of interest. The math is much easier and you can usually implement it faster and smaller than a full transform, even on small CPUs.
Fourier in Review
Just to review, a typical system that uses a Fourier transform will read a number of samples at some clock frequency using an analog to digital converter. For example, you might read 1024 samples at 1 MHz. Each sample will be 1 microsecond after the previous sample, and you’ll get just over a millisecond’s worth of data.
The transform will also produce 1024 results, often known as buckets. Each bucket will represent 1/1024 of the 1 MHz. Because of the Nyquist limit, the top half of the results are not very useful, so realistically, the first 512 buckets will represent 0 Hz to 500 Hz.
Here’s a spreadsheet (after all, this is the DSP Spreadsheet series of articles) that uses the XLMiner Analysis Addon to do Fourier transforms. The sample frequency is 1024Hz so each bucket (column G) is worth 1 Hz. Column B generates two sine waves together at 100 Hz and 20 Hz. We talked about generating signals like this last year.
The FFT data in column C are complex numbers and they do not update live. You have to use the XLMiner add on to recompute it if you change frequencies. The magnitude of the numbers appear in column F and you’ll see that the buckets that correspond to the two input frequencies will be much higher than the other buckets.
One Bucket of Many
If our goal was to find out if the signal contained the 100 Hz signal, we could just compute cell F102, if we knew how. That’s exactly what the Goertzel does. Here’s the actual formula for each bucket in a discrete Fourier transform (courtesy of Wikipedia):
If you apply Euler’s formula you get:
Mathematically, N is the block size, k is the bucket number of interest, xn is the nth sample. Since n varies from 0 to N-1, you can see this expression uses all of the data to compute each bucket’s value.
The trick, then, is to find the minimum number of steps we need to compute just one bucket. Keep in mind that for a full transform, you have to do this math for each bucket. Here we only need to do it for one. However, you do have to sum up everything.
Cut to the Chase
If you do a bit of manipulation, you can reduce it to the following steps:
- Precompute the frequency ratio in radians which is 2Πk/N
- Compute the cosine of the value you found in step 1 (call this C); also compute the sine (call this S)
- For each bucket in sequence, take the value of the previous bucket (P) and the bucket before that (R) and compute bucket=2*C*P-R and then add the current input value
- Now the last bucket will be set up to nearly give you the value for the kth bucket. To find the real part, multiply the final bucket by C and subtract P. For the imaginary part, multiply S times the final bucket value
- Now you want to find the magnitude of the complex number formed by the real and imaginary parts. That’s easy to do if you convert to polar notation (find the hypotenuse of the triangle formed by the real and imaginary part).
Of course, there’s a spreadsheet for this, too. With the limitations of a spreadsheet, it does a lot of extra calculations, so to make up for it, column D hides all the rows that don’t matter. In a real program you would simply not do those calculations (columns E-G) except for the last time through the loop (that is, on each row). You need a few samples to prime the pump, so the first two rows are grayed out.
There are other optimizations you’d make, too. For example, the cosine term is doubled each time through the loop but not on the final calculation. In a real program, you’d do that calculation once, but for a small spreadsheet, there’s not much difference.
Column B is the input data, of course. Column C is the running total of the buckets. Columns E-G on the last row represents the final computations to get the bucket k’s complex result. The real place to look for the results in around L11. That part of row 11 picks up the one visible part of columns E-G and also squares the magnitude since we are looking for power. There’s also a threshold and a cell that tells you if the frequency (at L1) is present or not.
Try changing the generated frequencies and you’ll see that it does a good job of detecting frequencies near 20 Hz. You can adjust the threshold to get a tighter or broader response, but the difference between 20 Hz and 19 or 21 Hz is pretty marked.
Real Code
Since a spreadsheet doesn’t really give you a feel for how you’d really do this in a conventional programming language here is some pseudo-code:
omega_r = cos(2.0*pi*k/N);
omega_r2=2*omega_r;
omega_i = sin(2.0*pi*k/N);
P = 0.0;
S = 0.0;
for (n=0; n<N; ++n)
{
y = x(n) + omega_r2*P -S;
S = P;
P = y;
}
realy = omega_r*P - S; // note that hear P=last result and S is the one before
imagy=omega_i*P;
magy=abs(complex(realy,imagy)); // mag=sqrt(r^2+i^2);
Since you can precompute everything outside of the loop, this is pretty easy to formulate in even a simple processor. You could even skip the square root if you are only making comparisons and want to square the result anyway.
Tradeoffs
Don’t forget, everything is a tradeoff. The sampling rate and the number of samples determine how long it takes to read each sample and also how broad the frequency range for the bucket is. If you are having to probe for lots of frequencies, you might be better off with a well-implemented Fourier transform library.
There are other ways to optimize both algorithms. For example, you can limit the amount of data you process by selecting the right sample rate and block size. Or, instead of using a Fourier transform, consider a Hartley transform which doesn’t use complex numbers.
In Practice
I doubt very seriously anyone needs to do this kind of logic in a spreadsheet. But like the other DSP spreadsheet models, it is very easy to examine and experiment with the algorithm laid out this way. Developing intuition about algorithms is at least as important as the rote recollection of them — maybe it is even more important.
The Goertzel is an easy way to sidestep a more complicated algorithm for many common cases. It is a tool you should have in your box of tricks.
Intel has a new, interesting vulnerability fix in their processors, 
