先交代一个真实场景我之前处理过一批机械振动数据时间域的波形密密麻麻肉眼根本看不出异常。随手用MATLAB画了个频谱图故障特征一下就跳出来了。也是从那个时候起我意识到频谱和功率谱不是科研专用工具而是每个做信号相关工作的人都该掌握的基础语言。但问题来了——身边不少朋友拿着MATLAB却画不出一张规范的功率谱图网上搜到的代码要么只有片段要么画完图连坐标轴单位都不敢写横轴到底是Hz还是rad/s纵轴幅度该不该乘2用fft算出来的和用pwelch算出来的差了好几个数量级到底哪个对这篇内容我就把一套完整可用的MATLAB频谱和功率谱画图程序拆开来从手算FFT到调用pwelch每一步都讲清楚为什么这样做再把我在实际调试中踩过的坑直接列出来。无论是处理声音、振动、电力谐波还是做噪声测试你都可以把这里的代码当成一个工具箱起点按需修改。内容比较多建议收藏后跟着敲一遍。1. 为什么总要把时域信号翻到频域里看1.1 一段波形藏不住的信息我们拿到的原始信号绝大多数是时间域的采样序列横轴是时间纵轴是幅值。时域波形只告诉你信号在什么时候变化多剧烈却不直接告诉你这个变化是由哪些频率组成的。比如一个由50Hz和120Hz正弦波叠加出来的信号在时域里看就是一条扭曲的曲线单凭肉眼很难拆开成分而在频域里看就是两个尖锐的谱峰清清楚楚。傅里叶变换做的事本质上就是换一副眼镜看信号——把信号看成无数个不同频率、不同幅度、不同初相位的正弦波的叠加。你用fft、pwelch画频谱图本质都是在对这副眼镜输出结果做针对性处理。1.2 频谱、功率谱和功率谱密度别把三个名字当成一回事不少新手把频谱和功率谱混着叫实际工程项目里这三者内涵差很多选错了会让结果无法解释。幅值谱是FFT结果取模后除以点数N得到的纵轴单位与原始信号幅值一致。比如输入信号的幅值是5幅值谱在对应频率处的峰值就接近5。它保留了幅度信息适合直观判断某个频率成分有多强。功率谱这里取的是信号功率的概念常用abs(fft(y)).^2/N计算纵轴单位是信号幅值的平方。它不保相位更关注能量分布。**功率谱密度PSD**再进一步把功率谱除以频率分辨率等效于除以Fs得到单位频率带宽内的功率纵轴单位是信号幅值平方除以Hz例如振动测试里常用的g^2/Hz。PSD是连续谱分析的标准选择也是pwelch默认输出的东西。三者之间有明确换算关系画图前必须想清楚自己到底要哪一类。否则别人问你这个峰的能量是多少你的纵轴单位答案都可能说错。2. 手写一版FFT频谱频率轴、幅值修正和单双边谱2.1 MATLAB里最容易被误解的fft函数MATLAB的fft(y)返回的是复数序列长度为N。很多新手直接abs(fft(y))之后画图发现峰值不对、纵轴量级巨大这就是没做归一化。fft本身是离散傅里叶变换的实现未做任何幅值归一化。对长度为N的序列某个频率分量的复振幅在FFT结果里大致等于原始振幅的N/2倍双边谱形式因此要恢复真实幅值必须除以N。一个很常见的错误是看到别人代码里除以N/2就盲目照抄。实际上单边谱除数和双边谱处理方式不同这牵涉到谱的对称性问题下面专门说。2.2 频率轴到底怎么构造假设采样率Fs、采样点数NFFT结果对应的数字频率间隔是Fs/N。MATLAB的fft输出顺序是从直流0Hz开始一直到Fs实际上是Fs - Fs/N后半部分与前半部分关于奈奎斯特频率Fs/2共轭对称。如果信号是实信号我们只看0 ~ Fs/2这一段就够了这叫单边谱。取Y(1:N/21)即可。频率轴常用N length(y); f (0:N/2) * Fs / N; % 单边谱频率轴这里N/2要按整除情况处理N通常取偶数。如果N是奇数取(N1)/2边界条件略有不同工程上我习惯让N为2的幂次省心得多。2.3 幅值修正为什么有个2实信号的双边谱中正频率和负频率各占一半幅度。我们只看单边谱相当于只取了正频率部分那么能量会减半所以要把除直流分量和奈奎斯特频率之外的点幅值乘以2。代码实现是Y fft(y); Y_single Y(1:N/21); % 取单边 mag abs(Y_single) / N; % 先归一化 mag(2:end-1) 2 * mag(2:end-1); % 去掉直流和奈奎斯特后再乘2直流分量是信号的平均值它不会成对出现奈奎斯特频率点也只有一个。因此这两个位置不需要乘2。乘2之前务必搞清楚位置索引否则结果会有固定偏差。2.4 一个可运行的示例下面这段代码可以直接复制运行合成一个50Hz加120Hz再叠加噪声的信号做一次完整的单边幅值谱画图Fs 1000; % 采样率 1000Hz t (0:Fs-1)/Fs; % 采样1秒 y 5*sin(2*pi*50*t) 3*sin(2*pi*120*t) 1.5*randn(size(t)); N length(y); Y fft(y); Y_single Y(1:N/21); mag abs(Y_single)/N; mag(2:end-1) 2*mag(2:end-1); f (0:N/2)*Fs/N; plot(f, mag); xlabel(频率 (Hz)); ylabel(幅值); title(单边幅值谱); grid on;运行后你会看到50Hz处峰值接近5120Hz处接近3这就是正确的物理幅值。背景的随机噪声在基线上方铺了一层不规则的小幅值幅值很小且分散在频谱各处。3. 功率谱估计的正路周期图、自相关法和Welch方法3.1 周期图与它的毛病把幅值谱的纵轴取平方再除以采样率就得到功率谱密度。最朴素的方式是直接对整段数据做一次FFT然后取模平方归一化这就是周期图法Pxx abs(fft(y)).^2 / (N * Fs); Pxx_single(2:end-1) 2 * Pxx_single(2:end-1);周期图实现简单但很脆。数据只有一段随机噪声没有机会平均谱线起伏剧烈如果信号被截断频谱泄漏也会被完整地保留下来。实际处理时用周期图得到的结果往往像一把乱草几乎无法解释。我遇到过一位A同学用周期图画加速度传感器数据出来的谱峰到处都是换了不同数据段结果还完全不同。后来换成分段平均图一下就干净了——问题的本质就是单个周期图方差太大。3.2 Welch法怎么改善方差Welch方法是把长信号分成长度相同的若干段每段可能重叠分别计算周期图然后对多段求平均。平均能压低随机噪声带来的方差代价是频率分辨率下降因为分段后每段点数变少了。重叠的目的也很直观如果信号是平稳的相邻分段信息高度相关完全不用重叠也可以但实际信号常有时变性有一定重叠比如50%或75%能保留更多样本用于平均估计结果更稳。MATLAB的pwelch就是Welch方法的成熟实现[pxx_w, f_w] pwelch(y, hamming(256), 128, 1024, Fs);第一行返回的pxx_w是单边功率谱密度频率范围0到Fs/2。直接plot就能画出平滑的PSD曲线。3.3 pwelch的四个核心参数用pwelch时绝大多数人不知道参数怎么调。我建议把四个参数拆开理解window每个分段的窗函数及长度。比如hamming(256)表示分段长度256点加汉明窗。窗函数抑制截断泄漏但会增大主瓣宽度。noverlap相邻分段重叠点数。一般取分段长度的50%~75%。设得过高计算量大且分段并非独立设得过小平均次数少方差降不下来。nfftF.T点数建议大于等于分段长度但通常取2的幂次以便加速。nfft越大频率网格越密但物理分辨率不由它决定。Fs采样率只影响频率轴刻度和PSD的量纲换算不影响频谱形状。它们的典型组合如下表场景分段长度窗函数重叠率nfft长稳态信号1024hann/hamming50%1024或2048噪声测试4096hann50%4096瞬态信号512hann75%20484. 一份可以直接抄作业的完整画图程序4.1 封装成函数为了不重复造轮子我把常用三种谱幅值谱、功率谱、功率谱密度封装成一个函数。你只需要把数据和采样率传进去选择方法就能得到一个标准的单边谱绘图结果。function [f, spec] plot_spectrum(x, Fs, method, window, noverlap, nfft) % 功能绘制单边幅值谱/功率谱/功率谱密度 % 输入 % x - 输入信号列向量 % Fs - 采样率 Hz % method - amp 幅值谱pow 功率谱psd 功率谱密度 % window - 窗函数向量或点数可选 % noverlap - 重叠点数可选 % nfft - FFT点数可选 % 输出 % f - 频率轴 % spec - 对应的谱值 x x(:); N length(x); switch lower(method) case amp if nargin 4 || isempty(window) Y fft(x); f (0:floor(N/2)) * Fs / N; spec abs(Y(1:length(f))) / N; spec(2:end-1) 2*spec(2:end-1); else xw x .* window(:); Y fft(xw, nfft); f (0:floor(N/2)) * Fs / N; spec abs(Y(1:length(f))) / sum(window); spec(2:end-1) 2*spec(2:end-1); end case pow if nargin 4 || isempty(window) Y fft(x); f (0:floor(N/2)) * Fs / N; spec (abs(Y(1:length(f))).^2) / N^2; spec(2:end-1) 2*spec(2:end-1); else xw x .* window(:); Y fft(xw, nfft); f (0:floor(N/2)) * Fs / N; spec (abs(Y(1:length(f))).^2) / sum(window)^2; spec(2:end-1) 2*spec(2:end-1); end case psd if nargin 4 || isempty(window) Y fft(x); f (0:floor(N/2)) * Fs / N; spec (abs(Y(1:length(f))).^2) / (N * Fs); spec(2:end-1) 2*spec(2:end-1); else [spec, f] pwelch(x, window, noverlap, nfft, Fs); end otherwise error(method 仅支持 amp / pow / psd); end使用时比如Fs 1000; t (0:Fs*2-1)/Fs; y 5*sin(2*pi*50*t) 3*sin(2*pi*120*t) 1.5*randn(size(t)); % 幅值谱 [f, spec_amp] plot_spectrum(y, Fs, amp); subplot(2,1,1); plot(f, spec_amp); title(幅值谱); % Welch PSD [f_psd, spec_psd] plot_spectrum(y, Fs, psd, hamming(256), 128, 1024); subplot(2,1,2); plot(f_psd, 10*log10(spec_psd)); title(PSD (dB/Hz));这里PSD用dB显示方便观察从低频到高频的宽动态范围。是实测数据时我几乎都会把纵轴设置为10*log10(spec)因为线性坐标下微弱成分会被淹没。4.2 输出到底用幅值谱、功率谱还是PSD根据使用场景不同选谱的方法也不一样如果只想知道信号里有哪些频率幅值谱最直观能看到每个频率处的真实幅值大小。如果关心噪声环境的能量水平用功率谱或PSD更好因为它能比较不同频带的能量差异。如果需要对接国际标准或测试报告比如振动ISO标准、声学倍频程分析几乎一律用PSD单位明确且可与其他系统比对。4.3 多信号叠加对比怎么画实际中经常要对比两种工况比如机器正常时和故障时的振动谱。此时可以用同一个频率轴把两条PSD曲线画在同一张图上并用颜色和线型区分[f1, psd1] pwelch(y_normal, hamming(1024), 512, 1024, Fs); [f2, psd2] pwelch(y_fault, hamming(1024), 512, 1024, Fs); plot(f1, 10*log10(psd1), b-, LineWidth, 1.2); hold on; plot(f2, 10*log10(psd2), r--, LineWidth, 1.2); legend(正常,故障);两张图叠加时一定要保证分段参数一致否则频率分辨率不同谱峰高低无法直接比较。这条经验让我的对比实验少走好多弯路。5. 画图阶段最值得注意的几个坑5.1 窗函数是买一赠一的双刃剑窗函数能抑制频谱泄漏原信号被截断相当于乘了一个矩形窗矩形窗的频谱旁瓣很高会污染邻近频率。加汉宁窗或汉明窗后主瓣变宽旁瓣大幅下降。可窗函数也带来了一个副作用——你的信号实际上被调制了如果直接用abs(fft(x.*window))除以N幅值会变低。因此在做幅值谱时不能简单除以N而是要除以窗函数的平均值或和。上面函数里我用了sum(window)做归一化这对幅值谱是近似正确的。严格讲加窗后还要考虑窗函数对能量修正的影响所以如果想做精确的幅值分析最简单的办法是不加窗直接FFT并做频谱泄露修正如果必须加窗就明确自己是在做功率/能量分析而不是精确幅值测量。我自己做噪声测试时一律用汉宁窗加pwelch因为最终指标是PSD不对具体谱线的峰高做绝对化要求。一旦有人拿幅值谱跟理论值对比我会让他先把窗去掉。5.2 分辨率与补零别被补出来的细节骗了频率分辨率由数据长度决定Δf Fs / N。如果你有一秒数据采样率1000Hz则分辨率是1Hz。想在谱上区分两个相差0.5Hz的频率成分靠的是一秒以上的数据而不是把nfft设大。把nfft从1024改成8192只是在原有谱之间插入更多插值点让曲线看起来更平滑但它不会把原本重叠的两个峰分开。这个区分必须从数据采集端解决——拉长采样时长或者用更高分辨率的估计方法比如子空间法。我见过有人把补零当成提高分辨率的法宝结果在无法分开的峰之间画出了双峰特效然后拿着假峰去诊断。这事务必记住补零只是给图美容不能无中生有创造物理信息。5.3 单位换算与坐标轴习惯画完图横轴纵轴的物理意义必须明确。常见错误包括横轴用角频率rad/s与Hz相差2π倍忘了标注导致误读。纵轴是幅值谱却写成了dB。dB本身是比值对数应用在幅值谱时需要指定参考值否则没有意义。pwelch默认输出功率谱密度其单位为幅值单位^2/Hz有些软件叫PSD有些叫auto-spectrum不要弄混。建议画图时用xlabel(频率 (Hz))PSD纵轴用ylabel(PSD (g^2/Hz))或dB/Hz幅值谱用ylabel(幅值)。单位写清楚后面报告评审阶段能省很多沟通成本。5.4 用已知信号做结果验证拿到任何新谱估计程序我建议第一步先用一个自己构造的理想信号验证而不是直接上实测数据。比如构造t (0:2047)/Fs; y 1.0*sin(2*pi*100*t) 0.5*sin(2*pi*300*t);用幅值谱程序画出来100Hz处峰值应接近1300Hz处接近0.5噪声底接近于0。如果峰值不对多半是幅值修正写错如果频率点偏移多半是频率轴构造错误如果谱线乱跳说明fft点数或重叠参数没写对。这个方法我每次写新代码都会做比对着任何教程都管用。6. 那么最终我给了自己一个什么工作习惯把整套程序跑通之后我养成了一个固定的三步走习惯先构造合成信号验证算法正确再对实测数据做FFT幅值谱快速摸底最后用pwelch精细估计PSD并保存结果。摸底阶段让我能快速找到关心的大致频带PSD阶段则负责给出可供报告和分析使用的平滑曲线。此外我会把谱估计参数记录在代码注释里比如% pwelch: window1024, overlap512, nfft2048。因为同样的数据用不同参数画出来的谱会有差异参数不写清楚三个月后自己翻代码都一头雾水。这个习惯帮我避免过好几次重复调参的尴尬。如果你现在正被频谱图画不准、单位搞不清的问题卡住建议先别急着抄复杂程序。用我上面最简单的FFT幅值谱跑通你的数据再逐步换到pwelch。记住先看谱峰的相对位置再看绝对幅值最后才谈平滑估计。等你把这一步走顺后面的进阶分析就水到渠成了。