企业官网建设流程全解析

MATLAB FFT 单边幅度谱：归一化、奇偶长度与频谱泄漏的完整实验

对一个幅度为 2.4 的余弦做 FFT，直接画abs(fft(x))，峰值通常不会显示成 2.4。除以样本数以后，正频率一侧又只剩一半。把整张谱乘以 2，看似修好了余弦，直流和 Nyquist 点却跟着错了。

这几个问题可以放进同一份小实验里检查：先用已知幅度的实信号验证单边谱，再把记录内的周期数从 20 改成 20.5，比较矩形窗与手写 Hann 窗。所有输入都在代码中生成，不需要音频文件、开发板或外部数据。

1. 先运行，再看结果

将第 3 节的完整代码保存为example01_fft.m，文件名必须与第一行函数名一致。在 MATLAB 中切换到该文件所在目录，执行：

example01_fft

也可以在已经配置好 MATLAB 命令行路径的终端中执行：

matlab -batch "example01_fft"

代码只使用基础 MATLAB，Hann 窗由余弦公式生成，没有调用hann、periodogram等工具箱函数。本次实际运行版本为 MATLAB25.2.0.2998904 (R2025b)；其他版本未在本次测试中验证。

默认生成output01/fft_amplitude_leakage.png和output01/fft_metrics.csv。图在后台导出，不会弹出窗口。同一输出目录再次运行会覆盖这些示例文件；需要保留多次结果时，可以改用example01_fft('output01_run2',false)。

2026 年 10 月 14 日的正向 batch 退出码为 0，11 项断言通过。主要结果如下。

这里最后一行很重要：相干增益校准后，离开 FFT 栅点的余弦峰值仍没有恢复到 2。

2. 单边谱究竟哪些点乘以 2

MATLAB 的fft计算未按输入长度归一化的 DFT。对长度为 N、恰好落在 DFT 栅点且不位于 DC/Nyquist 的实余弦，正、负频率各分到一半幅度。因此先计算abs(fft(x))/N，再合并共轭频率对，才能得到常见的单边幅度表示。函数定义和官方单边谱示例见 fft 文档。

• N 为偶数：保留 k=0 到 N/2，只有中间的正频率点乘以 2
• N 为奇数：保留 k=0 到 floor(N/2)，除 DC 外都乘以 2

差别出在末尾。偶数长度的 k=N/2 正好是 Nyquist 点，正负频率在离散域中重合；奇数长度没有恰好位于 Fs/2 的 DFT 栅点，最后一个正频率点仍有对应的负频率点。

代码中的 MATLAB 下标从 1 开始，所以数学上的 k=9 对应p(10)。图的横轴采用 DFT bin，不是 Hz；已知采样率 Fs 时，频率换算为f = k*Fs/N。两个不同 N 的案例是在各自的 k=9 上检查幅度，并没有声称它们在相同采样率下具有相同物理频率。

3. 完整可运行代码

下面保留了实际运行的整份函数及两个局部函数。negativeMode默认为 false；先按正常模式运行即可。

