4.1 快速傅里叶变换

4.1 快速傅里叶变换

本节摘要：快速傅里叶变换（FFT）是把一段时域采样信号拆解成频率成分的高效算法，scipy.fft 是 SciPy 官方实现的 FFT 模块。它相对 numpy.fft 多了三样东西：可切换的计算后端、完整的实数变换家族（rfft、hfft 等）、以及配套的频率轴与移频工具。本节先跑一段能用的最小频谱分析代码，再讲采样定理的直觉、窗函数与频谱泄漏的来源，最后说清后端切换什么时候值得做。

阅读收获

阅读完本节，你应当能够：

1. 用 rfft 与 rfftfreq 完成一段信号的频谱分析，正确读出峰值频率与幅值；
2. 说清 scipy.fft 相对 numpy.fft 的三个优势，并解释后端切换的含义；
3. 用采样定理判断一组采样参数能否无混叠地测量目标频率；
4. 解释窗函数为什么能抑制频谱泄漏，以及它带来的幅值代价；
5. 区分频率分辨率与补零的关系，避免把光滑误当成精确。

一、问题与直觉

场景是设备振动监测。传感器以 1000 赫兹采样，录了 1 秒数据。我们知道信号里应该有两个振动源，一个大约 50 赫兹，一个大约 120 赫兹，但两条正弦波叠在一起，波形缠成一团，肉眼根本分不清谁是谁。这时候唯一靠谱的办法是把信号拆成"频率成分"——这正是 FFT 干的事。

先上能跑的最小代码：

import numpy as np
from scipy.fft import rfft, rfftfreq

fs = 1000                                  # 采样率，单位赫兹
t = np.arange(0, 1, 1/fs)                  # 1 秒时间轴
x = 0.7*np.sin(2*np.pi*50*t) + 0.3*np.sin(2*np.pi*120*t)

X = rfft(x)
freq = rfftfreq(len(x), 1/fs)
idx = np.argmax(np.abs(X[1:])) + 1         # 跳过直流分量再取最大
print("主频:", freq[idx], "赫兹")
​

运行后输出主频 50 赫兹，把 idx 附近几个峰值打出来还能看到 120 赫兹。这就是频谱分析的最小闭环：rfft 把信号变到频域，rfftfreq 给出每个频点对应的物理频率，取幅值最大的位置就是主频。在你自己的数据上跑通这一步，比读十页理论都管用——先把工具用起来，再回头理解它在算什么。

你可能注意到 NumPy 里也有一个 numpy.fft，两边函数名几乎一样。历史原因：numpy.fft 是早年从 SciPy 的 fftpack 迁移过去的，功能停留在旧时代；scipy.fft 是后来重写的现代版本，统一了接口、补齐了实数变换与后端机制。所以用 SciPy 做科学计算时，优先从 scipy.fft 导入，别在两边各写一套。

二、核心原理

2.1 FFT 在算什么

傅里叶变换的朴素理解是匹配：把信号与一组不同频率的参考正弦波逐一做内积，内积大的频率就是信号里真实存在的成分。直接按这个定义算，N 个点要做 N 次内积、每次 N 次乘加，总复杂度是 N 的平方；FFT 利用旋转因子的对称性，把复杂度压到 N 乘 logN。数据量一大，这个差距是数量级的——1024 个点，直接算是一百多万次操作，FFT 只要一万次出头。这也是为什么"快速"两个字值得写进算法名字里。

FFT 的输出是复数，实部与虚部合起来同时携带幅度和相位信息：模长代表该频率的强度，辐角代表该频率的初始相位。多数工程分析只看幅值谱，但相位谱在系统辨识、解调任务里同样重要。另外，正逆变换是一对组合，用 ifft 把频谱变回时域，误差应该在 1e-15 量级——这个自检几乎不花成本，怀疑自己的频谱处理出了问题时，先跑一遍。

2.2 scipy.fft 相对 numpy.fft 的优势

两个模块都能算 FFT，数值结果几乎一致，毕竟底层数学是一样的。但 scipy.fft 是面向科学计算的正主，优势集中在三处。

一是后端切换。scipy.fft.set_backend 允许把底层计算换到其他 FFT 实现，比如 FFTW。默认后端 pocketfft 已经很快，一般不必折腾；但当你面对超大数组、或者团队规定必须统一用某套数值库时，这个开关是 numpy.fft 给不了的。二是实数变换家族完整。时域采样信号几乎都是实数，实数序列的频谱有共轭对称性，后半段全是前半段的镜像，rfft 只算一半频点，输出长度是 N/2 加 1，计算量接近减半。scipy.fft 还提供 hfft、ihfft 以及 rfft2、rfftn 等多维版本，numpy.fft 在这块是残缺的。三是配套工具齐全。fftfreq、rfftfreq 给出频率轴，fftshift、ifftshift 负责把零频移到中心，next_fast_len 帮忙找快速长度，fftconvolve 直接做频域卷积——这些函数在频谱分析里天天用，散落在 numpy 各处反而麻烦。

2.3 采样定理的直觉

采样率 fs 决定你能看到的最高频率：奈奎斯特频率是 fs 的一半，超过它的成分会伪装成低频信号混进来，这就是混叠。直觉版本：采样就像给信号拍快照，一个周期内至少拍两张照片，才能分清它在往哪边转。你隔一天才给月亮拍一张照片，当然看不出它绕地球的轨迹。所以 44.1 千赫的音频采样率能覆盖 20 千赫的人耳上限，是算了账的。反过来，如果你怀疑信号里有 500 赫兹的成分，采样率至少得 1000 赫兹，否则那个成分会以假身份出现在频谱里，误导后续所有分析。举个具体的数：采样率 1000 赫兹时，700 赫兹的成分会伪装成 300 赫兹出现，因为两者在采样点上无法区分。这提醒我们，采样的第一步就要确认信号的最高频率，必要时在采样前加抗混叠低通——硬件上常有这层保护，软件里却经常被忽略。

