简介一份基于MATLAB的快速傅里叶变换FFT频谱分析仿真例程面向信号处理初学者和需要快速上手频域分析的技术人员用于解决时域信号到频谱的可视化转换问题。该例程将完整的频谱分析流程封装为简洁脚本通过读取信号数据、调用fft函数执行变换再结合abs与angle提取幅度和相位信息最后用plot绘制频谱图直观展示信号在频域中的分布帮助识别特征频率、噪声与谐波成分。压缩包内共2个文件包含一个MATLAB主脚本及其自动备份文件.asv包体仅951B轻量易读方便在此基础上修改扩展。已有310人学习既适合课程实验、毕业设计或自学入门也可作为教师授课演示素材。借助该例程读者能深入理解FFT作为离散傅里叶变换高效算法的原理掌握频谱分析的标准流程为通信、音频处理、振动分析等领域进一步应用打下基础。1. 从 fft.rar 说起MATLAB 频谱仿真到底在仿真什么解压开一个名为 fft.rar 的压缩包里面往往是一堆 .m 文件和一张看不出参数的频谱截图。真正有价值的不是那张图而是把「时域信号 → fft → 频谱」这条链在 MATLAB 里跑通的方法。所谓频谱仿真就是先在时域构造一段成分已知的信号幅值、频率、相位都自己定再做 FFT 得到频谱用谱峰去反推和验证采样率、采样点数、窗函数这些参数是否合理。它的核心价值在于在接触真实采集数据之前先用一段完全已知的信号把算法边界摸清楚。无论是 STM32F4 上采集音频做实时频谱分析还是在 Vivado 里调用 FFT IP 核做频域处理上位机用 MATLAB 先仿真一遍都是最省事的对拍手段。所以下面按「原理 → 脚本 → 排坑 → 实测数据」的顺序把参数怎么定、代码怎么写、结果怎么验一次讲透。2. 写 MATLAB FFT 频谱仿真前必须先定的三个参数fs、N 与频率分辨率2.1 采样率 fs 与奈奎斯特频率先定采样条件再谈频谱FFT 不是把任意一串数变成频谱的魔法。对一段以 fs 采样得到的实信号FFT 输出的频率范围只有 0 到 fs/2 是可信的fs/2 就是奈奎斯特频率。采样定理要求 fs 必须大于信号最高频率的两倍否则高于 fs/2 的成分会折叠回低频段造成混叠。仿真时容易忽略的一点是信号是你自己生成的最高频率已知所以 fs 应取最高频率的 5 到 10 倍给后续可能的抗混叠处理留出过渡带余量而不是卡着 2 倍去选。代码里时间轴是这样建立的fs 1000; % 采样率 1000 Hz N 2048; % 采样点数 t (0:N-1) / fs; % 时间轴0 到 (N-1)/fs 秒 x sin(2*pi*50*t); % 50 Hz 正弦信号t 的间隔固定为 1/fs 秒这个间隔决定了后面所有频率计算。若把 t 写成 linspace(0, 1, N)最后一个采样点落在 1 秒整比正确时间轴多出近一个采样周期FFT 后的频率会整体偏移。这是仿真脚本里最常见的隐性错误之一看起来波形正常频谱峰值却对不上设定频率。2.2 频率分辨率 df fs/N两个相邻峰值能不能分开由它决定FFT 输出相邻两条谱线之间的频率间隔是 df fs/N也就是常说的频率分辨率。两个频率相差小于 df 时它们在频谱上会合成一个宽峰无法区分。所以仿真里能不能看清两个峰不取决于 fft 函数本身而取决于 fs 和 N 的选择。很多人以为换一个更高级的算法就能提升分辨率实际在同样的 N 和 fs 下任何离散傅里叶变换实现给出的谱线间隔都是一样的。采样率 fs (Hz)采样点数 N频率分辨率 df (Hz)可分辨的最小频差100010001.0约 1 Hz 以上100020480.488约 0.5 Hz800081920.977约 1 Hz8000655360.122约 0.12 Hz需要区分两个相差 Δf 的频率时先算 N fs / Δf。例如要分辨 50 Hz 与 50.5 HzΔf 0.5 HzN 至少要 2000取 2048 勉强可用取 4096 更稳。反过来如果只需要 1 Hz 级别的谱线N 没必要堆得太大否则计算量和内存白花幅值谱也不会更准——多出来的时间只会让信号包含更多不确定性。2.3 fft 最小骨架时间轴、频谱轴与幅值谱的一次完整对应一段能直接运行的仿真骨架长这样fs 1000; % 采样率 N 2048; % fft 点数 t (0:N-1) / fs; x sin(2*pi*50*t); % 幅值 1 的 50 Hz 信号 Y fft(x); % 复频谱长度 N f (0:N-1) * fs / N; % 频率轴第 k 条谱线对应 k*fs/N P abs(Y) / N; % 双侧幅值谱 plot(f, P); xlim([0 fs/2]); % 只看可信的 0 ~ fs/2 段 xlabel(频率 (Hz)); ylabel(幅值);这里 f (0:N-1) * fs / N 把谱线序号 k 映射到实际频率第 0 条是直流第 N/2 条正好是 fs/2。abs(Y)/N 除以 N 是把 FFT 的累加结果归一化成每个频率成分贡献了多少幅值但注意它给出的是双侧谱50 Hz 的真实幅值 1 被拆成两半各 0.5分别落在正负 50 Hz 上。运行后先确认峰值位置对不对再谈幅值。2.4 单边频谱的幅值修正2/N 从哪里来什么时候不乘对于实信号负频率半段只是正频率半段的镜像没有额外信息画图时通常只保留 0 到 fs/2并把正频率处的幅值乘 2 还原真实幅值这就是单边幅值谱P1 P(1:N/21); % 取单边含直流与奈奎斯特点 P1(2:end-1) 2 * P1(2:end-1); % 中间所有点乘 2 f1 (0:N/2) * fs / N; % 单边频率轴两个例外必须记住直流分量第一个点本来就没有镜像不乘 2当 N 为偶数时频率轴末尾的 fs/2 点奈奎斯特频率点也只出现一次同样不乘 2。代码里 P1(2:end-1) 恰好把两端排除在外。常见错误是整条 P1 都乘 2结果直流和奈奎斯特点幅值偏高整整一倍在含直流偏置的仿真里一眼就能看出来。提示用 fftshift 观察对称频谱时频率轴要同步平移为 (-N/2 : N/2-1) * fs / N此时 0 频在中间横轴不再是 0 到 fs两种画法不要混用。3. 用 MATLAB 构造多频叠加信号的可复现频谱仿真脚本3.1 时域信号合成把幅值、频率、初相写进一条表达式仿真第一步是让信号成分已知。把多个正弦叠加是最常用的构造方式fs 8000; % 采样率 8 kHz N 8192; % 点数df 0.9766 Hz t (0:N-1) / fs; f1 100; A1 1.0; % 主峰 f2 250; A2 0.5; % 半幅值 f3 1000; A3 0.25; % 四分之一幅值带初相 x A1 * sin(2*pi*f1*t) ... A2 * sin(2*pi*f2*t) ... A3 * sin(2*pi*f3*t pi/4);表达式中 A 是要在频域验证的幅值f 是谱峰位置初相 pi/4 供验证相位谱时使用——abs 画图看不到相位用 angle 才能看见。选择 100、250、1000 Hz 有三个用意都在 0 到 fs/2 可信区间内三者间隔远大于 df 0.977 Hz不会互相干扰频率跨度足够大能检验频率轴是否为线性映射。如果想顺便观察直流分量信号尾部再加一个 0.5 常数偏置FFT 的第 0 条谱线就会显示 0.5。3.2 绘制单边频谱图从 fft 结果到横轴 Hz 纵轴幅值的完整映射把 3.1 节信号做成可直接复制运行的完整脚本% spectrum_demo.m fs 8000; N 8192; t (0:N-1) / fs; x 1.0*sin(2*pi*100*t) 0.5*sin(2*pi*250*t) 0.25*sin(2*pi*1000*t pi/4); Y fft(x); % 变换 P2 abs(Y) / N; % 双侧幅值谱 P1 P2(1:N/21); % 截取单边 P1(2:end-1) 2 * P1(2:end-1); % 幅值修正 f (0:N/2) * fs / N; plot(f, P1, LineWidth, 1); xlabel(频率 (Hz)); ylabel(幅值); title(单边频谱fs8000 Hz, N8192); grid on;注意三处顺序不能乱先 abs(Y)/N 后乘 2先截取 P2(1:N/21) 再对中间段修正频率轴用 0:N/2 而不是 0:N。plot 线性坐标足以看清三个主峰当需要观察低于 0.05 的小幅值成分时把纵轴换成 semilogy噪底和毛刺会显示得更清楚这是实际振动频谱分析的常用画法。3.3 频谱仿真脚本的必调参数速查表参数含义典型取值与调整方向fs采样率取最高信号频率的 5~10 倍加大可推远混叠边界NFFT 点数由 df fs/N 反推辨频需求 Δf 时 N ≥ fs/ΔfNFFT实际 FFT 长度可大于 N尾部补零建议取 2 的幂窗函数抑制频谱泄漏稳态周期信号用 hann强旁瓣抑制用 blackman幅值系数单边修正2/N直流与奈奎斯特点除外去直流预处理只关心交流成分时先减 mean(x)这张表是排错入口峰值位置不对先查 fs 和频率轴映射幅值不对先查修正系数和窗函数归一化峰变宽再查泄漏和补零方式。仿真时每改一个参数先在脑子里预期峰值怎么变再运行——预期与结果不一致的瞬间就是理解最深的时刻。3.4 仿真结果怎么自检峰值频率、幅值与理论值对账运行上一节脚本图上应出现三条谱线100 Hz 处幅值 1.0250 Hz 处 0.51000 Hz 处 0.25。用数据而不是肉眼对账[val, idx] findpeaks(P1, f, MinPeakHeight, 0.05); table(idx, val, VariableNames, {Freq_Hz, Amp})findpeaks 返回峰频率和峰幅值。如果 val 恰好等于理论幅值说明 fs、N、频率轴、修正全部正确如果 val 是理论值的 2 倍说明某段谱线重复乘了 2是 0.5 倍说明漏了修正频率对不上则先怀疑 t 是否用 linspace 构造再查频率轴公式。对账通过后把这段频谱计算封装成函数后续换窗、换采样率、换实测数据都复用同一套逻辑。4. FFT 频谱仿真绕不开的四个坑频谱泄漏、栅栏效应、噪声与混叠4.1 频谱泄漏当峰值频率不是 fs/N 的整数倍时幅值会塌把 3.1 节里的 100 Hz 改成 100.3 Hz 再运行峰值幅值会从 1.0 掉到 0.9 左右峰脚下还出现一串衰减旁瓣。原因是信号周期与 N/fs 的观测窗长度不匹配能量从主瓣漏到相邻谱线——这就是频谱泄漏。它让幅值读数变差也让小幅值信号被旁瓣淹没。常见做法是加窗把时域信号逐点乘一个两端衰减到 0 的窗函数强制观测区间两端连续N 8192; fs 8000; t (0:N-1)/fs; x sin(2*pi*100.3*t); w hann(N, periodic); % 周期汉宁窗 xw x(:) .* w(:); % 逐点相乘 Yw fft(xw); Pw abs(Yw(1:N/21)) / sum(w); % 归一化改为 sum(w) Pw(2:end-1) 2 * Pw(2:end-1); fw (0:N/2) * fs / N; plot(fw, Pw);归一化从 /N 改成 /sum(w) 是关键窗函数改变了信号的累积能量用窗长总和修正后幅值才能还原。加汉宁窗后 100.3 Hz 处的峰幅值从 0.9 回升到接近 1.0旁瓣明显压低。代价是主瓣变宽两个间隔很近的频率更难分开加窗始终是分辨率和幅值精度之间的折中。4.2 栅栏效应与补零NFFT 补零只插值、不提高分辨率FFT 只能在 df 的整数倍位置输出谱线峰值恰好落在两条谱线之间时读数偏低这就是栅栏效应。不少人误以为把 NFFT 加大就能解决NFFT 8 * N; % 补零到 8 倍长度 Ypad fft(x, NFFT); % fft 自动在尾部补零 fpad (0:NFFT/2) * fs / NFFT; Ppad abs(Ypad(1:NFFT/21)) / N; % 归一化仍用原始 N Ppad(2:end-1) 2 * Ppad(2:end-1);补零后谱线更密、峰形更圆但两个真实频率仍然分不开因为区分频率靠的是实际观测时长 N/fs而不是 FFT 长度。df fs/N 里的 N 是实际采样点数补零只是在相邻谱线之间插值不产生新信息。判断标准很简单只关心峰位置的读数精度补零到 4 或 8 倍足够关心两个接近频率能否分开回到 2.2 节把原始采样点数加长这才是唯一的正路。4.3 含噪信号的频谱仿真SNR 与窗函数的配合真实采集的信号都带噪仿真里用 randn 加高斯白噪声最直接x 1.0 * sin(2*pi*100*t) 0.1 * randn(size(t)); % SNR 约 17 dB加噪后噪底整体抬高100 Hz 峰仍清晰可见但把幅值从 1.0 降到 0.05峰就会被噪底吞掉。此时换任何窗都救不回来——窗函数改变的是泄漏旁瓣不是白噪声的统计功率。有效手段只有三类加长 N 摊薄噪底均值、多段频谱取平均降低方差、时域先滤波再 FFT。仿真阶段建议把 SNR 当参数扫一遍从 40 dB 一路降到 5 dB记录峰值读数的抖动范围提前知道算法在什么信噪比下失效比拿到真实数据后手忙脚乱要好得多。4.4 采样率不足的混叠仿真600 Hz 在 1000 Hz 采样下出现在 400 Hz这个坑最隐蔽。仿真里信号频率可以随便写但若 fs 不满足奈奎斯特定理FFT 不会报错而是把频率折叠到别处fs 1000; N 1024; t (0:N-1)/fs; x sin(2*pi*600*t); % 600 Hz fs/2 500 Hz Y fft(x); f (0:N/2) * fs / N; plot(f, 2*abs(Y(1:N/21))/N);运行后峰值出现在 400 Hz 而不是 600 Hz。混叠频率按 fa |f - k*fs| 折叠600 - 1000 -400取绝对值 400。这类问题在两类场景最常翻车采集卡前置抗混叠滤波器没配好仿真时随意加大信号频率却不同步提高 fs。判断方法只有一个——记录最高频率核对 fs/2 是否大于它必要时在脚本开头加一行 assertassert(fs 2*maxf, 采样率不足发生混叠);提示仿真阶段遇到混叠是好事。真实系统里混叠一旦发生就无法在数字域还原只能重采或加硬件滤波器而仿真只需改 fs 就能完整看到折叠过程值得专门写一个混叠 demo 存进自己的频谱工具箱。5. 从仿真到实测把 CSV 数据导入 MATLAB 做 FFT 频谱分析5.1 不知道采样率时怎么求频率频谱从时间列反推 fs仿真脚本验证完后最常见需求是把实测 CSV 数据拿来做同样分析比如供水管网噪声记录仪的频带划分、振动传感器的时程数据。readmatrix 可以直接读入raw readmatrix(vibration.csv); % 自动识别数值忽略表头 t_col raw(:,1); % 第一列时间戳 x_col raw(:,2); % 第二列信号幅值 fs 1 / mean(diff(t_col)); % 时间戳间隔的倒数即采样率 x x_col - mean(x_col); % 去直流避免 0 Hz 大峰盖住细节这里用 mean(diff(t_col)) 取平均时间间隔比直接用 t_col(2)-t_col(1) 更抗抖动实测采集的时间戳常有轻微不均。如果 CSV 只有信号列没有时间列常见做法是先输出相对频率单位 cycles/sample横轴按 k/N 表示再乘估计的 fs 得到物理频率。这个估计没有捷径采集时长 T 已知时 fs ≈ N/TT 也不知道就去采集软件里查采样率设置页面不要随手编一个数代入。5.2 用 findpeaks 自动提取峰值频率与幅值避免肉眼读数实测频谱峰不止一个手动找峰不可复现。findpeaks 需要给出一套明确的找峰规则Y fft(x); N length(x); P2 abs(Y) / N; P1 P2(1:N/21); P1(2:end-1) 2 * P1(2:end-1); f_axis (0:N/2) * fs / N; [peaks, locs] findpeaks(P1, f_axis, ... MinPeakHeight, 0.05, ... % 低于 0.05 的峰不认 MinPeakDistance, 2); % 两个峰至少相距 2 Hz plot(f_axis, P1); hold on; plot(locs, peaks, ro);MinPeakHeight 过滤噪底毛刺MinPeakDistance 按 df 和经验值压制伪峰。两个参数要结合谱图来回调峰太多就抬高 MinPeakHeight漏了小峰值就调低。locs 和 peaks 可以直接对接后续需求——做振动频谱图分析时常要按频带统计能量把分频带边界写成数组对 P1 做 trapz 积分即可findpeaks 只负责定位不负责计量。5.3 分段频谱图一条兼顾分辨率与噪声的验证技巧实测数据很长时一次性做整段 FFT 分辨率虽高所有随机噪声也都保留在谱里。分段平均是把长数据切段、分别加窗做 FFT、再把幅值谱取平均seg 4096; % 每段长度 nseg floor(length(x) / seg); S zeros(seg/21, 1); for k 1:nseg xk x((k-1)*seg1 : k*seg); Yk fft(xk .* hann(seg, periodic)); S S abs(Yk(1:seg/21)); end S S / nseg; fseg (0:seg/2) * fs / seg;分段数 nseg 越多噪底越平稳但每段变短分辨率随之下降。调试时先取 seg 为总长的 1/10 画一张再取 1/100 画一张两张对比通常就能找到峰不糊、底不毛的折中段长——这个选择习惯比记住任何经验公式都可靠。仿真阶段就值得练一练给信号随机加噪分别画整段频谱与分段频谱观察噪底曲线的粗糙度差异你就知道为什么商业振动分析软件几乎都采用平均谱作为默认显示。本文还有配套的精品资源点击获取