functionexample01_fft(outDir,negativeMode)% Real-valued signals; base MATLAB only. No external input.ifnargin<1,outDir=fullfile(fileparts(mfilename('fullpath')),'output01');endifnargin<2,negativeMode=false;endif~isfolder(outDir),mkdir(outDir);endfprintf('MATLAB_VERSION=%s\nNEGATIVE_MODE=%d\n',version,negativeMode);tol=1e-10;fig=figure('Visible','off','Color','w','Position',[1001001100760]);cleanup=onCleanup(@()close(fig));tiledlayout(2,2,'TileSpacing','compact');forN=[128127]n=0:N-1;x=1.25+2.4*cos(2*pi*9*n/N);ifmod(N,2)==0,x=x+0.6*cos(pi*n);else,x=x+0.7*cos(2*pi*floor(N/2)*n/N);endp=oneSided(x);ifnegativeMode% Wrong recipe: double all bins, including DC and Nyquist.wrong=2*abs(fft(x))/N;p=wrong(1:floor(N/2)+1);endck(abs(p(1)-1.25)<tol,'DC','DC must not be doubled');ck(abs(p(10)-2.4)<tol,'Tone','Coherent cosine amplitude');ifmod(N,2)==0ck(abs(p(end)-0.6)<tol,'Nyquist','Even Nyquist must not be doubled');elseck(abs(p(end)-0.7)<tol,'OddLast','Odd final bin must be doubled');endfprintf('N=%d DC=%.12g TONE=%.12g LAST=%.12g\n',N,p(1),p(10),p(end));nexttile;stem(0:floor(N/2),p,'.');grid on;title(sprintf('N=%d: correct one-sided amplitudes',N));xlabel('DFT bin k');ylabel('Amplitude');endN=256;n=0:N-1;k=0:N-1;xc=2*cos(2*pi*20*n/N);xl=2*cos(2*pi*20.5*n/N);w=0.5-0.5*cos(2*pi*n/N);% PERIODIC Hann, implemented without hann().cg=mean(w);pc=oneSided(xc);pr=oneSided(xl);phc=oneSided(xc.*w)/cg;ph=oneSided(xl.*w)/cg;Yc=fft(xc);Yr=fft(xl);Yh=fft(xl.*w);offCoherent=~(k==20|k==N-20);coherentLeak=sum(abs(Yc(offCoherent)).^2)/sum(abs(Yc).^2);% Same guard bands for both windows; mirrored positive/negative main lobes.outside=~(abs(k-20.5)<=2.5|abs(k-(N-20.5))<=2.5);rectLeak=sum(abs(Yr(outside)).^2)/sum(abs(Yr).^2);hannLeak=sum(abs(Yh(outside)).^2)/sum(abs(Yh).^2);ck(coherentLeak<1e-20,'CoherentLeak','Integer-period leakage near roundoff');ck(rectLeak>0.01,'FractionalLeak','Noninteger-period leakage exists');ck(hannLeak<rectLeak/10,'WindowLeak','Hann reduces far-out leakage in this case');ck(abs(cg-0.5)<tol&&abs(phc(21)-2)<tol,'CoherentGain','Hann coherent gain calibration');ck(max(ph)<1.9,'Scalloping','Coherent gain does not fix half-bin scalloping');fprintf('COHERENT_LEAK=%.12g RECT_LEAK=%.12g HANN_LEAK=%.12g\n',coherentLeak,rectLeak,hannLeak);fprintf('HANN_CG=%.12g COHERENT_PEAK=%.12g HALF_BIN_PEAK=%.12g\n',cg,phc(21),max(ph));fprintf('LIMIT: isolated coherent tone away from DC/Nyquist; not PSD or broadband calibration.\n');nexttile;plot(0:N/2,pc,'-o',0:N/2,pr,'-s',0:N/2,ph,'-d','MarkerSize',3);xlim([1427]);grid on;xlabel('DFT bin k');ylabel('Amplitude');legend('Integer, rectangular','Half-bin, rectangular','Half-bin, Hann / CG','Location','best');title('Hann reduces leakage but retains scalloping');nexttile;semilogy(0:N/2,max(pr,1e-14),0:N/2,max(ph,1e-14),'LineWidth',1.2);xlim([070]);ylim([1e-73]);grid on;xlabel('DFT bin k');ylabel('Amplitude (log scale)');legend('Rectangular','Hann / CG');title('Far-out leakage; main lobe is broader');exportgraphics(fig,fullfile(outDir,'fft_amplitude_leakage.png'),'Resolution',150);metrics=table(coherentLeak,rectLeak,hannLeak,cg,max(ph),'VariableNames',...{'coherent_leak','rect_leak','hann_leak','hann_CG','half_bin_peak'});writetable(metrics,fullfile(outDir,'fft_metrics.csv'));fprintf('RESULT=PASS example01_fft\n');endfunctionp=oneSided(x)N=numel(x);q=abs(fft(x))/N;p=q(1:floor(N/2)+1);ifmod(N,2)==0,p(2:end-1)=2*p(2:end-1);else,p(2:end)=2*p(2:end);endendfunctionck(ok,id,msg)assert(ok,['lessonFFT:'id],msg);fprintf('ASSERT_PASS=%s\n',id);end

4. 对照实际导出的图

上排两图检查奇偶长度。两幅图的直流都为 1.25，k=9 的余弦都为 2.4。右端分别为 0.6 和 0.7，但两者采用不同的倍增规则。这正是只测偶数长度时容易漏掉的问题。

