Search code examples
pythonnumpyfft

extracting phase information using numpy fft

I am trying to use a fast fourier transform to extract the phase shift of a single sinusoidal function. I know that on paper, If we denote the transform of our function as T, then we have the following relations:

enter image description here

However, I am finding that while I am able to accurately capture the frequency of my cosine wave, the phase is inaccurate unless I sample at an extremely high rate. For example:

import numpy as np
import pylab as pl

num_t = 100000
t = np.linspace(0,1,num_t)
dt = 1.0/num_t
w = 2.0*np.pi*30.0
phase = np.pi/2.0

amp = np.fft.rfft(np.cos(w*t+phase))
freqs = np.fft.rfftfreq(t.shape[-1],dt)

print (np.arctan2(amp.imag,amp.real))[30]

pl.subplot(211)
pl.plot(freqs[:60],np.sqrt(amp.real**2+amp.imag**2)[:60])
pl.subplot(212)
pl.plot(freqs[:60],(np.arctan2(amp.imag,amp.real))[:60])
pl.show()

Using num=100000 points I get a phase of 1.57173880459.

Using num=10000 points I get a phase of 1.58022110476.

Using num=1000 points I get a phase of 1.6650441064.

What's going wrong? Even with 1000 points I have 33 points per cycle, which should be enough to resolve it. Is there maybe a way to increase the number of computed frequency points? Is there any way to do this with a "low" number of points?

EDIT: from further experimentation it seems that I need ~1000 points per cycle in order to accurately extract a phase. Why?!

EDIT 2: further experiments indicate that accuracy is related to number of points per cycle, rather than absolute numbers. Increasing the number of sampled points per cycle makes phase more accurate, but if both signal frequency and number of sampled points are increased by the same factor, the accuracy stays the same.

Solution

  • Your points are not distributed equally over the interval, you have the point at the end doubled: 0 is the same point as 1. This gets less important the more points you take, obviusly, but still gives some error. You can avoid it totally, the linspace has a flag for this. Also it has a flag to return you the dt directly along with the array.

    Do 

    t, dt = np.linspace(0, 1, num_t, endpoint=False, retstep=True)

    instead of

    t = np.linspace(0,1,num_t)
dt = 1.0/num_t

    then it works :)