First Order Goertzel Algorithm
Based on the equation of the Fourier transformation:
The first order Goertzel algorithm is based on the convolution of x[n] and an additional signal h[n] and ends, after a hell of a complicated explanation, in the simple formula:
with:
and:
I first found the implementation of this into a Fourier transformation in a very interesting book about Matlab. It basically looked like:
a1 = -exp(1i*2*pi*k/N);
y = x(1);
for n = 2:N
y = x(n) - a1*y;
end
y = -a1*y;
OK, we don’t analyse this syntax now. The same algorithm in C looks like this:
for (k = 0; k < N; k++)
{
c[k].real = y[0].real;
c[k].imag = y[0].imag;
w.real = -Math.Cos((double)(2.0 * Math.PI * (double)(k) / (double)(N)));
w.imag = Math.Sin((double)(2.0 * Math.PI * (double)(k) / (double)(N)));
for (n = 1; n <= N; n++)
c[k] = kdiff(y[n], kprod(c[k], w));
c[k] = kprod(c[k], w);
c[k].real = -c[k].real / (double)(N) * 2.0;
c[k].imag = -c[k].imag / (double)(N) * 2.0;
}
With the complex operations:
public TKomplex kdiff(TKomplex a, TKomplex b)
{
TKomplex res;
res.real = a.real - b.real;
res.imag = a.imag - b.imag;
return (res);
}
public TKomplex kprod(TKomplex a, TKomplex b)
{
TKomplex res;
res.real = a.real * b.real - a.imag * b.imag;
res.imag = a.real * b.imag + a.imag * b.real;
return (res);
}
A First Step to the Approach
Somehow, I had the feeling that I have seen something like
already somewhere else… There is this thing about the polynomial
y = a1*x + a2*x^2 + a3*x^3 + a4*x^4…
and its calculation according to Horners role for polynomial evaluation.
y = x * (a1 + x*(a2 + x*(a3 + x*(a4 …))))
If I implement this in a small loop, it almost looks like the formula above already. But a Fourier transformation does not look like this polynomial. The formula of the Fourier transformation is:
OK. But there is a way to bring this into a polynomial. The expression
is a complex vector of the length 1 that spins around the unit circle. In the complex calculation, we can say
And the Fourier transformation can be written like:
And that can be written as a polynomial.
If I put this into a loop, it starts like:
and goes on like:
and:
and:
Till we are at f(0):
And theoretically:
to finish this polynomial.
In the special case of the Fourier components, this last term can be skipped as:
Implemented in C, this looks like:
for (k = 0; k < N; k++)
{
c[k].real = y[N].real;
c[k].imag = y[N].imag;
w.real = Math.Cos((double)(2.0 * Math.PI * (double)(k) / (double)(N)));
w.imag = -Math.Sin((double)(2.0 * Math.PI * (double)(k) / (double)(N)));
for (n = N - 1; n > 0; n--)
c[k] = ksum(y[n], kprod(c[k], w));
c[k].real = c[k].real / (double)(N) * 2.0;
c[k].imag = -c[k].imag / (double)(N) * 2.0;
}
Getting Closer
That looks quite similar to the Goerzel algorithm already. But there are three points that differ:
- w is of opposite sign.
- in the inner loop, the indexing n starts at N-1 instead of 1. That means we work from the backside.
- In the loop, we have an addition instead of a subtraction.
Now the big question is: What does that mean?
If I run into the other direction in my loop, that’s no difference for the addressing of the samples y[n]. They are still the same. But w is a complex vector. If we start at 0 with n and increase it till N, w will spin counter clockwise. If we run into the other direction starting at N-1 and decreasing n till 1, w will spin clockwise. and that has some impact. Additionally, we do not start at the same value of w. We start at:
Instead of:
This is basically the conjugate-complex value of w. Fortunately with this starting value, w will spin clockwise automatically as it has a negative imaginary part.
And finally: As we run into the other direction, the loop does not end at e0 and we cannot just skip the last term. With this modification, the algorithm becomes like this:
for (k = 0; k < N; k++)
{
c[k].real = y[0].real;
c[k].imag = y[0].imag;
w.real = Math.Cos((double)(2.0 * Math.PI * (double)(k) / (double)(N)));
w.imag = Math.Sin((double)(2.0 * Math.PI * (double)(k) / (double)(N)));
for (n = 1; n <= N; n++)
c[k] = ksum(y[n], kprod(c[k], w));
c[k] = kprod(c[k], w);
c[k].real = c[k].real / (double)(N) * 2.0;
c[k].imag = -c[k].imag / (double)(N) * 2.0;
}
And Closer
But still: Goertzel algorithm does not work with the conjugate-complex of e-jx. It uses the negative value of it and subtracts this instead of adding it. Because of this, finally the real part of the Fourier components becomes negative and has to be negated at the end of the loop. Now the algorithm looks like this:
for (k = 0; k < N; k++)
{
c[k].real = y[0].real;
c[k].imag = y[0].imag;
w.real = -Math.Cos((double)(2.0 * Math.PI * (double)(k) / (double)(N)));
w.imag = Math.Sin((double)(2.0 * Math.PI * (double)(k) / (double)(N)));
for (n = 1; n <= N; n++)
c[k] = kdiff(y[n], kprod(c[k], w));
c[k] = kprod(c[k], w);
c[k].real = -c[k].real / (double)(N) * 2.0;
c[k].imag = -c[k].imag / (double)(N) * 2.0;
}
And this is the Goertzel algorithm. Derived from a polynomial calculation with some small modification, we got the same algorithm.
Advantage and Disadvantage of the Goerzel Algorithm
The main goal of the Goertzel algorithm is the calculation speed. It has to calculate much less sine and cosine values than the standard implementation of the DFT and therefore it works quite a bit faster. But in my article about the quick and lean DFT, I found a way to speed up the DFT algorithm even a bit more.
And I see a possible disadvantage: If we have a big amount of samples N will be big as well and due to that w will have a very small angle that has to be multiplied many, many times. A very small imaginary part has to be processed with a, compared to the imaginary part, big real part many times. That can cause a loss of accuracy of this calculation.
Conclusion
This approach shows quite an easy to be understood way to the first order Goertzel algorithm and shows the algorithm is more or less a polynomial evaluation by Horner running into the other direction. :-)
But of course: the second order Goertzel is quite another story.
