python - 为什么我的代码无法获得 FFT 相位谱？

Question

我一直在尝试对与或多或少的正弦函数相对应的一组谨慎的值进行傅立叶分析。频谱给出了很大的值，但是当我尝试做相位频谱时，一切都出错了：频谱中的点到处都是，并且确实具有任何物理或数学意义。谁能告诉我我的代码有什么问题？我尝试过使用两种方法来实现相位计算（一种只是较少的评论）。

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)```

标签： pythonnumpyfftphase