Q10112 "Non-standard" termination equations for Goertzel Algorithm

The equations shown on the Microstar Web site for finishing a Goertzel Filter calculation of a DFT term are different from the "classic" equations shown by other references. Who is correct?

As it turns out: both the "classic" equations and the modified termination
equations shown in the Web page article [1] are
correct — *depending on how you implement the algorithm!*

To show just how tricky this is, let us first describe the obvious
way to implement the filtering loop and classical termination equations
for the Goertzel Filter operation. A DFT calculation would start with a
zero accumulator, so the filter is started with zero values in its two
state variables. There are `M`

data terms to process, so the
filter loop can be executed `M`

times. The termination equation
is applied to finish.

The following listing shows the calculation using a Matlab/Octave
script notation [2]. Verify that this is
completely consistent with the description just given. (Keep in mind
that *index 1* in this script language corresponds to
*subscript 0* in a pure mathematical notation.)

% Configuration parameters k = 3 % desired DFT term M = 10 % number of terms in DFT data block W = exp(-j*2*pi*k/M); % Data to transform x = [ define M term data vector here ] % Initial state of filtering loop d1 = 0.0; d2 = 0.0; % Classic Goertzel filter calculations with M passes for n=1:M y(n) = x(n)+2*real(W)*d1-d2; d2 = d1; d1 = y(n); end goertx = d1 - W*d2

Unfortunately, the final result is **wrong!** You can verify
that the result is wrong numerically by comparing to a normal DFT
calculation.

You can prove that this implementation is flawed, by inspection of
the algorithm code and the expanded form of the Goertzel filter
equations (see reference 1 for example).
Observe that in the final pass of the algorithm, the last term
`x(M)`

has never been multiplied by any `W`

term.
If you compare to the fully expanded Goertzel filter calculations, you
can see that this last term should be multiplied by the factor
`W`

exactly one time. This eventually leads to the explanation
for what has happened: a final multiplication by factor `W`

is missing.

There is more than one way to correct this problem. An obvious way
is to muliply by the `W`

factor after the filtering sequence
has terminated. This is perfectly valid, but it requires an extra
complex-valued multiply operation. This undermines the intent of the
Goertzel algorithm to avoid as many of these complex-valued multiplications
as possible.

A less obvious way is to place an extra `0.0`

"padding value" at the end of the input data block, resulting in a data
sequence of length `M+1`

instead of `M`

. Then,
execute the filter `M+1`

passes instead of the obvious M
passes. *This modified algorithm produces a correct result*
— and it uses the "classical" termination equation. While this
works fine, the peculiar padded input sequence and extra filtering
pass make this seem like an inelegant solution.

However, there is a third approach. The valid `M+1`

filter
pass solution can be considered an `M`

-pass solution followed
by one special, artificial pass that uses the padding value 0 as input.
What happens during that added pass?

- Since the
`x`

padding term is zero, it has no effect on the filter output`y`

. Only the state terms d1 and d2 are involved during the extra pass. - There is an algebraic combination of the d1 and d2 terms to produce the final output value, where states d1 and d2 carry forward from the preceding pass.
- We can skip the updates of the d1 and d2 states following the final
`y`

calculation, since those values will never be used.

The separate operation avoids both the extra complex-valued multiply and the overhead of the extra algorithm pass. However, it remains a necessary calculation that must be performed before applying the "classical" Goertzel termination equations.

The supplementary calculation is just an algebraic combination of
the filter `d1`

and `d2`

states. But we can
observe that the "classic" Goertzel filter termination is a similar
algebraic combination of the filter `d1`

and `d2`

states. Thus, these two algebraic steps can be combined analytically.
What results is the modified termination equations as shown in the Web
article [1]. For reference, here is a correct
algorithm using `M`

filter passes and using the modified
termination conditions.

% Modified Goertzel filter calculations with M passes (correct!) for n=1:M y(n) = x(n)+2*real(W)*d1-d2; d2 = d1; d1 = y(n); end goert2 = conj(W)*d1 - d2

L-----

[1] View the article discussing applications of the Goertzel algorithm, in combinations with windowing techniques, at: Detecting A Single Frequency Efficiently.

[2] A Matlab/Octave script file can be run to compare the results of the basic DFT, classic Goertzel, and modified Goertzel implementations.