Extract phase information using numpy fft

Question

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

However, I find that although I can accurately capture the frequency of my cosine wave, the phase is inaccurate if I have not tried at an extremely high speed. For instance:

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 phase 1.57173880459.

Using num = 10000 points, I get phase 1.58022110476.

Using num = 1000 points, I get phase 1.6650441064.

What will go wrong? Even with 1000 points, I have 33 points per cycle, which should be enough to resolve this. Maybe there is a way to increase the number of calculated frequency points? Is there a way to do this with a "low" number of points?

EDIT: from further experiments, it seems that I need ~ 1000 points per cycle to accurately extract the phase. Why?

EDIT 2: Further experiments show that accuracy is related to the number of points per cycle, and not to absolute numbers. An increase in the number of sample points per cycle makes the phase more accurate, but if the signal frequency and the number of sample points increase by the same factor, the accuracy remains unchanged.

Answer 1

Your points are not distributed equally over the interval, you have a point at the end doubled: 0 - this is the same point as 1 . It becomes less important the more points you take, but it still gives some error. You can completely avoid this, linspace has a flag for this. He also has a flag that takes you back dt directly with the array.

Do

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

instead

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

then it works :)

Answer 2

The phase value in the resulting hopper of a non-rotating FFT is valid only if the input signal is exactly integer, periodic along the length of the FFT. Your test signal is missing, so the FFT measures something partially related to the phase difference of the signal break between the end points of the test sine wave. A higher sampling rate will create a slightly different last endpoint from the sine wave and, therefore, a possibly smaller gap.

If you want to reduce this FFT phase measurement error, create a test signal so that your test phase is tied to the exact center (sample N / 2) of the test vector (and not to sample 1), and then do fftshift (rotation by N / 2) so that there is no signal break between the 1st and last point of your resulting FFT vector of length N.