Problem with phase in FFT analysis

I have such problem. I am generating signal which is sinusoid with given parameters (such as amplitude, frequency and phase). Next I am calculating amplitude spectrum and phase spectrum in this way:

fs=1000; %sampling frequency

t=0:1/fs:1; %time

A=1; %amplitude

f=100; %frequency

fi=30*pi/180; %phase (calculated form deg. to rad.)

x = A*cos(2*pi*f*t+fi); %signal

X=fft(x); %FFT

df=fs/length(X); %frequency resolution

f=(0:1:length(X)/2)*df; %frequency axis

subplot(2,1,1);

M=abs(X)/length(x)*2; %amplitude spectrum

plot(f,M(1:length(f)));

subplot(2,1,2);

P=angle(X)*180/pi; %phase spectrum (in deg.)

plot(f,P(1:length(f)));

With amplitude spectrum is everything OK (for frequency 100Hz I have peak equal to 1 what is equal to amplitude) but phase spectrum is not OK (for 100Hz I have 50 deg). Please tell me if my code is OK or maybe I am not interpreting result in right way or maybe I should make some correction of phase results (as thought in case of amplitude spectrum). What is interesting the value of phase is changing due to change of signal frequency or signal sampling what in theory should not happened.

Regards

[Walter Roberson]

It might make a bit of difference to use

t = linspace(0,1,fs+1);

to avoid accumulated round-off errors that 0:1/fs:1 would have.

Answer 1

[David Young]

The problem is that there is not an exact whole number of cycles of the sine wave in the original time series. If you look at the values of x(1) and x(end), you'll find they are the same, because the final point is at time t=1, when a new cycle is just starting.

Unless you have an exactly whole number of cycles of a sine wave present in the initial series, there won't be a sample point in the frequency domain lying at exactly the frequency of the wave. In your case, the sample at X(101) (with its complex conjugate at X(902)) is close to, but not exactly at, the frequency f=100. This slight shift in frequency doesn't affect the amplitude too much, but it has a big effect on the phase.

A minimal change to get the answer you expect is to reduce the length of the input by 1, for example like this:

t=0 : 1/fs : 1-0.5/fs; %time

which will give a whole number of cycles provided f is an integer. (As Walter Roberson points out, linspace may be preferable - but rounding errors aren't in fact the problem in this case.) This change gives the phase as expected, as you can see by printing the value of P(101).

There's a general lesson: if you have real-world data, you are very unlikely to have an exactly whole number of cycles in the sequence, and so you can't get the phase easily from the FFT.