Goertzel algorithm

Numerical stability[edit]

It can be observed that the poles of the filter's Z transform are located at ${\displaystyle e^{+j\omega _{0}}}$ and ${\displaystyle e^{-j\omega _{0}}}$, on a circle of unit radius centered on the origin of the complex Z-transform plane. This property indicates that the filter process is marginally stable and vulnerable to numerical-error accumulation when computed using low-precision arithmetic and long input sequences.^[6] A numerically stable version was proposed by Christian Reinsch.^[7]

For the important case of computing a DFT term, the following special restrictions are applied.


Making these substitutions into equation (6) and observing that the term ${\displaystyle e^{+j2\pi k}=1}$, equation (6) then takes the following form:


We can observe that the right side of equation (9) is extremely similar to the defining formula for DFT term ${\displaystyle X[k]}$, the DFT term for index number ${\displaystyle k}$, but not exactly the same. The summation shown in equation (9) requires ${\displaystyle N+1}$ input terms, but only ${\displaystyle N}$ input terms are available when evaluating a DFT. A simple but inelegant expedient is to extend the input sequence ${\displaystyle x[n]}$ with one more artificial value ${\displaystyle x[N]=0}$.^[8] We can see from equation (9) that the mathematical effect on the final result is the same as removing term ${\displaystyle x[N]}$ from the summation, thus delivering the intended DFT value.


However, there is a more elegant approach that avoids the extra filter pass. From equation (1), we can note that when the extended input term ${\displaystyle x[N]=0}$ is used in the final step,


Thus, the algorithm can be completed as follows:


The last two mathematical operations are simplified by combining them algebraically:


Note that stopping the filter updates at term ${\displaystyle N-1}$ and immediately applying equation (2) rather than equation (11) misses the final filter state updates, yielding a result with incorrect phase.^[9]


The particular filtering structure chosen for the Goertzel algorithm is the key to its efficient DFT calculations. We can observe that only one output value ${\displaystyle y[N]}$ is used for calculating the DFT, so calculations for all the other output terms are omitted. Since the FIR filter is not calculated, the IIR stage calculations ${\displaystyle s[0],s[1]}$, etc. can be discarded immediately after updating the first stage's internal state.


This seems to leave a paradox: to complete the algorithm, the FIR filter stage must be evaluated once using the final two outputs from the IIR filter stage, while for computational efficiency the IIR filter iteration discards its output values. This is where the properties of the direct-form filter structure are applied. The two internal state variables of the IIR filter provide the last two values of the IIR filter output, which are the terms required to evaluate the FIR filter stage.