1. 为什么要从傅里叶聊到 DFT计算机的“数字化思维”1.1 连续世界的傅里叶很美但计算机不认信号处理这行干久了你会发现一个很有意思的现象刚入行的时候大家都在学傅里叶变换公式写得飞起连续信号、无穷积分、冲激函数一套一套的。可真到了写代码、跑实验、把算法部署到嵌入式板子上的时候几乎所有人用的都是 DFT离散傅里叶变换。为什么因为傅里叶变换的教科书版本是给“连续世界”设计的。它假设你有一个时间上连续、幅度上连续、理论上无限长的信号然后对它做积分变换。这个数学模型非常优美各种性质也很规整但它有个致命的前提你得能用一个数学表达式完整描述这个信号。可现实中的信号是这样的传感器每隔几毫秒采一个点网卡抓包抓到的是一段二进制流声音文件在硬盘上就是一堆 PCM 采样值。这些东西统统是离散的、有限长的。你没有办法把它塞进连续傅里叶变换的积分式里因为计算机根本没有“无限”和“连续”这两个概念。它只认识一个又一个的数。所以 DFT 的本质就是把傅里叶变换从“解析数学”搬进“数值计算”的一次妥协或者说是一次完美的落地。它不跟你谈无穷积分它只做一件事给你 N 个点我给你 N 个复数每一个复数代表这个信号在某个特定频率上的“含量”。就这一个功能支撑起了数字信号处理这座大厦的几乎所有楼层。1.2 DFT 到底在算什么从“无限积分”到“有限求和”如果只看结论DFT 的公式长得并不吓人。但初学者最容易卡住的地方恰恰是它和连续傅里叶变换之间那层窗户纸为什么积分变成了求和为什么指数上多了个 1/N为什么结果是一堆复数我的理解方式是这样的连续傅里叶变换是在一个“无限长的时域”里让信号去和一个频率为 f 的复指数做内积积分结果就是这个频率的分量有多强。DFT 把这件事做了两个改动。第一把积分换成求和。因为离散信号的时域本来就是一串点没有连续的积分域所以自然退化成对 N 个采样点逐个相乘再累加。第二把频域也离散化。连续傅里叶变换得到的是连续频谱频率可以取任意实数但计算机给不出无限多个频率点。DFT 只计算 N 个均匀分布的频率这 N 个频率恰好对应着信号里可能存在的“最基本信息单位”。所以 DFT 那句人话版本是我手里有一段长度固定、离散采样的信号我把它拆成 N 个不同频率的复指数信号分别看看原始信号和每一个频率“像不像”。像的程度用幅值表示像的时候差了多少相位用角度表示。这就是为什么 DFT 的输出是复数实部虚部合在一起既告诉你强度也告诉你时间偏移。这一节如果你能想通后面所有关于频谱、滤波、加窗、补零的内容都是在这个基础上盖楼。2. DFT 公式拆开看每一个输出点都是一个“匹配滤波器”2.1 公式里每个符号的含义与计算逻辑先把标准 DFT 公式摆出来X[k] Σ_{n0}^{N-1} x[n] · e^{-j(2πkn/N)}别被这个式子吓到逐个字母拆就很清楚x[n]时域信号的第 n 个采样点n 从 0 到 N-1一共 N 个点。X[k]频域输出k 从 0 到 N-1对应第 k 个频率点。这个频率点代表的实际物理频率是 f_k k·fs/Nfs 是采样率。e^{-j(2πkn/N)}一个旋转因子本质上是一个以特定速度旋转的复指数信号。这个复指数在 n 从 0 走到 N-1 的过程中恰好转了 k 圈。整个求和式将 x[n] 和这个旋转因子逐点相乘再累加得到的结果 X[k] 是一个复数。它的模长代表这个频率下的强度辐角代表相位。这个计算过程你可以把它想象成“对暗号”。x[n] 是你要检验的信号e^{-j(2πkn/N)} 是你准备好的一个模板。你把信号和模板逐点相乘再相加如果它们频率一致乘积在大多数点上是正数累加结果就很大如果频率不一致乘积正负交替累加结果趋近于零。我第一次学的时候觉得这个说法太玄后来发现它在数学上完全站得住脚DFT 的每一个输出点本质上就是一个匹配滤波器专门用来检测某个特定频率在信号里是否存在。做 N 个频率点的检测就得到 N 个输出也就是整个 DFT 向量。2.2 为什么除以 N从“求和放大”看归一化DFT 公式里其实没有除以 N但实际用的时候几乎都要除很多人就在这一步栽跟头。原因很简单。你想想如果信号里有一个频率分量幅度是 A那么它和同频率的模板逐点相乘每个点大约得到 A 附近的值N 个点加起来就成了 N·A。所以 DFT 输出的模长天然比真实幅度大了 N 倍。也就是说X[k] 算出来之后你还得除以 N 才能得到“真实幅度”。不过在工程上很少有人直接除以 N 就完事因为 DFT 默认算的是双边谱负频率那边还露着一半幅度。所以当你只想看正频率的时候常规做法是除 N 之后除了直流分量 k0 之外把其余频率点的幅值都乘以 2这样才能把负频率那半边的能量也算回来。一个具体的数值例子能让你感受更直观我生成了一个幅度为 1.0、频率为 50 Hz 的正弦波采样率 1000 Hz采样 1000 点。做完 DFT找到 k50 那个输出点的模长裸算出来大约是 500因为 N1000。然后我除以 N得到 0.5再乘以 2得到 1.0这才和原始信号幅度对上了。如果你漏了乘以 2 这个步骤频谱图上所有峰的高度都会比真实值少一半幅度校准怎么都对不上。2.3 参数决定视野点数 N、采样率 fs 与频率分辨率 ΔfDFT 能不能看清你想看的频率取决于三个参数怎么选。第一采样率 fs。采样率决定的是频率轴的最大范围也就是奈奎斯特频率 fs/2。高于这个频率的成分会被“伪装”成低频信号这就是混叠。第二节会专门展开这里先记住一个判断标准如果你关心信号里最高频率是 200 Hz采样率至少要到 400 Hz实际工程里通常会留 1.5 到 2 倍余量。第二点数 N。N 越大频率分辨率越细。这个关系用式子表达是 Δf fs/N。这个式子要把它刻在脑子里。它告诉你如果你用 1000 Hz 采样率采集了 1 秒钟数据也就是 N1000那你能分辨的最低频率间隔就是 1 Hz。两个频率相差不到 1 Hz 的正弦波在频谱图上会糊成一团你根本分不清是几个峰。第三时间长度。本质上 N 就是采样时间乘以采样率。想提高分辨率最直接的办法是“多采样一段时间”而不是盲目提高采样率。我见过很多初学者为了“看得更清楚”把采样率从 1000 提到 10000结果点数没变频率范围倒是变大了分辨率反而更差。这就是没想清楚参数之间的关系。参数选择这一段我的个人习惯是先定频率范围和分辨率需求再由分辨率需求倒推采样时长最后用采样时长乘带宽需求定采样率。顺序反过来后面调参一定别扭。3. DFT 的三大“坑”泄漏、混叠、栅栏效应3.1 频谱泄露截断带来的“假频率”DFT 处理的是有限长信号这就有一个隐藏前提你在时域上把一个无限长的信号“截断”了只留下中间的 N 个点。这个截断过程在频域上等于原始频谱和一个矩形窗的频谱做了卷积结果就是能量从真实的频率点向两边“溢”出去形成一圈旁瓣。这就是频谱泄露。你做个实验就能感受到生成一个 50 Hz 正弦波采样率 1000 Hz恰好采完整数周期做个 DFT频谱干净得像教科书只在 50 Hz 处有一个尖峰。但如果采样点数不是周期的整数倍比如你截取了 1.3 个周期频谱图就“糊”了50 Hz 旁边冒出一堆小峰主峰幅度也从 1.0 掉到了 0.8 左右。解决办法很成熟加窗。常用的有汉宁窗、海明窗、布莱克曼窗。窗函数的作用是让截断边界处的信号平滑地衰减到零而不是生硬切断。这样旁瓣会大幅降低但代价是主瓣变宽也就是说频率分辨率会受影响。这又是一个典型的工程权衡你要“干净”就得接受“粗”要“精细”就得接受“脏”。实际做频谱分析的时候我的经验是如果能控制采样点数恰好包含整数周期那就别加窗结果最准如果信号频率未知、周期无法对齐推荐汉宁窗起步它在旁瓣抑制和主瓣宽度之间比较平衡。加窗之后幅度记得要乘以归一化系数否则测出的幅值会偏小这一点最容易被忽略。3.2 混叠采样不够快的代价混叠这个问题学过奈奎斯特定理的人都知道名字但真正理解透彻的不多。我用一句话解释当采样率不够高时高于 fs/2 的频率成分会“冒充”一个低于 fs/2 的频率出现在频谱里而且你根本区分不出来。我给了个学生活的类比车轮转得比相机快门还快的时候视频里车轮看起来在倒转。这就是时域采样里的混叠视觉版的。信号处理世界里的混叠一模一样只是藏在数字数据里肉眼不可见。那么怎么防两个思路。第一提高采样率让所有值得一提的频率都低于奈奎斯特频率。第二在 ADC 之前加模拟低通滤波器也就是抗混叠滤波器把高于奈奎斯特频率的成分提前滤掉。工程上这两件事经常同时做。如果你已经采完了一组数据才发现有混叠那就没办法了信息已经在采样过程中永久丢失了神仙也救不回来。这点我强调很多次因为我在实际项目里看到过的混叠错误案例几乎都是同一个套路采样率随便取滤波器省掉等看到频谱里出现一个诡异的小峰时已经晚了只能重新采集。3.3 栅栏效应与补零是修图不是增分辨率栅栏效应这个词有点冷门但它描述的痛点每个人都会遇到。DFT 只计算 N 个固定频率点就好像隔着栅栏看待频谱你能看到的只是栅栏缝隙里的那些频率真正的峰值如果恰好落在两个缝隙之间你观察到的幅度就会比真实值低位置也会偏。好多人的第一反应是那我补零把 N 变大看不就更细了吗这个方法有效但出发点要想明白补零的真正作用是做一个“频谱插值”让频谱曲线看起来更平滑峰值定位更准但它不会增加任何新的信息量。它的本质是改变了 DFT 的观察网格并没有改变你手中的数据量。举个例子你有 1000 个采样点fft 之后除以 1000得到 1000 个频域点。如果你在后面补 9000 个零再做 10000 点 FFT频谱上每个点之间的间隔缩小到原来的十分之一曲线确实更细腻了。但注意基波位置还是由那 1000 个真实数据点决定的你能把峰值定位得更精确但两个真实靠得很近的频率该分不开还是分不开。想真正分辨出靠得很近的两个峰只有一条路采集更长时间的信号让原始数据点数增多也就是提高真实的分辨率。所以我的实用建议是补零用来做峰值的精细定位非常有效尤其在高精度测频场景里但如果你需要分辨两个相差 1 Hz 的信号请老老实实多采几秒数据不要指望补零能变魔术。4. 实操流程用 Python 从 0 到 1 做一次 DFT 分析4.1 准备工作与工具选型DFT 上手最快的路径就是 Python。它不用你手写 DFT因为有一个算法叫 FFT快速傅里叶变换能把 DFT 的计算复杂度从 O(N²) 降到 O(N log N)这也是为什么 DFT 在工程上能真正用得起来。标准库 numpy 里的 np.fft.fft 就是最常用的入口。我建议一开始不要跳进 numpy 以外的库先把 numpy 自带的 fft 玩明白因为它足够快、接口简单、文档齐全。matplotlib 用来画图。两者配合已经能完成 90% 的频谱分析需求。后文所有代码基于 Python 3 和 numpy 1.2x、matplotlib 3.x。你需要一个能跑 jupyter notebook 或者任意一个 Python 脚本运行环境就行。4.2 写代码生成信号、加窗、FFT、频谱绘制先准备一个混合信号。我造一个包含 50 Hz 和 120 Hz 两个正弦波幅度分别为 1.0 和 0.8再叠加一点随机噪声采样率 fs1000 Hz。代码里我特意注释清楚每个细节。import numpy as np import matplotlib.pyplot as plt # 参数设置 fs 1000 # 采样率 1000 Hz N 1000 # 采样点数采样时长 1 秒 n np.arange(N) # 样本序号 0~N-1 # 生成混合信号50Hz 幅度1.0120Hz 幅度0.8加噪声 t n / fs x 1.0 * np.sin(2 * np.pi * 50 * t) 0.8 * np.sin(2 * np.pi * 120 * t) x 0.05 * np.random.randn(N) # 生成频率轴单位 Hz freqs np.fft.fftfreq(N, 1/fs) # 做 FFT X np.fft.fft(x) # 幅度谱归一化先除以 N mag np.abs(X) / N # 单边谱除了直流分量外其余频率点乘 2 mag_single mag.copy() mag_single[1:] mag_single[1:] * 2 # 画图 plt.figure(figsize(12, 5)) plt.subplot(1, 2, 1) plt.plot(t, x) plt.title(时间域波形) plt.xlabel(时间 (s)) plt.ylabel(幅度) plt.grid(True) plt.subplot(1, 2, 2) plt.plot(freqs[:N//2], mag_single[:N//2]) plt.title(幅度谱单边) plt.xlabel(频率 (Hz)) plt.ylabel(幅度) plt.grid(True) plt.tight_layout() plt.show()跑完之后你会看到频谱图上 50 Hz 和 120 Hz 附近各有一个尖峰峰值高度分别接近 1.0 和 0.8和原始设置吻合。我这里加了噪声所以峰值周围会有一点底噪正常。有个容易被忽略又影响很大的环节噪声的随机种子没固定多次跑结果会有细微抖动。做演示无所谓做科研和工程测试时最好 np.random.seed(0) 固定种子保证结果可复现。这个习惯我建议从第一天就养成。4.3 结果怎么读幅度、相位、直流与镜像频率幅度谱读法很简单横轴是频率纵轴是幅度峰的位置对应信号里的频率峰的高度对应这个频率的强度。初学者容易出错的是这三个地方。第一直流分量。k0 那个点代表信号的直流成分也就是平均值。如果信号均值不为零直流分量会非常高可能把其他频率的峰“压扁”。处理完均值不为零的信号建议先减均值再 FFT或者画图时直接用 ylim 限制纵轴范围。第二种做法更方便。第二相位谱。大多数人做频谱分析只关心幅度谱但相位信息同样有用。相位谱是 np.angle(X)它告诉你每个频率成分的时间偏移。读取相位的时候要注意它会对频率和噪声非常敏感信号一有噪声相位值就抖动得很厉害所以实际工程中非必要不看相位谱。第三镜像频率。DFT 输出的后半段k N/2对应的是负频率对实信号来说负频率部分和正频率部分是共轭对称的信息重复。所以画图时我永远只画 freq[:N//2] 那半。你如果直接画全部 N 个点会看到左右对称的两个峰新手容易误以为信号里有两对频率其实是同一对的镜像。顺带说一句np.fft.fftfreq(N, 1/fs) 生成的频率轴前一半是正频率后一半是负频率顺序和 fft 输出一一对应。图里我只取前半段就是为了避开镜像。5. 常见问题与排查技巧实录5.1 结果全是 NaN 或 Inf这个问题在做高采样率、超长序列时尤其容易出现。主要原因几乎都是数据里有 NaN 或 Inf一旦出现这两个值FFT 会把它传染给所有频率点。排查方法很简单先检查输入数据。print(np.isnan(x).any(), np.isinf(x).any())如果返回 True就找到源头了。在信号预处理阶段我会统一做一次清洗去掉异常点或者用插值替换。千万别带着脏数据直接 FFT结果会让你怀疑人生。另一个不太常见但确实存在的坑是输入数据点数为非整数次幂有些旧的 FFT 实现会退化到很慢的朴算法。numpy 里不存在这个问题但如果你自己实现过 FFT或者用了某些老库建议把数据长度补到 2 的幂次再做。这不算 bug但能避免莫名奇怪的性能问题。5.2 频率峰值偏了一个点怎么精细估计频率你测一个 50 Hz 的信号fft 之后发现最高峰在 49.1 Hz 或者 50.9 Hz但直觉告诉你就是 50 Hz。这其实是栅栏效应的典型表现峰值真实位置在两个离散频率点之间。最省事的办法是补零把 N 扩充到 10 倍或 20 倍再做 FFT峰的位置会收敛到真实频率附近。补零不能增加真实分辨率但它可以给你一个更密的频率网格定位精度可以提高不少。如果你连补零都嫌麻烦还有一个经典的经验公式抛物线插值。取峰值点以及它左右两个点用二次抛物线拟合抛物线顶点的横坐标就是更精确的频率估计。这个方法在大多数场景下已经够用而且计算量几乎为零。我在实际项目里测频精度要求到 0.01 Hz 的时候用的就是“补零到 10 倍 抛物线插值”的组合拳。实测下来稳定性和精度都很好比单纯增加采样时长更经济。5.3 直流分量巨大把其他频谱“压扁”了如果你采集的信号本身有直流偏置比如传感器输出是 2.5V 基准上叠加小信号FFT 之后 0 Hz 处会有一个巨峰旁边的信号可能被它衬得看不见。这不代表小信号不存在只是画图比例被直流压住了。处理办法就一句话减均值。x x - np.mean(x)做完这一步直流分量基本归零频谱图的比例一下就正常了。这里有个细微差别减均值只是减掉常数分量不会影响你关心的交流频率成分所以完全可以放心用。还有另一种情况信号本身就是一个大直流上加了很小的变化量你既要看直流也要看小信号。这时候别减均值直接做 FFT然后把纵轴改成分贝刻度也就是 20*log10(mag)这样大动态范围的信号能同时显示出来。分贝刻度是看这类频谱的好朋友。5.4 加了窗幅度反而变小了加窗之后信号的能量有一部分被窗函数在边界处“吃掉”了所以算出来的幅度比真实值小。拿汉宁窗举例它的相干增益约 0.5也就是说加窗后的幅值大约只有真实值的一半。很多人在这里愣住以为是代码错了。解决办法是根据窗函数做幅度修正。标准做法是加窗之后幅值除以窗函数的均值。汉宁窗的均值是 0.5所以只要在幅度谱上乘以 2。更通用的做法在代码里应该是from numpy import hanning w hanning(N) xw x * w Xw np.fft.fft(xw) mag_w np.abs(Xw) / N / np.mean(w)关键点就是那个 np.mean(w)。它把窗函数造成的能量衰减“补回来”。用这个写法无论你换什么窗幅度都能校准到接近真实值。不过要提醒一句加窗永远是有代价的。主瓣宽度变宽可能让原本能分开的两个邻近频率开始粘连。所以在工程上如果要精确测幅值尽量用矩形窗不加窗或者确保采样周期是整数周期如果不确定别盯着幅度精度先保证频谱干净幅度再做修正。5.5 FFT 点数与采样点数不一致该补零还是截断有些人会为了性能或分辨率把长度凑到 65536 之类的数这时候你手里的数据可能只有 3000 点。处理办法分情况。如果你还是有原始数据的获取条件建议重新采集把点数和 FFT 点数对齐。这是最干净的做法。如果数据已经固定只能从手里的 3000 点出发那么两种选择补零到 65536或者直接 N3000。这两者到底有什么区别简单说补零不改变分辨率只改变频谱图的细致程度直接 3000 点虽然频率网格粗一点但每个频率点的能量是真实计算出来的没有“虚假插入”。碰到这种情况我个人的习惯是先做 N3000 的 FFT 看整体轮廓如果某个峰的位置需要精细定位再补零到更大的点数做插值。不会一上来就盲目补零那样只会让图看起来精致但并不会让数据更可信。从我开始用 DFT 到现在踩过的坑比想象中多得多。有一段时间我迷信所有“看起来更精细”的数据处理手段后来才发现DFT 最核心的素养不是会调用 fft 函数而是知道每一个参数控制了哪些频域行为每一个操作背后付出了什么代价。搞清楚这些你才算真正读懂了“计算机眼里的傅里叶”。