左下图改用 N=256。记录中包含 20 个周期时，幅度为 2 的余弦落在 k=20；改成 20.5 个周期后，频率落在两个栅点之间，多个频点都有可见幅度。有限记录相当于截断信号，窗函数的频谱会影响观测到的谱形。Hann 能压低远处旁瓣，同时也会扩宽主瓣，这是窗选择的取舍。MathWorks 的频谱估计说明介绍了这个机制。

右下图用对数纵轴展开旁瓣。它画的仍是幅度，既没有平方，也没有换算成 dB，更不是 PSD。代码中的max(...,1e-14)只用于避免对数坐标绘制零值，不参与断言和能量占比计算。

4.1 为什么校准要除以 0.5

本例使用周期型 Hann：

w[n] = 0.5 - 0.5*cos(2*pi*n/N)，其中n=0,...,N-1。

它的相干增益为mean(w)=0.5。加窗后的oneSided已经除以 N，再除以相干增益，相当于用窗系数之和归一化。对本例远离 DC/Nyquist 的孤立整频点余弦，校准后检查值为 2。

但当真实频率位于 k=20.5 时，FFT 只在整数 k 上采样谱线。最高采样点未必正好落在主瓣顶端，所以实测最大值约为 1.697662。这种栅栏损失不会被一个固定的相干增益系数自动消除。DC 附近、Nyquist 附近、多频率相互重叠以及宽带输入，都需要重新检查校准适用条件。

4.2 泄漏比例怎样计算

“整周期泄漏比例”和“矩形窗/Hann 泄漏比例”的区域定义不同，不能混在一起比较。

整周期案例只排除 k=20 与 k=256-20 两个目标频点，其余平方模之和除以整张双边谱的平方模之和。实测约 1.13e-29，处于浮点舍入量级。

半频点案例则围绕 k=20.5 与 k=256-20.5，各保留半宽 2.5 bin 的保护区；矩形窗和 Hann 都使用同一保护区，然后比较区域外能量占各自总频域能量的比例。得到约 6.7124% 与 0.01138%。这是当前信号、窗和保护区下的结果，不能据此给所有窗排一个通用名次，也不能把它理解成加窗前后总能量保持不变。

5. 让错误公式真的失败一次

执行：

example01_fft('output01_negative',true)

该模式故意把所有频点乘以 2，包括 DC。输入的直流为 1.25，错误算法给出 2.5，第一项断言立即失败。

本次负向 batch 的退出码为 1，捕获的错误标识是lessonFFT:DC。这个失败说明断言能抓住指定错误；如果实际遇到的是路径错误、函数不存在或图形环境异常，不能算作这个负对照通过。

6. 拿去分析真实采样前还要检查什么

• 这里处理实信号。复数 IQ 通常不满足相同的共轭对称条件，不能直接丢掉负频率再倍增
• DC 显示的是偏置的绝对值。Nyquist 测试使用0.6*cos(pi*n)；采样本身会丢失某些相位信息，不能推广成任意相位的连续 Nyquist 余弦都可恢复原幅度
• 输入必须有明确的采样时间关系。乱序、丢样和不均匀采样不能只靠一行 FFT 自动修正
• 补零能加密频谱采样，不能增加已采集记录的信息。本文源码没有运行补零实验；若扩展为 NFFT 点 FFT，单边端点规则应按 NFFT 的奇偶判断，幅度归一化仍要区分原始记录长度与窗和。可参阅 Amplitude Estimation and Zero Padding
• 这些数值是余弦峰幅，不是 RMS、功率或功率谱密度。需要 PSD 时应采用相应的定义、单位和归一化

7. 参考来源与延伸阅读

选题参考了 MathWorks 教学项目 Fourier Analysis 的固定版本 README。该版本的 LICENSE.md为 BSD 3-Clause；完整课程列有 Symbolic Math Toolbox 依赖。本文实验独立编写，没有复制仓库代码或图片，也不需要运行整套课程。

如果数据来自 FPGA，可以接着看 ILA 在线调试：Vivado 抓取片内信号与定位下板问题。采样时钟、有效数据条件和数值解释先确认清楚，后续幅度谱才有可靠的输入。

最后留一个具体问题：把整周期案例的 N 改成奇数，却仍保留p(2:end-1)的倍增范围，哪个频点会低一半？可以先根据端点规则判断，再观察OddLast断言。