# Goertzel algorithm produces incorrect phase

I've implemented Goertzel algorithm according to the Wikipedia page (https://en.wikipedia.org/wiki/Goertzel_algorithm) and another page (http://www.mstarlabs.com/dsp/goertzel/goertzel.html), which are consistent with each other. I am then testing it using a naive implementation for a single DFT frequency (https://en.wikipedia.org/wiki/Discrete_Fourier_transform). I am getting correct amplitude, but incorrect phase. Please help me find the mistake in my code, I've spent a couple of days on it now, and can't see the problem...

/**
* Complex imaginary number i.
*/
template<typename T>
constexpr auto COMPLEX_i = std::complex<T> (
static_cast<T> (0.0), static_cast<T> (1.0));

/**
* Convert frequency in Hertz to Omega (radian per sample).
*/
template<typename T>
T
hertz_to_omega (T hertz, T sample_rate)
{
return hertz / sample_rate * static_cast<T> (2.0) * M_PI;
}

/**
* Simple implementation of the Goertzel algorithm
* for detecting a given omega frequency from a range
* of raw signal values.
*/
template<typename T, typename Range>
std::complex<T>
goertzel (const Range &range, T omega)
{
const auto coeff = static_cast<T> (2.0) * std::cos (omega);

T s1 = 0.0;
T s2 = 0.0;
for (const auto value : range)
{
const auto s0 = coeff * s1 - s2 + value;
s2 = s1;
s1 = s0;
}

return std::complex<T> (
std::cos (omega) * s1 - s2,
std::sin (omega) * s1) /
static_cast<T> (range.size ());
}

BOOST_AUTO_TEST_CASE( test_goertzel )
{
using T = double;

const T FREQ = 543.2;
const T SAMPLE_RATE = 5678.9;
const T OMEGA = hertz_to_omega (FREQ, SAMPLE_RATE);
const T START_RADIANS = M_PI / 3.0;
const size_t LENGTH = 20000;

std::complex<T> dft_sum (static_cast <T> (0.0), static_cast<T> (0.0));
std::vector<T> data;
{
for (size_t sample_index = 0; sample_index < LENGTH; ++sample_index)
{
const T value = std::cos (radians);
dft_sum += value * std::exp (- COMPLEX_i<T>
* OMEGA * static_cast<T> (sample_index));
data.push_back (value);
}
}

const auto result1 = dft_sum / static_cast<T> (LENGTH);
const auto result2 = goertzel<T> (data, OMEGA);

std::cout << "DFT: " << result1 << ", ABS: " << std::abs (result1) << std::endl;
std::cout << "GOE: " << result2 << ", ABS: " << std::abs (result2) << std::endl;
}


Running the test above produces the following output:

Running 1 test case...
DFT: (0.250009,0.433004), ABS: 0.499997
GOE: (0.114572,0.486693), ABS: 0.499997


The magnitude is detected correctly (it should be 0.5), but the phases of the naive DFT implementation and the Goertzel algorithm output don't match. Please help me find the mistake in either my single-bin DFT or Goertzel implementation.

EDIT:

It seems that DFT calculation is detecting the phase correctly, as atan(0.433/0.25)/180 = 0.333 (and I specified the phase of pi / 3 for the generated signal on which I am testing), so its my implementation of Goertzel algorithm which is incorrect.

The short answer is that when Goertzel algorithm is used for omega, which doesn't correspond to of the finite set of frequencies computed by DFT (which correspond to integer k), then an additional correction factor should be applied. I've got an answer from Peter K after describing my question more clearly here: What's wrong with my Goertzel algorithm implementation?

If I understood your code correctly, you need a final multiplication for s2, that is s2*exp(-j2*pi*k/N) before returning it.

• According to formula (11) here en.wikipedia.org/wiki/Goertzel_algorithm, the final result is computed as exp(j 2 pi k / N) * s1 - s2. In my code, OMEGA = 2 pi k / N. So it's equivalent to exp(j omega) * s1 - s2 which, when represented as real and complex parts, is equivalent to real: cos(omega) - s2, imag: sin(omega) * s1, which is what I am returning. How did you reach what you offer in your suggestion?
– akuz
Aug 17 '18 at 9:16
• I am just looking at the code I've been using and that I have extensively tested. I am not expert of c++, might be that I have overlooked something. This is the code I have used for matlab bitbucket.org/signalteam/goertzel/src/… Aug 17 '18 at 9:19
• I see you are using formula (10) (from Wikipedia) by pushing an additional value of 0.0, and then using formula (2) to compute the final result -- this is also correct. I am using formula (11) which combines (2) and (10) for the last step, without explicitly using an additional value of 0.0. The result theoretically should be the same. But I will try your approach a little bit later, and let you know how it goes (if it works, I will accept your answer). Thank you! (Perhaps it's my naive implementation of DFT bin is incorrect?)
– akuz
Aug 17 '18 at 9:45
• It seems that DFT is detecting the phase correctly, as atan(0.433/0.25)/180 = 0.333 (and I specified the phase of pi / 3 for the generated signal on which I am testing), so its my implementation of Goertzel algorithm which is incorrect.
– akuz
Aug 17 '18 at 9:58
• I've replaced the return value of my goerstzel() function to match your code: const auto s0 = coeff * s1 - s2; return (s0 - s1 * std::exp (- COMPLEX_i<T> * omega)) / static_cast<T> (range.size ()); -- and I get exactly the same answer I was getting before. Weird.
– akuz
Aug 17 '18 at 11:58