CODEWITHSUNDEEP
pythonscipyfft
Rohit Pandey
Rohit Pandey

Reputation: 2681

Scipy FFT - how to get phase angle

I'm having trouble getting the phase of a simple sine curve using the scipy fft module in python. I followed this tutorial closely and converted the matlab code to python. However, no matter what phase I use for the input, the graph always shows 3. What am I missing?

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack
import cmath
A=10
fc = 10
phase=60
fs=32#Sampling frequency with oversampling factor of 32

t = np.arange(0,2,1/fs)

#Convert the phase shift to radians from degrees.
phi = phase*np.pi/180

x=A*np.cos(2*np.pi*fc*t+phi)

N=256
X = scipy.fftpack.fftshift(scipy.fftpack.fft(x,N))/N

df=fs/N #Frequency resolution.
sampleindex = np.arange(-N/2,N/2,1) #Ordered index for FFT plot.
f = sampleindex*df #x-axis index continued to ordered frequencies

raw_phases = np.angle(X)

X2=np.copy(X)#Store the FFT results in another array.
#Detect very small numbers and ignore them.
tau = max(abs(X))/10
X2[abs(X)<tau]=0

phase=[cmath.phase(i) for i in X2]
plt.plot(f,phase)
plt.show()

enter image description here

EDIT: Here is some simpler code. Still can't seem to get the phase.

y = 28*np.sin(2*np.pi/7*x)
yf = scipy.fftpack.fft(y)
xf = np.linspace(0.0, 1.0/(2.0*T), N/2)
phase = np.angle(yf)
yf = np.abs(yf[:N//2])
phase = phase[:N//2]
xf = xf[1:]
yf = yf[1:]
phase = phase[1:]
yf = yf-np.mean(yf)
#The frequencies seem to always be scaled by 0.1433, not sure why.
c = 2*np.pi/7/0.143301
freqs = xf[yf>np.std(yf)]*c
phases = phase[yf>np.std(yf)]

The frequencies I get are clustered around 2*np.pi/7. But the phases I get are:

array([-0.217795  , -0.22007488, -0.22226087,  2.91723935,  2.91524011,
    2.91333367])

While there should be no phase at all.

Upvotes: 5

Views: 16431

Answers (2)

Cris Luengo
Cris Luengo

Reputation: 60504

This is the simplest code to show how to get the angles.

Note that I've created the signal y such that there is an integer number of periods in it (as I suggested in a comment, and @hotpaw2 suggested in their answer).

This is the problem with OP's code.

I used linspace to create the time axis, using the endpoint=False option. This is important, if t were to contain the number 10, then we no longer have an exact integer number of periods. When working with the discrete Fourier transform, it helps to think of what happens when you repeat the signal. Simply take y, and concatenate a copy of itself: np.hstack((y,y)). Is this new signal still a sampling of a single sine wave, or did we create something more complex? What happens at that point where the two copies meet?

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

phase = np.pi / 4
t = np.linspace(0, 10, num=200, 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
plt.plot(f, p)
plt.show()

plot of the phase of the Frequency spectrum of a pure sine wave

Upvotes: 3

hotpaw2
hotpaw2

Reputation: 70693

An FFT measures circular phase, referenced to both the very beginning and very end of the input data window. If your input sine wave isn't exactly integer periodic in the FFT aperture, then there will be a discontinuity between the phase at the beginning and end of the window, thus the FFT phase measurement won't be what you might expect.

Instead, I recommend one begin (or reference) the input sinewave at the desired phase in the middle of your data window (at N/2), not the beginning. Then do an circular fftshift before the FFT. The resulting FFT phase measurement will then represent phase with respect to the middle of your original data window, where there is no discontinuity; and atan2() will represent the ratio between the even and odd components of your continuous input waveform as is more typically expected.

Upvotes: 2

Related Questions