Phase extraction from Fourier transform

Question

Is it possible in principle to correctly extract the phase from Fourier transform?

I just tried to do so using Python, here some attempts:

# attempt 1, phase=pi/2
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

fs = 1000
T = 1/fs
t = np.linspace(0,N*T, N)
N=len(t)

phase = np.pi/2
signal = 4*np.sin(2*np.pi*2*t-phase)
FFT = np.fft.fft(signal)
freqs = np.fft.fftfreq(len(signal), T)

plt.plot(freqs[0:N//2], 2/N*np.abs(FFT[0:N//2]), label='amplitude')
plt.plot(freqs[0:N//2], np.angle(FFT[0:N//2]), label='phase')

plt.axhline(phase, color='k', linestyle='--', label='phase value')

plt.legend()
plt.grid()

which outputs:

zoomed-in:

The first question that comes to my mind - why does the phase so blurred (why it is not just one point as for the amplitude spectrum?)

if I try another phase:

fs = 1000
T = 1/fs
t = np.linspace(0,N*T, N)
N=len(t)

phase = -np.pi/8
signal = 4*np.sin(2*np.pi*2*t-phase)
FFT = np.fft.fft(signal)
freqs = np.fft.fftfreq(len(signal), T)

plt.plot(freqs[0:N//2], 2/N*np.abs(FFT[0:N//2]), label='amplitude')
plt.plot(freqs[0:N//2], np.angle(FFT[0:N//2]), label='phase')

plt.axhline(phase, color='k', linestyle='--', label='phase value')

plt.legend()
plt.grid()

it outputs:

The second question that I have - why does the phase exist in the whole frequency range, not only the frequency I'm interested in?

According to these examples I can't extract exact phase of original signal, is it correct?

Answer 1

You need to evaluate the DFT at the frequency of the input signal.

I recommend doing this experiment with the DTFT instead of the DFT. You can implement it in just a few lines of code, and unless the signal is very long, it will still be calculated very quickly on modern computers.

Then, if the input signal is $x[n] = \cos(2\pi f_0 kT_s + \phi)$ and its DTFT is $X(f)$, then $\angle X(f_0) = \phi$.

Example code, in Julia, that illustrates estimating the phase of a single sinusoid using the DTFT:

fs = 1000
ts = 1/fs
t = range(0, 1-ts, step=ts)
# signal parameters
A = 4       # amplitude
f0 = 10     # frequency
phi = -pi/8 # phase
signal = A*cos.(2*pi*f0*t .+ phi)
# DTFT evaluated at f = f0
dtft_tone = exp.(-2im*pi*f0*t)
dtft_f0 = (1/length(t))*sum(signal .* dtft_tone)
# estimate parameters
A_est = 2*abs(dtft_f0)   # 4.0
phi_est = angle(dtft_f0) # -0.3927 == -pi/8

Now, change the analysis frequency:

dtft_tone = exp.(-2im*pi*20*t)

Here we're calculating the amplitude and phase of the signal at frequency 20. Both should be zero. However, we get an amplitude of 5.770404116019941e-16 and a phase of 1.1063743245570448. The amplitude is not zero due to numerical inaccuracy in the calculation, which is essentially unavoidable. In consequence, we also get a phase measurement, but it is meaningless.

One way out of it is to "zero" all DTFT (or DFT) measurements with an amplitude less than a threshold, say 1e-6. In your code, you would go through every element of FFT and force to zero those that are too small.