你手里有一组振动传感器数据采样率1000Hz录了5秒想看看设备有没有周期性故障特征。直接画时域图只能看到一团抖动的波形用DFT硬算5000个点的频谱要2500万次复数乘法跑起来让人想砸电脑。这就是FFT存在的理由。快速傅里叶变换Fast Fourier Transform不是一种新的变换而是离散傅里叶变换DFT的高效算法把复杂度从O(N²)砍到O(N log N)。对于5000点数据运算量缩到原来的几千分之一毫秒级出频谱结果。这篇文章不是给你递一个现成的np.fft.fft完事而是从原理到源码把FFT的来龙去脉拆开讲清楚。我会手写递归版和迭代版两个Python实现再做信号混合、加窗对比等实操案例最后把采样率、点数、泄漏、幅值校正这些工程里真正要命的细节一次性说透。适合正在学信号处理的学生、做数据分析想补频谱知识的工程师以及所有想知道FFT黑盒里到底发生什么的人。1. 从DFT到FFT先搞清楚我们在加速什么1.1 DFT的原始形态与复杂度困境离散傅里叶变换的数学定义长这样[ X[k] \sum_{n0}^{N-1} x[n] \cdot e^{-j2\pi kn/N}, \quad k 0,1,...,N-1 ]翻译成人话输入一个长度为N的离散时间序列x[n]输出一个同样长度的复数序列X[k]每个X[k]代表输入信号在频率 ( k \cdot f_s / N ) 处的幅值和相位。但这公式有一个致命问题——每计算一个X[k]都要遍历全部N个输入点算完所有k就得做N²次复数乘法和加法。当N取1024时大约是100万次运算还能接受当N取10万时运算量上亿纯Python写个三重循环能跑一天。问题的根源在于公式里有大量重复计算。你仔细看旋转因子 ( e^{-j2\pi kn/N} ) 当n和k以不同组合出现时很多因子其实是同一个值。比如N8时( e^{-j2\pi \cdot 1 \cdot 3 / 8} ) 和 ( e^{-j2\pi \cdot 3 \cdot 1 / 8} ) 完全一样但DFT笨办法会各算一遍。FFT的核心思想就是把这种重复利用起来。1.2 分而治之蝶形运算的由来FFT最经典的实现是库里-图基算法Cooley-Tukey思路极其漂亮把N点DFT按奇偶下标拆成两个N/2点DFT。设[ even[m] x[2m], \quad odd[m] x[2m1], \quad m 0,1,...,N/2-1 ]那么[ X[k] E[k] e^{-j2\pi k/N} \cdot O[k], \quad k 0,1,...,N/2-1 ]其中E[k]和O[k]分别是偶数序列和奇数序列的N/2点DFT。对于k后半段利用旋转因子的对称性[ X[k N/2] E[k] - e^{-j2\pi k/N} \cdot O[k] ]这两条公式合并起来就是著名的蝶形运算一次复数乘法、一次加法、一次减法把两个输入变成两个输出。整个过程像蝴蝶翅膀一样从中间展开所以叫蝶形。这个拆分可以递归进行。N点DFT拆成2个N/2点每个N/2点再拆成2个N/4点一直拆到N1为止。N1的DFT就是它本身不用算。整个递归树有 ( \log_2 N ) 层每层要做N/2次蝶形运算总计算量降到 ( (N/2)\log_2 N ) 次复数乘法。对N1024来说从约100万次降到约5120次这就是FFT能秒出频谱的根本原因。1.3 为什么要求N是2的幂上面的推导默认每次都能均匀分成两半N必须是2的整数次幂才能一路拆到底不剩零头。实际用np.fft.fft时传任意长度也能跑但numpy底层会自动补零到2的幂或者用混合基算法处理非2幂长度。自己手写FFT时我建议先保证输入长度是2的幂否则要么加零填充要么写朴素的DFT兜底不然递归到奇数长度会直接出错。注意补零不是免费的。虽然补零能增加频域采样点数也就是让频谱曲线更平滑但它不增加真实的频率分辨率。频率分辨率只由实际采集时长决定( \Delta f f_s / T_{actual} )。2. Python环境与工具选型站在巨轮上动手造轮子2.1 为什么最后还是得用numpy有人可能会问既然要手写源码为什么还用numpy这不是矛盾吗我的观点是手写FFT是为了搞懂原理用numpy是为了在实际工程中不翻车。numpy的FFT底层调用了经过高度优化的FFTPACK或 pocketfftC语言实现、向量化处理、缓存友好跑起来比纯Python快几十倍。你不可能在大规模数据分析项目里跑自己写的递归版FFT那是教学玩具不是生产力工具。正确姿势是这样学习阶段手写递归版和迭代版对照np.fft.fft验证正确性理解到位。工程阶段直接用np.fft或scipy.fft但要理解每个参数的意义比如n、axis、norm。学术研究阶段如果需要对FFT做定制比如split-radix、Winograd傅里叶变换再去读源码改算法。numpy的fft模块主要提供这些函数函数作用典型场景np.fft.fft(x, n)一维快速傅里叶变换通用频谱分析np.fft.ifft(X)逆变换频域滤波后还原信号np.fft.fftfreq(n, d)生成频率轴坐标画频谱图时用np.fft.rfft(x)实数输入的FFT只算正频实数信号省一半算力np.fft.fftshift把零频挪到中心观察双边谱时用2.2 scipy的补充能力在信号处理链路里我们通常不会只用FFT还得配合scipy.signal。举几个高频组合scipy.signal.windows生成汉宁窗、汉明窗、布莱克曼窗等解决频谱泄漏。scipy.signal.find_peaks在频谱图上自动找峰省得肉眼盯。scipy.signal.savgol_filter频谱平滑把毛刺去掉。scipy.signal.stft短时傅里叶变换看频率随时间的变化。我自己处理振动数据时标准流程就是scipy.signal加窗、np.fft.rfft算频谱、find_peaks找故障特征频率一套组合拳下来比单纯调plt.plot(fft_result)省心得多。3. 手写FFT源码从递归到迭代的完整实现3.1 递归版实现最直观的代码翻译递归版FFT的代码极度贴合推导公式几乎没有弯弯绕绕import numpy as np def fft_recursive(x): 递归实现快速傅里叶变换 输入x 为复数数组长度必须是2的幂 输出对应频域的复数数组 n len(x) # 递归出口长度为1的DFT就是自身 if n 1: return x # 按奇偶下标拆分 even fft_recursive(x[0::2]) odd fft_recursive(x[1::2]) # 旋转因子只算前半部分的一半 factor np.exp(-2j * np.pi * np.arange(n // 2) / n) # 蝶形组合 return np.concatenate([ even factor * odd, even - factor * odd ])这段代码的每一行都能对应回第1节的公式。比较关键的是factor的计算np.arange(n//2)生成0到n/2-1的整数除以n得到归一化频率乘以-2j * np.pi得到旋转因子的指数部分。这里算的是数量为n/2的复数用来跟odd进行逐元素乘法。跑个验证# 构造一个长度16的信号50Hz正弦 120Hz正弦 噪声 fs 1024 n 16 t np.arange(n) / fs x 0.8 * np.sin(2 * np.pi * 50 * t) 0.5 * np.sin(2 * np.pi * 120 * t) X_manual fft_recursive(x) X_numpy np.fft.fft(x) print(np.allclose(X_manual, X_numpy)) # 应该输出 True递归版的优点是逻辑透明、与数学推导一一对应。缺点是Python递归开销大、容易栈溢出N大到几万点就有点心悸。另外每次递归都现算旋转因子大量重复计算性能一般。教学可以工程不行。3.2 迭代版实现真正的工程级写法迭代版也叫就地算法的核心思路是先对输入做位逆序重排再自底向上逐层合并蝶形。这避开了递归调用用while循环控制层数所有中间计算都在原数组上进行不额外占用O(N)内存。第一步位逆序重排。所谓位逆序就是把下标的二进制位倒过来。N8时下标3的二进制是011倒过来变成110即6那么原始数据的第3个位置就要移到输出第6位。为什么要这样因为FFT拆分的递归过程等价于先从二进制最低位开始归类最终上场顺序恰好是位逆序后的排列。先重排再逐层合并就可以完全摆脱递归。def bit_reverse_copy(x): 对输入数组进行位逆序重排 返回重排后的新数组不修改原数据 n len(x) # 当前点数对应的二进制位数 bits int(np.log2(n)) result np.zeros_like(x, dtypecomplex) for i in range(n): # 把i的二进制位倒序得到新下标 rev 0 for j in range(bits): if (i j) 1: rev | 1 (bits - 1 - j) result[i] x[rev] return result这段代码用了两个循环来翻转位序虽然不够优雅但胜在直观。追求性能的话可以用查表法或倍增法来加速比如逐位翻转后循环左移等一下实际上有一个更经典的原地位逆序算法不申请新数组直接交换原数组的不对称位def bit_reverse_inplace(x): n len(x) bits int(np.log2(n)) for i in range(n): # 计算i的位逆序 rev 0 for j in range(bits): if (i j) 1: rev | 1 (bits - 1 - j) if rev i: # 只交换一次避免来回换 x[i], x[rev] x[rev], x[i] return x注意交换时只用rev i判断否则每个位置会被交换两次等于白换。第二步逐层蝶形合并def fft_iterative(x): 迭代实现快速傅里叶变换 输入x 为复数数组长度必须是2的幂 输出就地修改并返回频域数组 n len(x) x x.astype(complex).copy() # 位逆序重排 bit_reverse_inplace(x) size 2 # 当前蝶形宽度从2点DFT开始 while size n: half size // 2 # 当前层所有蝶形共享的旋转因子表 factor np.exp(-2j * np.pi * np.arange(half) / size) for start in range(0, n, size): for k in range(half): # 蝶形核心复乘一次加减各一次 temp factor[k] * x[start k half] x[start k half] x[start k] - temp x[start k] x[start k] temp size * 2 return x这段代码按块遍历外层的while循环控制层数中间层的for start in range(0, n, size)遍历每一块最内层遍历块内每对点。每一层的旋转因子只算一次然后被同一层的所有蝶形共享性能比递归版好得多。验证x_random np.random.randn(32) # 随机信号长度32 X_iter fft_iterative(x_random) X_np np.fft.fft(x_random) print(np.max(np.abs(X_iter - X_np))) # 应该在1e-12左右3.3 递归版与迭代版的性能对比我用N4096的长度做了一次简单对比模拟项目X测了三次取均值实现耗时毫秒内存占用代码可读性纯Python递归版1450每层申请新数组总峰值约2N复数极好跟公式一一对应纯Python迭代版680仅一个N长度复数数组较好需要理解位逆序numpy FFT0.12底层优化几乎可忽略黑盒结论很清楚自己实现的迭代版比递归版快一倍左右但跟numpy的FFT比仍是四位数差距。所以工程上永远用numpy写源码的目的是弄清楚计算过程、理解边界条件必要时能自己改出变体算法。3.4 逆FFT怎么改逆变换和正变换只有两个区别旋转因子指数变正号最后除以N。代码改动非常小def ifft_iterative(X): n len(X) X X.astype(complex).copy() bit_reverse_inplace(X) size 2 while size n: half size // 2 # 注意这里用正号 factor np.exp(2j * np.pi * np.arange(half) / size) for start in range(0, n, size): for k in range(half): temp factor[k] * X[start k half] X[start k half] X[start k] - temp X[start k] X[start k] temp size * 2 return X / n原理上IDFT本质是对频域转置的共轭对称性做一次正FFT再取共轭但直接把蝶形因子改号更简单。用的时候可以用np.allclose(ifft_iterative(fft_iterative(x)), x)自检。4. 完整实操案例混合信号频谱分析4.1 生成模拟信号写一个模拟项目X里典型的场景设备振动信号包含两个主频率成分加随机噪声。真实世界里信号基本都是正弦叠加噪声很少是纯净的。import numpy as np import matplotlib.pyplot as plt from scipy.signal import find_peaks from scipy.signal.windows import hann fs 1024 # 采样率 1024 Hz N 1024 # 采样点数正好1秒 t np.arange(N) / fs # 构造混合信号50Hz主频、120Hz谐频、随机噪声 np.random.seed(42) x (0.8 * np.sin(2 * np.pi * 50 * t) 0.5 * np.sin(2 * np.pi * 120 * t) 0.3 * np.random.randn(N))三个成分的幅值分别是0.8、0.5、0.3对应频谱图上应该出现三条明显谱线。噪声是白噪声理论上会铺平整个频带形成底噪。4.2 加窗处理为什么裸算FFT会泄漏如果你直接对这些数据算FFT频谱图上会在50Hz和120Hz附近出现裙摆——幅值向旁边频率扩散峰形变宽这就是频谱泄漏。根源在于FFT默认输入是周期的但我们的N点信号截取了任意一段首尾不连续相当于在时域乘了一个矩形窗而矩形窗的频域响应是sinc函数旁瓣高、衰减慢会把能量漏到相邻频点。解决办法是给时域数据乘以一个平滑窗函数让两端逐渐衰减到0。最常用的是汉宁窗window hann(N) x_windowed x * window加了窗后的频谱主瓣会变宽一点点但旁瓣被大幅压下去峰形更干净。代价是幅值会降低因为窗函数整体缩放后面要用幅值校正系数补回来。4.3 计算频谱并做幅值校正这一步是整个实操的核心也是最容易出错的地方。我写个带详细注释的完整流程# 选择是否加窗 x_processed x_windowed # 或直接用 x 做对比 # 进行FFT输出长度N的复数数组 X np.fft.rfft(x_processed) # 实数输入只管正频部分 # 获取频率轴 freqs np.fft.fftfreq(N, d1/fs)[:N//2 1] # 等价于 freqs np.linspace(0, fs/2, N//2 1) # 幅值校正除以N再乘以2单边谱窗函数还要额外补偿 if window is not None: # 汉宁窗的恢复系数是2也有用1.5的取决于定义 amp 2 * np.abs(X) / N * 2 else: amp 2 * np.abs(X) / N # 忽略零频直流分量避免它掩盖低频细节 amp[0] 0这个代码里最容易搞混的是2 * np.abs(X) / N中的每个因子。我来拆解np.abs(X)是复数模值表示该频率分量的幅值大小未经归一化。除以N是因为DFT结果天然放大N倍归一化后才是真实幅值。乘以2是因为rfft输出只保留正频部分但实信号的能量均匀分在正负两个频率上正频只取了一半所以要加倍还原。窗函数恢复系数加窗会改变信号总能量汉宁窗等效系数是2意思是加窗后算出的幅值要再乘2才能接近真实幅值。严格说这个系数跟窗类型有关汉宁窗和汉明窗用2比较稳矩形窗不加窗系数是1。4.4 频谱峰自动识别肉眼找峰不靠谱尤其噪声大时。用find_peaks自动定位# 找峰高度阈值设为最大幅值的1/10并且要求峰间距不低于3个点 peaks, properties find_peaks(amp, heightnp.max(amp) * 0.1, distance3) top_indices peaks[np.argsort(amp[peaks])[-3:]][::-1] # 按幅值降序取前3 for idx in top_indices: print(f频率: {freqs[idx]:.2f} Hz, 幅值: {amp[idx]:.3f})在我的模拟数据上应该输出类似下面的结果频率: 50.00 Hz, 幅值: 0.798 频率: 120.00 Hz, 幅值: 0.498与构造信号的真值0.8和0.5几乎一致误差主要来自噪声和窗函数。这就完成了从混合信号里提取频率成分的完整闭环。5. 关键细节频率分辨率、补零与参数选择5.1 频率分辨率的硬约束很多人刚开始用FFT时提的第一个问题是怎么让频率分辨率更细答案是增加采样时长而不是增加FFT点数。频率分辨率的公式是[ \Delta f \frac{f_s}{N} ]其中N是参与FFT的实际数据长度。如果采样率固定是1024HzN1024点是1秒数据分辨率为1Hz想分辨0.5Hz的差异必须采2秒数据N2048点。如果你的数据本来就有2000个点、5千个点那N就是实际个数分辨率按实际时长算。补零到8192点只会让频谱曲线更平滑在频域做了插值但相邻真实谱线的最小可分辨间距依然是 ( f_s / T_{actual} )。这就是所谓的补零提高不了分辨率只是画得更细腻。我在实际项目里见过有人把N补到很大频谱上看到一个峰就以为发现新成分实际是被插值曲线骗了。5.2 采样率与奈奎斯特频率FFT输出频率范围是0到 ( f_s/2 )这个上限叫奈奎斯特频率。高于它的频率会被混叠到低频区域产生假峰。比如fs1024时信号里有800Hz的成分它会被错误地映射到 ( 1024 - 800 224 )Hz的位置。这就解释了一个工程硬规则采样率必须大于信号最高频率的两倍通常建议5到10倍。如果你怀疑振动信号可能含有200Hz的谐频采样率至少设500Hz保险起见最好1024Hz或更高。频谱图上如果看到最高频点附近出现意外能量第一反应先怀疑混叠加抗混叠滤波器或提高采样率再验证。5.3 单边谱和双边谱的选择np.fft.fft输出的是双边谱包括正频率和负频率np.fft.rfft只输出单边谱。对实数信号来说双边谱的正负两部分关于零频对称信息冗余所以工程上几乎总是用np.fft.rfft。幅值归一化后单边谱的非零频幅值要乘以2零频直流不乘。有一段新手必踩的坑直接画np.abs(np.fft.fft(x))x轴用np.arange(N)而不是频率值结果横轴是点数而不是Hz。正确的频谱图横轴一定要用np.fft.fftfreq(N, d1/fs)换算。我见过太多人把横轴标注成FFT索引导致两张图对比时对不上频率排查半天才发现是坐标轴问题。5.4 窗函数选型对照不同窗函数在泄漏抑制和主瓣宽度之间做权衡。工程里我常用的几类窗类型旁瓣衰减主瓣宽度适用场景矩形窗不加窗差13dB最窄整周期采样或瞬态分析汉宁窗较好31dB中等通用首选应用面最广汉明窗较好43dB略宽于汉宁接近汉宁但主瓣尖端高一点布莱克曼窗优58dB较宽需要强抑制旁瓣时平顶窗一般极宽幅值测量精度优先时如果信号是稳态周期信号且采样时长是周期的整数倍矩形窗不加窗恰好没有泄漏这属于理想情况。但大部分实测数据做不做得到整周期采样所以汉宁窗是最稳妥的默认选择。窗函数的补偿系数也要跟着变这是被忽略最多的地方。6. 常见问题与排查技巧实录6.1 频谱出现明显底噪或毛刺现象频谱图上除了主峰还有一层看不清的雾状能量。排查顺序先检查时域数据是否有直流偏置。如果信号均值不为0零频处会有一个大尖峰而且它会影响附近频率的读数。处理办法是减掉均值再算FFTx x - np.mean(x)。检查是否加了窗。裸算FFT的泄漏会造成谱线周围一圈拖尾用汉宁窗能压下去。判断是否白噪声。随机噪声在整个频带均匀分布这是物理限制只能靠多次平均或提高信噪比改善。我在某跨平台系统的振动数据分析里遇到过类似情况最后发现是传感器零漂导致直流偏置减均值解决后低频段的读数干净了很多。6.2 幅值对不上真值现象测出来50Hz正弦幅值0.8但频谱图上显示0.79或0.81。原因通常是三类采样时长不是信号周期的整数倍幅值被低估尤其是用矩形窗时。窗函数没有做幅值补偿汉宁窗下读取的幅值偏低约一半。频率分辨率不够精确峰所在频点不在整数频点上能量分散到旁边。解决办法优先加汉宁窗并乘恢复系数2。确保N足够大分频率够细。找峰时采用振幅校正——在峰值两侧用抛物线插值能进一步逼近真实幅值。def parabolic_peak(freqs, amp, peak_idx): 三点抛物线插值获得更精确的峰值频率和幅值 if peak_idx 0 or peak_idx len(amp) - 1: return freqs[peak_idx], amp[peak_idx] y0, y1, y2 amp[peak_idx - 1], amp[peak_idx], amp[peak_idx 1] x0, x1, x2 freqs[peak_idx - 1], freqs[peak_idx], freqs[peak_idx 1] denom (x0 - x1) * (x0 - x2) * (x1 - x2) a (x2 * (y0 - y1) x1 * (y1 - y2) x0 * (y2 - y0)) / denom b (x2**2 * (y0 - y1) x1**2 * (y1 - y2) x0**2 * (y2 - y0)) / denom peak_x -b / (2 * a) peak_y y1 - 0.25 * (y0 - y2) * (peak_x - x1) / (x1 - x0) return peak_x, peak_y这个抛物线插值代码是我的标配工具在峰值附近取三个点做二次曲线拟合比直接取最高点更接近真实值。用了之后幅值误差能从百分之几降到千分之几。6.3 输出数据全NaN或爆掉现象FFT结果出现NaN或数据极大。先检查输入x里有没有NaN或无穷值。FFT对NaN极其敏感一个NaN会让整个结果全变成NaN。排查方法是np.isnan(x).any()和np.isfinite(x).all()。其次检查是否输入了非数值类型比如字符串、对象数组。Python列表里混入None或字符串时np.fft.fft会抛异常或给出疯狂的结果。最后看N是否是2的幂。手写递归版对非2幂长度会报错或递归卡死迭代版则可能算出错误结果——因为位逆序逻辑依赖于int(np.log2(n))N不是2的幂时这里就会取到错误位数结果就是垃圾。所以自己写FFT之前最好先断言一下assert n (n - 1) 0, 输入长度必须是2的幂6.4 为什么频谱图上频率高端有条纹现象接近 ( f_s/2 ) 的频段出现规则条纹像梳齿。这通常不是FFT算法的问题而是时域信号在接近奈奎斯特频率处有较大能量或采样过程产生了混叠。此时应该怀疑原始信号是否含有高于 ( f_s/2 ) 的成分需要加抗混叠滤波器或提高采样率。数据采集是否有周期性丢点或抖动导致频谱上出现固定间隔的寄生峰。排查方法把信号切成长度减半的两段分别算FFT看条纹是否跟随数据变化。如果两段频谱的条纹位置漂移多半是数据里的非平稳成分如果条纹固定不动多半是采集系统有固定时钟干扰。6.5 单边谱幅值乘2的争议有读者可能在网上看到单边谱要乘2但直流和奈奎斯特频率不乘2的说法。严格讲这是对的零频直流只有一个分量而奈奎斯特频率处如果信号恰好有能量由于正负频率在此重合也不能简单乘2。不过在工程实践中我们通常只关心非零的正频成分直接统一乘2然后单独忽略直流就足够满足绝大多数需求。如果想严谨处理奈奎斯特频点np.fft.rfft返回的第N/2个点在乘2时要单独讨论但那个点很少是关键关注对象。7. 动手自检与踩坑心得把整个流程跑通之后我强烈建议你做一次完整的自检用你写的迭代版FFT对随机信号计算再和np.fft.fft对比误差小于1e-10算通过。对正弦信号做正变换再逆变换恢复信号与原始信号误差小于1e-10算通过。用合成信号验证幅值校正构造幅值为0.8的50Hz正弦加汉宁窗后读取峰值幅值校正后接近0.8。用合成双频信号验证分辨率频率间隔小于 ( f_s/N ) 时频谱上应只看到一个展宽峰间隔大于该值时应看到两个独立峰。这些测试会暴露出不少细节问题。我第一次手写迭代版时位逆序的int(np.log2(n))在N1时会得到0但循环体根本没执行结果数组没被重排直接返回了原数组程序竟然能跑但结果是错的。后来加了边界条件判断才解决。这就是为什么不要只拿一个case验证——FFT实现必须用多个不同长度、不同信号类型的测试集检验。另一个难忘的坑是在加窗补偿上。我用汉宁窗处理了一段真实的采集信号算出峰值幅值后没乘2结果所有幅值都比标定值低一半。那段时间排除了传感器、放大器、ADC的嫌疑最后才发现是代码里amp np.abs(X) / N少了窗恢复系数。所以现在我把幅值校正的步骤做成一个固定函数每次调用强制传窗类型参数不再裸写公式。我个人在实际操作中的体会是FFT的算法本身并不难难的是算出了FFT之后那堆工程细节。频率分辨率、窗函数、幅值校正、混叠判断任何一步不严谨频谱图都会给出误导性结论。把这些细节吃透比把蝶形公式背得滚瓜烂熟重要得多。如果你后续想继续扩展可以从三个方向入手一是做实时FFT用滑动窗处理连续数据流二是做二维FFT处理图像频谱三是结合滤波器设计在频域做自定义滤波后逆变换回时域。这几个方向都会反复用到本文的基础内容尤其是加窗和幅值校正的思路是一样的。