Python中快速傅里叶变换（FFT）实现及脚本正确性验证

你的FFT脚本存在几处关键错误，下面逐一说明并给出修正方案：

1. 样本数量N的获取错误 #

N = np.shape(x) 会返回一个元组（比如这里x是长度40的数组，返回(40,)），后续计算时间步长时会触发类型错误。应该用整数类型的样本数，改为：

N = len(x)  # 或者 N = x.size

2. 时间步长计算逻辑错误 #

从你的x数组来看，相邻元素间隔是5（比如-25到-20差5），假设这个间隔对应的实际时间是5 * 10e-12秒（你原代码里的10e-12应为时间单位系数），正确的时间步长Δt应该是相邻x值的差乘以时间系数：

dt = (x[1] - x[0]) * 10e-12  # 得到5e-11秒

原代码里(200/N)*10e-12的逻辑不成立：200并不是x的总范围（x总范围是170 - (-25) = 195），且时间步长是单个采样间隔，不是总范围除以样本数。

3. fftfreq参数错误 #

fftfreq的第二个参数是采样间隔Δt（每个样本之间的时间差），而非1/T。正确调用方式是传入计算好的dt：

xf = np.fft.fftfreq(N, d=dt)[:N//2]

4. 细节优化 #

• 需明确fft和fftfreq是numpy的函数，加上np.fft.前缀避免命名冲突；
• 纵坐标标签f(t)不准确，建议改为Amplitude或|Y(f)|，因为FFT结果是频域幅值。

修正后的完整脚本 #

import numpy as np
import matplotlib.pyplot as plt

# 原始数据
y = np.array([9.88706879e-05, -1.80853647e-05, 2.42572582e-05, 
              1.12215205e-04, 1.32105126e-04, 1.13424614e-05, -1.58262175e-04, 
              -2.62013276e-04, -2.58070932e-04, -1.53975865e-04, -8.19357356e-05, 
              -1.55734157e-04, -2.90791620e-04, -3.70294471e-04, -3.46855608e-04, 
              -2.23495910e-04, -1.35441615e-04, -2.11411786e-04, -4.21891416e-04, 
              -6.77753516e-04, -8.09657243e-04, -6.97948704e-04, -5.01935670e-04, 
              -4.20075723e-04, -5.28464040e-04, -8.14942203e-04, -1.03669983e-03, 
              -9.76604755e-04, -7.50889655e-04, -5.34882634e-04, -4.06928662e-04, 
              -3.96093220e-04, -4.31306957e-04, -4.25399844e-04, -3.26933980e-04, 
              -1.32440493e-04, 5.40550849e-06, -4.87299567e-05, -2.04672372e-04, 
              -3.15870097e-04])

x = np.array([-25, -20, -15, -10,  -5,   0,   5,  10,  15,  20,  25,  30, 35,
              40,  45,  50,  55,  60,  65,  70,  75,  80,  85,  90,  95, 100,
              105, 110, 115, 120, 125, 130, 135, 140, 145, 150, 155, 160, 165,
              170])

# FFT计算
N = len(x)
dt = (x[1] - x[0]) * 10e-12  # 采样间隔，可根据实际物理意义调整
xf = np.fft.fftfreq(N, d=dt)[:N//2]
yf = np.fft.fft(y)

# 绘图
plt.figure()
plt.plot(xf, 2.0/N * np.abs(yf[0:N//2]))
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
plt.title('FFT Result')
plt.show()

额外说明 #

• 若x数组的单位不是对应10e-12秒，需根据实际物理意义调整dt的计算；
• 2.0/N * np.abs(yf[0:N//2])这部分是正确的，它将FFT结果幅值归一化，仅保留正频率部分的有效值。

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