1. 为什么我要自己实现 FFT补零和幅值修正这两个细节决定成败最近在做一个振动特征提取的预研采集器每帧只给我 250 个采样点但后端的硬件 FFT 模块要求输入长度必须是 2 的幂。这就免不了要把数据补零到 256 点再送进处理器。你可能会说直接调用 MATLAB 的 fft() 不就行了吗问题是自编的 FFT 在嵌入式侧验证时更清晰因为你能一步步看到补零发生在哪、每级蝶形的旋转因子怎么引入、最后幅值修正该除以哪个 N。折腾下来我发现真正影响频谱质量的并不是 FFT 算法本身而是补零和幅值修正这两个容易被忽略的环节。这个场景适合几类人看正在把 MATLAB 里的 FFT 算法改写成 C/MEX 的工程师需要理解补零为什么能平滑频谱的科研人员以及一直被频谱幅值对不上实际信号幅值困扰的初学者。我会把自编 radix-2 FFT 的完整思路、补零的取舍、频谱图和相位图的绘制判读一起串起来讲代码可以直接复制到 MATLAB 里跑。1.1 数据不足位补零到底补的是什么所谓“不足位补 0”在 MATLAB 里写成两行N_actual length(x); nfft 2^nextpow2(N_actual); x_padded [x(:); zeros(nfft - N_actual, 1)];第一行取原始长度第二行找一个不小于该长度的 2 的幂。这里要清楚补零不是给你凭空增加有效信息而是在频域做插值。原始 250 点的 DFT 结果有 250 条谱线补零到 256 点之后变成 256 条谱线。多出来的谱线并没有让原始信号的物理分辨率提高它们只是把原来谱线之间的空隙用 sinc 插值填上了。实际采集时长不变的情况下补零后的频谱更平滑通过搜索峰值也更容易读准频率。如果非要较真“补零能不能分辨两个频率很接近的信号”答案是否定的。物理频率分辨率由实际采样时长决定也就是 fs/N_actual。补零只能让谱线更细密但两个靠得很近的峰值该糊成一团还是糊成一团。1.2 幅值修正是自编程里最大的坑直接用 abs(fft(x)) 画频谱画出来的峰值往往跟你输入的振幅对不上。原因很简单MATLAB 自带的 fft() 不自动除以 N三角函数变换的结果默认分布在正负频率两个方向上。以幅度为 2 的余弦信号为例单侧频谱上峰值修正后应当是 2不做修正你会看到一个接近 1 的数值还要乘一个 2 才是真正的单边幅值。自编 FFT 也一样内部每个蝶形只是累加和旋转因子相乘没有任何归一化。所以最终一定要做两步mag_raw abs(X); mag_one_sided mag_raw / nfft; mag_one_sided(2:end-1) mag_one_sided(2:end-1) * 2;这么处理之后你才能在频谱图上直接读出信号的真实幅值。我在下面第 4 节会详细拆解公式这里先记住一个原则——幅值修正是“除以 FFT 点数、单边谱补乘 2、直流和奈奎斯特频率不乘 2”。2. 自编 FFT 的实现思路从 DFT 公式到位反转和蝶形运算FFT 本质上是把 DFT 的运算拆成一层层蝶形。DFT 的原始定义是X(k) sum_{n0}^{N-1} x(n) * exp(-j2pikn/N)直接实现要算 N^2 次复乘N 到 256 可能无所谓但工程里 N 经常是 1024 甚至 4096开销就大了。radix-2 FFT 利用旋转因子的周期性和对称性把计算量降到 N*log2(N) 的量级。2.1 我自己采用的递归 radix-2 实现自编 FFT 在 MATLAB 里有很多写法。为了可读性我选择了 DIT按时间抽取递归版本核心只有十行左右function Y my_fft_radix2(x) % 递归实现 radix-2 按时间抽取 FFT N length(x); if N 1 Y x; return; end x x(:); % 抽取偶数位和奇数位 x_even my_fft_radix2(x(1:2:end)); x_odd my_fft_radix2(x(2:2:end)); % 旋转因子 Wn exp(-1j * 2 * pi * (0:N/2-1) / N); % 蝶形运算 Y [x_even Wn .* x_odd; x_even - Wn .* x_odd]; end这段代码每次把序列一分为二偶数下标的序列和奇数下标的序列分别递归调用自己再用蝶形公式合并。旋转因子 Wn 就是单位圆上等间隔取的相位点。Wn .* x_odd是把奇数序列先做一次相旋再与偶数序列相加。2.2 外层补零封装函数实际工程里我不希望手写补零逻辑散落在各处所以又套了一层接口名字也叫 myfft和 MATLAB 内置 fft 保持相同调用习惯function [X, f] myfft(x, fs, nfft) % myfft: 自编 FFT不足 nfft 位自动补零 % 输入: % x: 时域信号 % fs: 采样率 % nfft: FFT 点数必须是 2 的幂若缺省则自动取 nextpow2 if nargin 3 nfft 2^nextpow2(length(x)); end x x(:); if nfft length(x) error(nfft 不能小于输入数据长度); end if mod(log2(nfft), 1) ~ 0 error(自编 radix-2 FFT 要求 nfft 为 2 的幂); end % 核心不足位补 0 x_pad [x; zeros(nfft - length(x), 1)]; % 调用 radix-2 递归 FFT X my_fft_radix2(x_pad); f (0:nfft-1) * fs / nfft; end这里有个测试细节如果 nfft 由2^nextpow2()自动生成肯定满足 2 的幂条件但封装函数里最好还是加一个幂判断防止外部传参时手误传了 100 或者 300 这样的数。递归实现会直接算错而且很难察觉。2.3 递归 vs 迭代到底选哪个递归版本代码短最适合讲解和验证原理但性能一般。实际嵌入式里通常用迭代 DIT。迭代 DIT 需要先把输入序列按二进制位反转重新排列然后一层层蝶形合并。同样的逻辑可以这么写function Y my_fft_iterative(x) N length(x); xd bitrevorder(x); % 也可以用自己写的位反转 Y xd; len 1; while len N half len; len len * 2; for k 0:half:N-1 for j 0:half-1 w exp(-1j * 2 * pi * j / len); u Y(k j 1); v Y(k j half 1) * w; Y(k j 1) u v; Y(k j half 1) u - v; end end end end位反转是迭代 DIT 里最麻烦的部分。MATLAB 有 bitrevorder 可以直接用但在嵌入式环境里要自己实现比如对下标做二进制逆序。这就是为什么我建议先用递归版本把算法流程吃透等移植到效率要求高的场景再换迭代版本。3. 补零的本质频谱插值而不是提高频率分辨率很多初学 FFT 的人都误以为补零到 1024 点频率分辨率就变成 fs/1024。这是不对的。频率分辨率的物理下限是实际有效信号长度决定的 fs/N_actual补零之后虽然谱线间距变小但谱线的包络还是原来的包络。我习惯把它类比成图像放大放大后的像素变多了你可以更精细地看轮廓但没法看清照片里原本就没有的细节。3.1 分辨率与谱线间距的对照还是用刚才的 250 点例子fs1000Hz实际分辨率 df_real 1000/250 4Hz补零到 256 点后谱线间距 df_padded 1000/256 3.90625Hz。这两者看起来很接近但含义完全不同。如果另一个信号的频率比你当前信号低 4Hz只有真实采集 250 点以上才能区分出来补零再多也只是让已有峰的形状更光顺不会把两团靠在一起的峰拆开。我可以再做一组对照实验fs1000Hz分别采集 0.25s 和 0.5s两个时长都补零到 1024 点。0.25s 的频谱虽然也画出了 1024 条线可一旦两个正弦间隔小于 4Hz它绝对区分不开而 0.5s 的频谱可以轻松分辨之前无法区分出的接近频率。这说明什么想要更好的频率分辨率请去延长采样时间而不是一味加大 nfft。3.2 nfft 到底应该设成多少实际项目中我的经验分两种如果仅是为了满足 radix-2 长度要求nfft 2^nextpow2(length(x)) 就够。如果是为了让频谱图线条更平滑、峰值搜索更精确可以在原始点数基础上乘 4 或 8但不要几百倍地补。因为补零并不是不花代价。补零太多会把峰值附近的旁瓣画得特别细人眼觉得“变好看了”但旁瓣泄漏也会被插值成连续起伏反而干扰峰值判读。而且嵌入式场景里 nfft 越大存储和计算量增长越明显。补零时还有一个不起眼的点如果原始信号的 DC 分量不为 0x_pad尾部突然从非零变成 0会在频域产生一个小冲击成分。所以必要时应先减掉均值再补零或者把均值单独换算到 0Hz 谱线中。4. 幅值修正与相位还原不能只看 abs(FFT) 就结束自编 FFT 也好MATLAB 内置 fft 也好默认输出都是复数。很多人画频谱图时只取了幅值把 phase 丢在一边结果到了需要判断时间延迟、滤波器群延时或相位噪声时又重新采数据。所以我在项目里固定要求频谱图、相位图一起画提前习惯就好。4.1 单边谱的幅值修正公式标准的修正流程四步nfft length(X); % 实际用的 FFT 点数包含补零后的点数 halfN floor(nfft/2) 1; % 单边谱有效索引范围 f_one f(1:halfN); mag_one abs(X(1:halfN)) / nfft; % 除了 DC(索引1) 和 Nyquist(索引 halfN) 外其余频率点要乘 2 mag_one(2:halfN-1) mag_one(2:halfN-1) * 2; phase_one angle(X(1:halfN));为什么是单边谱因为实信号经过傅里叶变换后正频率和负频率是对称的。我们把幅度按正负各半重新合并到单边所以除非直流或奈奎斯特频率都需要乘以 2。幅值修正的细节用一张表看更直观谱线位置修正方式原因DC0 Hzabs(X(1)) / nfft直流分量只出现在 0Hz不参与正负频率折叠中间频率2 * abs(X(k)) / nfft正负频各占一半需要合并Nyquistfs/2abs(X(halfN)) / nfft奈奎斯特频率同样只有一个4.2 相位图和 angle() 的坑相位直接用 angle(X) 得到的是主值范围在 [-pi, pi]。频谱图上会看到很多 2pi 跳变那不是相位突变而是角度被截断在了这个区间。要想看到线性相位趋势要用 unwrap() 展开phase_unwrap unwrap(angle(X(1:halfN))); plot(f_one, phase_unwrap);还有一个常见问题某些 bin 上信号能量接近 0 时幅值很小相位会变成一堆杂乱的角度看起来毫无规律。这是数值噪声导致的。判读相位时只关心幅值峰值附近的谱线别把整个相位图都当成有效信息。工程上我通常会先找谱峰再在当前峰值周围 ±2 根谱线范围内看相位一致性。5. 频谱图和相位图的绘制与判读把这个过程的完整脚本放在一起。为了演示构造一个 60Hz、幅值 2 的余弦信号采样率 1000Hz采集 250 点然后用自编 myfft 补零到 256 点fs 1000; N_actual 250; f0 60; A0 2; n 0:N_actual-1; x A0 * cos(2*pi*f0*n/fs); nfft 2^nextpow2(N_actual); [X, f] myfft(x, fs, nfft); halfN nfft/2 1; f_one f(1:halfN); mag_one abs(X(1:halfN)) / nfft; mag_one(2:halfN-1) mag_one(2:halfN-1) * 2; phase_one unwrap(angle(X(1:halfN))); figure; subplot(2,1,1); plot(f_one, mag_one, LineWidth, 1.2); xlabel(频率 (Hz)); ylabel(幅值); title(频谱图幅值已修正); grid on; subplot(2,1,2); plot(f_one, phase_one, LineWidth, 1.2); xlabel(频率 (Hz)); ylabel(相位 (rad)); title(相位图); grid on;运行之后会看到两个子图上图在 60Hz 附近出现一个明显峰值修正后峰值接近 2下图相位在峰值附近保持一条接近线性的曲线远离峰值的地方相位跳变得比较乱。这里需要说明一点由于 250 点采集 60Hz 并不是整个 256 点的整数周期所以在 60Hz 两侧会出现若干旁瓣这是频谱泄漏不是自编 FFT 算错了。如果想知道精确幅值在加窗前通常取峰值附近的包络或改用整周期采样把 f0 设计成 bin 的整数倍。下面把从图里读信息的几个小经验分享给你峰值的横坐标往往不是正好落在 f0。如果信号存在泄漏峰值位置会偏离真实频率这时用二次插值法parabolic interpolation在峰值附近拟合比直接读谱线更准。相位图只需要看峰位附近的局部。设计滤波器时用频率轴上所有有效带的相位响应查群延时。如果只在 0Hz 附近有突然抬起的幅值先检查是否忘记减均值而不是立刻怀疑补零或 FFT 代码。6. 实测中踩过的坑和工程建议自编 FFT 的坑和直接用 MATLAB ffT 的坑不完全一样因为你能控制的地方越多越容易在细节上翻车。我把最近项目中踩过的几个典型问题按“症状—原因—处理”列出来希望能帮大家少走点弯路。6.1 幅值修正后仍然偏小大多是频谱泄漏导致症状输入幅值 2 的余弦修正后峰值只有 1.7 左右。原因信号持续时长不是 nfft 的整数周期补零后时域信号存在截断能量被泄漏到旁瓣。处理先看峰值附近所有谱线能量的总和或者加窗再修正。加窗后不能直接沿用 2/N 这个系数因为窗函数本身改变了信号能量需要按窗增益同步修正。窗函数修正的经验公式加窗后单边幅值约为 2 * abs(X_win) / (nfft * sum(win)/nfft)其中 sum(win)/nfft 是窗的相干增益。对普通矩形窗这个值约等于 1对汉宁窗约等于 0.5所以幅值修正时要把 2 替换成 2/0.54 这样的调整系数。别看它只差一个系数工程里差一点后面的阈值判断都会出错。6.2 直接调用自编函数nfft 不是 2 的幂时静默出错递归 radix-2 版本我在前面强调过要加幂判断。因为递归实现里如果 N 不是 2 的幂函数还能继续拆分成奇偶子序列到某些深度的时候子序列长度不一致最后结果根本不是正确的频谱但不会直接报错。这个很坑。我在迭代版本里同样需要检查输入长度。更保险的做法是在 myfft 入口强制要求nfft 2^nextpow2(length(x))除非用户明确指定更大的值。业务逻辑里避免随手写nfft length(x)那等于完全没补零如果长度恰好不是 2 的幂自编 radix-2 FFT 根本无法运行。6.3 相位图上的脉冲跳变也许不是信号相位跳变信号本身相位连续但相位图在幅度接近 0 处会出现大量毛刺。这时候不要急着下结论说信号相位不稳定。可以先画 abs(X)对幅度小于阈值比如最大幅度的 1e-6的索引做掩膜把对应相位显示成 NaN画出相位图就不会被噪声刷屏。我现在常用的绘图代码mag_for_phase abs(X(1:halfN)); threshold max(mag_for_phase) * 1e-6; phase_clean phase_one; phase_clean(mag_for_phase threshold) NaN; plot(f_one, phase_clean);这样画出来的相位图只在有意义的峰值附近有值其他位置留白既能看清相位趋势又不会被噪声干扰。6.4 补零不等于窗函数更不能替代抗混叠滤波器最后一个容易混淆的概念补零只是在数字域把尾部填 0它不影响模拟域的抗混叠判断。如果采样前模拟通道混入了高于奈奎斯特频率的噪声补零再多频谱图上照样会有混叠成分。所以我处理采集数据时的顺序是先用抗混叠滤波器处理再决定是否加窗最后才考虑补零到 2 次幂。顺序反了后面所有操作都是白费。6.5 验证自编 FFT 是否正确的标准流程我每次写完自编 FFT不会直接上真实数据而是先用一组能构造出理论结果的信号验证构造两个频率、已知幅度和初相的正弦叠加再用内置 fft 和我自编 myfft 做对比最大误差应该小于 1e-12 量级。如果浮点误差到了 1e-9先去查 nfft 是否幂、位反转是否出错误差在 1e-12 以内才能认为算法本身没问题。具体的验证代码fs 1000; N 250; n 0:N-1; x 2*cos(2*pi*60*n/fs 0.3) 1.5*cos(2*pi*120*n/fs - 0.7); nfft 512; X1 fft(x, nfft); [X2, ~] myfft(x, fs, nfft); err max(abs(X1 - X2)); fprintf(最大误差 %e\n, err);我在实际运行中误差大约在 1e-14 到 1e-13 之间属于浮点舍入的正常范围。到这个精度你就可以放心地画频谱图和相位图了。说到底自编 FFT 不算什么高大上的工作它真正考验的是对补零、位反转、旋转因子、幅值修正这些细节的理解。把这些细节搞清楚后再看频谱图上的一个峰你会知道它是真实频率、是泄漏旁瓣还是混叠伪影再看相位图的一跳你会知道它是物理跳变还是 angle 主值被截断。能到这个程度FFT 在你的项目里才真正算工具而不是黑盒。 自编FFT详解：补零与幅值修正的工程实践
1. 为什么我要自己实现 FFT补零和幅值修正这两个细节决定成败最近在做一个振动特征提取的预研采集器每帧只给我 250 个采样点但后端的硬件 FFT 模块要求输入长度必须是 2 的幂。这就免不了要把数据补零到 256 点再送进处理器。你可能会说直接调用 MATLAB 的 fft() 不就行了吗问题是自编的 FFT 在嵌入式侧验证时更清晰因为你能一步步看到补零发生在哪、每级蝶形的旋转因子怎么引入、最后幅值修正该除以哪个 N。折腾下来我发现真正影响频谱质量的并不是 FFT 算法本身而是补零和幅值修正这两个容易被忽略的环节。这个场景适合几类人看正在把 MATLAB 里的 FFT 算法改写成 C/MEX 的工程师需要理解补零为什么能平滑频谱的科研人员以及一直被频谱幅值对不上实际信号幅值困扰的初学者。我会把自编 radix-2 FFT 的完整思路、补零的取舍、频谱图和相位图的绘制判读一起串起来讲代码可以直接复制到 MATLAB 里跑。1.1 数据不足位补零到底补的是什么所谓“不足位补 0”在 MATLAB 里写成两行N_actual length(x); nfft 2^nextpow2(N_actual); x_padded [x(:); zeros(nfft - N_actual, 1)];第一行取原始长度第二行找一个不小于该长度的 2 的幂。这里要清楚补零不是给你凭空增加有效信息而是在频域做插值。原始 250 点的 DFT 结果有 250 条谱线补零到 256 点之后变成 256 条谱线。多出来的谱线并没有让原始信号的物理分辨率提高它们只是把原来谱线之间的空隙用 sinc 插值填上了。实际采集时长不变的情况下补零后的频谱更平滑通过搜索峰值也更容易读准频率。如果非要较真“补零能不能分辨两个频率很接近的信号”答案是否定的。物理频率分辨率由实际采样时长决定也就是 fs/N_actual。补零只能让谱线更细密但两个靠得很近的峰值该糊成一团还是糊成一团。1.2 幅值修正是自编程里最大的坑直接用 abs(fft(x)) 画频谱画出来的峰值往往跟你输入的振幅对不上。原因很简单MATLAB 自带的 fft() 不自动除以 N三角函数变换的结果默认分布在正负频率两个方向上。以幅度为 2 的余弦信号为例单侧频谱上峰值修正后应当是 2不做修正你会看到一个接近 1 的数值还要乘一个 2 才是真正的单边幅值。自编 FFT 也一样内部每个蝶形只是累加和旋转因子相乘没有任何归一化。所以最终一定要做两步mag_raw abs(X); mag_one_sided mag_raw / nfft; mag_one_sided(2:end-1) mag_one_sided(2:end-1) * 2;这么处理之后你才能在频谱图上直接读出信号的真实幅值。我在下面第 4 节会详细拆解公式这里先记住一个原则——幅值修正是“除以 FFT 点数、单边谱补乘 2、直流和奈奎斯特频率不乘 2”。2. 自编 FFT 的实现思路从 DFT 公式到位反转和蝶形运算FFT 本质上是把 DFT 的运算拆成一层层蝶形。DFT 的原始定义是X(k) sum_{n0}^{N-1} x(n) * exp(-j2pikn/N)直接实现要算 N^2 次复乘N 到 256 可能无所谓但工程里 N 经常是 1024 甚至 4096开销就大了。radix-2 FFT 利用旋转因子的周期性和对称性把计算量降到 N*log2(N) 的量级。2.1 我自己采用的递归 radix-2 实现自编 FFT 在 MATLAB 里有很多写法。为了可读性我选择了 DIT按时间抽取递归版本核心只有十行左右function Y my_fft_radix2(x) % 递归实现 radix-2 按时间抽取 FFT N length(x); if N 1 Y x; return; end x x(:); % 抽取偶数位和奇数位 x_even my_fft_radix2(x(1:2:end)); x_odd my_fft_radix2(x(2:2:end)); % 旋转因子 Wn exp(-1j * 2 * pi * (0:N/2-1) / N); % 蝶形运算 Y [x_even Wn .* x_odd; x_even - Wn .* x_odd]; end这段代码每次把序列一分为二偶数下标的序列和奇数下标的序列分别递归调用自己再用蝶形公式合并。旋转因子 Wn 就是单位圆上等间隔取的相位点。Wn .* x_odd是把奇数序列先做一次相旋再与偶数序列相加。2.2 外层补零封装函数实际工程里我不希望手写补零逻辑散落在各处所以又套了一层接口名字也叫 myfft和 MATLAB 内置 fft 保持相同调用习惯function [X, f] myfft(x, fs, nfft) % myfft: 自编 FFT不足 nfft 位自动补零 % 输入: % x: 时域信号 % fs: 采样率 % nfft: FFT 点数必须是 2 的幂若缺省则自动取 nextpow2 if nargin 3 nfft 2^nextpow2(length(x)); end x x(:); if nfft length(x) error(nfft 不能小于输入数据长度); end if mod(log2(nfft), 1) ~ 0 error(自编 radix-2 FFT 要求 nfft 为 2 的幂); end % 核心不足位补 0 x_pad [x; zeros(nfft - length(x), 1)]; % 调用 radix-2 递归 FFT X my_fft_radix2(x_pad); f (0:nfft-1) * fs / nfft; end这里有个测试细节如果 nfft 由2^nextpow2()自动生成肯定满足 2 的幂条件但封装函数里最好还是加一个幂判断防止外部传参时手误传了 100 或者 300 这样的数。递归实现会直接算错而且很难察觉。2.3 递归 vs 迭代到底选哪个递归版本代码短最适合讲解和验证原理但性能一般。实际嵌入式里通常用迭代 DIT。迭代 DIT 需要先把输入序列按二进制位反转重新排列然后一层层蝶形合并。同样的逻辑可以这么写function Y my_fft_iterative(x) N length(x); xd bitrevorder(x); % 也可以用自己写的位反转 Y xd; len 1; while len N half len; len len * 2; for k 0:half:N-1 for j 0:half-1 w exp(-1j * 2 * pi * j / len); u Y(k j 1); v Y(k j half 1) * w; Y(k j 1) u v; Y(k j half 1) u - v; end end end end位反转是迭代 DIT 里最麻烦的部分。MATLAB 有 bitrevorder 可以直接用但在嵌入式环境里要自己实现比如对下标做二进制逆序。这就是为什么我建议先用递归版本把算法流程吃透等移植到效率要求高的场景再换迭代版本。3. 补零的本质频谱插值而不是提高频率分辨率很多初学 FFT 的人都误以为补零到 1024 点频率分辨率就变成 fs/1024。这是不对的。频率分辨率的物理下限是实际有效信号长度决定的 fs/N_actual补零之后虽然谱线间距变小但谱线的包络还是原来的包络。我习惯把它类比成图像放大放大后的像素变多了你可以更精细地看轮廓但没法看清照片里原本就没有的细节。3.1 分辨率与谱线间距的对照还是用刚才的 250 点例子fs1000Hz实际分辨率 df_real 1000/250 4Hz补零到 256 点后谱线间距 df_padded 1000/256 3.90625Hz。这两者看起来很接近但含义完全不同。如果另一个信号的频率比你当前信号低 4Hz只有真实采集 250 点以上才能区分出来补零再多也只是让已有峰的形状更光顺不会把两团靠在一起的峰拆开。我可以再做一组对照实验fs1000Hz分别采集 0.25s 和 0.5s两个时长都补零到 1024 点。0.25s 的频谱虽然也画出了 1024 条线可一旦两个正弦间隔小于 4Hz它绝对区分不开而 0.5s 的频谱可以轻松分辨之前无法区分出的接近频率。这说明什么想要更好的频率分辨率请去延长采样时间而不是一味加大 nfft。3.2 nfft 到底应该设成多少实际项目中我的经验分两种如果仅是为了满足 radix-2 长度要求nfft 2^nextpow2(length(x)) 就够。如果是为了让频谱图线条更平滑、峰值搜索更精确可以在原始点数基础上乘 4 或 8但不要几百倍地补。因为补零并不是不花代价。补零太多会把峰值附近的旁瓣画得特别细人眼觉得“变好看了”但旁瓣泄漏也会被插值成连续起伏反而干扰峰值判读。而且嵌入式场景里 nfft 越大存储和计算量增长越明显。补零时还有一个不起眼的点如果原始信号的 DC 分量不为 0x_pad尾部突然从非零变成 0会在频域产生一个小冲击成分。所以必要时应先减掉均值再补零或者把均值单独换算到 0Hz 谱线中。4. 幅值修正与相位还原不能只看 abs(FFT) 就结束自编 FFT 也好MATLAB 内置 fft 也好默认输出都是复数。很多人画频谱图时只取了幅值把 phase 丢在一边结果到了需要判断时间延迟、滤波器群延时或相位噪声时又重新采数据。所以我在项目里固定要求频谱图、相位图一起画提前习惯就好。4.1 单边谱的幅值修正公式标准的修正流程四步nfft length(X); % 实际用的 FFT 点数包含补零后的点数 halfN floor(nfft/2) 1; % 单边谱有效索引范围 f_one f(1:halfN); mag_one abs(X(1:halfN)) / nfft; % 除了 DC(索引1) 和 Nyquist(索引 halfN) 外其余频率点要乘 2 mag_one(2:halfN-1) mag_one(2:halfN-1) * 2; phase_one angle(X(1:halfN));为什么是单边谱因为实信号经过傅里叶变换后正频率和负频率是对称的。我们把幅度按正负各半重新合并到单边所以除非直流或奈奎斯特频率都需要乘以 2。幅值修正的细节用一张表看更直观谱线位置修正方式原因DC0 Hzabs(X(1)) / nfft直流分量只出现在 0Hz不参与正负频率折叠中间频率2 * abs(X(k)) / nfft正负频各占一半需要合并Nyquistfs/2abs(X(halfN)) / nfft奈奎斯特频率同样只有一个4.2 相位图和 angle() 的坑相位直接用 angle(X) 得到的是主值范围在 [-pi, pi]。频谱图上会看到很多 2pi 跳变那不是相位突变而是角度被截断在了这个区间。要想看到线性相位趋势要用 unwrap() 展开phase_unwrap unwrap(angle(X(1:halfN))); plot(f_one, phase_unwrap);还有一个常见问题某些 bin 上信号能量接近 0 时幅值很小相位会变成一堆杂乱的角度看起来毫无规律。这是数值噪声导致的。判读相位时只关心幅值峰值附近的谱线别把整个相位图都当成有效信息。工程上我通常会先找谱峰再在当前峰值周围 ±2 根谱线范围内看相位一致性。5. 频谱图和相位图的绘制与判读把这个过程的完整脚本放在一起。为了演示构造一个 60Hz、幅值 2 的余弦信号采样率 1000Hz采集 250 点然后用自编 myfft 补零到 256 点fs 1000; N_actual 250; f0 60; A0 2; n 0:N_actual-1; x A0 * cos(2*pi*f0*n/fs); nfft 2^nextpow2(N_actual); [X, f] myfft(x, fs, nfft); halfN nfft/2 1; f_one f(1:halfN); mag_one abs(X(1:halfN)) / nfft; mag_one(2:halfN-1) mag_one(2:halfN-1) * 2; phase_one unwrap(angle(X(1:halfN))); figure; subplot(2,1,1); plot(f_one, mag_one, LineWidth, 1.2); xlabel(频率 (Hz)); ylabel(幅值); title(频谱图幅值已修正); grid on; subplot(2,1,2); plot(f_one, phase_one, LineWidth, 1.2); xlabel(频率 (Hz)); ylabel(相位 (rad)); title(相位图); grid on;运行之后会看到两个子图上图在 60Hz 附近出现一个明显峰值修正后峰值接近 2下图相位在峰值附近保持一条接近线性的曲线远离峰值的地方相位跳变得比较乱。这里需要说明一点由于 250 点采集 60Hz 并不是整个 256 点的整数周期所以在 60Hz 两侧会出现若干旁瓣这是频谱泄漏不是自编 FFT 算错了。如果想知道精确幅值在加窗前通常取峰值附近的包络或改用整周期采样把 f0 设计成 bin 的整数倍。下面把从图里读信息的几个小经验分享给你峰值的横坐标往往不是正好落在 f0。如果信号存在泄漏峰值位置会偏离真实频率这时用二次插值法parabolic interpolation在峰值附近拟合比直接读谱线更准。相位图只需要看峰位附近的局部。设计滤波器时用频率轴上所有有效带的相位响应查群延时。如果只在 0Hz 附近有突然抬起的幅值先检查是否忘记减均值而不是立刻怀疑补零或 FFT 代码。6. 实测中踩过的坑和工程建议自编 FFT 的坑和直接用 MATLAB ffT 的坑不完全一样因为你能控制的地方越多越容易在细节上翻车。我把最近项目中踩过的几个典型问题按“症状—原因—处理”列出来希望能帮大家少走点弯路。6.1 幅值修正后仍然偏小大多是频谱泄漏导致症状输入幅值 2 的余弦修正后峰值只有 1.7 左右。原因信号持续时长不是 nfft 的整数周期补零后时域信号存在截断能量被泄漏到旁瓣。处理先看峰值附近所有谱线能量的总和或者加窗再修正。加窗后不能直接沿用 2/N 这个系数因为窗函数本身改变了信号能量需要按窗增益同步修正。窗函数修正的经验公式加窗后单边幅值约为 2 * abs(X_win) / (nfft * sum(win)/nfft)其中 sum(win)/nfft 是窗的相干增益。对普通矩形窗这个值约等于 1对汉宁窗约等于 0.5所以幅值修正时要把 2 替换成 2/0.54 这样的调整系数。别看它只差一个系数工程里差一点后面的阈值判断都会出错。6.2 直接调用自编函数nfft 不是 2 的幂时静默出错递归 radix-2 版本我在前面强调过要加幂判断。因为递归实现里如果 N 不是 2 的幂函数还能继续拆分成奇偶子序列到某些深度的时候子序列长度不一致最后结果根本不是正确的频谱但不会直接报错。这个很坑。我在迭代版本里同样需要检查输入长度。更保险的做法是在 myfft 入口强制要求nfft 2^nextpow2(length(x))除非用户明确指定更大的值。业务逻辑里避免随手写nfft length(x)那等于完全没补零如果长度恰好不是 2 的幂自编 radix-2 FFT 根本无法运行。6.3 相位图上的脉冲跳变也许不是信号相位跳变信号本身相位连续但相位图在幅度接近 0 处会出现大量毛刺。这时候不要急着下结论说信号相位不稳定。可以先画 abs(X)对幅度小于阈值比如最大幅度的 1e-6的索引做掩膜把对应相位显示成 NaN画出相位图就不会被噪声刷屏。我现在常用的绘图代码mag_for_phase abs(X(1:halfN)); threshold max(mag_for_phase) * 1e-6; phase_clean phase_one; phase_clean(mag_for_phase threshold) NaN; plot(f_one, phase_clean);这样画出来的相位图只在有意义的峰值附近有值其他位置留白既能看清相位趋势又不会被噪声干扰。6.4 补零不等于窗函数更不能替代抗混叠滤波器最后一个容易混淆的概念补零只是在数字域把尾部填 0它不影响模拟域的抗混叠判断。如果采样前模拟通道混入了高于奈奎斯特频率的噪声补零再多频谱图上照样会有混叠成分。所以我处理采集数据时的顺序是先用抗混叠滤波器处理再决定是否加窗最后才考虑补零到 2 次幂。顺序反了后面所有操作都是白费。6.5 验证自编 FFT 是否正确的标准流程我每次写完自编 FFT不会直接上真实数据而是先用一组能构造出理论结果的信号验证构造两个频率、已知幅度和初相的正弦叠加再用内置 fft 和我自编 myfft 做对比最大误差应该小于 1e-12 量级。如果浮点误差到了 1e-9先去查 nfft 是否幂、位反转是否出错误差在 1e-12 以内才能认为算法本身没问题。具体的验证代码fs 1000; N 250; n 0:N-1; x 2*cos(2*pi*60*n/fs 0.3) 1.5*cos(2*pi*120*n/fs - 0.7); nfft 512; X1 fft(x, nfft); [X2, ~] myfft(x, fs, nfft); err max(abs(X1 - X2)); fprintf(最大误差 %e\n, err);我在实际运行中误差大约在 1e-14 到 1e-13 之间属于浮点舍入的正常范围。到这个精度你就可以放心地画频谱图和相位图了。说到底自编 FFT 不算什么高大上的工作它真正考验的是对补零、位反转、旋转因子、幅值修正这些细节的理解。把这些细节搞清楚后再看频谱图上的一个峰你会知道它是真实频率、是泄漏旁瓣还是混叠伪影再看相位图的一跳你会知道它是物理跳变还是 angle 主值被截断。能到这个程度FFT 在你的项目里才真正算工具而不是黑盒。