# Phase extraction from Fourier transform

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?

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.

• ok, thank you for advice! but how can I interpret the result that I plotted? why does the phase have such a behaviour? May 16, 2022 at 14:07
• Probably because the amplitude of the FFT at other frequencies is not exactly zero, due to numerical inaccuracy in the floating point operations, and/or because you're not analyzing an integer number of periods of the sine wave (it's hard to tell from your code what the value of N is).
– MBaz
May 16, 2022 at 14:20
• I've added some code and further explanation to my answer. And you're right: the numerical error does not behave like typical noise, maybe contrary to intuition.
– MBaz
May 17, 2022 at 16:29
• Just make f0 vary between -fs/2 and fs/2 in whatever resolution (step size) you want.
– MBaz
May 17, 2022 at 18:42
• I tried different thresholds to "zero" data and found that it works as it should be from the values of ~ >5e-2, but the general tendency is clear for me know, thank you! May 24, 2022 at 13:14