Goertzel algorithm to get the phase?

Question

I am using Goertzel algorithm to get the amplitude of a certain frequency. I am trying now to get the phase from it, and I don't know how.

Can some one explain, and show me how to get the phase of a certain-f from this code?

Also, I am using it to 16khz, with sample rate 44.1. What's the smallest length of samples that I can run it on?

double AlgorithmGoertzel( int16_t *sample,int sampleRate, double Freq, int len )
{
    double realW = 2.0 * cos(2.0 * M_PI * Freq / sampleRate);
    double imagW = 2.0 * sin(2.0 * M_PI * Freq / sampleRate);
    double d1 = 0;
    double d2 = 0;
    double y;
    for (int i = 0; i < len; i++) {
        y=(double)(signed short)sample[i] +realW * d1 - d2;
        d2 = d1;
        d1 = y;
    }
    double rR = 0.5 * realW *d1-d2;
    double rI = 0.5 * imagW *d1-d2;

    return (sqrt(pow(rR, 2)+pow(rI,2)))/len;
}

Answer 1

I don't think the algorithm consists of multiplying the sequence by a constant, but by the complex signal exp(n*i*2pi*freq/samplerate); 0<=n<=length, and getting the average magnitude (or power of the signal).

As the complex output is R*exp(i theta), R gives the power at the given frequency and theta gives the phase. (theta == atan2 ( imag, real))

Answer 2

Do a rectangular to polar conversion. That will give you phase and magnitude.

magnitude = sqrt ((Vreal * Vreal) + (Vimag * Vimag))

phase = atan2 (Vimag, Vreal)

Answer 3

The number of samples you need to feed a Goertzel filter will be inversely proportional to your desired or required filter bandwidth. A Goertzel provides a Sinc shaped bandpass filter, with the main lobe width proportional to 2*Fs/N.

If you use a complex Goertzel, the resulting phase will be relative to some point in the filter's data window. You may thus have to calculate an offset to get phase relative to some other reference point in time.