Octave:Incorrect FFT phase spectrum

Question

A small program that I wrote in octave does not yield desired phase spectrum. The magnitude plot is perfect though.

 f = 200;
 fs = 1000;
 phase = pi/3; 
 t = 0: 1/fs: 1;
 sig = sin((2*pi*f*t) + phase);
 sig_fft = fft(sig);

 sig_fft_phase = angle(sig_fft) * 180/pi;

 sig_fft_phase(201)

sig_fft_phase(201) returns 5.998 (6 degrees) rather than 60 degrees. What am I doing wrong? Is my expectation incorrect?

Answer 1

In your example, if you generate the frequency axis: (sorry, I don’t have Octave here, so Python will have to do—I’m sure it’s the same in Octave):

faxis = (np.arange(0, t.size) / t.size) * fs

you’ll see that faxis[200] (Python is 0-indexed, equivalent to Octave’s 201 index) is 199.80019980019981. You think you’re asking for the phase at 200 Hz but you’re not, you’re asking for the phase of 199.8 Hz.

(This happens because your t vector includes 1.0—that one extra sample slightly decreases the spectral spacing! I don’t think the link @Sardar_Usama posted in their comment is correct—it has nothing to do with the fact that the sinusoid doesn’t end on a complete cycle, since this approach should work with incomplete cycles.)

The solution: zero-pad the 1001-long sig vector to 2000 samples. Then, with a new faxis frequency vector, faxis[400] (Octave’s 401st index) corresponds to exactly 200 Hz:

In [54]: sig_fft = fft.fft(sig, 2000);

In [55]: faxis = np.arange(0, sig_fft.size) / sig_fft.size * fs

In [56]: faxis[400]
Out[56]: 200.0

In [57]: np.angle(sig_fft[400]) * 180 / np.pi
Out[57]: -29.950454729683386

But oh no, what happened? This says the angle is -30°?

Well, recall that Euler’s formula says that sin(x) = (exp(i * x) - exp(-i * x)) / 2i. That i in the denominator means that the phase recovered by the FFT won’t be 60°, even though the input sine wave has phase of 60°. Instead, the FFT bin’s phase will be 60 - 90 degrees, since -90° = angle(1/i) = angle(-i). So this is actually the right answer! To recover the sine wave’s phase, you’d need to add 90° to the phase of the FFT bin.

So to summarize, you need to fix two things:

1. make sure you’re looking at the right frequency bin. For an N-point FFT (and no fftshift), the bins are [0 : N - 1] / N * fs. Above, we just used a N=2000 point FFT to ensure that 200 Hz be represented.
2. Understand that, although you have a sine wave, as far as the FFT is concerned, it gets two complex exponentials, at +200 and -200 Hz, and with amplitudes 1/(2i) and -1/(2i). That imaginary value in the denominators shifts the phase you expect by -90° and +90° respectively.
   • If you happened to have used cos, a cosine wave, for sig, you wouldn’t have run into this mathematical obstacle 😆, so pay attention to the difference between sin and cos in the future 💪!