如何用FFT将加速度计时间序列数据转换为功率谱？绘图异常求助

我用加速度计采集了某部件的振动数据，生成了包含accX、accY、accZ列的时间序列DataFrame。计划通过快速傅里叶变换（FFT）将时域数据转换为频域数据以得到功率谱，使用numpy库的fft函数尝试后，在放大查看accY的功率谱绘图时发现异常。

已尝试的代码 #

初始实现代码 #

# 对每列应用FFT
fft_result_X = np.fft.fft(accX)
fft_result_Y = np.fft.fft(accY)
fft_result_Z = np.fft.fft(accZ)

# 计算每列的功率谱幅值
power_spectrum_X = np.abs(fft_result_X) ** 2
power_spectrum_Y = np.abs(fft_result_Y) ** 2
power_spectrum_Z = np.abs(fft_result_Z) ** 2

# 生成频谱对应的频率轴
frequencies = np.fft.fftfreq(len(accX), 1/62.5)

单列绘图测试代码 #

data = dado['accX'].values
ps = np.abs(np.fft.fft(data))**2

time_step = 1 / 62.5  # 采样间隔，对应采样频率62.5Hz
freqs = np.fft.fftfreq(data.size, time_step)
idx = np.argsort(freqs)

plt.plot(freqs[idx], ps[idx])
plt.show()

该测试代码成功生成了accX的功率谱图，但accY的绘图放大后出现异常，以下是针对性的排查与解决方法：

异常排查与解决步骤 #

• 处理频谱对称性：numpy的FFT输出包含正、负频率对称分量，全范围绘图会导致重复数据，放大时易出现视觉混乱。只需保留正频率部分（0到Nyquist频率，即采样频率的一半）即可：

  # 计算Nyquist频率对应的索引
  nyquist_idx = len(frequencies) // 2
  # 提取正频率及对应功率谱
  positive_freqs = frequencies[:nyquist_idx]
  power_spectrum_Y_pos = power_spectrum_Y[:nyquist_idx]
  # 绘制正频率范围的功率谱
  plt.plot(positive_freqs, power_spectrum_Y_pos)
  plt.xlabel('频率 (Hz)')
  plt.ylabel('功率谱')
  plt.show()
  

• 去除直流分量：振动数据若存在非零均值（直流分量），会在0Hz处产生超大峰值，掩盖其他频段的细节。先对accY做去均值处理：

  # 去均值
  accY_centered = accY - np.mean(accY)
  # 重新计算FFT和功率谱
  fft_result_Y = np.fft.fft(accY_centered)
  power_spectrum_Y = np.abs(fft_result_Y) ** 2
  

• 验证数据有效性：检查accY列是否存在缺失值、异常值：

  # 检查缺失值数量
  print(accY.isna().sum())
  # 绘制箱线图查看异常值分布
  plt.boxplot(accY)
  plt.show()
  

  若存在缺失值可通过accY.fillna(method='ffill')填充，异常值可通过截断或插值处理。
• 聚焦关键频段：放大绘图时指定x轴范围，避免全范围绘图导致细节被淹没，比如聚焦0到30Hz（采样频率62.5Hz的Nyquist频率为31.25Hz）：

  plt.plot(positive_freqs, power_spectrum_Y_pos)
  plt.xlim(0, 30)
  plt.xlabel('频率 (Hz)')
  plt.ylabel('功率谱')
  plt.show()
  

• 确认频率轴正确性：确保np.fft.fftfreq的参数正确，d为采样间隔（即1/采样频率），你的代码中1/62.5是正确的，但需确认数据长度与采样频率匹配，避免频率轴错位。

内容的提问来源于stack exchange，提问作者Renan Saraiva dos Santos