12,881,824 members (29,027 online)
alternative version

#### Stats

19.1K views
21 bookmarked
Posted 7 Oct 2013

# Another Approach to the First Order Goertzel Algorithm

, 20 Oct 2013 CPOL
 Rate this:
This approach shows quite an easy to be understood way to the first order Goerzel algorithm.

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

1. w is of opposite sign.
2. in the inner loop, the indexing n starts at N-1 instead of 1. That means we work from the backside.
3. 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:

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.

## About the Author

 Software Developer (Senior) Güdel AG, Langenthal Switzerland
Computers are very straight... They always do exactly what we tell them to do... Only, much too often what we tell them to do is not really what we want them to do

Writing Software is one of the most creative tings one can do. I have been doing this for more than ten years now and still having a lot of fun with it. Besides doing software for HMI's on C# for business, I enjoy very much to implement interesting algorithms and analyse the mathematics they are based on in my leisure time

For more detailed descriptions and math visit me on my own page

www.mosismath.com

 Pro

## Comments and Discussions

 First Prev Next
 Comment FatCatProgrammer14-Dec-15 4:05 FatCatProgrammer 14-Dec-15 4:05
 My vote of 5 den2k8820-Nov-14 2:39 den2k88 20-Nov-14 2:39
 Vote of 5 Kenneth Haugland13-Sep-14 4:16 Kenneth Haugland 13-Sep-14 4:16
 My vote of 5 phil.o6-Jun-14 10:23 phil.o 6-Jun-14 10:23
 Re: My vote of 5 Mosi_6210-Jun-14 7:19 Mosi_62 10-Jun-14 7:19
 Good, maybe but... deebee++14-Oct-13 21:11 deebee++ 14-Oct-13 21:11
 Re: Good, maybe but... Mosi_6220-Oct-13 4:50 Mosi_62 20-Oct-13 4:50
 Re: Good, maybe but... H.Brydon11-Nov-13 3:26 H.Brydon 11-Nov-13 3:26
 Re: Good, maybe but... Nelek15-Sep-14 12:18 Nelek 15-Sep-14 12:18
 Re: Good, maybe but... BigTimber@home20-Oct-13 14:28 BigTimber@home 20-Oct-13 14:28
 Goertzel? roscler14-Oct-13 11:00 roscler 14-Oct-13 11:00
 Re: Goertzel? Mosi_6220-Oct-13 4:30 Mosi_62 20-Oct-13 4:30
 Some parts can be done faster JohnWallis427-Oct-13 13:19 JohnWallis42 7-Oct-13 13:19
 Last Visit: 31-Dec-99 18:00     Last Update: 23-Apr-17 6:50 Refresh 1

General    News    Suggestion    Question    Bug    Answer    Joke    Praise    Rant    Admin

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.