如何用简单示例验证快速傅里叶变换(FFT)结果符合理论预期

• 第一个问题：FFT计算函数与频率序列生成函数不匹配

  你调用的scipy.fftpack.fft()返回的是包含正负频率的完整双边谱结果，长度等于采样点数N；但scipy.fftpack.rfftfreq()是为实信号快速傅里叶变换rfft()设计的，仅返回正半部分的频率序列，长度为N//2 + 1。二者长度不一致，绘图时会自动截断FFT结果，导致频率和幅值对应关系完全错误。

• 第二个问题：单边谱幅值缩放规则错误

  实信号的能量均匀分布在正负频率轴上，如果你要得到和时域幅值对应的单边幅频谱，除了乘以采样间隔做能量归一化之外，还需要对正频率部分的幅值额外乘以2的系数。

提供两种可直接运行的正确实现，二选一即可：

方案1：使用实信号专用FFT接口rfft（推荐，计算量更小） #

import numpy as np
from scipy.fftpack import rfft, rfftfreq
import matplotlib.pyplot as plt

N = 100
# 加endpoint=False避免首尾采样点重复导致的频谱泄漏
x = np.linspace(0.0, 1, N, endpoint=False)
y = np.sin(5 * 2.0*np.pi*x) + 0.5*np.sin(2 * 2.0*np.pi*x)
# 用rfft对应rfftfreq，乘以2做单边谱幅值校正
yf = rfft(y) * (x[1]-x[0]) * 2
freq = rfftfreq(N, x[1]-x[0])

plt.plot(freq, np.abs(yf))
plt.xlabel('频率(Hz)')
plt.ylabel('幅值')
plt.show()

方案2：保留原fft接口，匹配对应频率序列 #

import numpy as np
from scipy.fftpack import fft, fftfreq
import matplotlib.pyplot as plt

N = 100
x = np.linspace(0.0, 1, N, endpoint=False)
y = np.sin(5 * 2.0*np.pi*x) + 0.5*np.sin(2 * 2.0*np.pi*x)
yf = fft(y)*(x[1]-x[0])
freq = fftfreq(N, x[1]-x[0])
# 只取正频率部分，同时幅值乘以2校正
pos_mask = freq >= 0
pos_freq = freq[pos_mask]
pos_yf = np.abs(yf[pos_mask]) * 2

plt.plot(pos_freq, pos_yf)
plt.xlabel('频率(Hz)')
plt.ylabel('幅值')
plt.show()

两种方案运行后都可以在2Hz处看到幅值0.5的峰值，5Hz处看到幅值1的峰值，和理论预期完全一致。

内容的提问来源于stack exchange，提问作者Jonny Joker