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