简介面向电子测量、信号处理方向学习者这份基于MATLAB GUI的频谱分析仪设计文档系统介绍如何借助虚拟仪器技术观测电信号频谱结构并测量失真度、调制度、谱纯度等参数适合通信、电子信息类学生完成课程设计或快速入门频谱分析原理。文档为doc格式压缩包共1个文件、约625KB正文按概述、技术路线、实现方法三大模块展开实现部分依次讲解GUI界面搭建、多路信号输入选择信号、声卡输入、读取wav文件、信号发生器、时域特征计算与基于FFT的频域分析同时给出Simulink仿真思路。读者可跟随文档逐步搭建简易虚拟频谱分析仪掌握界面布局、信号采集与频谱绘制的完整流程文档还包含作者梳理的现存问题、致谢与参考文献有助于进一步拓展和优化设计。目前已有867人学习可作为电子测量与虚拟仪器方向的一份实用参考资料。1. 频谱分析仪在 MATLAB 里是什么很多人以为把波形丢给 fft 就是频谱分析仪了其实那只是把时域数据变换到频域而已。仪器级的频谱分析还要解决怎么采样、怎么截断、怎么把谱线和真实频率对应起来。这个标题说的 matlab 频谱分析仪是把 MATLAB 当成一台能出数的仪器给定信号的采样率与分析时长加对窗做单边谱再把频率、幅度、相位从谱线里读出来。适合处理设备振动测试、电网谐波、控制器噪声、射频基带信号这类任务的人。如果你只看阶跃响应和时域指标频谱是另一双眼睛一旦开始盯谱线这套链路值得完整过一遍。这里不讨论按钮在哪而是把能带走的那套参数和脚本讲清楚。2. 从采样到 FFT搭起频谱分析的最小链路2.1 先定三个基本量采样率、分析时长、频率分辨率我把频谱分析拆成三个基本量采样率 fs、分析时长 T、FFT 点数 N。它们不是独立选的。采样率首先被信号带宽卡住按 Nyquist 至少两倍工程上一般取目标最高频率的 5 到 10 倍采样率太高不会带来更高频率分辨率只会让 N 变大内存和时间都上涨。分析时长 T 决定频率分辨率相邻两根谱线间隔 Δf 1/T写成 FFT 参数就是 fs/N。想分辨 1 Hz 的两根线至少要 1 秒数据想分辨 0.1 Hz至少 10 秒。这是绕不开的物理约束改 FFT 点数改变不了它。采样率 fs分析时长 TFFT 点数 N频率分辨率 Δf最高分析频率10000 Hz0.1 s100010 Hz5000 Hz10000 Hz1 s100001 Hz5000 Hz10000 Hz10 s1000000.1 Hz5000 Hz最高分析频率由 fs/2 决定与 N 无关。我一般会把目标谱线的最小间距除以 5 到 10再决定分析时长50 Hz 和 55 Hz 两条线至少要 1 到 2 秒数据Δf 才会明显小于它们间距。注意很多版本的采集软件里能设置“采样点数”它决定的是 N而不是分辨率分辨率看的是这一帧数据对应了多长的真实时间。若信号源自带时基漂移标称 fs 和实际采样率有偏差长数据的谱线位置会向右偏这是后话。2.2 用 fft 而不是手工 DFT最小频谱代码先跑通最小链路再讨论窗。假设你有一段 1 秒的仿真信号10 kHz 采样包含 123.4 Hz 和 2000 Hz 两个正弦fs 10000; % 采样率 10 kHz t (0:fs*1-1) / fs; % 1 秒时间向量 x 1.2*sin(2*pi*123.4*t 0.5) ... 0.8*sin(2*pi*2000*t) ... 0.05*randn(size(t)); % 加一点噪声 N length(x); % 10000 X fft(x); % 复频谱 plot(abs(X)); % 先看一眼原始输出fft 是 DFT 的快速实现输出长度和输入相同结果是复数包含幅度和相位。X(1) 对应 0 HzX(k) 对应 (k-1)*fs/N。这一版直接画 abs(X) 的问题很明显横轴是谱线索引不是频率纵轴幅度是错误放大的0 Hz 处还有一根冲天谱线。频率轴用类似 linspace 的函数构造也可以但我习惯用 (0:half-1)*fs/N下标从 0 开始算不容易漏点。2.3 画图前必须处理的三件事直流、单边、幅值修正继续上面的 X在画图之前我一般会做三步修正x x - mean(x); % 1) 去直流 N numel(x); X fft(x); half floor(N/2) 1; % 单边谱长度 mag abs(X(1:half)); mag(2:end-1) 2 * mag(2:end-1); % 2) 正负频率合并 f (0:half-1) * fs / N; % 3) 频率轴标定 plot(f, mag); xlim([0 2500]); grid on;去直流是最容易漏的一步。ADC 或传感器输出通常带一个直流偏置不去掉的话 0 Hz 处会有个大峰它的泄漏足以把几十 Hz 以内的小信号淹没。单边谱把正负频率的能量合回来幅度谱才读得出真实幅值但 0 Hz 和最高频率点没有镜像不能乘 2。频率轴标定后横轴才是能直接给别人看的 Hz。这段代码跑完124 Hz 附近的主峰应该接近 1.22000 Hz 那个峰接近 0.8而 123.4 Hz 峰的两侧会有一串拖尾这就是下一章要解决的泄漏。MATLAB 画图时如果低频段谱线太多糊在一起先用 xlim 限频带再用 xticks 给刻度不要试图通过改 N 来“变清晰”——那是两回事。3. 关键参数与窗函数频谱泄漏怎么压下去3.1 泄漏的来源截断等价于乘矩形窗FFT 在数学上把这 N 点数据当作周期信号的一个周期来处理要求首尾能无缝拼接。实际采样很少正好截到整数个周期这段信号里 123.4 Hz 在 1 秒内是 123.4 个周期边界必然不连续。截断本质上是给无限长信号乘了一个矩形窗频域相当于原信号谱和 sinc 卷积主瓣变宽旁瓣一串往下掉而 sinc 旁瓣衰减只有 -13 dB。主信号 1 V 时旁边 10 mV 的小信号会被旁瓣裙边盖住这就是泄漏。它不是你操作错而是截断本身带进来的解决手段是换窗函数。3.2 常见窗函数的选型参数表MATLAB 里用 window(winname, N) 或直接调 hann、hamming、blackman、flattopwin 都能生成窗。常用窗的参数如下窗函数MATLAB 函数主瓣宽度旁瓣衰减幅度精度常用场景矩形rectwin2 bins-13 dB高瞬态信号或谱线间距极近Hannhann4 bins-31 dB好通用首选Hamminghamming4 bins-43 dB好旁瓣敏感且主瓣要窄Blackmanblackman6 bins-58 dB中动态范围高的测量Flattopflattopwin8 bins很低最高只测幅值不测频率bin 是一根谱线。主瓣越宽两个相近频率越难分开旁瓣越低大信号对小信号的掩盖越轻。选窗就是在分辨力和动态范围之间权衡不知道选什么就用 Hann要分辨靠得近的谱线用矩形或 Hann单测幅值用 Flattop对旁瓣敏感比如测谐波时旁边有个很强的基波用 Blackman。表中的数值是近似值具体到每个窗变体会略有出入。3.3 加窗后的幅度修正与两种修正的区别加窗会把信号两端压小FFT 峰值幅度整体变低必须修正。幅度谱用的是平均增益修正w hann(N, periodic); % 频谱分析推荐周期窗 sc sum(w) / N; % 平均增益约 0.5 xw (x - mean(x)) .* w; Xw fft(xw); magw abs(Xw(1:floor(N/2)1)); magw(2:end-1) 2 * magw(2:end-1); magw magw / sc; % 恢复真实幅度sc 是窗函数的平均值不是最大值。Hann 窗的 sum(w) 接近 N/2sc 约 0.5所以恢复时要把幅度乘 2这一条经常有人漏掉。hann(N,periodic) 与默认的 symmetric 变体在高精度测量里有区别频谱分析我一般用 periodic。修正之后 magw 的峰值才接近真实幅值 1.2 和 0.8。注意这里是幅度谱修正如果后面要算功率谱密度得用 ENBW 而不是平均增益这张表放到第 5 章。4. 读谱线从频谱里提取频率、幅度和相位4.1 用 findpeaks 自动找峰频谱画出来之后手工读峰值坐标只能糊弄一两次。我一般直接在幅度谱上跑 findpeaks[pks, idx] findpeaks(magw, ... MinPeakHeight, 0.05, ... % 幅度阈值 MinPeakDistance, 5); % 至少间隔 5 根谱线 fpeaks (idx - 1) * fs / N;MinPeakHeight 按你要报的最小信号幅度设不确定时就先取 max(magw)*0.02避免把噪声峰当结果。MinPeakDistance 的单位是谱线根数用来把窗函数旁瓣产生的假峰挡掉和频率间隔的换算是 ceil(minFreqSpan / Δf)。如果背景噪声起伏很大MinPeakHeight 不好定改用 MinPeakProminence 会更稳定它要求峰比周围突出一定倍数。这里我没有把频率轴 f 传进 findpeaks先拿谱线索引再换算 Hz因为不同版本对第三参的“坐标单位”解释不完全一致先索引后换算最省心。下表是 findpeaks 在频谱分析里最常用的三个参数参数作用常见取值MinPeakHeight幅度阈值低于此值不认为有峰max(magw)*0.02MinPeakDistance两个峰之间的最小谱线根数去掉旁瓣假峰ceil(最小频率间隔/Δf)MinPeakProminence峰相对周围地形的突出度背景起伏大时更稳约目标幅度的 0.5~1 倍4.2 抛物线插值和 Goertzel 精测把频率读到亚谱线精度FFT 是离散网格真实峰值经常落在两根谱线之间。三点抛物线插值可以估计偏移量 d再把频率修正到亚谱线精度k idx(1); % 峰值索引 if k 1 k length(magw) d 0.5 * (magw(k-1) - magw(k1)) / ... (magw(k-1) - 2*magw(k) magw(k1)); f0 (k - 1 d) * fs / N; % 修正后的频率 end Xk goertzel(xw, round(f0 / (fs/N)) 1); % 单点 DFT A0 2 * abs(Xk) / (N * sc); % 精测幅度 phi0 angle(Xk); % 相位d 应该在 ±0.5 之间超出范围说明旁边有干扰峰插值结果不可信。goertzel 相当于单点 DFT在已知频率附近精测比整段 FFT 省内存幅度也不再受“峰值恰好落在两根谱线之间”的影响goertzel 的索引同样是 1-based第 1 点对应 0 Hz所以 round 之后要 1。如果目标频率附近有其他谱线干扰抛物线会偏这时才需要把窗函数频响建模后做拟合用优化工具箱的 lsqcurvefit 能解但一般场景三点插值已经足够。4.3 相位读取和两路信号的相位差无窗 FFT 在峰值处的 angle 就是初相但加窗之后频谱相位混入了窗的群延迟直接读 angle 会得到和频率成正比的大偏移。要还原绝对初相可以补偿 2pif0*(N-1)/2/fs符号要按窗的对称方式确定很容易错。我一般建议工程上改测相对相位两路信号用同样的窗、同一段起点对齐做 FFTangle(X1) - angle(X2) 会消掉公共群延迟得到的就是该频率点的相位差。如果两路采样有固定延迟记得先把时延标定掉否则相位差会随频率线性变化。5. 排错与底噪频率轴、补零和功率谱密度的坑5.1 频率轴标定和绘图显示还有原始数据的字节序最常见的错误是把 FFT 结果的横轴当成样本点。正确做法是用真实采样率换算f (0:floor(N/2)) * fs / N; plot(f, magw); xlim([0 500]); % 只看低频段 xticks(0:50:500);显示范围太大、谱线太密时低频部分会糊成一团这和采集时间无关是显示问题先限制频带再看。比如设备振动分析常常只关心 0~500 Hz直接 xlim 解决。如果是嵌入式采集卡拿回来的裸数据常见是 16 进制字符串或者 uint16/int16 整数进 FFT 之前必须按字节序转成有符号数并按满量程标定raw fread(fid, inf, uint16); x double(typecast(uint16(raw), int16)) / 32768;typecast 按本机字节序解释跨平台传输时要先确认文件是大端还是小端必要时用 swapbytes。漏掉这一步频谱里会冒出一条虚假直流和镜像频率排错时非常难发现。我一般会在去直流之前先打印 min/max确认数据范围正常再做后续处理。5.2 补零改变不了分辨率栅栏效应与 zoom FFT有时觉得谱线太粗就有人把数据补一堆零再 FFTN numel(x); xpad [x; zeros(7*N, 1)]; % 补到 8 倍长度 Xpad fft(xpad); fpad (0:floor(numel(xpad)/2)) / numel(xpad) * fs;补零后曲线变得圆滑看起来好看了但两个真实频率能否被分开仍然由原始数据时长决定补零只是插值没有增加任何信息。原始 1 秒数据对应的 Δf 还是 1 Hz只是谱线在 1 Hz 网格之间多画了几个点。真正想要更细的分辨率只有两种情况加长分析时长或者在窄带内做 zoom FFT——把目标频带搬到零频、低通滤波、再降采样然后 FFT。降采样 D 倍之后同样的输出点数覆盖的带宽变成 fs/D等效于把该频带展宽 D 倍但需要更长的原始记录才能保住谱线数。如果硬件资源够加长数据直接长 FFT 是更省事的办法zoom FFT 主要给数据量受限的场景。5.3 噪声底和功率谱密度幅度谱会骗你PSD 不会幅度谱的噪声底会随 FFT 点数变化点数变大噪声能量分散到更多谱线上底看起来更低容易给人“换了大 FFT 噪声变小”的错觉。要看噪声的真实水平用功率谱密度[psd, f] pwelch(x, hann(N/4, periodic), [], N/4, fs); plot(f, 10*log10(psd)); ylabel(Power/frequency (dB/Hz));pwelch 把数据分段做平均段数越多噪声曲线越平滑但频率分辨率变差。窗长取 N/4 表示四段平均具体取多少取决于你要平滑度还是要分辨力。PSD 归一化到每 Hz 带宽读数不再随 N 变化多个窗之间也能直接比较。比较时注意窗的 ENBW 不同窗函数ENBWbins矩形1.00Hann1.50Hamming1.36Blackman1.73Flattop3.77同一段噪声用 Hann 和矩形窗得到的 PSD 底就差 10*log10(1.5) ≈ 1.76 dB这是窗的有效带宽造成的正常差异不是仪器坏了。这个常数在对比底噪测试结果时经常被忽略。6. 把频谱分析仪沉淀成一个可复用函数6.1 一个封装了去直流、加窗、修正、找峰的函数骨架前面几节的步骤每次都要重打一遍我习惯收进一个函数后面任何数据进来一行调用function result specAna(x, fs, windowName, minPeakHeight) if nargin 3 || isempty(windowName), windowName hann; end if nargin 4, minPeakHeight []; end x x(:) - mean(x); % 去直流 N numel(x); w window(windowName, N); % 依赖 Signal Processing Toolbox w w(:); sc sum(w) / N; % 幅度修正系数 X fft(x .* w); half floor(N/2) 1; amp abs(X(1:half)); amp(2:end-1) 2 * amp(2:end-1); amp amp / sc; f (0:half-1) * fs / N; result.f f; result.amp amp; if ~isempty(minPeakHeight) [result.pks, idx] findpeaks(amp, ... MinPeakHeight, minPeakHeight, ... MinPeakDistance, 5); result.fpeaks (idx - 1) * fs / N; end endwindow 和 findpeaks 都来自信号处理工具箱没有工具箱时Hann 窗可以手动写 0.5*(1-cos(2pi(0:N-1)/(N-1)))。minPeakDistance 固定 5 根谱线要检测更近的峰时按 ceil(minFreqSpan/Δf) 改掉。返回的 result 是一个结构体批量处理时把它放进数组即可。6.2 两条使用路线交互核查和批量出报告交互式查看时直接开内置频谱分析仪spectrumAnalyzer(x);新版本对应 spectrumAnalyzer老版本里是 dsp.SpectrumAnalyzer界面里可以切窗函数、峰值标注和单位适合边调边看。批量出报告就走上面的函数r specAna(x, fs, hann, 0.1); T table(r.fpeaks, r.pks, VariableNames, {Freq_Hz, Amp}); writetable(T, peaks.csv); figure; plot(r.f, r.amp); xlim([0 500]); grid on; print(gcf, -depsc2, spectrum.eps); % 矢量图print 导出的是 eps 矢量图写报告缩放不糊比截图强。批量分析上百个文件时把 specAna 放进 for 循环逐帧把峰值表累积起来最后再统一绘图交互式核查交给 App批量出数交给这个骨架这套组合能覆盖从振动测试到谐波分析的多数频谱任务。本文还有配套的精品资源点击获取 
长河君 · 2026/9/17 16:12:15 · 阅读 1.2k MATLAB频谱分析仪实战：核心参数、窗函数与谱线读取
简介面向电子测量、信号处理方向学习者这份基于MATLAB GUI的频谱分析仪设计文档系统介绍如何借助虚拟仪器技术观测电信号频谱结构并测量失真度、调制度、谱纯度等参数适合通信、电子信息类学生完成课程设计或快速入门频谱分析原理。文档为doc格式压缩包共1个文件、约625KB正文按概述、技术路线、实现方法三大模块展开实现部分依次讲解GUI界面搭建、多路信号输入选择信号、声卡输入、读取wav文件、信号发生器、时域特征计算与基于FFT的频域分析同时给出Simulink仿真思路。读者可跟随文档逐步搭建简易虚拟频谱分析仪掌握界面布局、信号采集与频谱绘制的完整流程文档还包含作者梳理的现存问题、致谢与参考文献有助于进一步拓展和优化设计。目前已有867人学习可作为电子测量与虚拟仪器方向的一份实用参考资料。1. 频谱分析仪在 MATLAB 里是什么很多人以为把波形丢给 fft 就是频谱分析仪了其实那只是把时域数据变换到频域而已。仪器级的频谱分析还要解决怎么采样、怎么截断、怎么把谱线和真实频率对应起来。这个标题说的 matlab 频谱分析仪是把 MATLAB 当成一台能出数的仪器给定信号的采样率与分析时长加对窗做单边谱再把频率、幅度、相位从谱线里读出来。适合处理设备振动测试、电网谐波、控制器噪声、射频基带信号这类任务的人。如果你只看阶跃响应和时域指标频谱是另一双眼睛一旦开始盯谱线这套链路值得完整过一遍。这里不讨论按钮在哪而是把能带走的那套参数和脚本讲清楚。2. 从采样到 FFT搭起频谱分析的最小链路2.1 先定三个基本量采样率、分析时长、频率分辨率我把频谱分析拆成三个基本量采样率 fs、分析时长 T、FFT 点数 N。它们不是独立选的。采样率首先被信号带宽卡住按 Nyquist 至少两倍工程上一般取目标最高频率的 5 到 10 倍采样率太高不会带来更高频率分辨率只会让 N 变大内存和时间都上涨。分析时长 T 决定频率分辨率相邻两根谱线间隔 Δf 1/T写成 FFT 参数就是 fs/N。想分辨 1 Hz 的两根线至少要 1 秒数据想分辨 0.1 Hz至少 10 秒。这是绕不开的物理约束改 FFT 点数改变不了它。采样率 fs分析时长 TFFT 点数 N频率分辨率 Δf最高分析频率10000 Hz0.1 s100010 Hz5000 Hz10000 Hz1 s100001 Hz5000 Hz10000 Hz10 s1000000.1 Hz5000 Hz最高分析频率由 fs/2 决定与 N 无关。我一般会把目标谱线的最小间距除以 5 到 10再决定分析时长50 Hz 和 55 Hz 两条线至少要 1 到 2 秒数据Δf 才会明显小于它们间距。注意很多版本的采集软件里能设置“采样点数”它决定的是 N而不是分辨率分辨率看的是这一帧数据对应了多长的真实时间。若信号源自带时基漂移标称 fs 和实际采样率有偏差长数据的谱线位置会向右偏这是后话。2.2 用 fft 而不是手工 DFT最小频谱代码先跑通最小链路再讨论窗。假设你有一段 1 秒的仿真信号10 kHz 采样包含 123.4 Hz 和 2000 Hz 两个正弦fs 10000; % 采样率 10 kHz t (0:fs*1-1) / fs; % 1 秒时间向量 x 1.2*sin(2*pi*123.4*t 0.5) ... 0.8*sin(2*pi*2000*t) ... 0.05*randn(size(t)); % 加一点噪声 N length(x); % 10000 X fft(x); % 复频谱 plot(abs(X)); % 先看一眼原始输出fft 是 DFT 的快速实现输出长度和输入相同结果是复数包含幅度和相位。X(1) 对应 0 HzX(k) 对应 (k-1)*fs/N。这一版直接画 abs(X) 的问题很明显横轴是谱线索引不是频率纵轴幅度是错误放大的0 Hz 处还有一根冲天谱线。频率轴用类似 linspace 的函数构造也可以但我习惯用 (0:half-1)*fs/N下标从 0 开始算不容易漏点。2.3 画图前必须处理的三件事直流、单边、幅值修正继续上面的 X在画图之前我一般会做三步修正x x - mean(x); % 1) 去直流 N numel(x); X fft(x); half floor(N/2) 1; % 单边谱长度 mag abs(X(1:half)); mag(2:end-1) 2 * mag(2:end-1); % 2) 正负频率合并 f (0:half-1) * fs / N; % 3) 频率轴标定 plot(f, mag); xlim([0 2500]); grid on;去直流是最容易漏的一步。ADC 或传感器输出通常带一个直流偏置不去掉的话 0 Hz 处会有个大峰它的泄漏足以把几十 Hz 以内的小信号淹没。单边谱把正负频率的能量合回来幅度谱才读得出真实幅值但 0 Hz 和最高频率点没有镜像不能乘 2。频率轴标定后横轴才是能直接给别人看的 Hz。这段代码跑完124 Hz 附近的主峰应该接近 1.22000 Hz 那个峰接近 0.8而 123.4 Hz 峰的两侧会有一串拖尾这就是下一章要解决的泄漏。MATLAB 画图时如果低频段谱线太多糊在一起先用 xlim 限频带再用 xticks 给刻度不要试图通过改 N 来“变清晰”——那是两回事。3. 关键参数与窗函数频谱泄漏怎么压下去3.1 泄漏的来源截断等价于乘矩形窗FFT 在数学上把这 N 点数据当作周期信号的一个周期来处理要求首尾能无缝拼接。实际采样很少正好截到整数个周期这段信号里 123.4 Hz 在 1 秒内是 123.4 个周期边界必然不连续。截断本质上是给无限长信号乘了一个矩形窗频域相当于原信号谱和 sinc 卷积主瓣变宽旁瓣一串往下掉而 sinc 旁瓣衰减只有 -13 dB。主信号 1 V 时旁边 10 mV 的小信号会被旁瓣裙边盖住这就是泄漏。它不是你操作错而是截断本身带进来的解决手段是换窗函数。3.2 常见窗函数的选型参数表MATLAB 里用 window(winname, N) 或直接调 hann、hamming、blackman、flattopwin 都能生成窗。常用窗的参数如下窗函数MATLAB 函数主瓣宽度旁瓣衰减幅度精度常用场景矩形rectwin2 bins-13 dB高瞬态信号或谱线间距极近Hannhann4 bins-31 dB好通用首选Hamminghamming4 bins-43 dB好旁瓣敏感且主瓣要窄Blackmanblackman6 bins-58 dB中动态范围高的测量Flattopflattopwin8 bins很低最高只测幅值不测频率bin 是一根谱线。主瓣越宽两个相近频率越难分开旁瓣越低大信号对小信号的掩盖越轻。选窗就是在分辨力和动态范围之间权衡不知道选什么就用 Hann要分辨靠得近的谱线用矩形或 Hann单测幅值用 Flattop对旁瓣敏感比如测谐波时旁边有个很强的基波用 Blackman。表中的数值是近似值具体到每个窗变体会略有出入。3.3 加窗后的幅度修正与两种修正的区别加窗会把信号两端压小FFT 峰值幅度整体变低必须修正。幅度谱用的是平均增益修正w hann(N, periodic); % 频谱分析推荐周期窗 sc sum(w) / N; % 平均增益约 0.5 xw (x - mean(x)) .* w; Xw fft(xw); magw abs(Xw(1:floor(N/2)1)); magw(2:end-1) 2 * magw(2:end-1); magw magw / sc; % 恢复真实幅度sc 是窗函数的平均值不是最大值。Hann 窗的 sum(w) 接近 N/2sc 约 0.5所以恢复时要把幅度乘 2这一条经常有人漏掉。hann(N,periodic) 与默认的 symmetric 变体在高精度测量里有区别频谱分析我一般用 periodic。修正之后 magw 的峰值才接近真实幅值 1.2 和 0.8。注意这里是幅度谱修正如果后面要算功率谱密度得用 ENBW 而不是平均增益这张表放到第 5 章。4. 读谱线从频谱里提取频率、幅度和相位4.1 用 findpeaks 自动找峰频谱画出来之后手工读峰值坐标只能糊弄一两次。我一般直接在幅度谱上跑 findpeaks[pks, idx] findpeaks(magw, ... MinPeakHeight, 0.05, ... % 幅度阈值 MinPeakDistance, 5); % 至少间隔 5 根谱线 fpeaks (idx - 1) * fs / N;MinPeakHeight 按你要报的最小信号幅度设不确定时就先取 max(magw)*0.02避免把噪声峰当结果。MinPeakDistance 的单位是谱线根数用来把窗函数旁瓣产生的假峰挡掉和频率间隔的换算是 ceil(minFreqSpan / Δf)。如果背景噪声起伏很大MinPeakHeight 不好定改用 MinPeakProminence 会更稳定它要求峰比周围突出一定倍数。这里我没有把频率轴 f 传进 findpeaks先拿谱线索引再换算 Hz因为不同版本对第三参的“坐标单位”解释不完全一致先索引后换算最省心。下表是 findpeaks 在频谱分析里最常用的三个参数参数作用常见取值MinPeakHeight幅度阈值低于此值不认为有峰max(magw)*0.02MinPeakDistance两个峰之间的最小谱线根数去掉旁瓣假峰ceil(最小频率间隔/Δf)MinPeakProminence峰相对周围地形的突出度背景起伏大时更稳约目标幅度的 0.5~1 倍4.2 抛物线插值和 Goertzel 精测把频率读到亚谱线精度FFT 是离散网格真实峰值经常落在两根谱线之间。三点抛物线插值可以估计偏移量 d再把频率修正到亚谱线精度k idx(1); % 峰值索引 if k 1 k length(magw) d 0.5 * (magw(k-1) - magw(k1)) / ... (magw(k-1) - 2*magw(k) magw(k1)); f0 (k - 1 d) * fs / N; % 修正后的频率 end Xk goertzel(xw, round(f0 / (fs/N)) 1); % 单点 DFT A0 2 * abs(Xk) / (N * sc); % 精测幅度 phi0 angle(Xk); % 相位d 应该在 ±0.5 之间超出范围说明旁边有干扰峰插值结果不可信。goertzel 相当于单点 DFT在已知频率附近精测比整段 FFT 省内存幅度也不再受“峰值恰好落在两根谱线之间”的影响goertzel 的索引同样是 1-based第 1 点对应 0 Hz所以 round 之后要 1。如果目标频率附近有其他谱线干扰抛物线会偏这时才需要把窗函数频响建模后做拟合用优化工具箱的 lsqcurvefit 能解但一般场景三点插值已经足够。4.3 相位读取和两路信号的相位差无窗 FFT 在峰值处的 angle 就是初相但加窗之后频谱相位混入了窗的群延迟直接读 angle 会得到和频率成正比的大偏移。要还原绝对初相可以补偿 2pif0*(N-1)/2/fs符号要按窗的对称方式确定很容易错。我一般建议工程上改测相对相位两路信号用同样的窗、同一段起点对齐做 FFTangle(X1) - angle(X2) 会消掉公共群延迟得到的就是该频率点的相位差。如果两路采样有固定延迟记得先把时延标定掉否则相位差会随频率线性变化。5. 排错与底噪频率轴、补零和功率谱密度的坑5.1 频率轴标定和绘图显示还有原始数据的字节序最常见的错误是把 FFT 结果的横轴当成样本点。正确做法是用真实采样率换算f (0:floor(N/2)) * fs / N; plot(f, magw); xlim([0 500]); % 只看低频段 xticks(0:50:500);显示范围太大、谱线太密时低频部分会糊成一团这和采集时间无关是显示问题先限制频带再看。比如设备振动分析常常只关心 0~500 Hz直接 xlim 解决。如果是嵌入式采集卡拿回来的裸数据常见是 16 进制字符串或者 uint16/int16 整数进 FFT 之前必须按字节序转成有符号数并按满量程标定raw fread(fid, inf, uint16); x double(typecast(uint16(raw), int16)) / 32768;typecast 按本机字节序解释跨平台传输时要先确认文件是大端还是小端必要时用 swapbytes。漏掉这一步频谱里会冒出一条虚假直流和镜像频率排错时非常难发现。我一般会在去直流之前先打印 min/max确认数据范围正常再做后续处理。5.2 补零改变不了分辨率栅栏效应与 zoom FFT有时觉得谱线太粗就有人把数据补一堆零再 FFTN numel(x); xpad [x; zeros(7*N, 1)]; % 补到 8 倍长度 Xpad fft(xpad); fpad (0:floor(numel(xpad)/2)) / numel(xpad) * fs;补零后曲线变得圆滑看起来好看了但两个真实频率能否被分开仍然由原始数据时长决定补零只是插值没有增加任何信息。原始 1 秒数据对应的 Δf 还是 1 Hz只是谱线在 1 Hz 网格之间多画了几个点。真正想要更细的分辨率只有两种情况加长分析时长或者在窄带内做 zoom FFT——把目标频带搬到零频、低通滤波、再降采样然后 FFT。降采样 D 倍之后同样的输出点数覆盖的带宽变成 fs/D等效于把该频带展宽 D 倍但需要更长的原始记录才能保住谱线数。如果硬件资源够加长数据直接长 FFT 是更省事的办法zoom FFT 主要给数据量受限的场景。5.3 噪声底和功率谱密度幅度谱会骗你PSD 不会幅度谱的噪声底会随 FFT 点数变化点数变大噪声能量分散到更多谱线上底看起来更低容易给人“换了大 FFT 噪声变小”的错觉。要看噪声的真实水平用功率谱密度[psd, f] pwelch(x, hann(N/4, periodic), [], N/4, fs); plot(f, 10*log10(psd)); ylabel(Power/frequency (dB/Hz));pwelch 把数据分段做平均段数越多噪声曲线越平滑但频率分辨率变差。窗长取 N/4 表示四段平均具体取多少取决于你要平滑度还是要分辨力。PSD 归一化到每 Hz 带宽读数不再随 N 变化多个窗之间也能直接比较。比较时注意窗的 ENBW 不同窗函数ENBWbins矩形1.00Hann1.50Hamming1.36Blackman1.73Flattop3.77同一段噪声用 Hann 和矩形窗得到的 PSD 底就差 10*log10(1.5) ≈ 1.76 dB这是窗的有效带宽造成的正常差异不是仪器坏了。这个常数在对比底噪测试结果时经常被忽略。6. 把频谱分析仪沉淀成一个可复用函数6.1 一个封装了去直流、加窗、修正、找峰的函数骨架前面几节的步骤每次都要重打一遍我习惯收进一个函数后面任何数据进来一行调用function result specAna(x, fs, windowName, minPeakHeight) if nargin 3 || isempty(windowName), windowName hann; end if nargin 4, minPeakHeight []; end x x(:) - mean(x); % 去直流 N numel(x); w window(windowName, N); % 依赖 Signal Processing Toolbox w w(:); sc sum(w) / N; % 幅度修正系数 X fft(x .* w); half floor(N/2) 1; amp abs(X(1:half)); amp(2:end-1) 2 * amp(2:end-1); amp amp / sc; f (0:half-1) * fs / N; result.f f; result.amp amp; if ~isempty(minPeakHeight) [result.pks, idx] findpeaks(amp, ... MinPeakHeight, minPeakHeight, ... MinPeakDistance, 5); result.fpeaks (idx - 1) * fs / N; end endwindow 和 findpeaks 都来自信号处理工具箱没有工具箱时Hann 窗可以手动写 0.5*(1-cos(2pi(0:N-1)/(N-1)))。minPeakDistance 固定 5 根谱线要检测更近的峰时按 ceil(minFreqSpan/Δf) 改掉。返回的 result 是一个结构体批量处理时把它放进数组即可。6.2 两条使用路线交互核查和批量出报告交互式查看时直接开内置频谱分析仪spectrumAnalyzer(x);新版本对应 spectrumAnalyzer老版本里是 dsp.SpectrumAnalyzer界面里可以切窗函数、峰值标注和单位适合边调边看。批量出报告就走上面的函数r specAna(x, fs, hann, 0.1); T table(r.fpeaks, r.pks, VariableNames, {Freq_Hz, Amp}); writetable(T, peaks.csv); figure; plot(r.f, r.amp); xlim([0 500]); grid on; print(gcf, -depsc2, spectrum.eps); % 矢量图print 导出的是 eps 矢量图写报告缩放不糊比截图强。批量分析上百个文件时把 specAna 放进 for 循环逐帧把峰值表累积起来最后再统一绘图交互式核查交给 App批量出数交给这个骨架这套组合能覆盖从振动测试到谐波分析的多数频谱任务。本文还有配套的精品资源点击获取