recovering phase of sine signal from FFT

Question

I have a simple sine function as $sin(2\pi ft + \phi)$. I want to obtain the phase signal $\phi$. I tried to use FFT to calculate $\phi$. In matlab I do the following

f=200; %frequency of sine wave
overSampRate=30; %oversampling rate
fs=overSampRate*f; %sampling frequency
phase = 3/5*pi; %desired phase shift in radians
nCyl = 5; %to generate five cycles of sine wave

t=0:1/fs:nCyl*1/f; %time base

x=sin(2*pi*f*t+phase); %replace with cos if a cosine wave is desired

NFFT=1024; %NFFT-point DFT
X=fft(x,NFFT); %compute DFT using FFT
XX=2*abs(X(1:NFFT/2+1));
[tt ind]=max(XX);
phase_Estimate=angle(X(ind);

This result makes almost no sense to me. For example, when $\phi=0.523$, phase_Estimate is obtained $-0.98$.

Answer 1

It makes sense after all: your phase estimate is an estimate of the unknown phase angle $\phi$ minus $\pi/2$, because you're using a sine function. So the phase of the FFT at the appropriate frequency index corresponds to $\phi-\pi/2$. So in your example, $\phi\approx -0.98+\pi/2=0.59$, which is a lot closer to the actual phase. If you used a cosine function, your way of estimating the phase would be correct.

Answer 2

Using an non-interpolated FFT result phase only works if the period of sinusoid is exactly an integer submultiple of the FFT length. In your example, the sine wave isn't integer periodic in aperture.

If not, you will need to interpolate the phase to get a better estimate. Here's one method to get an better interpolated phase:

First fftshift (rotate by N/2) the data to move the zero phase reference point to the center of the window before doing the FFT. (This is needed to keep the phase from flipping/alternating between adjacent FFT result bins. * )

Then do the FFT and estimate the frequency of the sinusoid by parabolic or, better yet, Sinc interpolation.

Then use the estimated frequency to linearly interpolate the phase between the nearest two FFT result bin phases.

Then use the estimated frequency and phase at the center of the window to calculate the phase at some other point, such as the beginning of the FFT window.

Also note that the phase of a sine is different from the phase of a cosine wave by pi/2. atan(im,re) returns the cosine phase.

(* as an alternative to pre-fftshit-ing the data, one could also post-flip the phase of the odd FFT result bins.)