ARTICLE · INTELLIGENCE

战地情报 · 详情页

来自尧图项目组的一线实战观察与深度解析

DFT离散傅里叶变换详解:从公式到Python实现与频谱分析陷阱

DFT离散傅里叶变换详解:从公式到Python实现与频谱分析陷阱 你见过示波器上那条“纠缠不清”的时域波形吗十几个频率叠在一起看起来完全是一堆乱麻可同一段信号放到频谱图上每条谱线都清清楚楚地站在属于自己的位置上。想让计算机帮你完成这个“变乱为清”的过程靠的就是DFT——离散傅里叶变换。它本质上就是傅里叶变换在计算机世界里的“方言版”把连续积分改成离散求和把无限区间截成有限数组让只会做加减乘除的CPU也能听懂傅里叶在说什么。这篇文章是信号处理入门系列的第28篇我会从一个工程师的角度把DFT的公式含义、手算过程、Python实现以及几个最常见的坑一次讲透。适合刚学完傅里叶变换、却不知道怎么在代码里落地的同学也适合那些已经用过np.fft.fft、但对结果“知其然不知其所以然”的人。1. 傅里叶变换到计算机这里卡在哪一步1.1 教科书公式里的三个“计算机做不到”翻开任何一本信号处理教材傅里叶变换长这样X(f) ∫ x(t) e^(-j2πft) dt这个式子看起来优雅但你要是真打算拿它去写代码马上会撞上三个坎。第一时间变量t是连续的。计算机里没有“连续”这个概念内存里存的永远是有限个离散的数值点。你说“给我一段连续的信号”我拿到手的只能是示波器或ADC采样出来的一串数字。就好比让一个只会数豆子的人去估算一条河的水流量他没法把整条河搬回家只能一瓢一瓢地舀出来数。第二积分符号处理不了。CPU的指令集里没有“求积分”这条指令它只会做加法、减法、乘法、除法。数学上的积分在计算机里必须改写成“把无数个小碎块累加到一起”也就是离散求和。第三积分区间是负无穷到正无穷。这个更直接——计算机的内存是有限的别说无穷长时间就算是一段几分钟的高保真音频存储起来也要小心规划。所以连续傅里叶变换在计算机里是跑不起来的。这不是数学工具不好而是计算模型不匹配。1.2 计算机习惯的信号长什么样采样、截断、周期化要让傅里叶变换能在计算机里落地必须先把信号改造成计算机喜欢的样子一共三件事。第一件采样。用固定时间间隔T去读信号的瞬时值得到一串离散数字x[0], x[1], ..., x[N-1]。这一步通常由ADC完成采样间隔的倒数就是采样率fs。连续时间t变成了nT信号从x(t)变成了x[n]。第二件截断。你不可能存无穷多个采样点只能取一段有限长的数据比如取N个点。这N个点对应的时间长度是T0 N/fs也是你这次频谱分析能看到的“观察窗口”。第三件周期化。这一步很多人会忽略但它是理解DFT的关键DFT在处理你给的这N个点时会默认它们是一个周期信号的一个周期也就是说x[0]和x[N-1]是首尾相邻的信号从最后一个点跳回第一个点时被当作“本来就连着”。这个隐含假设会带来很多后续影响比如频谱泄漏后面我会专门讲。把这三件事做完连续傅里叶变换就顺理成章地变成了DFT。所以我说DFT是“计算机眼里的傅里叶”——它看到了相同的频谱本质但看的方式必须符合计算机的离散、有限、周期化的习惯。2. DFT公式里的n、k、N拆开看就不难了2.1 从连续积分到离散求和变量怎么对应DFT的标准公式长这样X[k] Σn0到N-1 x[n] · e^(-j2πkn/N)第一次看到这个式子的人往往被下标和指数吓住。其实它就是从连续傅里叶变换变形过来的每一项都有明确的来路。连续公式里时间t被替换成nTT的倒数会被频率变量的1/(NT)吸收掉最后归一化到每个采样点上就变成e的指数里那个kn/N。频率f被替换成k/(NT)其中NT就是总时长所以k越大代表的频率越高。积分符号变成从n0到N-1的累加dt被吸收进归一化系数里。换句话说n负责遍历每一个采样点时间维k负责遍历每一个待计算的频率点频率维N是总点数。你算出来的X[k]就是一个复数值既包含频率k这个分量的大小也包含它的相位。2.2 一个“频点桶”到底代表什么频率X[k]的k不是Hz它只是一个下标代表第k个“频点桶”。要把k换算成物理频率公式是f_k k · fs / N举个例子采样率fs 8000Hz点数N 1024那么k10对应的物理频率是10 × 8000 / 1024 ≈ 78.125Hz。相邻两个频点桶之间的间隔是Δf fs / N 1 / (N/fs) 1 / T0这里T0是截断后信号的总时长。注意这个关系频率分辨率Δf只跟采样时长有关跟采样率无关。你想分辨两个相差1Hz的频率观察窗口至少要1秒想分辨0.1Hz窗口就要10秒。这个不等式是写在物理规律里的后面补零那节还会再碰到它。2.3 旋转因子e^(-j2πkn/N)本质上是在绕圈公式里那个e的复数指数项初看很劝退。但只要你记得欧拉公式e^(jθ) cosθ j·sinθ就明白它其实是在复平面上画圆。随着n从0走到N-1指数项的角度2πkn/N一圈一圈地转。k越大同样走完N个点它转的圈数越多所以对应的频率越高。你可以把这件事想象成绕圆形操场跑步。一个人跑完一圈相当于最低频跑完十圈就是十倍频率。DFT干的事就是拿着一个带测速器的“标准跑者”——也就是这组旋转因子——去跟你的信号逐点做乘法、累加看你的信号里包含多少个“一圈选手”、多少个“十圈选手”。这正是傅里叶分析的精髓一个看似复杂的信号可以被分解成不同转速的“圆运动”的叠加。DFT用离散的方式完成了这个分解。3. 手算一个N4的DFT把公式踩实一点3.1 选一个特别简单的信号看出频谱怎么出来背十遍公式不如动手算一遍。我选一个最经典、也最能说明问题的信号x [1, 0, -1, 0]这N4个点其实就是cos(2πn/4)在n0,1,2,3处的取值也就是一个恰好完整采到一个周期的余弦波频率等于fs/4。选它的原因是频率正好落在DFT的整数倍频点上不会出现频谱泄漏结果非常干净。3.2 列出旋转因子表逐步算出X[0]到X[3]N4时旋转因子是W e^(-j2π/4) e^(-jπ/2) -j。它的各次幂很有规律幂次值W^01W^1-jW^2-1W^3jW^41接下来套公式。k0时e的指数是0所以X[0] 10(-1)0 0。这很合理因为x里没有直流分量。k1时 X[1] 1×W^0 0×W^1 (-1)×W^2 0×W^3 1 0 (-1)×(-1) 0 2k2时 X[2] 1×W^0 0×W^2 (-1)×W^4 0×W^6 1 0 - 1 0 0k3时 X[3] 1×W^0 0×W^3 (-1)×W^6 0×W^9 1 0 (-1)×(-1) 0 2所以DFT结果是 X [0, 2, 0, 2]。3.3 从手算结果反推物理意义这个结果里有几个非常有价值的细节。第一能量出现在k1和k3两个桶上恰好对应物理频率fs/4和-fs/4也就是镜像频率。这说明一个实信号在DFT里会被拆成正负两个频率分量这就是为什么后面要讨论单边谱。第二幅度还原的标定。X[1]的值是2但要还原原信号中这个余弦分量的真实幅度得用2/N乘以|X[k]|也就是2/4×21正好等于cos(2πn/4)的幅度1。注意直流分量X[0]不用乘2。第三X[1]和X[3]互为共轭2和2在这里恰好相等但如果信号带相位你会发现它们模相等、幅角相反。这个对称性是由实信号决定的不随N变化。我自己教新人时总让他们亲手算一遍N4的例子。算完这个DFT就不再是公式纸上的符号而是能看见的数学过程了。4. 从零写DFT并用FFT验证一组Python代码搞定4.1 最朴素的二重循环实现纸上得来终觉浅来写代码。先实现一个完全照公式来的DFT不优化、不加速纯粹为了验证思路import numpy as np def dft_manual(x): N len(x) X np.zeros(N, dtypecomplex) for k in range(N): for n in range(N): X[k] x[n] * np.exp(-2j * np.pi * k * n / N) return X这段代码就是公式的逐字翻译外层循环遍历所有k内层循环遍历所有n累加求和。缺点是复杂度是O(N²)N稍微大一点就慢得离谱。比如N4096需要算约1678万次复数乘加纯Python跑起来能等你好几秒。但它用来做教学验证非常直观。4.2 用NumPy的FFT验证我们的结果既然手写DFT慢那就拿它和NumPy内置的FFT对拍一下验证正确性N 64 fs 8000 t np.arange(N) / fs # 构造一个包含三个频率分量的测试信号 x (0.6 * np.sin(2 * np.pi * 1000 * t) 0.3 * np.sin(2 * np.pi * 2000 * t 0.5) 0.2 * np.sin(2 * np.pi * 3000 * t)) X1 dft_manual(x) X2 np.fft.fft(x) print(np.max(np.abs(X1 - X2)))实际跑出来最大误差在10的负十二次方量级基本就是浮点数本身的误差。这个结果说明两个事一是我们手写的DFT逻辑是对的二是FFT算出来的结果和DFT是同一个东西只是算法更快。4.3 完整频谱分析练习频率轴怎么画、幅度怎么标定验证完正确性就可以做正经的频谱分析了X np.fft.fft(x) freqs np.fft.fftfreq(N, 1 / fs) # 生成频率轴 mag np.abs(X) * 2 / N # 单边谱幅值标定 mag[0] / 2 # 直流分量不乘2 for f, m in zip(freqs[:N//2], mag[:N//2]): if m 0.01: print(f{f:.1f} Hz: {m:.3f})输出结果大概是1000 Hz0.6002000 Hz0.3003000 Hz0.200和构造信号时的幅度完全对得上。这里的要点是np.fft.fft的结果里实信号的能量是对称分布在整个0到fs频率区间上的画单边谱时取前半段并把幅度乘以2/N才能还原真实的余弦幅度。直流分量特殊它只出现一次不乘2。很多人第一次用FFT时画出来的幅度谱数值跟原信号对不上九成都是漏了2/N这个标定系数。5. 三个让新手翻车的DFT陷阱泄漏、栅栏、混叠5.1 频谱泄漏频率不在整数频点上的代价前面N4的例子之所以干净是因为信号频率恰好落在DFT频点桶上。但实际工程里这种情况少之又少。比如采样率1000Hz采样0.75秒那么Δf1/0.75≈1.333Hz。你给一个50Hz的正弦50并不是1.333的整数倍它的能量就会“漏”到附近一堆频点桶上频谱上出现一个鼓包而不是一根干净谱线。泄漏的根源是我在1.2节提过的“截断”取N个点相当于信号乘以一个矩形窗矩形窗的频谱是sinc函数有主瓣还有一串旁瓣旁瓣的能量叠加到信号频谱上就成了那条不干净的尾巴。处理办法最常用的是加窗。加窗就是把这N个点的两端平滑地压下去让截断没那么生硬。比如汉宁窗主瓣比矩形窗宽一点但旁瓣衰减快得多。窗函数怎么选是个权衡我做个小对比窗函数主瓣宽度旁瓣衰减适合场景矩形窗不加窗窄差频率精确落在谱线上的单频检测汉宁窗中好大多数通用频谱分析平顶窗宽很好更关心幅度精度而非频率精度布莱克曼窗宽更好需要更干净旁瓣的场合一句话频谱分析之前想清楚你是更在乎频率分辨能力还是更在乎幅度读数。没有免费的午餐。5.2 栅栏效应与补零的真实作用DFT只在有限个离散频率点上计算结果就好像你透过一道栅栏看风景只能看到栅栏缝隙里那几条风景缝隙之间的东西看不到。这叫栅栏效应。很多人听说“补零可以提高分辨率”于是把4096点信号补零到8192点再去算FFT以为能看到更精细的频谱。可惜不对。补零确实让频点间隔变小了看起来谱更“光滑”但它没有增加任何真实信息。真正的频率分辨率由信号本身的总时长决定Δf_real 1/T0。你补进去的零没有延长观测时间原信号里两个靠得非常近的频率补零后依旧分不开。补零真正的作用是“插值”——在已有频谱信息之间插出更多采样点让谱线更平滑方便你找峰值的精确位置。这就像把照片放大看但照片的清晰度并不会因为放大变高。我通常的做法是先看原N点的FFT确定大概频率范围如有必要再补零找准确峰位。别指望补零变出信号里本来就没有的分辨能力。5.3 混叠别把高频“算”成低频很多初学者做完频谱分析发现某个低频处莫名其妙出现一个大峰怎么查都查不出原因。我第一反应往往是采样率不够高频信号折叠下来了。根据奈奎斯特采样定理采样率必须大于信号最高频率的两倍否则高于fs/2的频率分量会被“折叠”回低频区。举个例子fs1000Hz时信号里有一个700Hz的成分它会被混叠到300Hz的位置。你在频谱上看到300Hz的峰其实根本不存在300Hz的信号这是采样造成的假象。这还不是最可怕的。700Hz折叠到300Hz还好发现如果采样率只是略低比如信号最高频率是520Hz折叠到480Hz看起来像是合法的低频信号极难察觉。所以工程上必须在ADC采样之前加一级抗混叠滤波器把高于fs/2的成分提前滤掉。这个滤波器是必须的不是可选的。我见过不少新手在仿真里用理想正弦测算法完全忽略混叠结果搭真实系统时频谱一塌糊涂。DFT终究算的是你给它的那些离散点离散点采歪了后面再怎么处理也救不回来。6. DFT与FFT的父子关系以及一次搜索资料时的撞名经历6.1 FFT只是DFT的一种快速算法很多人把DFT和FFT当成两个不同的东西其实不是。它们是同一个数学变换的两种实现方式DFT是定义FFT只是DFT的一种高效算法。直接按公式算DFT需要N²次复数乘加。当作FFT利用旋转因子的周期性和对称性把大DFT拆成若干小DFT计算量降到N·log2(N)这个量级。拿N1024来算直接DFT约104万次FFT约10240次差了约100倍。如果N16384差距会拉到上百倍。这也是为什么你在NumPy、MATLAB里调用的“傅里叶变换”函数内部几乎全是FFT算法。但要注意一点FFT算出来的结果跟按定义算的DFT在数学上是完全等价的不存在“FFT更准”或“FFT有误差”的说法。所差的只是速度和舍入误差的细微区别。6.2 芯片测试里的DFT完全不是一回事说到这里我必须提一段自己的经历。有次工作里需要临时查一下DFT在某种硬件加速器上的实现细节我打开搜索引擎输入“DFT”前几页出来的全是tessent dft、testmax dft、dft flow这类东西看得我一头雾水。愣了半分钟才反应过来在芯片设计领域DFT是Design for Test的缩写即可测性设计属于芯片生产测试环节的专用术语。这个撞名特别容易坑人。如果你搜资料是为了学信号处理记得在关键词里加上“信号处理”“离散傅里叶变换”或者“信号与系统”这样的限定词反过来如果你是搞芯片的搜DFT也要带上“可测性设计”才能精准命中。同一组字母两个领域的关注点完全不同一个是数学变换一个是工程方法别搞混了。文章写到这里该讲的都讲了。最后分享一个我自己的实操习惯。拿到一段数据准备做频谱分析之前我一定先做三件事第一确认采样率估算信号最高频率有没有超过fs/2第二看一眼截取时长T0算出频率分辨率Δf1/T0判断我想分辨的两个频率能不能分开第三如果信号频率不是整周期对齐主动加一个窗再跑FFT。这三件事做完基本不会再被上面那些陷阱坑到。DFT说到底不复杂难的是在动手之前把采样和截断这两个“输入条件”想清楚。希望这篇能帮你少走一点弯路。
RELATED READING

延伸阅读

更多一线实战笔记与深度复盘,助您持续精进