2.4 频率分辨率

频率分辨率等于 fs 除以 N，也等于采样总时长的倒数。想分清 50 赫兹和 50.5 赫兹，至少需要 2 秒的数据，因为 1 秒只能分辨 1 赫兹的间隔。补零可以把频谱曲线插得更密、看起来更光滑，但它不增加真实信息，分辨率不会因此变好——这是初学者最容易踩的误解：图变细了，不代表能分开两个挨着的频率。判断能不能分辨，只看采样时长，不看画图参数。

2.5 窗函数与频谱泄漏

FFT 假设信号是周期延拓的：它把截取的这一段当作无限循环。当截取长度不是信号周期的整数倍，边界处就会产生"接缝"，能量从真实频率漏向旁瓣，在频谱上拖出一条条逐渐衰减的裙带，这叫频谱泄漏。泄漏最麻烦的地方是会把小幅值的真实频率淹没在旁瓣里。

解决办法是加窗：把数据两端乘到接近零，让"接缝"消失。汉宁窗、汉明窗是常用选择，scipy.signal.windows 里有一整排。代价有两个：主瓣变宽，两个很近的频率更难分开；峰值幅度被压低，做定量分析时要按窗函数的相干增益修正回来。反过来，如果信号恰好是整周期截断（同步采样、锁相触发），不加窗反而更准，加了窗纯属多此一举。

加窗在代码里只有一行：x 乘以 signal.windows.hann(N)，把窗函数与信号逐点相乘，再做 rfft。注意窗函数同时压低数据两端，幅度谱的峰值会被削掉约一半，定量分析时要么乘回修正系数，要么干脆只做相对比较。另外，窗函数还有一个延伸用途：把长数据切成多段、每段加窗后再平均频谱，就是 Welch 功率谱估计的基本思路，第 5 章会用到。

三、工程实践要点

3.1 一组稳定的频谱分析模板

固定顺序：造时间轴 → 按需加窗 → rfft → rfftfreq → 取模。要不要加窗，取决于你的数据是不是整周期截断。实测数据几乎都不是，所以默认加汉宁窗是稳妥的；合成信号测试时，保证截断长度是周期整数倍，可以不加窗直接看干净的谱线。

3.2 fft 与 rfft 的选择

判断标准一句话：输入是实数组，就用 rfft，别犹豫。它不仅省一半计算，输出的频率轴还短一半，画图、找峰值都更清爽。

⚠️ 常见坑：画幅值谱时忘了修正。实数信号的单边谱，幅值要乘 2 再除以 N，直流分量除外。不修正的话，读出来的幅值跟时域波形对不上，定量分析直接翻车；修正时又把直流分量也乘了 2，零频处出现一个莫名其妙的凸起。
💡 关键直觉：频率轴永远用 rfftfreq 或 fftfreq 生成，不要手写 linspace 数格子。差一个点是常有的事，而频率轴错了，整个频谱图都白画。

3.3 参数细节：n、norm 与 workers

fft 的 n 参数可以补零或截断：补零能加密频点，截断能固定输出长度，两者都不改变真实分辨率。norm 参数控制归一化方式，默认 backward 是标准定义，做能量守恒分析时用 ortho 更顺手——同一份数据在不同归一化下幅值差一个系数，混着用会算错。workers 参数在多核机器上对超大数组有效，小数组上开了反而有调度开销。

补充一个常见疑问：到底要不要归一化？如果你的目标只是找峰值频率，归一化方式不影响结果；如果你要读绝对幅值、算能量，就必须把 norm 与 2/N 修正统一起来，自始至终用同一套约定，混用必错。

3.4 什么时候值得折腾后端

默认的 pocketfft 在绝大多数场景已经够快。只有当 FFT 是应用的绝对热点、数组尺寸固定且巨大、并且你确认瓶颈在变换本身时，才值得尝试切换后端或预计算计划。多数情况下，把时间花在窗函数选择和数据质量上，回报大得多。

3.5 二维场景：图像频谱

fft2 与 rfft2 处理图像：把图像变换到频域后，低频集中在角落，用 fftshift 搬到中心，就得到常见的"中心亮斑加十字"的频谱图。图像去噪的一种经典思路就是频域操作——把高频噪声对应的区域乘一个衰减系数，再逆变换回来。这个视角我们在 4.3 和 4.4 还会碰到，这里先记住：一维的一切结论，在二维都有对应版本。

本节速览

• 要点一：rfft 配 rfftfreq 是实数信号频谱分析的标准组合，输出长度是 N/2 加 1。
• 要点二：scipy.fft 相对 numpy.fft 的优势在后端切换、实数变换家族与配套频率轴工具。
• 要点三：奈奎斯特频率是采样率的一半，超过它的成分会混叠成假低频。
• 要点四：频率分辨率等于采样总时长的倒数，补零只光滑曲线、不提升分辨率。
• 要点五：非整周期截断会产生频谱泄漏，加窗能抑制旁瓣，但会压峰值、宽主瓣。
• 要点六：单边幅值谱要乘 2 除以 N 修正，直流分量除外。
• 要点七：默认后端已经够快，后端切换是锦上添花，不是必需品。

看清楚了频谱，下一步就是动手改：下一节我们用 scipy.signal 设计滤波器，把干扰和噪声按频率砍掉。