Goertzel Filtering in Python

The Goertzel Filter Equation

There’s a lot of ways to write the Goertzel filter equation. I’m going to work with the one that looks like a DFT term, since it’s the simplest in some ways and makes the analogy with the DFT clear. This is Equation (6) in the Wikipedia article.

\[ y[n] = e^{j w_0 n} \sum_{k=0}^n x[k] e^{-j w_0 k} \]

Some things to note here:

• This is an equation in the complex numbers. \(j\) is \(\sqrt{-1}\) (often written \(i\)).

• \(x[k]\) is the \(k\)-th input sample in the filter window — a real number. The filter window contains previous samples, numbered \(0 \ldots n\) inclusive (so the counts get a little weird) in reverse order. In particular, \(k=0\) corresponds to the current sample. Of course, the reverse of a sine wave is just the same sinewave except for phase, so for this application we don’t really care about forwards vs backwards.

• \(y[n]\) is the output filter value at the current time point: a complex number. The \(n\) is somewhat redundant, but it reflects a common DSP convention.

• \(w_0\) is the target frequency of the filter (the frequency you are trying to measure) in radians per second. To get this frequency, you must scale your desired frequency \(f\) in Hz (cycles per second) in two ways:

  • Multiply by \(2 \pi\) to convert cycles to radians.

  • Divide by the sample rate to scale the target frequency to the range \(0 \ldots 1\).

  If your target frequency is 2025 Hz and your sample rate is 48000 samples per second, for example, you find that

  \[\begin{eqnarray*} w_0 &=& \frac{2 \pi~2025}{48000} \\ &=& 0.265\ldots \end{eqnarray*}\]

• In our application you don’t care about the phase of the output, just the magnitude — a real number. So you’re really looking for

  \[ |y[n]| = |e^{j w_0 n} \sum_{k=0}^n x[k] e^{-j w_0 k}| \]

  where

  \[ |a + jb| = \sqrt{a^2 + b^2} \]

  [If you were feeling clever (I’m not) you could go back through the math and use Euler’s Formula

  \[ e^{jx} = cos(x) + j\sin(x) \]

  and some trigonometry and get rid of all the complex numbers.]

Implementing The Filter

Given the above, it’s pretty straightforward to implement the filter equation. You could write a Python function

    def goertzel(x, f):
    """x is an array of samples, f is the target frequency.
    Returns the output magnitude."""

To implement this function, you would like to be able to do complex arithmetic. In Python, the numpy package provides this: basic math on complex numbers works, numpy.exp() does complex exponentials, numpy.abs() computes the magnitude of a complex number.

If your language had no complex arithmetic support, you could rewrite the filter equation using Euler’s Formula and then do it all as real arithmetic.

Saving The Filter Coefficients

If you look at the filter equation

\[ |y[n]| = |e^{j w_0 n}n \sum_{k=0}^n x[k] e^{-j w_0 k}| \]

the sum term looks like a dot product. For a given filter size, the normalizing coefficient \(e^{j w_0 n}\) and filter coefficients \(e^{-j w_0 k}\) can be precomputed and reused.

In Python with numpy this would look something like

    norm = numpy.exp(1j * w0 * n)
    coeffs = numpy.array([-1j * w0 * k for k in range(n)])

and then

    y = numpy.abs(norm * numpy.dot(coeffs, xs))

This is fast and easy. I defined a Goertzel class that saved the norm and coefficients, then had a filter() method to apply the filter.