# Using Goertzel Algorithm in under-sampling

I plan to calculate a signal's phase using Goertzel Algorithm. I have 2 signals coming to microcontroller's ADC. Need to measure the phase difference between them. Signals are 15MHz sinusoids. Sample rate will be about 1MHz or lower. Is it possible to use Goertzel with under-sampling? It would be much easier than calculating FFT of both signals, then search for highest beams and count phase of them. Does anyone tried it? Here is nice article about Goertzel:

Banks K. - The Goertzel Algorithm

And here is the code from 'listing1':

Listing 1: a Goertzel implementation

#include <stdio.h>
#include <math.h>

#define FLOATING    float
#define SAMPLE  unsigned char

#define SAMPLING_RATE   8000.0  //8kHz
#define TARGET_FREQUENCY    941.0   //941 Hz
#define N   205 //Block size

FLOATING coeff;
FLOATING Q1;
FLOATING Q2;
FLOATING sine;
FLOATING cosine;

SAMPLE testData[N];

/* Call this routine before every "block" (size=N) of samples. */
void ResetGoertzel(void)
{
Q2 = 0;
Q1 = 0;
}

/* Call this once, to precompute the constants. */
void InitGoertzel(void)
{
int   k;
FLOATING  floatN;
FLOATING  omega;

floatN = (FLOATING) N;
k = (int) (0.5 + ((floatN * TARGET_FREQUENCY) / SAMPLING_RATE));
omega = (2.0 * PI * k) / floatN;
sine = sin(omega);
cosine = cos(omega);
coeff = 2.0 * cosine;

printf("For SAMPLING_RATE = %f", SAMPLING_RATE);
printf(" N = %d", N);
printf(" and FREQUENCY = %f,\n", TARGET_FREQUENCY);
printf("k = %d and coeff = %f\n\n", k, coeff);

ResetGoertzel();
}

/* Call this routine for every sample. */
void ProcessSample(SAMPLE sample)
{
FLOATING Q0;
Q0 = coeff * Q1 - Q2 + (FLOATING) sample;
Q2 = Q1;
Q1 = Q0;
}

/* Basic Goertzel */
/* Call this routine after every block to get the complex result. */
void GetRealImag(FLOATING *realPart, FLOATING *imagPart)
{
*realPart = (Q1 - Q2 * cosine);
*imagPart = (Q2 * sine);
}

/* Optimized Goertzel */
/* Call this after every block to get the RELATIVE magnitude squared. */
FLOATING GetMagnitudeSquared(void)
{
FLOATING result;

result = Q1 * Q1 + Q2 * Q2 - Q1 * Q2 * coeff;
return result;
}

/*** End of Goertzel-specific code, the remainder is test code. */

/* Synthesize some test data at a given frequency. */
void Generate(FLOATING frequency)
{
int   index;
FLOATING  step;

step = frequency * ((2.0 * PI) / SAMPLING_RATE);

/* Generate the test data */
for (index = 0; index < N; index++)
{
testData[index] = (SAMPLE) (100.0 * sin(index * step) + 100.0);
}
}

/* Demo 1 */
void GenerateAndTest(FLOATING frequency)
{
int   index;

FLOATING  magnitudeSquared;
FLOATING  magnitude;
FLOATING  real;
FLOATING  imag;

printf("For test frequency %f:\n", frequency);
Generate(frequency);

/* Process the samples */
for (index = 0; index < N; index++)
{
ProcessSample(testData[index]);
}

/* Do the "basic Goertzel" processing. */
GetRealImag(&real, &imag);

printf("real = %f imag = %f\n", real, imag);

magnitudeSquared = real*real + imag*imag;
printf("Relative magnitude squared = %f\n", magnitudeSquared);
magnitude = sqrt(magnitudeSquared);
printf("Relative magnitude = %f\n", magnitude);

/* Do the "optimized Goertzel" processing */
magnitudeSquared = GetMagnitudeSquared();
printf("Relative magnitude squared = %f\n", magnitudeSquared);
magnitude = sqrt(magnitudeSquared);
printf("Relative magnitude = %f\n\n", magnitude);

ResetGoertzel();
}

/* Demo 2 */
void GenerateAndTest2(FLOATING frequency)
{
int   index;

FLOATING  magnitudeSquared;
FLOATING  magnitude;
FLOATING  real;
FLOATING  imag;

printf("Freq=%7.1f   ", frequency);
Generate(frequency);

/* Process the samples. */
for (index = 0; index < N; index++)
{
ProcessSample(testData[index]);
}

/* Do the "standard Goertzel" processing. */
GetRealImag(&real, &imag);

magnitudeSquared = real*real + imag*imag;
printf("rel mag^2=%16.5f   ", magnitudeSquared);
magnitude = sqrt(magnitudeSquared);
printf("rel mag=%12.5f\n", magnitude);

ResetGoertzel();
}

