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