One Bin at a Time: The Goertzel Resonator

My homemade bat detector has a flaw I’m oddly fond of. The front-end transducer is a cheap 40 kHz ultrasonic ranging piezo — a resonant lump of ceramic that rings loudly at its favourite frequency and stays nearly deaf everywhere else. It’s a mechanical band-pass filter I got for free, and it means my “broadband” detector really only listens to a 2-3 kHz sliver around 40 kHz. A pipistrelle at 45 kHz sails past unheard.

The software version of that ringing ceramic is the Goertzel algorithm, and once you see the shape you can’t unsee it. A full FFT hands you every frequency bin at once — hundreds of them, most of which you’ll throw away. But if you already know which frequency you’re hunting for (a species’ peak call, a touch-tone digit, a carrier), computing the entire transform is wasteful. Goertzel evaluates a single DFT bin using a second-order recurrence: two multiply-adds per sample, two state variables, no complex arithmetic in the loop. It is a digital tuned circuit. You feed it samples; it accumulates energy only if the input contains its target frequency.

Here it is in Rust, fed 1024 samples of a synthetic 45 kHz tone at a 384 kHz sample rate:

fn goertzel(samples: &[f64], target: f64, sr: f64) -> f64 {
    let n = samples.len() as f64;
    let k = (n * target / sr).round();
    let coeff = 2.0 * (2.0 * std::f64::consts::PI * k / n).cos();
    let (mut s1, mut s2) = (0.0, 0.0);
    for &x in samples {
        let s0 = x + coeff * s1 - s2;
        s2 = s1;
        s1 = s0;
    }
    (s1 * s1 + s2 * s2 - coeff * s1 * s2).sqrt()
}

25000 Hz ->      0.0
45000 Hz ->    512.0
55000 Hz ->      0.0

The resonator lights up at 45 kHz and ignores its neighbours. The same recurrence in Haskell collapses the whole state machine into a single fold — the running pair (s1, s2) is the accumulator, and the loop body becomes the step function:

import Data.List (foldl')

goertzel :: Double -> Double -> [Double] -> Double
goertzel sr target xs = sqrt (s1*s1 + s2*s2 - coeff*s1*s2)
  where
    n        = fromIntegral (length xs)
    k        = fromIntegral (round (n * target / sr) :: Int)
    coeff    = 2 * cos (2 * pi * k / n)
    (s1, s2) = foldl' (\(p1, p2) x -> (x + coeff*p1 - p2, p1)) (0, 0) xs

Goertzel had its heyday in the 1980s inside DTMF decoders — the chips that turned your telephone keypad’s dual tones into digits. Eight target frequencies, eight tiny resonators running in parallel on hardware that couldn’t dream of a real-time FFT. That’s the constraint it was born for: many known frequencies, cheap silicon, no cycles to spare. My bat detector is the same problem inverted — instead of decoding a keypad, I’d point one resonator at 45 kHz, another at 55, a third at 82, and let whichever one rings tell me which bat just flew over. The ceramic transducer picks one bin in the analogue domain because physics gave it a resonance. Goertzel picks one bin in the digital domain because I told it a target and it kept only two numbers to find it.