int main(void)
{
FLOATING freq;

InitGoertzel();

/* Demo 1 */
GenerateAndTest(TARGET_FREQUENCY - 250);
GenerateAndTest(TARGET_FREQUENCY);
GenerateAndTest(TARGET_FREQUENCY + 250);

/* Demo 2 */
for (freq = TARGET_FREQUENCY - 300; freq <= TARGET_FREQUENCY + 300; freq += 15)
{
GenerateAndTest2(freq);
}

return 0;
}


What changes I should do in the code to play well with under-sampling?

Thanks! Mariusz

Goertzel algorithm is derived from DFT definition of one frequency bin:

• $X(k) = \displaystyle\sum_{m=0}^{N-1} x(n)e^{-j\cdot2\cdot\pi\cdot nk/N}$,

and also it can be rewritten:

• $X(k) = \displaystyle\sum_{m=0}^{N-1} x(n)e^{-j\cdot2\cdot\pi\cdot -k(N-n)/N}$, since $W_{N}^{-kN} = 1$.

Frome the property above recurrent equation can be derived (there is no sense to show all the conclusion here, I suppose):

• $y_r(k) = W_{N}^{-k} \cdot (s(r) + y_{r-1}(k))$,

where $s(r)$ is the sample of the input signal. It is also prooved that $y_r(k) = X(k)$ for $r = N-1$. So after $N$ iterations you have exact DFT bin value $X(k)$.

From the other side FFT also calulates the exact DFT bins. Both of the algorithms exploit DFT equation and DFT properties in the origin to compute frequency bins, there is no approximation. So for both algorithms you will get the equivalent results (the same input signal is assumed, of course). The only differences may occure due to impementation of 2d order IIR filter which is in the core of Goertzel algorithm.

Hope this helps.

• Thank you Serj. Now I know that I may use the Goertzel the same way as normal FFT when undersampling. Just to make me sure. If I have Fsyg= 14.318MHz and Fsample= 200kHz and N= 2000 (number of samples) then my wanted bin will have number 180 ? Or I've calculated it wrong? Jul 24 '14 at 16:57
• Maybe I do not understand something but I've made a simple model with your inputs and for real input signal wanted bins are 821 and 1181. I've just ploted FFT result and looked it through. Undersampling is downsampling, you may model it by yourself and check it. If your input signal contains noise, you should filter it out before performing downsampling otherwise your result would be totally corrupted. And yes, Goertzel and FFT lead to the same result.
– Serj
Jul 25 '14 at 2:55
• I made FFT aslo in matlab just now. And hmm.. I suppose we both are right. I received 821 and 1181 bin numbers as you did. However I received this (I checked on the graph after ploting) when on the x-axis I had bin numbers 0:N. When I used x-axis as: -N/2:N/2 (centered around 0 and also now the result of positive frequency represents bins from 0:N/2 as it should from what we know about spectrum) then my bin' number is 181. I made mistake above, forgotten to add +1 to my calculations. So, to sum up. The bin' number we want depends of how we show the spectrum on x-axis and we both are right :) Jul 25 '14 at 7:27
• My bin' calculations to receive 181: F0=143180, Fs=200k, N=2000, k=F0/Fs*N, k=143180, but it's to much, as our baseband is only N/2=1000 bins. so: floor(143180/baseband)=143 then: 143180-143*baseband = 180 our wanted bin = 180+1 =181 I don't know if it's commonly used calculation in searching bin number in undersampling but I thing it's working well(assuming our final spectrum is 0:N/2). If you use groetzel then bin' number to calculate equations (k=181 in this example) need to be integer and don't need to be power of 2 as in FFT. Jul 25 '14 at 7:27
• in above should be F0=14.318M. Sorry. Jul 25 '14 at 7:34

If you know the approximate frequency of your signal well enough to know which peak FFT result bin in which it will end up, a complex output Goertzel filter at that bin's frequency over the same length as the FFT will give you the same result for that bin (plus or minus slight numerical differences).

• I know that I can calculate a phase of undersampled signal using fft and then search for highest beam, assuming that I know the signal's frequency. My question here is if I'm able to do the same with Goertzel algorithm? I mean, does the wanted beam to calculate phase using Goertzel will be the same beam as when I would use normal FFT ? I know that there will be alliasing. Jul 23 '14 at 23:19