为何我的Python傅里叶分析代码无法获取正确的FFT相位谱

代码存在的核心错误 #

• 第一个错误：频谱数组浅拷贝导致原始数据被篡改
  你写的SP = sp属于numpy数组的引用赋值，修改SP的元素时原始的频谱数组sp也会被同步修改，后续你找峰值用的sp已经是被置0操作修改过的错误数组，需要改用深拷贝SP = sp.copy()。
• 第二个错误：相位谱绘制了无意义的噪声分量相位
  你对幅度低于阈值的频谱分量做了置0处理，但依然把所有频率点（包括置0的噪声点）的相位全部绘制了出来。这些噪声分量的相位本身就是随机无意义的，自然会显示为杂乱分布，你应该只绘制幅度大于阈值的有效分量的相位。
• 第三个错误：峰值索引查找逻辑错误
  你用np.amax(sp[1:n])找峰值，直接对复数数组取最大值，numpy中复数对比默认仅比较实部，根本无法找到实际幅度最大的峰值点，正确做法是对幅度取最大值np.amax(np.abs(sp[1:n]))。
• 第四个错误：缺失依赖库导入
  代码中使用了signal.detrend但没有导入对应模块，需要在开头添加from scipy import signal才能正常运行。

修正后的核心代码片段 #

# 首先补全导入
import numpy as np
import matplotlib.pyplot as plt
import math
from scipy import signal # 新增缺失的导入

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.copy() # 改为深拷贝，不修改原始sp数组
tau = max(abs(sp))/10000
SP[abs(sp) < tau] = 0
phase=np.arctan2(np.imag(SP),np.real(SP))
valid_idx = np.abs(SP) > 0 # 过滤掉置0的噪声点
plt.plot(freq[valid_idx], phase[valid_idx], '.', 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(np.abs(sp) == np.amax(np.abs(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)

额外优化建议 #

• 如果你处理的是实信号，还可以只绘制正频率部分的频谱和相位，避免重复显示共轭对称的负频率分量，结果会更直观。
• 阈值tau可以根据实际噪声水平调整，避免过滤掉有效小分量或者保留过多噪声。

内容的提问来源于stack exchange，提问作者Tomé Silva