81б. The Goertzel algorithm: finding a tone without an FFT

In section 81а the signal was filtered: pass some frequencies, suppress others. In section 117 you took the whole spectrum at once — with the fast Fourier transform (FFT). Often you need less: is there a 697 Hz tone in the signal, is the 2 kHz buzzer sounding, has a DTMF digit arrived. For a few frequencies the FFT computes hundreds of unneeded bins. The Goertzel algorithm computes one bin — with a second-order recursion and one multiplication per sample. This section covers how it works, how to choose the block length and the bin number, what it costs compared with the FFT, how to avoid overflow in q31 and how to build from it a dual-tone multi-frequency (DTMF) detector that does not trigger on music.

The basis: our own PC model in pure Python (_verify/goertzel81c/goertzel81c.py) — the algorithm was checked against a direct DFT calculation, and the detector was run on thousands of blocks with noise and detuning. The C code was built with armclang for STM32F407 and STM32G474 without warnings; the same file built with gcc on a PC produced the same q31 powers as the model, bit for bit. The listing of the inner loop comes from the fromelf disassembler. It was not run on a board, and cycles were not measured.

Terms used in the documentation

Term	Meaning
Goertzel algorithm	one DFT bin computed by a second-order recursion
DFT bin	the frequency k·fs/N
frequency resolution	the bin width fs/N
spectral leakage	a tone spreads into neighboring bins
Hamming window	a smooth taper applied to the block
DTMF (dual-tone multi-frequency)	tone dialing: a digit is two tones at once
twist	the level difference between the two tones of a digit
second harmonic	the doubled frequency
guard bits, headroom	spare bits for the number to grow into

81б.1 One frequency instead of a spectrum

The DFT of a block of N samples in bin k is the sum X[k] = Σ x[n]·e^(−j2πkn/N). Goertzel gets its magnitude without storing the block and without computing sines on every sample. Each sample goes through one second-order stage:

s[n] = x[n] + c·s[n−1] − s[n−2],     c = 2cos(2πk/N)

After the last sample of the block, the bin power is obtained from the two state values s1 = s[N−1] and s2 = s[N−2]:

P = s1² + s2² − c·s1·s2 = |X[k]|²

The amplitude of a tone that falls into the bin is 2√P / N. Per sample — one multiplication and two additions; the coefficient c is constant and is computed once.

The PC model checked this against the direct DFT sum on a random block of 205 samples: in all bins from 0 to 102 the powers agreed to 11 digits. The bin number does not have to be an integer: with k = 17.861 (the exact 697 Hz frequency) Goertzel gave the same power as the DFT sum at that frequency.

81б.2 N and k: resolution and leakage

The bin width is fs / N, the block length is N / fs. Example — DTMF dialing: a digit is one tone from the low group (697, 770, 852, 941 Hz) and one from the high group (1209, 1336, 1477, 1633 Hz). With fs = 8000 Hz and N = 205 the block lasts 25.6 ms and a bin is 39.02 Hz. The closest pair is 697 and 770 Hz: 73 Hz, 1.87 bins — they can be told apart.

The bin number k = f·N / fs is rounded to an integer, and the tone ends up off the center of the bin. The model gives the worst-case amplitude loss over 16 tone phases, without a window:

Why exactly 205: the model tried every N from 100 to 250 and for each one found the worst offset of the eight tones from the bin centers. The best is N = 205, 0.269 bin; next come 113 (0.292) and 185 (0.297). The 256 samples usual for an FFT give 0.360. Goertzel does not need a power of two — you choose N to suit the frequencies.

Leakage is the other side of the same thing. A tone falls not only into its own bin but also into the neighboring ones. Without a window the 770 Hz tone outweighs its neighbor in the group by only 16.2 dB. Detuning makes it worse: if the frequency drifts by 1.5 % (for 1633 Hz that is 24.5 Hz, 0.63 bin), the worst-case amplitude loss is −11.7 dB.

81б.3 The window and the cost compared with the FFT

The window. The block is multiplied by a Hamming window before the stage — just as a Hann window was applied before the FFT in section 117. The bin becomes wider and more tolerant of detuning: at ±1.5 % the worst-case loss is −4.4 dB instead of −11.7, and the tone's margin over its neighbor in the group is 27.3 dB instead of 16.2. The price is one multiplication per sample, shared by all frequencies, and a window table: 205 q15 numbers — 410 bytes. The window also reduces the tone itself: the window gain for a tone is 0.538, and the amplitude is divided by it before comparison with the threshold.

