Why can’t I get a FFT phase spectrum with my code?

  fft, numpy, phase, python

I’ve been trying to do a fourier analysis to some discreet set of values corresponding to more or less sinusoidal functions. The frequency spectrum is giving great values, but when I try to do the phase-frequency spectrum it all goes wrong: the points in the spectrum is all over the place, and does make any physical or mathematical sense. Can anyone tell me what’s wrong with my code? I have tried using implementing the phase calculation with two method (one is simply less as commentary).

import numpy as np
import matplotlib.pyplot as plt
import math

LX=64
LY=32
font = '12'
font2 = '14'
lwd = 0.5


arr=np.loadtxt("output/obstPos.dat",delimiter=' ')
t = arr[:,0]
x = arr[:,1]
plt.plot(t, x, linestyle="-", color = "r", linewidth=lwd)
x_detrended = signal.detrend(x)     #removes any linear trend (0 Hz)


plt.xlabel(r"$rm t$", fontsize=font2)
plt.ylabel(r'$rm x$',  fontsize=font2)
plt.yticks(fontsize=font)
plt.xticks(fontsize=font)
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0), useMathText=True)
plt.savefig('output/pos-time.png', dpi=600, bbox_inches='tight')
plt.clf()


#frequency measurement (Fourier)

time_step = 500 #step between the points
n = x.size
sp = np.fft.fft(x_detrended)
freq = np.fft.fftfreq(n, d=time_step)

#plot frequency spectrum
plt.plot(freq, np.abs(sp), '.', color='blue')

plt.ylabel(r"$rm vert N_{x}vert$", fontsize=font2)
plt.xlabel(r'$rm freq $',  fontsize=font2)
plt.yticks(fontsize=font)
plt.xticks(fontsize=font)
plt.rcParams['font.size']=font
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0), useMathText=True)
plt.xlim(-0.0005, 0.0005)
plt.ylim(0, 6000)
plt.savefig('output/frequencies.png', dpi=200, bbox_inches='tight')
plt.clf()

#plot frequency-phase spectrum
#phase = np.angle(sp)
SP = sp
tau = max(abs(sp))/10000
SP[abs(sp) < tau] = 0
phase=np.arctan2(np.imag(SP),np.real(SP))
plt.plot(freq, phase, '.', color='blue')

plt.ylabel(r"$rm phase$", fontsize=font2)
plt.xlabel(r'$rm freq $',  fontsize=font2)
plt.yticks(fontsize=font)
plt.xticks(fontsize=font)
plt.rcParams['font.size']=font
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0), useMathText=True)
plt.xlim(-0.0005, 0.0005)
plt.ylim(-3.2, 3.2)
plt.savefig('output/phases.png', dpi=200, bbox_inches='tight')
plt.clf()

#select the interval which is the peak
index = np.where(sp == np.amax(sp[1:n]))[0][0]
freqMax = np.abs(freq[index])
print("The oscillation frequency is: %.8f"%freqMax)
phaseMax = phase[index]
print("The oscillation phase is: %.8f"%phaseMax)```

Source: Python Questions

LEAVE A COMMENT