Goertzel Algorithm

Power Calculations

I have decided to drop power calculations because they cause all sorts of problems.

Just let power=amplitude*amplitude; if you really want to use it.

In the end a power calculation needs a load impedance to actually make any sense anyway.

Continuous Mode

What would be really nice is a form of Goertzel's algorithm that is continuous like the exponential smoothing function. Because of the windowing I don't think it is possible?

AlanX

Goertzel's Algorithm

Here is the core code that works as expected:

  // Generate samples
  signal=0.0;
  for (i=0;i<N;i++) {
    // samples[i]=20*((1-0.19)*exp(-i/75.0)+0.19*exp(-i/675.0));
    signal=1.0*(i%125==0?1.0-signal:signal);
    // signal=sin(2*M_PI*(i%125)/125); // 1kHz
    samples[i]=signal;
    // Window the data
    // data[i]=samples[i]*(0.5-0.35*cos(2*M_PI*i/N));  // Hanning Window
    // amplitudeFactor=2.0;
    // powerFactor=1.63;
    data[i]=samples[i]*(0.54-0.46*cos(2*M_PI*i/N)); // Hamming Window
    amplitudeFactor=1.85;
    powerFactor=1.59;
  }

  // Scan frequencies
  for (freq=fmin;freq<=fmax;freq+=1) {
    w=M_PI*2*freq/SampleFreq;
    c=2*cos(w);
    s0=0.0;
    s1=0.0;
    s2=0.0;
    for (i=0;i<N;i++) {
      // Goertzel
      s0=data[i]+c*s1-s2;
      s2=s1;
      s1=s0;
    }
    amplitude=amplitudeFactor*sqrt(s2*s2+s1*s1-c*s1*s2)/N;
    ...
}

The need for windowing

Here is Goertzel's algorithm without windowing:

While the peaks are correct what is that background noise? They call it spectral leakage. Basically Goertzel's algorithm see an impulse/sudden start-up/stop, not a continuous 1v/0v square wave. The purpose of windowing is the get rid of the startup noise.

Here is the Hanning window:

A big improvement but not perfect. Although not the best window to suppress impulse noise, it is the easiest to adjust to preserve the fourier transform magnitudes.

What? Basically the default Hanning window (0.5-0.25*cos(2*M_PI*i/N)) reduces the area of the input signal by 50% resulting in a 50% reduction in the fourier transform magnitude. It is easy to see that multiplying by 2 will fix this problem.

You could work out the adjustment factors for the other windows but not by inspection.

Update

Finally found a reference for this (thanks to Peter Schaldenbrand of Siemens):

https://community.sw.siemens.com/s/article/window-correction-factors

Here it says you can correct the amplitude or the energy (i.e. power) but not both at the same time.

Here is the Hamming window:

The other correction factors don't match my window functions.

The K Term Problem

  k/N=Fc/Fs

where:

  Fc  = Test (centre) Frequency

  Fs = Sample Frequency

  N = Number of samples

The problem is that many mathematical texts say that k must be an integer.

In practice this does not appear to be necessary. I use double Fc/Fs rather the integer k/N in my code.

Finally the Power Normalisation problem

It is surprising that so few program codes I have looked did not normalise the power/magnitude (i.e. divide by N^2):

power=(s2*s2+s1*s1-c*s1*s2)/N/N;
or
amplitude=sqrt(s2*s2+s1*s1-c*s1*s2)/N;

AlanX