baffled by fft phase spectrum!

Question

A very simple MATLAB experiment:

f = 200;  
fs = 1000;  
t = 0: 1/fs : 1;
x = cos(2*pi*f*t);  
plot(angle(fftshift(fft(x))));  

And here's the output:

Now, made a minor change to the above code snippet; reducing the time duration by just 1 sample, as follows:

f = 200;  
fs = 1000;  
t = 0: 1/fs : 1 - 1/fs;
x = cos(2*pi*f*t);  
plot(angle(fftshift(fft(x)))); 

And the phase spectrum goes totally crazy:

Questions:

1. In the first plot, I was hoping to see a zero phase at bin 700 which corresponds to the positive frequency of 200 in this example. That doesn't seem to be the case. Secondly, I don't understand the linear parts of the graph in plot 1. I do appreciate phase components that could potentially exist due to the so-called numerical noise, but then how could that noise be so 'linear' in phase?

2. In the second plot, why removing just one sample would have such a drastic impact on the phase plot?

3. Am I doing something fundamentally wrong here?

Answer 1

You aren't doing anything wrong, but you also aren't thinking carefully about what you should expect to see, which is why you're surprised at the result. For question 1, your conjecture is close, but you actually have things backwards; it's numerical noise that's plaguing your second one, not your first one.

Pictures may help. Here's plots of the magnitude and phase for the first trial:

x = Cos[2.0 \[Pi] 200 Range[0, 1, 1/1000]];
fx = Fourier[x];
ListLinePlot[Abs[fx], PlotRange -> All]
ListLinePlot[Arg[fx], PlotRange -> All]

And the second one:

x = Cos[2.0 \[Pi] 200 Range[0, 1 - 1/1000, 1/1000]];
fx = Fourier[x];
ListLinePlot[Abs[fx], PlotRange -> All]
ListLinePlot[Arg[fx], PlotRange -> All]

So what's going on here? The second one is easiest to explain. First, the FFT of the second one has zero amplitude everywhere except at the two peaks visible in the magnitude spectrum; this follows since the definition of the FFT using 1,000 datapoints returns frequencies of the form $k/1000$ for $0\leq k\leq 999$, and hence your signal falls exactly on a frequency bin. As a result, at the other 998 points your signal is entirely due to machine noise due to floating point errors, and hence your phase spectrum is nonsense, because it's literally the phase of pseudorandom numbers.

For the first one, however, the definition of the FFT encompasses frequencies of the form $k/1001$ for $0\leq k \leq 1000$, whereas your signal frequency is $200/1000$, which is not of the form $k/1001$. As a result, your signal is broadened by spectral leakage, and will be nonzero almost everywhere. I will not comment on the physical form of the phase plot, but I will say that it admits a closed analytic form.

In general, I think graphs of phase angle alone are a really bad idea for conveying information, precisely for this reason; first, you can't tell if you're looking at the phase of low-amplitude garbage or actual signal, and second, it's not translation-invariant, and it's easy to get utterly bewildering graphs for simple inputs. Much better, if you're still looking for something which conveys phase information, is a graph which simultaneously portrays phase and amplitude information in the same visual manner, such as a plot where phase is encoded as hue and magnitude is encoded as brightness.

ADDENDUM: Here's a couple pictures from Mathematica which illustrate the principle I stated in the preceding paragraph:

hue = Compile[{{z, _Complex}}, {(1.0 Arg[-z] + \[Pi])/(2 \[Pi]), 
Exp[1 - Max[Abs[z], 1]], Min[Abs[z], 1]}, 
CompilationTarget -> "C", RuntimeAttributes -> {Listable}];
L = 500;
data = Table[Boole[x <= 11 && y <= 11], {x, L}, {y, L}];
Image[hue@
RotateRight[
10 RotateRight[Fourier[RotateRight[data, {-5, -5}]], {L/2, L/2}]], 
ColorSpace -> Hue, Magnification -> 1]
Image[hue@
RotateRight[
10 RotateRight[Fourier[RotateRight[data, {-4, -4}]], {L/2, L/2}]], 
ColorSpace -> Hue, Magnification -> 1]
Image[hue@
RotateRight[
10 RotateRight[Fourier[RotateRight[data, {0, 0}]], {L/2, L/2}]], 
ColorSpace -> Hue, Magnification -> 1]

All three images are 2D Fourier transforms of the same input signal (an $11\times 11$ square of 1's padded with zeros to a length of $500\times 500$), but the inputs have been cyclically rotated by 5, 4, and 0, and 200 datapoints. The magnitude spectra (encoded by pixel brightness) are identical, but the phase spectra are completely different! Phase encoding is done so that 1 maps to red, $i$ maps to green, $-1$ maps to cyan, and $-i$ maps to purple. This is what I mean when I say that phase spectra are non-shift invariant, and thus not amenable to human visual understanding. For example, with a cyclic shift of 200 datapoints, it's completely impossible to tell what's going on in the phase, since it just looks like static, but the input signal is no more complicated than the other input cases.

Answer 2

If you want to vary the frequency of a signal or the FFT length so that the signal varies between exactly periodic and not exactly periodic in the FFT aperture, and don't want to see the phase of the peak magnitude bin for that signal change, one can reference the initial phase of the signal to the center if the FFT aperture instead of the beginning (for a generated sin(t), put t=0 in the center of the FFT array).

Answer 3

The Gaussian Waves site details the part about the phase and its random-like behaviour: just a matter of floating-point error as DumpsterDoofus said