fft - Is there an algorithm to compute the phase for a single frequecy?

If you have a function $f(t)=A \cdot \sin(\omega t+\phi)$, and reference sin wave $\sin(\omega x)$ what would be a fast algorithm to compute $\phi$?

I was looking at Goertzel algorithm, but it doesn't seem to deal with phase?

Answer

Use a DFT at the specific frequency. Then compute amplitude and phase from the real/imag parts. It gives you the phase referenced to the start of the sampling time.

In a 'normal' FFT (or a DFT computed for all N harmonics), you typically compute frequency with f = k*(sample_rate)/N, where k is an integer. Although it may seem sacrilegious (especially to members of the Church of the Wholly Integer), you can actually use non-integer values of k when doing a single DFT.

For instance, suppose you've generated (or obtained) N = 256 points of a sine wave of 27 Hz. (let's say, sample_rate = 200). Your 'normal' frequencies for a 256 point FFT (or N point DFT) would correspond to: f = k*(sample_rate)/N = k*(200)/256, where k is an integer. But a non-integer 'k' of 34.56 would correspond to a frequency of 27 Hz., using the parameters listed above. It's like creating a DFT 'bin' that is exactly centered at the frequency of interest (27 Hz.). Some C++ code (DevC++ compiler) might look as follows:

#include 
#include 
#include 
#include 

using namespace std;

// arguments in main needed for Dev-C++ I/O
int main (int nNumberofArgs, char* pszArgs[ ] ) {
const long N = 256 ;
double sample_rate = 200., amp, phase, t, C, S, twopi = 6.2831853071795865; 
double  r[N] = {0.}, i[N] = {0.}, R = 0., I = 0. ;
long n ;

// k need not be integer

double k = 34.56;

// generate real points
for (n = 0; n < N; n++) {
    t =  n/sample_rate;
    r[n] = 10.*cos(twopi*27.*t - twopi/4.);
}  // end for

// compute one DFT
for (n = 0; n < N; n++) {

    C = cos(twopi*n*k/N); S = sin(twopi*n*k/N);
    R = R + r[n]*C + i[n]*S;
    I = I + i[n]*C - r[n]*S;
} // end for

cout<<"\n\ndft results for N = " << N << "\n";
cout<<"\nindex k     real          imaginary       amplitude         phase\n";

amp = 2*sqrt( (R/N)*(R/N) + (I/N)*(I/N) ) ;
phase = atan2( I, R ) ;

// printed R and I are scaled
printf("%4.2f\t%11.8f\t%11.8f\t%11.8f\t%11.8f\n",k,R/N,I/N,amp,phase);

cout << "\n\n";
system ("PAUSE");
return 0;
} // end main

//**** end program

(PS: I hope the above above translates well to stackoverflow – some of it might wrap around)

The result of the above is a phase of -twopi/4, as shown in the generated real points (and amp is doubled to reflect the pos/neg frequency).

A few things to note – I use cosine to generate the test waveform and interpret results – you have to be careful about that – phase is referenced to time = 0, which is when you started sampling (ie: when you collected r[0]), and cosine is the correct interpretation).

The above code is neither elegant nor efficient (eg: use a look-up tables for the sin/cos values, etc.).

Your results will get more accurate as you use larger N, and there's a little bit of error due to the fact that the sample rate and N above are not multiples of each other.

Of course, if you want to change your sample rate, N, or f, you'd have to change the code and the value of k. You can plunk down a DFT bin anywhere on the continuous frequency line – just make sure that you're using a value of k that corresponds to the frequency of interest.