5

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.

Rohit Pandey
  • 2,443
  • 7
  • 31
  • 54
  • You should compute the FFT of the full signal, not of the first 256 samples. Your experiment will only work if you have an integer number of periods of the sine wave in your signal. The `N` parameter to `scipy.fftpack.fft` causes the signal to be trimmed (or zero-padded) to `N` samples. Leave it out. – Cris Luengo Jan 31 '19 at 07:01
  • The length of the vector changes from 256 to 64 if I remove the N. Need to look into why. – Rohit Pandey Jan 31 '19 at 07:23
  • Yes, indeed, you make a vector t from 0 to 2 in steps of 1/32. Makes sense. The FFT always has as many samples as the input signal. – Cris Luengo Jan 31 '19 at 07:34
  • @CrisLuengo Added some more code that uses plain, vanilla fft. I notice that for this one, the very first component of the frequency spectrum is always very large. So, I have to ignore it. But still, when I try to get the phases, I get unexpected results. See my edit. – Rohit Pandey Jan 31 '19 at 07:39
  • The zero frequency contains the sum of all sample values. Because your signal does not sample the sine wave over an integer number of periods, its sum is not zero. And this is also the reason you are not getting the expected nice peaks at the frequency of the sine wave. See my answer for some simple code that shows what you want to see. – Cris Luengo Jan 31 '19 at 18:35

2 Answers2

3

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

Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
  • Thanks a lot. Quick follow up - in the real world, we probably won't get data that is an integer multiple of the period. Moreover, the signal could be composed of two sine waves with different periods and phases, making it impossible to have an integer number of periods for both. In such cases, is it impossible to get the phase information? – Rohit Pandey Feb 01 '19 at 07:40
  • 1
    @RohitPandey: You can, but each sine wave will have more than a single peak (like you've noticed). You need to identify the peak in the magnitude plot, and look only at the phase there. You should always use a suitable [windowing function](https://en.wikipedia.org/wiki/Window_function). – Cris Luengo Feb 01 '19 at 07:54
2

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.

hotpaw2
  • 70,107
  • 14
  • 90
  • 153