extracting phase information using numpy fft

user_2521423 0
Viewed 8861

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.

fft numpy python
3 Answers

This snippet of code might help:

def reconstruct_ifft(data):
"""
In this function, we take in a signal, find its fft, retain the dominant modes and reconstruct the signal from that

Parameters
----------
data : Signal to do the fft, ifft

Returns 
-------
reconstructed_signal : the reconstructed signal

"""
N = data.size
yf = rfft(data)

amp_yf = np.abs(yf) #amplitude

yf = yf*(amp_yf>(THRESHOLD*np.amax(amp_yf)))

reconstructed_signal = irfft(yf)

return reconstructed_signal

The 0.01 is the threshold of amplitudes of the fft that you would want to retain. Making the THRESHOLD greater(more than 1 does not make any sense), will give fewer modes and cause higher rms error but ensures higher frequency selectivity. (Please adjust the TABS for the python code)

user_159281340
Related