Difficulty understanding the phase calculated by the FFT. Short matlab demo to illustrate

Question

I'm testing the phase output of an fft of a sin signal and a cos signal. The script below creates the signals and performs an FFT on them. Bins who's amplitude is below a threshold are zeroed for the phase spectrum because I am only interested in the phase of the signals.

% 10khz 10 second long time interval
t = 0:1 / 10000:10;

%1khz cos
c = cos(2 * pi * 1000 .* t);
%1khz sin
s = sin(2 * pi * 1000 .* t);

%ffts
C = fft(c)/length(c);
S = fft(s)/length(s);

%magnitude and phases of ffts
CA = abs(C); %cos magnitude
SA = abs(S); %sin magnitude

Cthresh = max(CA) * 0.5;
Sthresh = max(SA) * 0.5;

%find all indeces below the threshold
Crange = find(CA < Cthresh);
Srange = find(SA < Sthresh);

%set the indeces below the threshold to 0 - phase will be meaningless for
%noise values
CP = angle(C);
CP(Crange) = 0;
SP = angle(S);
SP(Srange) = 0;

If you plot CP - the phase of the cos - you will get a phase of 0.3142 in the bins corresponding to the frequency of the cos signal and zeros elsewhere. This is pi/10. I'm expecting to get pi. Why is this?

If you plot SP you get values of 1.2566. I'm expecting to get pi/2 or 1.5708. 80% of the expected value. What is causing these errors?

Answer 1

If your input signal is not perfectly periodic in the FFT aperture length (an exact integer number of full periods), the sinusoids will be discontinuous across the ends of the FFT aperture. Thus you will get a phase that is the average of the two different phases at both ends of the FFT input vector.

If you want a more sensible phase, reference the phase of your sinusoids to the center of the FFT input vector, and do an fft shift before the FFT. This will result in a continuous sinusoid at the zero phase reference position, with a single phase instead of a weird average value.

Also note that matlab may reference the phase to the second point in a sampled sinusoid, e.g. vectorElement[i=1], not the first, vectorElement[i=0]. This will have a phase of pi/10 for a sinusoid of period = 20 samples.

Answer 2

The issue you have is exactly what hotpaw2 has stated. You have 100001 samples in t, so you do not have a perfectly periodic signal, therefore you have leakage. That means that you've got a sin()/sin() function from your implicit rectangular window convolved with your solution. That's what changes your phase.

If instead you try the following:

t = 0:1 / 10000:9.9999;
c = cos(2 * pi * 1000 .* t);
%1khz sin
s = sin(2 * pi * 1000 .* t);
%ffts
C = fft(c)/length(c);
S = fft(s)/length(s);

you would find that the phase of the cosine is zero (which is what you would expect) and that the phase of sine is pi/2.

Performing a linear shift in the time domain (using fftshift) will merely introduce a linear phase term in the frequency domain, and will not resolve the original problem.

In practice, rather than trying to set the length of the sequence precisely to match the period of the signal, windowing should be applied if the signal is to be examined in the frequency domain. In that case you really should make sure that your signals are appropriately aligned so that the window attenuates the end points, thus smoothing out the discontinuity. This has the effect of broadening the main lobe of the FFT, but it also controls the leakage.