Multiplications. Goertzel for one frequency is N multiplications in the loop and 3 at the end (the last term is conveniently computed as s2·(s2 − c·s1)). A simple radix-2 FFT is N/2·log2N complex multiplications, each of them four real ones:

Per frequency Goertzel stores three numbers: s1, s2 and the coefficient — 12 bytes. There is no need to store the block of samples: each sample goes into the stages immediately, even from the ADC interrupt. The FFT needs the whole block in memory plus twiddle tables — in CMSIS-DSP for 256 samples that is 2920 bytes of Flash (section 117).

At the same N = 256, Goertzel for one frequency is 259 multiplications, and the FFT becomes more efficient from about 16 frequencies (2·log2N). For the eight DTMF tones Goertzel is cheaper — and it works with N = 205, which a radix-2 FFT cannot take.

81б.4 Fixed point: headroom for growth

The stage is a resonator: at its own frequency the state grows from sample to sample. By how much depends on k. The bound on growth over a block is the sum of the absolute values of the stage's impulse response; the model computed it for all 16 bins: the largest at k = 18 — 249.0 (almost 2⁸), the smallest — 131.2. A full-scale sine exactly at the bin frequency drives the state up to 195.6 — the bound has margin.

In the code the q15 input is shifted into the int32 state by 7 bits. The worst case is the product c·s1, where c is almost 2: 32768·2^7·249·2 = 2 088 763 392 — less than 2³¹ − 1 = 2147483647. With a shift of 8 it comes to 4177526784 — it does not fit. This condition is in the code as a _Static_assert: the build fails if someone increases the shift.

The model checked full scale on bin 18 with an exact copy of the q31 code:

Accuracy. A tone at different levels (697 and 941 Hz, 8 phases each) was compared with the same calculation in double on the same q15 samples. The second variant is a state in int16: for growth of 249 times to fit in 16 bits, the input has to be shifted right by 9, and of the 16 bits of the sample only 7 remain.

q31 with a 64-bit product is accurate over the whole range of a q15 sample; a q15 state is only good for loud tones. Float on Cortex-M4 has a 24-bit mantissa, and on a weak bin next to a strong tone it is worse than q31: on bin 20 with full scale on bin 18 (117.8 dB lower) float on the PC was off by a factor of 2.6, q31 by 0.4 %. For DTMF thresholds this does not matter; for measuring a weak harmonic next to a strong one it does.

81б.5 C code and what the compiler does

The block in q31. The coefficient is stored as cos·2³¹ (2cos does not fit in q31), and the product is shifted by 30 instead of 31 — that is exactly the multiplication by 2cos:

int64_t goertzel_q31(const int16_t *x, uint32_t n, int32_t cq)
{
    int32_t s1 = 0, s2 = 0;
    GZ_ONE_BY_ONE
    for (uint32_t i = 0; i < n; i++) {
        int32_t s0 = (int32_t)(((int64_t)cq * s1) >> 30) - s2 + x[i] * (1 << GZ_SHIFT);
        s2 = s1;
        s1 = s0;
    }
    int32_t t = (int32_t)(((int64_t)cq * s1) >> 30);       /* 2cos·s1 */
    return (int64_t)s1 * s1 + (int64_t)s2 * s2 - (int64_t)t * s2;
}

The power is in units of 2⁴⁴ (twice the 15 bits of the sample plus the 7-bit shift): P = result / 2⁴⁴. Next to it in the file are a float variant, a table of 16 bin numbers and coefficients (tones and second harmonics) and two _Static_assert checks on the headroom.

F407 and G474 have the same Cortex-M4 core with an FPU (in the pack headers both have core_cm4.h and __FPU_PRESENT 1), and armclang with -O2 produced identical code for them. On its own it unrolls the loop by four samples; GZ_ONE_BY_ONE is a pragma that disables this so that the listing shows one sample. The inner loop of goertzel_q31 from fromelf:

