1
$\begingroup$

I've implemented Goertzel algorithm, now very simply in Python. But I cannot obtain the correct answer it's supposed to produce for a single DTF frequency. here is my code in Python:

import numpy as np

# omega frequency: radians/sample
omega = 0.23

# length of the data to generate
length = 200

# updated parameters
dft = 0.0
s1 = 0.0
s2 = 0.0

# initial phase
radians = 0.0
for i in range(length):

    # generate signal value
    value = np.cos(radians)

    # goertzel algorithm update
    s0 = value + 2.0 * np.cos(omega) * s1  - s2
    s2 = s1
    s1 = s0

    # dft sum update
    dft += value * np.exp(-1j * omega * i)

    # next value
    radians += omega

# final goertzel algorithm update
s0 = 2.0 * s1 * np.cos(omega) - s2
goe = s0 - np.exp(-1j * omega) * s1

# print results
print "DFT: {:.4f}".format(dft / length)
print "GOE: {:.4f}".format(goe / length)

Output:

DFT: 0.4979-0.0097j
GOE: -0.2065+0.4531j

The correct answer (with infinite data) is 0.5 + 0.0j.

DFT is therefore providing the correct result.

What's wrong with my Goertzel algorithm implementation?

https://en.wikipedia.org/wiki/Goertzel_algorithm

$\endgroup$
  • $\begingroup$ Here are two implementations for MatLab - mathworks.com/matlabcentral/fileexchange/… $\endgroup$ – Juha P Aug 20 '18 at 12:58
  • $\begingroup$ Omg! All this time spent puzzling over this, only now I realise I’m trying to use Goertzel algorithm for non integer frequency bin! Please post your comment as answer and I will accept it! $\endgroup$ – akuz Aug 20 '18 at 13:02
  • $\begingroup$ Interestingly, in that article they propose to apply correction np.exp (-1j * omega * (length - 1)), whereas the correct one seems to be np.exp (-1j * omega * length) as proposed in Peter K's answer. $\endgroup$ – akuz Aug 20 '18 at 13:36
1
$\begingroup$

Why do you think the two should be equal?

Check out the final line of equation (6) in the article you link to:

Wikipedia equation

I believe your misunderstanding is to not take account of the $e^{j\omega_0n}$ term.

$\endgroup$
  • $\begingroup$ Indeed, if I multiply my Goertzel algorithm output by np.exp( -1j * omega * length ), I get exactly the same output as the DFT calculation. Equation (6) you're quoting has to do with unrolling the frequency back into the time domain, from which it is obvious that the correction should be applied. But why then the section "DFT computations" doesn't mention this correction, so that the output of Goertzel algorithm produces the correct phase? Is it because they only consider integer k? Please provide a little more detail, and I will accept your answer. $\endgroup$ – akuz Aug 20 '18 at 13:33
  • $\begingroup$ Sorry, my misunderstanding, y[n] is still the frequency! I guess that term is not present in section "DFT computations" because it equals 1? $\endgroup$ – akuz Aug 20 '18 at 13:44
  • $\begingroup$ Answering my own last question :) The correction term equals to 1 for frequency bins computed by DFT, but not for any omega (in which case, the additional term specified by Peter K should be applied). Thank you!!! $\endgroup$ – akuz Aug 20 '18 at 14:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.