从公式到代码：手写DFT彻底搞懂FFT与频谱分析-北京尧图网络

1. 从公式到代码DFT到底在算什么搞信号处理的人十有八九都经历过这样一个阶段教材翻到离散傅里叶变换那一章公式看了无数遍笔记抄了厚厚一摞考试也能拿高分但一问到DFT的代码到底怎么写或者用Python跑一个DFT需要几步就会有点卡壳。这不是基础不扎实而是教学和实操之间始终隔着一层窗户纸——公式是抽象的语言代码是具体的工具两边得对上账。先把这个账算清楚。离散傅里叶变换的教科书定义是长这样的[ X[k] \sum_{n0}^{N-1} x[n] \cdot e^{-j 2\pi k n / N} ]这里面(x[n])是时域信号的第n个采样点(N)是总采样点数(X[k])是频域的第k个频率分量。看着挺简洁但你仔细品一品这里头藏着一个迭代套迭代的过程对于频域里的每一个k你都要把所有的n都遍历一遍把每个时域点乘以一个旋转因子(e^{-j 2\pi k n / N})然后累加。所以DFT的暴力计算复杂度是O(N²)N一旦过万跑起来就有点吃力了这也是为什么后来有人搞出了FFT快速傅里叶变换——它本质上是DFT的一种分治加速算法但很多人在学FFT之前连最基础的DFT都没有亲手实现过。我见过不少同学论文里明明是使用FFT对信号进行分析结果跑完一看频谱横坐标根本对应不上信号的物理频率波形画出来乱七八糟。问题往往就出在他对DFT内部的索引和归一化没有概念直接用别人封装好的函数出了错也找不着方向。所以这篇博文的目的很明确把DFT从公式到代码一行一行拆给你看。你不需要是数学专业的只要能看懂简单的循环跟着我一起把这段代码从零写出来你就能彻底搞懂DFT在每个索引、每个参数、每个归一化系数里到底干了什么。等你再回头看Matlab的fft、Python的numpy.fft会发现它们都只是一层壳子壳子底下就是这篇博文里讲的东西。这篇内容适合人群包括正在学《数字信号处理》的本科生、刚入门音频算法或者振动分析的工程师、以及所有会用fft但说不清fft为什么能这么算的实践派玩家。2. 为什么我们要手写一遍DFT你可能第一反应是这都什么年代了numpy一行fft不香吗为什么还要自己写一遍DFT这里我不是让你项目里舍近求远、非去手写O(N²)的代码而是强烈建议你在学习阶段至少完整地手写一次。理由很简单你只有亲手实现了DFT才能真正理解FFT的优化动机和频谱图中的种种细节。这就像你见过别人开车但自己真正上路开一趟之后才明白油门、刹车、方向盘在转弯时是怎么配合的。函数库是黑盒写代码时你就是白盒。具体来说手写一遍DFT能帮你解决这几个平时问得最多的问题频域的X[k]到底代表什么k1和k2之间的间隔是多少Hz为什么采样率fs决定了频谱的最大频率范围频谱图上为什么要取一半共轭对称是怎么回事为什么好多代码里最后要除以N不除会怎样DC分量就是频率为0的那一项为什么是个巨大的数这些问题你如果不看DFT的代码现场只拿着公式去推也能推出来但那叫学会了数学你把代码写出来、数据打出来亲眼看着某个k对应的复数和信号频率对上号这才叫学会了工程。另外一个很实际的原因是在做嵌入式开发或者FPGA实现时很多时候你没有现成的FFT库可以用甚至连浮点运算能力都有限。这时候你需要根据DFT的定义去实现一个够用就行的变换逻辑。比如在某些故障诊断场景里你只需要算特定几个频率点的幅值根本不需要算整个频谱这时候直接按DFT公式计算那几条特定谱线反而比FFT全家桶要快得多。这个优化的前提条件是什么就是你深刻地理解DFT公式里每个符号对应的代码含义。所以我建议这篇文章读完以后你先别急着用现成的numpy.fft而是先按照我的代码把DFT从零到一撸一遍。哪怕你只是照着敲了一遍收获也远比用过一百次fft要大得多。3. DFT代码的骨架循环、复数与索引接下来进入正题。我们先从最本原的算法结构说起然后一步一步写出可运行的Python代码。3.1 算法骨架三重循环还是双重循环先看公式[ X[k] \sum_{n0}^{N-1} x[n] \cdot e^{-j 2\pi k n / N} ]这其实是一个双重循环结构。外层循环遍历频率索引k从0到N-1内层循环遍历时间索引n从0到N-1。每个元素的计算都是一次复数乘法再加上一次复数加法。写成伪代码就是for k in range(N): sum 0 for n in range(N): sum x[n] * exp(-1j * 2 * pi * k * n / N) X[k] sum这里三重循环是个口语化的笼统说法严格来说就是嵌套双重循环。有些教材里会把旋转因子(W_N^{kn} e^{-j2\pi kn/N})单列出来看起来好像多了一重循环其实那只是为了写公式方便实际计算时它的复杂度还是O(N²)。3.2 旋转因子的几何直觉很多人一看到(e^{-j\theta})就头大其实在代码里它就是cos(theta) - 1j * sin(theta)本质上是复数平面上的一个单位圆上的点。当我第一次自己写DFT代码的时候让我真正开窍的是这个角度每个k值对应的旋转因子其实就是在复数平面上以不同步长绕圈。当k0时旋转因子是(e^01)所以X[0]就是把所有x[n]加起来平均值的N倍也就是DC分量。当k1时旋转因子按步长(2\pi/N)绕完整一圈对应频率为(\frac{1}{N \cdot T_s})其中(T_s)是采样间隔即对应频率(f_s/N)。当k2时步长变成(2 \times 2\pi/N)绕两圈对应频率(2f_s/N)。换句话说在N个点的序列里k对应的频率就是[ f_k \frac{k \cdot f_s}{N} ]这一点是后续一切频域分析的基石。你在代码里看到的所有频率坐标轴、频谱图横轴label它的底层来源就是这一条公式。3.3 复数表示与共轭对称代码里我们用的exp(-1j * 2 * pi * k * n / N)这个表达式里1j在Python中就是虚数单位。这里符号是负的-1j这是DFT的正向定义如果你用的是1j那就是逆变换方向或者你自定义的变体会造成频谱方向相反实信号倒还好但复数信号的频谱会整体翻转容易出问题。另外当输入信号是实信号时X[k]和X[N-k]是共轭对称的。这个对称性意味着频谱只有前一半N/2是有信息的后一半是镜像。这也是为什么几乎所有频谱图都只画前一半。不了解这个对称性的人很容易画出两条对称的山峰然后一脸困惑地问为什么我的频谱有两个峰答案很简单一个是正频率一个是负频率对实信号来说它们互为镜像。后面我会在示例里再演示一遍。4. 完整的DFT代码实现与逐步解析这里我用Python写一个最基础、最接近公式定义版本的DFT函数保证思路清晰不依赖第三方库。import math def dft_basic(x): 最基础的离散傅里叶变换实现 x: 一维实数或复数序列长度为N 返回: 长度为N的复数频谱序列 N len(x) X [0j] * N # 初始化复数数组 for k in range(N): sum_real 0.0 sum_imag 0.0 for n in range(N): angle -2.0 * math.pi * k * n / N # 欧拉公式: e^{j*angle} cos(angle) j*sin(angle) # 这里 angle 本身是负数所以写成加法正好对应 e^{-j*2*pi*k*n/N} sum_real x[n] * math.cos(angle) sum_imag x[n] * math.sin(angle) X[k] complex(sum_real, sum_imag) return X这段代码里我故意没有用cmath.exp而是直接展开成cos和sin。为什么要这么做一方面是因为在基础版代码里实部和虚部分开计算你更能看清复数乘法的本质另一方面在后续做性能优化或者移植到C语言的时候cos和sin拆分形式往往比直接调用复数指数函数更直白也方便做定点化处理。4.1 测试一个简单的正弦信号写完了函数得验证一下对不对。我们生成一个采样率(f_s1000Hz)、频率(f_0100Hz)、时长1秒的标准正弦波然后做DFT。import numpy as np import matplotlib.pyplot as plt fs 1000 # 采样率 1000 Hz T 1.0 # 时长 1 秒 N int(fs * T) # 采样点数 1000 t np.linspace(0, T, N, endpointFalse) # 时间轴 x np.sin(2 * np.pi * 100 * t) # 100 Hz 正弦波 X dft_basic(x)跑完以后我们可以把频谱幅值画出来mag np.abs(X) freqs np.arange(N) * fs / N plt.plot(freqs, mag) plt.xlabel(Frequency (Hz)) plt.ylabel(Magnitude)理论上这个频谱应该在(f100Hz)和(f900Hz)两个位置有尖峰。为什么900Hz也有因为1000Hz采样率下奈奎斯特频率是500Hz针对实信号来说100Hz对应负频率-100Hz也就是900Hz的镜像。如果你只画前一半频率0到500Hz那就只会看到100Hz这一个峰。这就是DFT代码里第一个经典新手惊讶点。记住实信号频谱是共轭对称的后一半是前一半的镜像。4.2 幅值的数值验证与归一化继续看数值。如果直接打印np.abs(X[100])和np.abs(X[900])你会发现它们的值大约是500。但信号的幅度明明是1为什么频谱幅度是500答案是DFT的定义里只有累加没有除以N。N1000一个完整正弦波的幅度为1经过DFT累加得到的峰值幅度就是(N/2 500)。如果你想画出来让人觉得物理上说得通就得对幅值做归一化单边频谱只看正频率就乘以(2/N)双边频谱正负频率都看就乘以(1/N)这个归一化的坑几乎所有画频谱图的人都踩过。有人画出来的频谱峰值动不动上几百有人画出来的DC分量巨大都是没搞明白归一化导致的。我自己做振动信号分析时对幅值的准确度要求很高这时候就特别需要弄清自己用的是单边幅值还是双边幅值否则后续所有的物理量单位都会不对。4.3 用numpy验证结果正确性有了手写版以后最好再和标准库做一次对照这样能确认代码逻辑没有隐藏bug。X_np np.fft.fft(x) # 对比前10个点的误差 diff np.max(np.abs(np.array(X) - X_np)) print(fMax difference: {diff})如果你跟我一样用的是上面那段dft_basicdiff应该会非常小大概在1e-13级别。之所以不是严格为0是因为浮点数的计算顺序差异以及numpy.fft内部可能用了FFT算法和不同精度的三角函数近似。看到这个量级的误差就可以放心了你的代码逻辑没问题频率和幅值都对应得上。5. 代码里的隐藏细节索引、频率分辨率和泄露写代码只是第一步真正让DFT从玩具变成工具的是理解这些代码细节背后的工程含义。这一部分我们聊几个一写代码就会碰到的隐藏问题。5.1 频率分辨率由什么决定看代码里的这一句freqs np.arange(N) * fs / N这说明频率轴的间隔是[ \Delta f \frac{f_s}{N} ]这就是频率分辨率。分母是N所以采样点数越多频率分辨率越高而采样率越高频率分辨率反而越粗因为总频率范围变宽了。在采样点数固定的前提下你没法同时提高最高频率和频率分辨率这是一个硬约束和算法实现无关。举个例子如果你采样率是1000Hz采样点数N100那么(\Delta f 10Hz)。这意味着两个频率相差小于10Hz的正弦信号在频谱上会糊成一团无法分辨。想提高分辨率要么加长采样时间增大N要么用更复杂的算法比如zoom FFT或者参数化频谱估计但这是另一篇文章的故事了。5.2 为什么频率轴只取前一半在DFT结果里对实信号来说X[k]和X[N-k]是共轭对称的。所以如果你把0到N-1全部画出来会看到一条以N/2为对称轴的镜像曲线。大多数频谱图只画0到N/2-1或者0到奈奎斯特频率(f_s/2)。代码处理时也很常见half N // 2 plt.plot(freqs[:half], mag[:half])或者转换为单边幅值谱mag_one_sided mag[:half] * 2 / N这里要注意DC分量和奈奎斯特频率处的那两个点不要乘以2因为它们在双边谱里没有镜像项。我见过很多人在这个细节上翻车明明第0个点就是DC结果还乘了2/N把直流幅值放大了两倍。5.3 频率泄露DFT代码跑出来不准的元凶另外一个常常在代码里体现出来的现象是频谱泄露。看我们前面的测试为什么用100Hz正弦信号频谱会在别的频率上也有非零幅值理论上一个完美的单频信号DFT频谱应该只在100Hz处有一个非零值。但实际跑出来你会发现100Hz两侧有一些小尾巴。原因在于当我们把信号截取为有限长度N时相当于给无限长的正弦信号乘上了一个矩形窗。这个矩形窗在频域里对应着sinc函数所以频谱上就出现了主瓣和旁瓣能量由主瓣向旁瓣泄露。从DFT代码角度看只有当信号频率正好等于(\frac{k f_s}{N})的某一个整数k时频谱才会是干净的单线如果频率不是整数格点比如100.7Hz那么它的能量就会铺开在多个k上这就叫泄露。实际工程处理中可以通过加窗汉宁窗、海明窗、布莱克曼窗等来抑制旁瓣但主瓣会变宽这是分辨率和动态范围之间的一个权衡。这些内容要在代码上做实验也非常直观同样的正弦数据加窗之后再做DFT旁瓣明显就降低了。所以学会裸写DFT之后你后续去理解窗函数、滤波器设计都会顺利得多。6. 从暴力DFT到FFT代码优化背后的思路前文提到暴力DFT的复杂度是O(N²)那么当N1024时大约需要100万次复数运算当N16384时运算量已经超过2.6亿次。在纯Python里跑这个速度慢到你几乎想砸键盘。这个时候我们自然会去想能不能把计算过程变得更高效答案是能那就是FFT。但FFT不是另外一套理论它只是利用了旋转因子的周期性和对称性把一个大DFT拆成若干个小的DFT递推合并最终让复杂度降到O(NlogN)。我不是要在这一篇文章里完整推导FFT蝶形运算但我想提醒你一点如果你手工实现过DFT再回头去看FFT的蝶形图会事半功倍。因为你会明白每一个蝶形运算里的加减乘除都是在算某几个特定k、n组合下的部分累加而不是什么神秘的黑魔法。6.1 理解旋转因子的周期性有助于优化旋转因子(W_N^{kn})有一个关键性质[ W_N^{kN} W_N^{k}, \quad W_N^{kN/2} -W_N^{k} ]前者是周期N后者是半周期反相。基于这两个性质很多重复的复数乘法可以共用计算值这就是FFT里查表法加速的根源。在实际工程中用Python实现FFT时我们通常都会预先算好一张旋转因子表# 预计算旋转因子表 twi_table [cmath.exp(-2j * math.pi * k / N) for k in range(N // 2)]这样在算法内层循环里就不需要反复调用三角函数可以省下大量时间。这个预处理思想你理解了DFT之后再去优化会觉得非常自然。6.2 什么时候该手写DFT、什么时候该用FFT这是一个很实际的问题。我的经验是如果你的N很小比如几十、几百暴力DFT和FFT差距并不明显手写DFT完全可以胜任。如果你只需要计算少数几个频率点比如工频50Hz、二倍频100Hz、三倍频150Hz可以直接按DFT公式只算这几个k复杂度从O(N²)降到O(MN)M是目标频率数非常划算。如果你需要完整的频谱且N比较大几千以上一定要用FFT库或者优化过的自研FFT否则运行时间会让你怀疑人生。在嵌入式系统里有时候片上内存有限FFT蝶形运算的倒位序重排比较麻烦而DFT的暴力循环结构天然适合顺序访问配合查表法也可以获得不错的性能。所以DFT已过时这种说法只存在于考试卷子上工程实际里DFT依然有它的位置。7. 常见问题与排查技巧实录最后这部分我把自己多年来看别人写DFT代码、以及自己调试DFT过程中遇到的问题整理成一份速查表按出现频率排个序方便你对号入座。7.1 频谱峰位置不对频率轴标错了这个问题出现的频次最高。很多人用某个现成库算出结果然后直接把k当成Hz画图结果100Hz的峰跑到了100的位置看起来也对但一个坐标是点号k一个是物理频率Hz两者相差一个比例因子(f_s/N)。如果你忘了乘这个因子或者乘错了整个频谱图都会偏移。排查方法先找一个已知频率的信号比如工频50Hz跑一遍DFT然后定位最大值对应的k值再验算(k \cdot f_s/N)是否等于50。不相等就是频率轴写错了。7.2 幅值不对归一化没搞清楚很多人画频谱图纵轴直接是abs(X)结果看到峰值几百几千一脸迷茫。这个我在4.2里已经解释过核心原因是DFT求和没有除以N。如果你是物理量分析通常需要的是单边幅值谱也就是乘以2/N如果你是功率谱密度分析那要用的又是另一种归一化方式。用错归一化后续所有读数都是错的。7.3 信号末端出现突变未使用正确的采样时间轴生成测试信号时经常有人这么写t np.linspace(0, T, N) # 错误包含了 T 这个点这看似没什么问题实际上最后一个点(tT)和第一个点(t0)的相位不连续相当于在信号末端人为制造了一个跳变在频谱上会引入杂散频率成分。正确做法是t np.linspace(0, T, N, endpointFalse) # 不包含右端点或者等价地写成t np.arange(N) / fs这个细节非常坑尤其在做正弦波扫频测试时稍有疏忽频谱里就会多出一堆莫名其妙的峰。7.4 用复数信号测试时有峰消失如果你输入的是复数指数信号(e^{j2\pi f_0 t})那么DFT结果应该是只有一个正频率峰。但如果你用实信号的正弦波去测试频谱会有两个对称峰。很多人拿着复数信号的期望去套实信号结果看到两个峰就以为代码写错了。请记住这两个峰是正常现象因为实信号同时包含正负频率分量。7.5 频率分辨率不够无法区分两个临近频率两个频率分别为100Hz和102Hz的正弦信号若采样率是1000Hz采样点数N100那么(\Delta f10Hz)这两个频率在频谱上会混在一起完全区分不开。这种时候不要怀疑DFT代码有bug要怀疑参数设计有问题。解决办法是增加采样时间、增大N或者换用更高分辨率的谱估计算法。7.6 信号里有直流偏置导致0频处巨大峰值很多传感器采集的信号都会有直流偏移比如加速度计的输出可能是2.5V±0.1V这时候DFT结果里X[0]会非常大把其他频率分量都衬托得看不清。处理方法是先做去直流处理x x - np.mean(x)或者在绘图时用对数坐标/只画非零以后的部分。这个操作在很多工程博客里被叫做去趋势项本质就是把公式里的k0项从结果里拿掉。7.7 纯Python手写DFT太慢怎么加快如果你需要验证小规模数据上的算法逻辑纯Python没问题但如果你用来处理工程信号建议至少使用numba加速或者直接改用numpy.fft。举个例子N1024时纯Python的双循环大概要跑几十毫秒numba版DFT可以微秒级完成numpy.fft更快。纯Python版本适合学习不适合生产环境这是很重要的一条工程经验。8. 一个完整的实验从时域到频谱图输出文章最后我放一个完整的Python脚本把前面讲到的知识点全部串起来。你可以直接复制运行也可以在此基础上做各种改动实验。import numpy as np import matplotlib.pyplot as plt def dft_basic(x): N len(x) X [0j] * N for k in range(N): sum_real 0.0 sum_imag 0.0 for n in range(N): angle -2.0 * np.pi * k * n / N sum_real x[n] * np.cos(angle) sum_imag x[n] * np.sin(angle) X[k] complex(sum_real, sum_imag) return np.array(X) fs 1000 T 1 N fs * T t np.arange(N) / fs # 构造多频信号50Hz正弦 120Hz正弦直流偏置 x 1.5 2.0 * np.sin(2 * np.pi * 50 * t) 0.8 * np.sin(2 * np.pi * 120 * t) # 计算DFT X dft_basic(x) mag np.abs(X) freqs np.arange(N) * fs / N # 绘制频谱只画前一半 half N // 2 plt.figure(figsize(10, 4)) plt.plot(freqs[:half], mag[:half]) plt.xlabel(Frequency (Hz)) plt.ylabel(Magnitude) plt.title(DFT Result (one-sided, unnormalized)) plt.grid(True) plt.show() # 单边归一化幅值 mag_norm mag[:half] * 2 / N mag_norm[0] mag[0] / N # DC 不乘2 plt.figure(figsize(10, 4)) plt.plot(freqs[:half], mag_norm) plt.xlabel(Frequency (Hz)) plt.ylabel(Amplitude) plt.title(DFT Result (one-sided, normalized amplitude)) plt.grid(True) plt.show()运行这个脚本你会看到第一个频谱图上50Hz和120Hz峰值都在几百的量级DC分量更是冲上天第二个频谱图里50Hz和120Hz的幅值分别在2.0和0.8附近和构造信号时的幅度完全吻合。这就是归一化的意义让频谱图上的读数直接对应时域信号的物理幅值。我实际调试这套代码时经常做的小实验是把120Hz改成非整数频率比如123.4Hz再看看频谱发生了什么。你会发现峰变矮了、变胖了周围长出一圈旁瓣。那就是频率泄露在作怪。再用汉宁窗去处理一遍观察谱线的变化你对DFT的理解又会深一层。DFT代码解析这件事情说难不难说简单也绝不简单。它考验的不是数学推导能力而是你能不能把公式里的上下标、求和范围、指数符号一个个翻译成循环、数组下标和实数加减。一旦你完成了这个翻译过程再把FFT、频谱分析、滤波器设计这些后续知识捡起来会发现所有东西都变得具体起来。希望这篇手把手的拆解能帮你迈过这道坎。