цикл:   MOV    lr,r12
        SMULL  r4,r12,r12,r2
        LSRS   r4,r4,#30
        ORR    r12,r4,r12,LSL #2
        LDRSH  r4,[r0,#2]!
        SUB    r3,r12,r3
        ADD    r12,r3,r4,LSL #7
        SUBS   r1,#1
        MOV    r3,lr
        BNE    цикл

Ten instructions per sample. Six of them do useful work: SMULL — the 64-bit product cos·s1, LSRS and ORR — the shift by 30 from the two halves, LDRSH — the sample, SUB and ADD — subtract s2 and add the sample shifted by 7. The two MOVs move the state, SUBS and BNE are the counter. In float it is nine instructions: VLDR, VMUL, VADD, VSUB and five housekeeping ones. The compiler itself rewrote the end of the block as s1² + s2·(s2 − c·s1). For 16 bins and a block of 205 samples that is 32800 loop instructions per 25.6 ms.

The check on the PC: the same file built with gcc and -DGZ_HOST_TEST processed 19 blocks — 16 digits, a weak tone, full scale and a sawtooth — over 16 bins. All 304 q31 powers matched the model bit for bit.

81б.6 The DTMF detector and false triggers

Eight stages give eight amplitudes. The model's rules:

• in each group the strongest tone is taken; both are above a threshold of 0.05 of full scale;
• the tone outweighs the second strongest in its group by at least 10 dB;
• the twist between the tones of a digit is no more than 8 dB;
• optionally: the second harmonic of a tone is 15 dB below it, or the two tones carry at least half of the total block power.

Runs of 256 digits each: exact frequencies; frequencies with a deviation of up to ±1.5 %; white noise with a signal-to-noise ratio of 10 dB and 0 dB. Tones at 0.3 of full scale, twist up to ±4 dB, random phases. False triggers — a sawtooth from 60 to 500 Hz in 1 Hz steps (441 signals) with harmonics up to 3900 Hz: a signal rich in harmonics, like voice or music.

What this shows:

• The window cures detuning: 256 digits out of 256 at ±1.5 % instead of 229. But the wider bin also lets in more foreign harmonics — 43 false triggers instead of 12.
• Checking the second harmonic of both tones rejects valid digits even without noise: 240 out of 256. The harmonics of the low tones, 1394, 1540 and 1704 Hz, lie 1.77, 1.15 and 2.15 bins from the high tones 1336, 1477 and 1633 Hz — the harmonic's bin picks up leakage from the real high tone. Checking only the high tone does not suffer from this and cuts the 43 false triggers down to 17.
• The power share removes false triggers completely, but is stricter with noise: at 0 dB — 53 digits out of 256 versus 180.
• Wrong digits — not a single one in all runs: the detector would rather stay silent than make a mistake.

The detector in q31 with a q15 window made the same decisions on 768 digits as the one in double.

81б.7 Pitfalls

The tone is there, but the detector is silent. The frequency has drifted from the bin center. Without a window a 1.5 % drift costs up to 11.7 dB, with a Hamming window 4.4 dB.

N = 256 "like for the FFT". The DTMF tones land up to 0.36 bin from the centers instead of 0.27 with N = 205. Goertzel does not need a power of two.

The amplitude threshold "shifted" after adding the window. The window reduces the tone: gain 0.538. Divide the amplitude by it or compute the threshold with the window.

Music dials digits. The level and margin rules alone let the sawtooth through: 43 false triggers out of 441 with a window. You need the second harmonic or the power share.

The harmonic check rejects valid digits. The 1540 Hz harmonic of a low tone is 1.15 bins from the high 1477 Hz tone: its bin sees the leakage of the high tone. Check the harmonic of the high tone only, or use the power share.

The q31 power is garbage on a loud signal. The input shift is greater than 7: the state grows up to 249 times, and c·s1 leaves int32 — an error of 812 % with a shift of 8. Guard the headroom with _Static_assert.

Weak tones get lost in q15. An int16 state requires shifting the input by 9: at −60 dB the power error is 20.9 dB. Keep the state in int32 and the product in 64 bits.

The power grows from block to block. s1 and s2 were not zeroed before a new block: the stage continues the previous block.

81б.8 Summary

• Goertzel is one DFT bin: s[n] = x[n] + 2cos(2πk/N)·s[n−1] − s[n−2], and at the end P = s1² + s2² − c·s1·s2. Checked against the DFT sum to 11 digits; k can be fractional.
• DTMF at 8 kHz: N = 205, block 25.6 ms, bin 39.02 Hz, tones in bins 18…42 no further than 0.27 bin from the centers — the best N from 100 to 250.
• Hamming window: with ±1.5 % detuning the loss is 4.4 dB instead of 11.7, the margin over the neighbor 27.3 dB instead of 16.2.
• Cost: 8 tones — 1664 multiplications per block and 96 bytes for the state and coefficients; an FFT of 256 — 4096 multiplications and 2048 bytes of buffers. Break-even is around 16 frequencies.
• q31: growth up to 249 times — an input shift of 7 bits; with a shift of 8 the error is 812 %; a q15 state loses 20.9 dB on a −60 dB tone. The loop is 10 instructions per sample (float — 9).
• Detector: window + power share — 0 false triggers on 441 sawtooth signals and 256 of 256 digits at ±1.5 %; checking the harmonic of the low tone rejects valid digits.

Full sources

goertzel81c.c — the Goertzel algorithm: a block in q31 and in float.

/* Раздел 81б. Алгоритм Гёрцеля: блок в q31 и во float.
 *
 * Для STM32F407 и STM32G474 (Cortex-M4F) собирается armclang без предупреждений;
 * с -DGZ_HOST_TEST — gcc на ПК: читает блоки отсчётов из stdin и печатает мощности,
 * проверка сверяет их с моделью goertzel81c.py. */
#include <stdint.h>

#define GZ_N      205u      /* отсчётов в блоке: 25,6 мс при 8 кГц */
#define GZ_BINS   16u       /* 8 тонов DTMF и 8 вторых гармоник */
#define GZ_SHIFT  7         /* вход q15 сдвигается в состояние на 7 бит */
#define GZ_GROWTH 250       /* рост состояния за блок, с запасом: модель даёт 249,0 при k = 18 */

/* цикл по отсчёту без разворота: так листинг показывает ровно один отсчёт;
 * -DGZ_UNROLL оставляет armclang его обычный разворот */
#if defined(__clang__) && !defined(GZ_UNROLL)
#define GZ_ONE_BY_ONE _Pragma("clang loop unroll(disable)")
#else
#define GZ_ONE_BY_ONE
#endif

/* 2cos·s1 должно поместиться в int32 — иначе переполнение знакового, неопределённое поведение */
_Static_assert((int64_t)32768 * (1 << GZ_SHIFT) * GZ_GROWTH * 2 <= INT32_MAX, "запас по переполнению");
_Static_assert((int64_t)32768 * (1 << (GZ_SHIFT + 1)) * GZ_GROWTH * 2 > INT32_MAX, "сдвиг 8 уже не помещается");

/* номера бинов k и cos(2πk/N)·2^31: тоны 697…1633 Гц, затем их вторые гармоники */
const uint8_t gz_k[GZ_BINS] = {18, 20, 22, 24, 31, 34, 38, 42, 36, 39, 44, 48, 62, 68, 76, 84};
const int32_t gz_cos_q31[GZ_BINS] = {
    1828878034, 1756489927, 1677503679, 1592215995, 1248909980, 1083228028, 848335877, 600708818,
    967599315, 787480650, 473275393, 213560429, -694828920, -1054685608, -1477235229, -1811414887,
};

int64_t goertzel_q31(const int16_t *x, uint32_t n, int32_t cq)
{
    int32_t s1 = 0, s2 = 0;
    GZ_ONE_BY_ONE
    for (uint32_t i = 0; i < n; i++) {
        int32_t s0 = (int32_t)(((int64_t)cq * s1) >> 30) - s2 + x[i] * (1 << GZ_SHIFT);
        s2 = s1;
        s1 = s0;
    }
    int32_t t = (int32_t)(((int64_t)cq * s1) >> 30);       /* 2cos·s1 */
    return (int64_t)s1 * s1 + (int64_t)s2 * s2 - (int64_t)t * s2;
}

float goertzel_f32(const float *x, uint32_t n, float c)
{
    float s1 = 0.0f, s2 = 0.0f;
    GZ_ONE_BY_ONE
    for (uint32_t i = 0; i < n; i++) {
        float s0 = x[i] + c * s1 - s2;
        s2 = s1;
        s1 = s0;
    }
    return s1 * s1 + s2 * s2 - c * s1 * s2;
}

#ifdef GZ_HOST_TEST
#include <math.h>
#include <stdio.h>

int main(void)
{
    int16_t q[GZ_N];
    float f[GZ_N];
    char name[32];
    int blocks = 0;
    while (scanf("%31s", name) == 1) {
        for (uint32_t i = 0; i < GZ_N; i++) {
            int v;
            if (scanf("%d", &v) != 1)
                return 1;
            q[i] = (int16_t)v;
            f[i] = (float)v / 32768.0f;
        }
        for (uint32_t b = 0; b < GZ_BINS; b++) {
            float c = (float)(2.0 * cos(2.0 * 3.14159265358979323846 * gz_k[b] / GZ_N));
            printf("%s %u q31 %lld f32 %.9e\n", name, (unsigned)gz_k[b], (long long)goertzel_q31(q, GZ_N, gz_cos_q31[b]),
                   (double)goertzel_f32(f, GZ_N, c));
        }
        blocks++;
    }
    printf("блоков %d\n", blocks);
    return 0;
}
#endif

Next: section 82 "Indication on a prototype: PWM, RGB, shift registers".