## 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.