ARTICLE DETAIL

建站实战干货

来自一线的建站与推广经验沉淀,每一条都经过真实交付验证。

DFT——计算机眼里的傅里叶:从公式到频谱分析实战

2026/9/16 22:59:35 拓冰建站 浏览量
DFT——计算机眼里的傅里叶:从公式到频谱分析实战 但凡上手做过信号处理的人早晚都会碰到一个知识裂谷教科书上的傅里叶变换是一个漂亮的积分公式变量是连续的、积分范围是正负无穷可一旦你打开 Python、MATLAB 或者 C 语言准备把傅里叶变换四个字变成代码能调用的函数却是fft、DFT输入的数组是离散的、长度是有限的输出也是一长串离散的数字。你会隐约觉得这俩是同一个东西但公式怎么对不上数组下标和频率到底怎么换算为什么算出来的频谱有一半是对称的这篇文章就来把这个裂谷填上。标题叫DFT——计算机眼里的傅里叶它解决的就是那个核心痛点为什么计算机只能理解离散且有限的信号以及它眼里的傅里叶跟你手算的傅里叶到底差在哪一步。这篇系列入门文章适合信号处理初学者、准备做频谱分析但没系统学过数字信号处理的人也适合那些已经跑过fft但完全不知道输出每个点代表什么、想真正看懂频谱图的读者。我会从连续傅里叶变换的不可计算性讲起一直讲到 DFT 公式的每一个符号、频谱图的横纵轴、泄漏与补零的真相最后带一段可以照抄的 Python 实战。1. 为什么连续傅里叶变换在计算机里根本算不出来1.1 教科书公式的三重不可能随便翻开一本信号与系统教材傅里叶变换长这样[ X(f) \int_{-\infty}^{\infty} x(t) e^{-j2\pi ft} dt ]这个公式在数学上极其优雅它告诉我们任意信号都能分解成不同频率的正弦波的叠加。但你要是真打算把它写成程序会立刻撞上三堵墙。第一堵墙是时间无限。积分下限负无穷到正无穷意味着你得拿到一个从宇宙诞生到宇宙毁灭都在持续的信号才有资格算精确傅里叶变换。可我们手里的数据永远是某个时间窗口内采集的一段一帧语音、一段加速度计读数、一秒钟的雷达回波。数据进来那一刻时间就已经被切断了。第二堵墙是变量连续。物理世界里的信号在时间轴上是连续的但计算机内存里存不了连续的东西。无论你采样率多高最终落进内存的只能是一串离散的采样点就像用拍立得给一条连续曲线拍照无论相机多高级照片上只有有限个像素点。第三堵墙是积分本身。计算机最擅长的事情是加法和乘法让机器去算一个真正的积分它只能通过数值近似来骗自己。数值积分本身就是把连续区域切成一小块一小块求和你算得再精也不过是近似。所以结论很直接连续傅里叶变换是个完美的数学玩具但计算机根本没法直接运行它。要让傅里叶分析在数字世界里落地必须做一套离散化改造。1.2 三管齐下的离散化改造要让傅里叶变换能进计算机你必须在三个维度上做妥协这套妥协方案可以说是数字信号处理最早的地基时域采样把连续时间信号x(t)按固定间隔Ts取值得到离散序列x[n]其中n是整数下标。这个过程对应的是 ADC 采样代价是引入最高频率上限奈奎斯特频率fs/2。时间截断序列不可能无限长只取其中N个点。这在数学上相当于拿一个长度为 N 的矩形窗去切信号代价就是后面要细说的频谱泄漏。频域采样把原本连续的频率轴f也离散化只算 N 个等间隔的频率点。这样频率轴就从连续变成了梳子齿每根齿之间的距离叫频率分辨率。有意思的是这三步妥协互相咬合在一起时域采样决定了频域会周期延拓时间截断决定了频域不再是一条条干净的谱线而是有了裙边频域采样则把频谱变成了一组有限的离散值。最终得到的就是 DFT——离散傅里叶变换。1.3 计算机真正运行的公式长什么样经过上述离散化改造计算机真正执行的 DFT 公式是[ X[k] \sum_{n0}^{N-1} x[n] \cdot e^{-j 2\pi k n / N}, \quad k 0, 1, \dots, N-1 ]这个式子看着复杂拆开其实非常简单。它做的其实就是一件事拿一个频率为k/N归一化频率的复指数信号与你原始信号的每一个采样点相乘然后求和。如果原始信号里刚好含有这个频率的成分两者会共振求和结果很大如果没有这个频率成分求和结果趋于零。你把k从 0 到 N-1 各算一遍就得到了 N 个数每个数对应一个频率点上的相关程度。注意这里X[k]是复数——它的模长告诉你这个频率成分的强度它的幅角告诉你这个频率成分的相位。这就是 DFT 的全部秘密它本质上是 N 个频率匹配滤波器或者说是一个由 N 个模板构成的匹配检测器分别去原信号里找每个频率。这段理解非常重要因为它解释了为什么 DFT 的输出是复数。很多初学者只看频谱图的模值把相位信息丢掉了这在某些场景下比如需要精确重建信号是会出大问题的。2. 从公式到频谱图DFT 的每个输出点到底代表什么2.1 横轴下标 k 怎么换算成赫兹DFT 算完你得到的是一个长度为 N 的复数数组X[k]问题来了数组下标 k 和物理频率 Hz 是什么关系答案是[ f_k k \cdot \frac{f_s}{N} ]其中fs是采样率。简单说相邻两个下标之间的频率间隔是fs/N这个值通常叫频率分辨率bin size。它意味着如果你的信号里有两个频率相差不到fs/N的分量DFT 基本没法把它们区分开。举个例子采样率fs 1000 Hz做N 1000点的 DFT那么下标 0 对应 0 Hz下标 1 对应 1 Hz下标 2 对应 2 Hz……一直到下标 999 对应 999 Hz。频率分辨率就是 1 Hz。这里有个几乎所有初学者都会问的问题下标 k 到底能取到多大最高对应的频率是多少理论上 k 能取到 N-1也就是说最高频率能算到fs - fs/N。但真实信号通常只关心到fs/2原因看下一节。2.2 对称性为什么频谱图只看一半如果你对一个实信号做 DFT然后画出|X[k]|幅度谱会发现一个规律整个序列关于N/2对称。X[1]和X[N-1]的幅度几乎一样X[2]和X[N-2]的幅度几乎一样以此类推。这不是巧合而是因为实信号的频谱具有共轭对称性X[k] X*[N-k]。换句话说前一半0 到 N/2-1是正频率部分对应 0 到fs/2奈奎斯特频率的物理频率后一半N/2 到 N-1是负频率部分对应-fs/2到 0 的频率。为什么实信号会有负频率你可以这样理解一个实正弦波cos(2πft)用欧拉公式展开是(e^{j2πft} e^{-j2πft})/2那一项e^{-j2πft}就是负频率分量。数学上负频率是复指数正交基的必要组成部分不是一个可以随意丢掉的另一半。但好消息是对于纯实信号负频率部分不携带任何额外信息完全由正频率部分决定。所以我们画频谱图时通常只画k 0到k N/2这一段也就是从 0 到fs/2的单边谱。如果你直接用numpy.fft.fft的裸结果画全谱会发现右半边像个镜像一样那不是 bug是 DFT 的天然属性。2.3 纵轴怎么从复数模值还原真实振幅幅度谱的纵轴也是很多人的迷惑点。X[k]是一个复数它的模|X[k]|算出来往往是一个很大的数和你心里预期的这个正弦波的振幅是多少伏完全对不上。原因在于 DFT 求和公式里把所有采样点都加了一遍累加规模大约正比于 N 点之和。系数校正其实很简单对于一个幅度为A的实正弦信号在完全对齐频率的情况下即正弦频率正好落在某个 k 上且没有泄漏|X[k]| A * N/2。所以还原真实振幅要乘以2/N [ A \approx \frac{2}{N}|X[k]| ]直流分量k0是个例外它的值是A_dc * N还原直流振幅只需要除以 N。我第一次做频谱分析时算出来的峰值振幅比预期大了近 100 倍后来才反应过来漏掉了这个2/N的系数。这个细节在教材里一句归一化就带过去了但实际工程里是花了很多时间才排查出来的。另外X[k]的幅角angle(X[k])给出了该频率分量的初相位如果你想从频谱里恢复某个正弦的精确参数幅度、频率、相位这三个量缺一不可。3. 频谱泄漏、栅栏效应与补零绕不开的三个高频翻车点3.1 泄漏的本质矩形窗在频域惹的祸上一篇说时间截断是 DFT 的必要条件但它的代价就是频谱泄漏。当你只取一段 N 个点的数据时相当于把一个无限长的信号乘以一个矩形窗[ x_{\text{截断}}[n] x[n] \cdot w[n] ]时域相乘频域卷积。矩形窗的频谱是 sinc 函数——主瓣很窄但旁瓣很高、衰减很慢。这个旁瓣就是泄漏的来源原本应该只在某个频率上出现的一根干净谱线会被 sinc 的旁瓣涂抹到附近的很多频率点上。更麻烦的是如果你的正弦频率不在某个离散的 k 上大多数真实场景都是如此泄漏会变得更加明显。教科书里那种漂亮的单根谱线只存在于正弦频率严格落在某个 k 值上的运气好情况。实际测出的频谱往往带一个山包主峰附近的点都明显偏大初看起来像是信号有好几个频率成分其实都是同一个正弦泄漏造成的假象。3.2 加窗在分辨率和泄漏之间做取舍抑制泄漏最常用的手段是加窗。所谓加窗就是对原始信号逐点乘一个非矩形的窗函数让数据两端逐渐衰减到零而不是被矩形窗那样一刀切。常见窗函数有窗函数主瓣宽度点数最高旁瓣电平适用场景矩形窗2-13 dB瞬态信号、频率分辨率优先汉宁窗4-31 dB一般频谱分析均衡性好海明窗4-41 dB窄带信号旁瓣抑制更好布莱克曼窗6-57 dB强调旁瓣抑制弱化分辨率加窗的代价也很明确主瓣变宽频率分辨率变差。本来两个频率差一点点就能分开加窗后会糊成一团。所以加窗不是免费的午餐它是一个分辨率和泄漏之间的折中选择。我自己的经验是如果你不确定选什么窗先试汉宁窗它基本是万金油如果信号里有非常强的干扰分量再考虑布莱克曼。还有一个必须提醒的细节加窗之后信号能量变少了因为两端被乘上了小于 1 的系数所以振幅还原公式要修正。比如汉宁窗恢复振幅的系数大约要乘 2 倍左右严格要按窗函数的相干增益算。具体系数等于N / sum(w)即用窗函数所有点的求和值去归一化。3.3 补零的真相插值不等于鱼目混珠补零zero padding是另一个被严重误解的操作。做 FFT 时如果你在原始 N 个点后面补M个零然后对NM个点做 DFT会得到更细腻的频谱看起来好像分辨率变高了。但真相是补零没有增加任何信息它只是在原始的 DFT 结果之间做了插值让频谱曲线更平滑。原来两个 k 之间的频率点现在有了一些插值出来的点但真正的物理分辨率还是由原始数据长度决定的仍然是fs/N而不是fs/(NM)。怎么区分插值和真实分辨率提升一个简单测试假设原始数据是 1000 点采样率 1000 Hz两个正弦频率分别为 100 Hz 和 110 Hz距离 10 Hz大于fs/N 1 Hz可以分辨。如果把数据裁剪成 100 点分辨率为 10 Hz这两个正弦就很难分开了——这时候你补 900 个零频谱会显得很光滑但 100 Hz 和 110 Hz 还是糊在一起分不开。补零解决不了数据本身不够长的问题。不过补零依然有实用价值尤其是在做人耳听觉分析或粗略找峰时插值出来的光滑曲线更容易读数。它不能提升分辨率但可以让你更精准地看出峰值位置相当于从粗尺换成了游标卡尺。4. 用 Python 做一次完整的 DFT 频谱分析4.1 构造一个故意不太好算的测试信号光讲理论很难建立直觉下面用一段可以直接复制的 Python 代码走一遍完整流程。先构造一个包含两个正弦分量和一个直流偏置的合成信号import numpy as np import matplotlib.pyplot as plt fs 1000 N 1000 t np.arange(N) / fs # 直流分量 2.0 50 Hz 振幅 1.5 120 Hz 振幅 0.8 x 2.0 1.5 * np.cos(2 * np.pi * 50 * t) 0.8 * np.cos(2 * np.pi * 120 * t)这里我故意选择了与fs/N 1 Hz能整除的频率50 Hz 和 120 Hz 都是整数这样频谱不会泄漏结果更干净方便对比验证。等你跑通后再把频率改成50.5和120.3看看泄漏是什么样的体会会更深。4.2 一次 FFT 以及横轴的正确生成方式直接用numpy.fft.fft做变换然后用numpy.fft.fftfreq生成对应的频率轴X np.fft.fft(x) freqs np.fft.fftfreq(N, d1/fs) # 单边谱只取前一半 half N // 2 freqs_1sided freqs[:half] mag_1sided np.abs(X[:half]) * 2 / N # 直流分量单独修正 mag_1sided[0] np.abs(X[0]) / N plt.stem(freqs_1sided, mag_1sided, basefmt ) plt.xlabel(Frequency (Hz)) plt.ylabel(Amplitude) plt.xlim(0, 200) plt.grid(True, alpha0.3) plt.show()看到结果后你应该会发现频谱在 50 Hz 和 120 Hz 处各有一个干净的尖峰幅度分别约为 1.5 和 0.80 Hz 处为 2.0。这里有个容易犯的错np.fft.fftfreq返回的频率轴是从 0 到 fs 的完整序列直接用plot(freqs, abs(X))会把负频率部分也画出来。我建议一套固定流程就是初始化fftfreq→ 切片到前一半 → 同时切片 X 的前一半这样横坐标和纵坐标永远不会错位。我还想特别强调np.fft.fftshift的用法。当你对频谱做fftshift(X)后零频会跑到数组正中间横轴变成了-fs/2到fs/2这种双边谱图形在一些论文里常见。但我个人的习惯是如果只是做常规频谱分析就用单边谱横轴直观、不用解释只有真正需要观察负频率行为比如做调制解调分析时才用fftshift。4.3 噪声场景下的加窗实战对比为了看清加窗的实际效果再生成一个带噪声、频率不在整数 bin 上的测试信号rng np.random.default_rng(42) x_noisy 1.0 * np.cos(2 * np.pi * 50.5 * t) 0.3 * rng.standard_normal(N) X_raw np.fft.fft(x_noisy) mag_raw np.abs(X_raw[:N//2]) * 2 / N win np.hanning(N) x_win x_noisy * win coherent_gain win.sum() / N X_win np.fft.fft(x_win) mag_win np.abs(X_win[:N//2]) * 2 / (N * coherent_gain) freqs_half np.fft.fftfreq(N, d1/fs)[:N//2] plt.figure(figsize(10, 4)) plt.plot(freqs_half, mag_raw, alpha0.7, labelno window) plt.plot(freqs_half, mag_win, alpha0.7, labelhanning) plt.xlabel(Frequency (Hz)) plt.ylabel(Amplitude) plt.xlim(30, 70) plt.legend() plt.grid(True, alpha0.3) plt.show()运行这段代码你会看到不加窗的时候50 Hz 附近的尖峰底部有一层裙边旁瓣整个谱线看起来像涂了一圈粉底加汉宁窗后旁瓣明显压低但主峰变宽了峰顶略低一些。这就是加窗的直观效果。注意我代码里那个coherent_gain——因为加窗改变了信号的能量振幅还原必须除以窗函数平均增益否则算出来的振幅显著偏小。这是很多工程代码里容易漏掉的一环。5. 实战中几个容易翻车的细节与自查习惯5.1 三个让我花费过大量时间的错误第一个是忘记频率轴是物理频率还是归一化频率。numpy的 FFT 本身不关心你的采样率它只按下标工作。所有关于这个峰对应多少 Hz的换算都必须由你通过fs/N来手动完成。如果你直接把np.fft.fft的结果拿去做频率定位很容易把采样率或分辨率算错。我初期调试频谱时峰值总是跟预期差 2 倍后来发现是采样率填错了fs写成了1/dt但单位没换算。第二个是谱线计数错误。单边谱只取前N//2个点但很多人的代码里取的是int(N/2)当 N 是奇数时会丢一个点还有人直接从下标 0 画到 N/2把中点那个负频率交界点算进去了。这个细节看着小一旦生成 MATLAB 风格的专业报告图上多出来的毛刺往往就是这个原因。第三个是DC 分量淹没其他频率分量。如果你的信号有一个很大的直流偏置比如 ADC 数据带了 3.3V 的共模电平做 FFT 后 0 Hz 处会有一个巨大的尖峰旁瓣可能瞬间把附近几个 bin 的真实信号盖住。一个很实用的预处理方法是先去掉信号的均值再做 FFT这相当于只分析交流分量或者在画图时用plt.ylim把 DC 峰截掉只看局部两种方式各有适用场景但都要有意识去处理 DC而不是任它污染整张图。5.2 计算前必看的清单最终总结一套我自己每次做 DFT 前都会过的自查清单希望能帮你省去和我一样多的排查时间确认采样率fs是多少奈奎斯特频率是否大于信号中你感兴趣的频率上限确认数据长度 N 与频率分辨率fs/N你想区分的两个频率差是否大于这个值确认 FFT 结果使用时是否生成了正确的频率轴fftfreq用对没有是否取单边谱单边谱的边界是否处理好了需要精确振幅时是否做了2/N放大加窗场景还要除以窗的相干增益信号里有没有 DC 分量、有没有较强工频或机械振动等干扰源它们会不会通过泄漏影响你关心的频段如果做峰值频率的精确定位是否考虑了补零插值的影响是否明白插值不能替代真实分辨率保存结果时是否把复数频谱的模和相位都留下来了后面如果还需要重建时能不能恢复出原始信息。我自己的体会是DFT 本身只是一个工具它不会骗人但工具用错了场景结果就是一次次的回头排查。信号处理在实践中最花费时间的问题几乎都是参数不匹配和单位不统一这类基础问题而不是算法本身的原理。上面这张清单能挡住大部分低级的坑。