# What's wrong with my Goertzel algorithm implementation?

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 = 0.23

# length of the data to generate
length = 200

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

# initial phase
for i in range(length):

# generate signal value

# 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

# 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

• Here are two implementations for MatLab - mathworks.com/matlabcentral/fileexchange/… Aug 20 '18 at 12:58
• 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!
– akuz
Aug 20 '18 at 13:02
• 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.
– akuz
Aug 20 '18 at 13:36

Why do you think the two should be equal?

Check out the final line of equation (6) in the article you link to: I believe your misunderstanding is to not take account of the $e^{j\omega_0n}$ term.

• 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.
– akuz
Aug 20 '18 at 13:33
• Sorry, my misunderstanding, y[n] is still the frequency! I guess that term is not present in section "DFT computations" because it equals 1?
– akuz
Aug 20 '18 at 13:44
• 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!!!
– akuz
Aug 20 '18 at 14:04