MATLAB 离散傅里叶变换DFT、逆离散傅里叶变换IDFT、快速傅里叶变换FFT的实现做信号处理这几年我越发觉得一个道理特别实在别人写好的fft()函数和你真正懂傅里叶变换完全是两码事。早些年我自己也偷懒遇到频谱分析就一行Y fft(x)结果一旦碰到频谱不对、相位怪异的bug整个人就懵了——因为我不知道函数背后到底做了什么。后来花了一整周把离散傅里叶变换DFT、逆离散傅里叶变换IDFT、以及快速傅里叶变换FFT在MATLAB里用笨办法从头实现了一遍才算真正打通了这条路的任督二脉。这篇东西就围绕这三件事展开先不讲调库而是用MATLAB一行一行写出DFT和IDFT再自己实现一个简化版的FFT彻底搞懂它们的数学原理和代码映射关系。之后再回头对比MATLAB自带的fft和ifft聊聊工程实践中绕不开的细节和坑。不管你是刚接触信号与系统的学生还是工作中偶尔要处理波形、振动、音频数据的工程师只要想搞明白频谱背后的机理这篇文章应该能给你省下不少弯路。1. 整体设计思路为什么非要自己写一遍DFT在动手写代码之前我必须先讲清楚一个很多教材里忽略的问题DFT到底在干什么你随便找一本信号与系统的书都会看到这样一个公式X(k) Σ_{n0}^{N-1} x(n) * e^{-j2πk*n/N}简单来说DFT做的就是一件事拿着一根固定频率的标准尺子去测量原始信号中有多少这个频率的成分。每一个k值就对应一个离散频率点把所有k都试一遍就得到了信号的频谱。这就像你在一个混合果汁里想要分别测出苹果、橙子、西瓜各占多少比例。你把每一种纯果汁当作标准模板拿它和混合果汁比对相似度高说明含量多相似度低就说明含量少。DFT的道理就是这个只不过标准模板换成了不同频率的单位复指数信号。1.1 系统需求和参数定义我在实现之前先把需求和参数定清楚不然写代码容易写到一半逻辑混乱输入信号一个有限长度的离散序列 x(n)长度记作 N。输出结果长度为 N 的复数序列 X(k)表示不同频率分量的幅度和相位。需要实现的功能正向变换从时间域 x(n) 计算频率域 X(k)。逆变换从频率域 X(k) 精确恢复 x(n)。快速算法用分治思想加速正向变换。使用工具MATLAB R2021a 以上版本无需额外工具箱。这里我想特别强调一个容易被新人忽略的问题MATLAB下标是从1开始的但DFT数学公式里的n和k都是从0开始的。这是所有初学者第一个莫名其妙的错误来源。很多人照着教材公式写代码一上来就下标越界或者结果错位。后面我的代码里都会刻意处理好这个映射关系。1.2 方案选型背后的考量搞DFT实现其实有几种路线可以选择第一种是直接调用MATLAB内置的fft()函数。这当然是最快的但得不到任何底层理解出了问题也排查不了。第二种是矩阵乘法实现。把DFT看成是一个N×N的复数矩阵乘以一个N×1的向量这在MATLAB里极其简洁一行就能写完X W * x其中W是DFT矩阵。这种写法的好处是代码酷炫、向量化程度高运行起来也不慢缺点是矩阵元素太多生成W的过程本身就等于在变相套公式。第三种是老老实实套循环。这是最笨的写法两层循环嵌套外层遍历所有k内层遍历所有n累加乘加。但恰恰是这种笨办法最直观地反映了DFT的原始定义部分和过程一目了然最适合用来理解原理。我最后选择的是第三种为主第二种为辅的方案。先用循环写一遍让数学公式和代码逐行对应再用矩阵方式重写一遍作为性能对照和验证。这样既保证了能看懂又兼顾了MATLAB作为一个矩阵实验室的天然优势。提示如果你是为了完成课程作业直接调用fft()当然没问题。但如果是为了真正掌握DFT请务必、务必先手写一遍循环版本。这不是效率问题是我踩了无数次坑后换来的体会。2. DFT核心原理与MATLAB逐步实现这一节的内容是整个文章的地基。我会从公式出发一步一步拆解然后在MATLAB里手写出来。每一个关键细节我都会说明为什么这样做。2.1 从连续傅里叶变换到离散傅里叶变换连续时间信号的傅里叶变换公式是X(f) ∫_{-∞}^{∞} x(t) * e^{-j2πft} dt但是在计算机里你没法处理连续信号只能对采样后的离散序列做变换。把积分换成求和把连续频率f换成离散频率点k就得到了DFT的公式。这里有一个关键直觉DFT的输入其实默认是一个周期的周期信号。你拿到一段长度为N的序列DFT会在心里默认为它循环重复。这也直接解释了后面要讲的频谱泄漏问题——当你取到的信号截断时周期延拓会在边界处引入突变频谱就会被抹脏。2.2 手写DFT的完整代码这是我认为全文最重要的代码段之一请逐行看function X my_dft(x) % MY_DFT 手写离散傅里叶变换 % 输入 x: 长度为N的行向量或列向量 % 输出 X: 对应DFT结果的复数向量与x同长同形状 N length(x); X zeros(size(x)); % 预分配保持维度统一 % 确保按列向量处理方便循环逻辑 if isrow(x) x x.; end for k 0:N-1 sum_val 0; for n 0:N-1 sum_val sum_val x(n1) * exp(-1j * 2 * pi * k * n / N); end X(k1) sum_val; end % 如果输入是行向量输出转置为行向量 if isrow(x) 0 size(X,1) 1 X X.; end end让我快速解读几个容易踩坑的地方第一个关键是下标偏移。数学里写的是x(n)n从0开始MATLAB里却是x(n1)因为下标1才是第一个元素。输出X(k1)同理。这个偏移虽然只是加1的问题但在实际写代码时极其容易出错。第二个关键是预分配。X zeros(size(x))这一步看起来可有可无但如果你在MATLAB循环里不断扩展数组数据量大的时候运行速度会急剧下降。养成预分配的好习惯是在MATLAB里写循环的必备意识。第三个关键是复数指数直接用exp(-1j*2*pi*k*n/N)。MATLAB里1j就是虚数单位别把它定义成变量。另外注意公式里是负号这个负号决定了你是在做正变换还是逆变换一旦写反结果将是错误的。2.3 验证手写DFT的正确性写完不能就算完必须验证。我个人的习惯是先用已知结果的场景做测试而不是直接上真实信号。比如说我构造一个纯正弦信号频率恰好落在DFT的某个频率点上。理论上它的频谱应该只有一个峰值其余都是0附近的小数。% 构造测试信号 fs 1000; % 采样率 1000 Hz N 64; % 采样点数 n 0:N-1; f0 125; % 信号频率 125 Hz x sin(2*pi*f0*n/fs); % 125Hz正弦波 % 计算DFT X1 my_dft(x); X2 fft(x); % MATLAB内置函数作为对照 % 对比误差 err max(abs(X1 - X2)); fprintf(最大误差: %e\n, err);这段代码跑下来理论上误差应该在10的负12次方量级这完全来自于浮点数的舍入误差。如果误差很大说明你的手写实现有问题。这个对比思路是调试的黄金方法——永远用内置函数当作验证参照物因为它经过万千用户的检验基本可以信任。2.4 计算复杂度分析与性能痛点如果你写过上面的代码就一定感受过那个酸爽的等待过程。当N64时全无压力但当你把N提高到1024时两层循环的时间就开始明显变长N4096时肉眼可见地慢N16384时你已经开始怀疑是不是死循环了。为什么会这样因为DFT直接实现的时间复杂度是O(N²)。每一对(k, n)都要做一次复数乘法总计算量是N²次。N4096的话就是1677万次复数乘加在MATLAB纯循环环境里这真的够喝一壶的。这就是FFT存在的根本原因。FFT把计算量从O(N²)降到了O(N log₂N)。N4096时DFT需要1677万次乘法FFT只需要大约4096×12≈5万次这是300多倍的差距而且N越大差距越恐怖。3. IDFT与信号重建从频域回到时间域前面把信号从时间域变换到频率域之后经常还需要把它变回来。逆离散傅里叶变换IDFT干的就是这个活。3.1 IDFT的数学原理IDFT的公式长这样x(n) (1/N) * Σ_{k0}^{N-1} X(k) * e^{j2πk*n/N}注意三个关键点第一复指数上的符号变成了正号。正变换是负号逆变换是正号这是数学上的定义不能搞混。第二前面多了一个1/N。这个系数起到归一化的作用。你可以这样理解正变换把信号的能量分布到N个频率分量上逆变换再把它们合起来如果不除以N恢复出来的信号幅度会放大N倍。第三求和范围和正变换一样都是从0到N-1。我第一次写IDFT的时候就因为在公式里忘了乘1/N导致恢复出来的波形幅度全是错的。这种错误非常隐蔽因为谱形状看起来几乎一样只有幅度比例不对。3.2 手写IDFT的MATLAB代码function x my_idft(X) % MY_IDFT 手写逆离散傅里叶变换 % 输入 X: 频域复数向量 % 输出 x: 恢复的时间域信号 N length(X); x zeros(size(X)); xvec X(:); % 强制变为列向量 for n 0:N-1 sum_val 0; for k 0:N-1 sum_val sum_val xvec(k1) * exp(1j * 2 * pi * k * n / N); end x(n1) sum_val / N; end x reshape(x, size(X)); % 恢复形状 end这段代码和DFT那段的区别只在指数符号和多了1/N。所以写代码的时候有一种很好的检查方式把IDFT代码里的指数符号改回去、去掉1/N它就应该能当DFT用。3.3 用仿真信号验证变换-逆变换的往返一致性做信号处理的人常说一句话变换一致性的金标准是往返测试。什么意思你把x做DFT得到X再用X做IDFT得到x_recovered理想情况下x_recovered必须和x一模一样。如果中间有任何差错这个测试立刻暴露。% 构造一个复杂一点的测试信号两个频率叠加加噪声 fs 4096; N 128; t (0:N-1)/fs; x 1.5*sin(2*pi*300*t) 0.8*cos(2*pi*1200*t) 0.2*randn(size(t)); % 正向变换 X my_dft(x); % 逆向重建 x_rec my_idft(X); % 计算误差 max_abs_err max(abs(x - x_rec)); rms_err sqrt(mean((x - x_rec).^2)); fprintf(最大绝对误差: %e\n, max_abs_err); fprintf(均方根误差: %e\n, rms_err); % 可视化对比 figure; subplot(2,1,1); plot(t, x, b-, LineWidth, 1.5); hold on; plot(t, x_rec, r--, LineWidth, 1.5); legend(原始信号, IDFT重建信号); title(原始信号与IDFT重建信号的对比); grid on;正常情况下最大误差应该在10的负13次方量级左右。这个误差纯粹来源于浮点数舍入完全在可接受范围内。如果误差到了10的负3次方甚至更大那就要回去查代码了。3.4 归一化的习惯与信号能量验证除了直接做往返测试以外还有一个更数学化的验证方式帕塞瓦尔定理。这个定理说的是信号在时间域的总能量等于在频率域的总能量具体到离散情况就是Σ|x(n)|² (1/N) Σ|X(k)|²用代码验证就是energy_time sum(abs(x).^2); energy_freq sum(abs(X).^2) / N; fprintf(时域能量: %f\n, energy_time); fprintf(频域能量: %f\n, energy_freq);这两个数值应当几乎完全相等。如果不相等最常见的原因就是正变换和逆变换的系数分配不一致。注意MATLAB内置的fft和ifft这对函数其实内部已经帮你处理好了归一化——ifft(fft(x))严格等于x。但如果你自己实现就必须记住1/N这一项。我建议大家在代码里用注释标出系数1/N这个位置防止以后回来读代码时犯迷糊。4. FFT的原理剖析与代码实现理解了DFT和IDFT是怎么用最笨的办法实现之后我们现在来谈FFT。FFT不是一种新的变换而是DFT的一种高效算法它和DFT的数学结果完全一致只是计算速度更快。4.1 FFT的核心思想分治与大而化小FFT最基本的版本是库利-图基算法它的核心思想可以浓缩成八个字奇偶分解递归分治。如果你有一个长度为N假设N是2的幂的序列x(n)你可以把它拆成偶数下标的部分x(0), x(2), x(4)...和奇数下标的部分x(1), x(3), x(5)...。然后用数学推导可以证明原来长度为N的DFT可以被表示成两个长度为N/2的DFT组合出来的结果。具体公式是这样的X(k) X_even(k) W_N^k * X_odd(k) X(k N/2) X_even(k) - W_N^k * X_odd(k)这里的X_even和X_odd分别是偶数序列和奇数序列的N/2点DFT而W_N exp(-j*2π/N)就是所谓的旋转因子或twiddle factor。这个过程可以一直递归下去直到序列长度变成1。长度为1的DFT就是它本身不需要任何计算。这也就是为什么FFT对序列长度要求是2的幂——因为只有这样才能不断二分下去直到长度为1。我打一个比方你有一堆积木拼成了一个复杂城堡你想数清楚里面有几块积木。笨办法是从头到尾一块一块数聪明的办法是先把城堡分成左半边和右半边每半边再分成上下两层每一次细分都只要把两个部分的数量简单相加。FFT的分解就是这个感觉。4.2 MATLAB内置FFT的使用与频率轴建立MATLAB的fft函数在工程中使用频率极高。最基础用法是这样x randn(1, 1024); X fft(x);如果x的长度不是2的幂MATLAB也不会报错。它会自动在尾部补零凑到一个合适的长度算法效率会相应变化。但如果你自己要把手写FFT和内置对比建议用2的幂长度否则内置函数可能是在用混合基算法你的递归实现没法一一对应。建立频率轴是另一个高频踩坑点。很多初学的人拿到X fft(x)之后直接用plot(abs(X))画图得到的是一个从0到N-1的频率索引横轴。但工程中咱们需要真正的物理频率单位Hz就需要知道采样率fs然后建立对应的频率矢量fs 1000; N 1024; f (0:N-1) * fs / N; % 频率轴单位Hz这段代码的本质在于DFT输出的N个点均匀分布在0到fs之间不包含fs这一点每个相邻频率点之间的间隔是fs/N。这个频率间隔我们专门有个名字叫做频率分辨率。频率分辨率越细意味着你越能分辨两个靠得很近的频率成分。但有一个细节你得特别小心真实的频谱其实只有前N/2个点是有用的因为后面的一半是前面一半关于fs/2的镜像。对于实信号后半段频谱和前半段是共轭对称的没有额外的信息量。所以你画幅频图时一般只画前N/21个点f_plot f(1:N/21); mag_plot abs(X(1:N/21)); plot(f_plot, mag_plot);如果你习惯把频率轴的范围设为从-fs/2到fs/2画出双边谱来观察对称性那就需要用到fftshift函数。这个函数的作用是把数组的左右两半对调让零频跑到中心f_shift (-N/2:N/2-1) * fs / N; X_shift fftshift(X); plot(f_shift, abs(X_shift));4.3 自己动手实现一个递归版FFT写递归版FFT不仅是验证理解的绝佳方式而且递归结构的代码和算法思想本身高度对应。以下是我写的版本function X my_fft(x) % MY_FFT 手写快速傅里叶变换库利-图基算法递归版 % 要求输入长度为2的幂 N length(x); % 递归基长度为1DFT就是本身 if N 1 X x; return; end % 保证N是2的幂 if mod(log2(N), 1) ~ 0 error(输入长度必须为2的幂); end % 偶数和奇数下标分组 % 注意MATLAB下标从1开始数学中下标从0开始 x_even x(1:2:end); % 数学下标 0,2,4,... x_odd x(2:2:end); % 数学下标 1,3,5,... % 递归计算两个子DFT X_even my_fft(x_even); X_odd my_fft(x_odd); % 组合阶段 X zeros(1, N); for k 0:N/2-1 W exp(-1j * 2 * pi * k / N); % 旋转因子 X(k1) X_even(k1) W * X_odd(k1); X(k1N/2) X_even(k1) - W * X_odd(k1); end end这个实现有几点值得说明。第一递归基是N1。长度为1的序列的DFT就是它自身这是所有递归的出口。没有这个条件程序会无限递归崩溃。第二分组用的是x(1:2:end)和x(2:2:end)。这里直接体现了MATLAB下标从1开始的实际效果——数学中的偶数下标0、2、4对应MATLAB下标1、3、5奇数同理。第三组合循环只需要从0到N/2-1。因为FFT利用的是DFT结果的周期性X_even(k N/2) X_even(k)旋转因子W_N^(kN/2) -W_N^k。这个对称性让一个长度为N/2的DFT的每个值都能被用到两次。4.4 验证手写FFT和内置fft及手写DFT的一致性写完了当然要验证。我用三种方法计算同一个序列的DFT对比它们的结果N 512; rng(42); % 固定随机种子可复现 x randn(1, N) 1j * randn(1, N); % 复信号更全面 % 三种方法 X_loop my_dft(x); % 原始DFT逐循环 X_rec my_fft(x); % 手写FFT X_built fft(x); % 内置FFT % 输出误差 err1 max(abs(X_loop - X_built)); err2 max(abs(X_rec - X_built)); fprintf(DFT循环 vs 内置fft: %e\n, err1); fprintf(手写FFT vs 内置fft: %e\n, err2); % 顺便计时比较 % 对于N4096这种规模差距会非常明显 N_big 4096; x_big randn(1, N_big); tic; X1 my_dft(x_big); t1 toc; tic; X2 my_fft(x_big); t2 toc; tic; X3 fft(x_big); t3 toc; fprintf(手写DFT耗时: %.4f秒\n, t1); fprintf(手写FFT耗时: %.4f秒\n, t2); fprintf(内置FFT耗时: %.4f秒\n, t3);实际跑下来当N4096时手写DFT可能需要几秒而手写FFT通常在几十毫秒内就能完成内置FFT更是快到毫秒以下。这种速度差距不是代码技巧能弥补的纯粹是算法复杂度决定的。4.5 FFT的工程细节补零、窗函数与频谱泄漏既然讲到FFT的工程应用了这里必须提几个绕不开的实战问题。补零是一种常见的操作。有时候你的采样点数只有500个但你想让频谱看起来更平滑可以在后面补零到1024个点。但必须搞清楚一个重要概念补零不会提高频率分辨率。它只是让FFT输出点更密看起来更平滑但真实的分辨率仍然由原始数据长度决定即fs/N。如果你想真正提高频率分辨率只有增加采样时长这是物理事实不是数学技巧。频谱泄漏是另一个经典问题。当你截取一段有限长度的信号时如果截断窗口内包含的信号周期数不是整数那么信号在周期延拓时就会在边界处发生跳变导致频谱中出现不该有的成分。这就是频谱泄漏。举个例子一个频率为100Hz的正弦波采样率1000Hz如果你恰好采样10个点就是完整一个周期。但如果你采样9个点或11个点边界处就会断裂频谱就不再是干净的一根线。应对泄漏的标准做法是加窗。给原始信号乘上一个窗函数比如汉宁窗、海明窗让信号两端平滑地衰减到零这样边界跳变就消失了泄漏被大幅抑制N 1024; n (0:N-1); x sin(2*pi*100*n/1000); % 100Hz正弦波但N非整数周期 w hanning(N); % 汉宁窗 xw x .* w; F abs(fft(xw));加窗也有代价主瓣变宽了频率分辨率进一步变差而且信号的幅度会被压制因为窗的加权平均不等于1。工程上这是一个经典的权衡问题。5. 三个变换的横向对比与工程选型写到这里我把DFT、IDFT、FFT三者放在一起做个系统的横向对比帮助你在实际应用里快速决策该用哪个。对比维度DFT手写循环IDFT手写循环FFT快速算法数学本质时域到频域的正变换频域到时域的逆变换DFT的高效算法时间复杂度O(N²)O(N²)O(N log₂N)输入长度要求任意任意通常为2的幂代码复杂度简单直观简单直观递归/迭代稍难理解数值误差与内置fft约10⁻¹²量级与内置ifft约10⁻¹²量级与内置fft约10⁻¹²量级使用场景教学、原理验证教学、时域重建实际工程中的首选这个表格可以读出几个关键结论。先说应用场景的取舍。在工程项目里你永远不会用我的my_dft去处理几万点的数据计算效率完全不能接受。但在教学、调试和验证核心逻辑时手写DFT的价值无可替代——它让你能够清晰地看到每一步中间结果定位问题。再说精度问题。很多人担心FFT因为采用了递归拆分会比直接DFT的精度差。实际上两者在数学上是完全等价的只存在浮点数舍入误差的细微差异。对于绝大多数工程应用来说这个差异完全可以忽略。最后说一下内存占用。直接DFT如果用矩阵乘法实现需要一个N×N的复数矩阵当N增长时内存占用呈平方级暴涨。而FFT利用蝶形运算只需要O(N)的存储空间。这在嵌入式等资源受限的场景下是生死攸关的区别。我自己的习惯是教学演示用手写循环工程开发用内置fft/ifft性能敏感且数据量极大时考虑C扩展或GPU加速。这个决策链条基本覆盖了90%以上的场景。5.1 各维度细微差异与选择建议可能有读者会问既然都推荐使用内置fft了那我为什么还要费劲自己写我要说这个费劲是有回报的。当你真正手动实现过一遍FFT你再去看频域分析的结果脑子里会出现一个清晰的计算过程图景哪些点对应哪个频率为什么会有对称性为什么要除以N补零对频谱做了什么——这些疑问会像拼图一样快速拼整。而在实际的MATLAB开发里我认为应该灵活使用两者如果你在做频谱分析、滤波等常规任务直接内置fft不要自己造轮子。如果你在做教学、做算法验证或者需要修改频域处理流程的某个环节那你绝对需要手写版本当作参考基准。如果你的数据量极大到了几百万点级别建议考虑分块处理或者专用硬件加速纯MATLAB里跑fft虽然也很快但内存开销不可忽视。6. 手写DFT/FFT过程中的常见错误与排查技巧最后这一部分我专门整理我在反复实践中遇见过的典型错误和排查心得。这些内容在教科书里通常不会写但实操中几乎人人都会踩。6.1 典型错误速查表错误现象可能原因排查方法输出结果整体偏移频谱左右颠倒下标从0还是1没处理好检查循环中n和k的起始值核对MATLAB的x(n1)映射逆变换恢复的信号幅度放大N倍IDFT忘了除以N确认逆变换公式是否包含1/N系数频谱峰值不在预期频率处频率轴建立错误采样率和点数没对应检查f (0:N-1)*fs/N 是否正确频谱两边对称且包含镜像实信号的固有性质不是bug只取前N/2点或使用fftshift手写FFT结果和内置fft完全不同数组分组用错MATLAB的奇偶分组方式重点检查x(1:2:end)和x(2:2:end)频谱整体被抹脏主峰附近有许多杂散频率频谱泄漏截断边界突变加窗函数抑制或调整采样长度为整数周期计算速度慢到无法忍受使用了O(N²)的循环DFT改用FFT或矩阵向量化6.2 我的调试心法构造解析信号对照排查DFT问题最有效的方法是构造一个你事前知道答案的信号去测试。比如构造一个纯正弦信号频率正好落在整数倍的基频上。这样理论上DFT结果应该只有一两个非零点其余全部为零。这是最直观的判别式测试。当初我写第一个DFT时就是因为期望的峰值位置和实际输出差了整整一个格子查出是下标映射的问题。这种用已知信号做测试的方式比拿真实复杂信号去人肉判断看起来对不对要可靠得多。6.3 用内置函数对照定位错误内置的fft函数是你的最佳参考答案。当手写结果和内置不一致时我建议你用这样的步骤来定位第一取一个最小的测试向量比如N4手动计算一遍DFT出来。然后分别跑my_dft和fft逐元素对比看是哪个点出错了。4点DFT的数学计算完全可以手推这样每一步都能对照。第二对比过程中如果发现某个k点的值不对就单独把那个k点对应的循环拉出来打印中间的累加值一步步检查和手算是否一致。第三如果手写FFT出错以迭代式调试为主——写一个小的递归深度打印函数观察分组和组合阶段的中间结果。内存和时间开销的想法我再多提一嘴遇到性能瓶颈最优先检查数据类型。如果数据是单精度singleMATLAB内置FFT会自动优化对应的算法路径。如果你的代码里不小心把复数数组转成了double性能也会受到细微影响。6.4 关于代码风格与可复现的几个建议最后分享几个文本代码层面的小习惯。变量命名要有语义。my_dft、my_fft、my_idft这样的名字清楚明白。函数头部写注释标明输入输出的格式以及算法来源。每个关键公式旁边写注释引用对应的数学表达式。重要脚本里用rng固定随机种子。这样每次跑出来结果一致调试和复现都方便。我记得有次排查一个随机噪声干扰问题因为没固定种子每次实验结果都不一样排查花了好几天后来加上rng(42)问题立刻暴露了。单步调试时善用断点。在循环内部设置断点可以检查n、k、exp指数等中间变量的数值是否合理。尤其是复数部分的实部和虚部分开看常常能快速发现问题。7. 从原理到工具箱完成一次完整的频谱分析小任务作为一个综合练习我把上面所有内容串成一个完整的小任务。这也是我认为读者可以从零复现的最有价值的一段流程。场景是这样的我们有一组振动传感器的采样数据采样率fs为2000Hz一共采集了2048个点。传感器信号中包含一个150Hz的机械振动分量和一个大约350Hz的电气干扰分量同时叠加了随机噪声。我们希望通过频谱分析把两个主要频率成分找出来并且测量它们的幅度。7.1 生成模拟信号fs 2000; % 采样率 2000 Hz N 2048; % 采样点数 t (0:N-1)/fs; % 时间向量 % 真实信号150Hz幅度1.2350Hz幅度0.8叠加白噪声 rng(7); x 1.2*sin(2*pi*150*t) 0.8*sin(2*pi*350*t pi/4) 0.15*randn(1,N);7.2 使用手写FFT进行频谱分析X my_fft(x); % 使用我们手写的FFT % 频率轴 f (0:N-1)*fs/N; % 取正频率部分 half N/2; f_plot f(1:half1); mag_plot abs(X(1:half1)) * 2 / N; % 幅度归一化 % 找到峰值 [pks, locs] findpeaks(mag_plot, MinPeakHeight, 0.5); fprintf(检测到的频率分量:\n); for i 1:length(locs) fprintf(频率 %.2f Hz幅度 %.4f\n, f_plot(locs(i)), pks(i)); end这段代码里有几个地方特别需要注意。第一幅度归一化。FFT输出的幅度和信号本身不是直接相等的。对于单频正弦信号来说它的峰值幅度的FFT值大约是幅度乘以N/2。所以要恢复原始幅度信息需要乘以2再除以N。这是我见过初学者问得最多的问题之一。第二频率分辨率是fs/N ≈ 0.9766Hz这个分辨率下150Hz和350Hz的频率点都能被准确测量只要它们不在同一个频率点上。7.3 结果的解释与检查实际跑完这段代码输出的频率和幅度会非常接近真实值。150Hz的幅度读出来大概是1.19几分350Hz的大概0.79几分和设定的1.2、0.8相当接近。误差主要来自噪声和频谱泄漏的残余影响。如果我想进一步抑制噪声影响提高谱线质量我可以在做FFT之前对信号加窗。加窗后主峰处的幅度会被窗函数拉低因此需要额外的幅度修正因子。比如汉宁窗的幅度恢复因子大约是2.0如果用了汉宁窗就要把幅度再乘以2。我在实际项目里的操作习惯是这样先用不加窗的原始FFT快速看个大概找到主要频率区域如果需要精确定量某几个成分的幅度和相位再用加窗处理并做相应的幅度修正。8. 写在最后的经验分享如果你从头到尾把代码自己敲了一遍我相信你对DFT、IDFT、FFT的理解已经和只会调用fft的人拉开了实质性的距离。我在实际使用中的体会是手写变换的价值不在于替代内置函数而在于获得一种掌控感。当频谱出现奇怪的形状、相位出现异常偏移时你能很快判断出是信号的物理问题还是算法的处理问题而不至于一头雾水。这套代码我后来在几个模拟项目中反复使用过无论是机械振动频谱分析、电网谐波检测还是语音信号的短时频域处理对底层计算的透彻理解都帮了大忙。即使最终交付的代码全部换成了高效的内置实现调试过程中手写版本的参考价值依然无可替代。如果你有精力我建议你再把这个练习向前推一步尝试把递归版FFT改成迭代版基2蝶形算法画出蝶形运算图观察数据重新排列的规律。这个过程会让你对FFT的底层有个更深入的体感。我第一次画出数据重排的表时有种原来如此的通透感。动手试一遍比读十遍公式都管用。