Click here to Skip to main content
11,641,222 members (61,550 online)
Click here to Skip to main content

Tagged as

Another Approach to the First Order Goertzel Algorithm

, 20 Oct 2013 CPOL 14.4K 268 19
Rate this:
Please Sign up or sign in to vote.
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&hellip; 

and its calculation according to Horners role for polynomial evaluation.

y = x * (a1 + x*(a2 + x*(a3 + x*(a4 &hellip;))))

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:

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. Smile | :)

But of course: the second order Goertzel is quite another story.

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)

Share

About the Author

Mosi_62
Software Developer (Senior) Aebi Maschinenfabrik, Burgdorf
Switzerland 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 Smile | :)

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 farming Tractors for business, I enjoy very much to implement interesting algorithms and analyse the mathematics they are based on in my leisure time Smile | :)

You may also be interested in...

Comments and Discussions

 
GeneralMy vote of 5 Pin
den2k8820-Nov-14 2:39
memberden2k8820-Nov-14 2:39 
QuestionVote of 5 Pin
Kenneth Haugland13-Sep-14 4:16
professionalKenneth Haugland13-Sep-14 4:16 
GeneralMy vote of 5 Pin
phil.o6-Jun-14 10:23
professionalphil.o6-Jun-14 10:23 
GeneralRe: My vote of 5 Pin
Mosi_6210-Jun-14 7:19
professionalMosi_6210-Jun-14 7:19 
QuestionGood, maybe but... Pin
deebee++14-Oct-13 21:11
memberdeebee++14-Oct-13 21:11 
AnswerRe: Good, maybe but... Pin
Mosi_6220-Oct-13 4:50
memberMosi_6220-Oct-13 4:50 
GeneralRe: Good, maybe but... Pin
H.Brydon11-Nov-13 3:26
memberH.Brydon11-Nov-13 3:26 
GeneralRe: Good, maybe but... Pin
Nelek15-Sep-14 12:18
memberNelek15-Sep-14 12:18 
AnswerRe: Good, maybe but... Pin
BigTimber@home20-Oct-13 14:28
memberBigTimber@home20-Oct-13 14:28 
QuestionGoertzel? Pin
roscler14-Oct-13 11:00
professionalroscler14-Oct-13 11:00 
AnswerRe: Goertzel? Pin
Mosi_6220-Oct-13 4:30
memberMosi_6220-Oct-13 4:30 
QuestionSome parts can be done faster Pin
JohnWallis427-Oct-13 13:19
memberJohnWallis427-Oct-13 13:19 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Rant Rant    Admin Admin   

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

| Advertise | Privacy | Terms of Use | Mobile
Web02 | 2.8.150731.1 | Last Updated 20 Oct 2013
Article Copyright 2013 by Mosi_62
Everything else Copyright © CodeProject, 1999-2015
Layout: fixed | fluid