Plotting and extracting fft phase

Question

Here is a code that compares fft phase plotting with 2 different methods :

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack

phase = np.pi / 4
f = 1
fs = f*20
dur=10
t = np.linspace(0, dur, num=fs*dur, endpoint=False)
y = np.cos(2 * np.pi * t + phase)
Y = scipy.fftpack.fftshift(scipy.fftpack.fft(y))
f = scipy.fftpack.fftshift(scipy.fftpack.fftfreq(len(t)))

p = np.angle(Y)
p[np.abs(Y) < 1] = 0

fig, ax = plt.subplots(2, 1)
ax[0].plot(t, y)
ax[1].plot(f*fs, p, label='from fft')
ax[1].phase_spectrum(y, fs, window=None, label='from phase_spectrum')
plt.legend()
plt.show()

here is the result :

Here is result when signal number of period is not an integer :

I have several questions :

• why are phase plot with phase_spectrum or using fft and the angle so different? Using fft and then np.angle produces good results, but how can we explain magnitude_spectrum result?
• Here we are in a very simple case where we have sine periodic signal with N periods For if I have a widepand signal and I want to extract phase at f how can I do? For example here with both methods presented in example, I'm not sure I can extract a precise phase. With phase_spectrum, at f = 1 I cannot find back pi/4. And with fft and then np.angle, in order to extract the good phase I need to be sure signal number of period is an integer.

Answer 1

Before the answer, just a small note:
Remove the p[np.abs(Y) < 1] = 0 line. Most of your spectrum has magnitude below 1, which is why, with this line, your spectrum looks mostly like a flat line at zero.

Now to the answer:
phase_spectrum does three things different than you:

• It applies phase unwrapping.
  • If you want to apply phase unwrapping in your code, just do np.unwrap(np.angle(Y)).
  • If you want matplotlib to plot the spectrum without unwrapping, use angle_spectrum instead.
• It applies a windowing function to the data before computing the spectrum.
  • I know you passed a window=None, but, for some reason, matplotlib decided that window=None means "use a hanning window, please" (see the docs).
  • If you do not want matplotlib to apply the window, one solution is to pass window=lambda x: x.
    • The docs actually suggest passing window=matplotlib.mlab.window_none, but the source for it is just a def window_none(x): return x.
• It computes the one-sided version of your spectrum. The docs say this is the default whenever the input is real, not complex.
  • To get the normal two-sided version, pass sides='twosided' to the phase_spectrum call.

---

Now, about getting the phase at a frequency f:

To do this, you must use the phase without unwrapping.

You are right that you can't directly extract the phase of the single tone signal if you do not have an integer number of cycles. That is because the singal's frequency does not fall exactly on top of any frequency bin in the FFT. You can get an approximation with the phase of the nearest bin, though. You could also do a sinc interpolation of the spectrum to get its value at the frequency you want.

If you only care about the phase of a single frequency f, then you shouldn't use FFT at all. The FFT computes the phase and amplitue at all frequencies. If you only care about a single frequency, just do Y_at_f = y @ np.exp(2j * np.pi * f * t) and get that phase by np.angle(Y_at_f).

Answer 2

You can extract the phase referenced to the center of your data window by doing an fftshift (circular rotate by N/2) before the FFT. This is because, after an fftshift, atan2() always is related to the ratio of oddness to eveness of the data around its center (as decomposed to an odd function plus an even function).

So calculate the phase of the signal in the middle of your window during its generation, and use that instead of the phase at the beginning.