Data Acquisition (DAQ) and Control from Microstar Laboratories

Knowledge Base: Processing

  • See the list of other topics for category Processing

Q10112 "Non-standard" termination equations for Goertzel Algorithm

Tags: Help, processing, custom command modules, DSP, Goertzel Algorithm, mathematics

Applies to: custom command programming, DFT calculations

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?

  1. 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.
  2. 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.
  3. 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-----

References:

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