从公式到代码：手写DFT彻底搞懂FFT与频谱分析

相关推荐

Automatisch Disqus 集成：New Comments 与 New Flagged Comments 触发器的配置与实现原理

Kali Linux 安装全攻略：虚拟机与物理机实战指南

用 GitHub Actions 自动发布 Scalar Docs 项目：基于 @scalar/cli 的 CI 发布实践

更多资讯

手搓教程：工程师的确定性防线与AI协同方法论

高并发排行榜方案：Redis ZSET + 快照缓存 + 双TTL

工学椅选型核心：坐深、腰靠落点、坐高膝高比三大生物力学指标

2026年AI论文检测与智能降重工具全解析

WordPress主题CoreNext免授权版安装配置与安全检查指南

技术内容标题设计：原则、技巧与实战案例

今日更新

51单片机开环控制磁阻传感器的硬件匹配与代码实现

AI Agent安全操作数据库：SQL注入防御与参数化查询实战

Android逆向实战：还原激励类App签到接口的签名算法

本周精选

AI驱动的自动化专利分析平台技术解析

LeetCode-Go 题解：1252. Cells with Odd Values in a Matrix（奇数值单元格计数）双解法剖析

在线答疑系统Java毕设实战：Spring Boot前后端分离与状态流转

本月热门

基于YOLO的焊缝缺陷检测：从数据标注到工业部署全流程实践

分布式服务的决策记录

创业产品的首版范围

读完文章，想聊聊您的网站？