很多做故障诊断、状态监测的朋友应该都有过这种经历设备明明有异常你抓了一段振动信号直接做FFT肉眼可见就是一堆“糊在一起”的谱线低频和高频成分混在一起根本分不清哪个是哪个。尤其是滚动轴承早期故障、齿轮点蚀这类信号故障特征往往藏在某个高频共振频段的调制信号里普通的傅里叶频谱压根看不出来。这时候就该小波变换和小波包变换上场了。它俩不是替代FFT的而是补FFT的短板的——FFT只能告诉你“这个信号里有什么频率”小波还能告诉你“这个频率成分发生在什么时刻、落在哪个频段”。再加上包络谱这一套组合拳就能把藏在调制信号里的故障特征频率扒出来。这一篇我就把自己常用的MATLAB小波分析流程完整拆开讲一遍。不是教科书式地讲一堆数学公式再给个demo而是按我实际做项目的顺序怎么做分解、怎么做重构、怎么选小波基、怎么算包络谱每一步都直接给出能跑的代码。你把这个流程摸透了再遇到“分解与重构信号、提取特征频率、画频谱图和包络谱图”这类需求基本就是套模板的事。1. 从“FFT看不清故障频率”开始小波变换到底解决什么问题先聊点实在的。FFT在稳态、平稳信号面前表现确实不错但工程上的实测信号几乎都是非平稳的启动过程中的变转速振动、冲击引起的瞬态响应、早期故障产生的周期性微弱冲击这些信号要么频率随时间变化要么能量又小又短。你用FFT算出来的结果只是这些信息在时间轴上的一个“平均值”时间位置完全丢失了。调幅信号是个特别典型的例子。轴承外圈出现局部缺陷后每个滚动体滚过缺陷位置时会产生一个冲击这个冲击激励起轴承系统的高频固有振动所以传感器上测到的信号不是一个单纯的高频正弦波而是一个高频分量被一个低频分量调制后的“调幅波”。这个低频分量就是故障特征频率我们真正关心的就是它。直接把这种信号做FFT你会在高频段看到一大圈“鼓包”但这个鼓包的细节结构被噪声和邻近频率成分淹没低频调制频率也往往不明显。你需要先把目标频段从混合信号里抠出来再做一个解调处理也就是求包络才能看清故障频率。小波变换在这里的核心贡献有两个多分辨率分析低频段用窄频带、长窗口去看高频段用宽频带、短窗口去看。这种“变焦距”能力能在保留时间信息的同时把信号一层一层拆开变成不同频段的分量。按频带重构分解出来的每一层分量都对应原始信号的某个频段你可以把不关心频段丢到一边把关心频段重组出来再对这个重构信号做包络谱分析。小波包变换则更进一步。小波变换离散小波的分解思路是高精度往下拆低频高频部分一拆完就扔在一边不继续处理了。如果你关心的特征频率恰好落在高频段小波变换就无能为力。小波包则对高频段也做同样的二分拆解最后得到的是一个完整的频带二叉树低频和高频都有精细分辨率。所以这两兄弟的适用场景非常清晰信号特征在低频段用小波变换就够特征频率落在中高频段或者你想用多频带能量做特征向量做模式识别那就上小波包。2. 小波基和分解层数怎么选这一步决定了结果能不能看标题里既然点名了小波变换和小波包变换的MATLAB程序这一步我就展开讲讲。选错了小波基或者分解层数后面画出来的图就一个字——糊。2.1 小波基的“匹配”逻辑小波基函数要跟信号的形状越像越好。MATLAB里经常用到的有这么几类小波族代表特点适用场景Daubechiesdb3、db8正交、紧支撑与冲击成分相似度较高轴承/齿轮冲击类信号工程实用最多Symletssym4、sym8近似对称相位畸变小需要尽量保留波形形态时Coifletscoif2、coif4对称性更好但支撑长度更长信号比较平滑时Biorthogonalbior3.5、bior4.4对称且可精确重构图像类、需要线性相位时我自己在机械故障特征提取里用得最多的还是db系列和sym系列尤其希望波形保真度高一点时会选sym。正交性对小波重构来说很重要混淆和冗余会少一些。如果你只是做一个特征频率的包络谱分析不是特别在意相位db6到db9都挺稳。有个误区是“小波阶数越高越好”。实际并不是。db阶数越高滤波器长度越长边界效应越明显计算量也越大。我一般先在db4到db8之间试观察重构信号的时域波形和频谱是否干净再定具体阶数。2.2 分解层数怎么算不是随便拍脑袋分解层数直接决定了频带宽度。比如采样频率是fsN层分解后最高分解层近似分量的频率范围0 ~ fs/(2^(N1))第1层细节分量的频率范围fs/4 ~ fs/2第2层细节分量的频率范围fs/8 ~ fs/4如果你已经知道目标特征频率落在哪个范围就用这个表反推需要几层。公式化一点你想把信号拆到最低频段的带宽是Δf那么分解层数N需要满足 fs/(2^(N1)) ≤ Δf。举个例子采样频率fs10000 Hz你的故障特征频率在100 Hz左右高频共振频带在2500 Hz附近。这时候N4比较合适。为什么因为4层分解后各个频带范围是这样的分量频带范围Hz备注A40 ~ 312.5轴转频、低频故障频率都在这里D4312.5 ~ 625D3625 ~ 1250D21250 ~ 2500D12500 ~ 50002500 Hz共振频带落在这里这样目标频段就被自然隔离开了。分解层数太少低频部分还是混频层数太多某一层信号长度变短边界效应影响增大重构质量下降。经验上除非分析超低频趋势项一般不建议超过8层。3. MATLAB小波分解与重构实战从wavedec到指定频带信号重建这部分直接上代码。我用MATLAB做的流程基本固定为构造/采集信号 → wavedec做多层小波分解 → 用appcoef/detcoef提取各层系数 → 用wrcoef把指定系数重构回时域信号。下面用一个仿真信号演示。3.1 构造一个带噪的调制测试信号clear; clc; close all; % ---------- 生成一个模拟振动信号 ---------- fs 10000; % 采样频率 10 kHz t (0:fs*2-1)/fs; % 2 秒信号 f_rot 25; % 转频 25 Hz f_fault 91.5; % 故障特征频率例如外圈BPFO f_resonance 2500; % 系统高频共振频率 % 低频轴频成分 signal 1.0 * cos(2*pi*f_rot*t); % 高频调制成分故障冲击激励起共振 modulation 1 0.85 * cos(2*pi*f_fault*t); signal signal 0.7 * modulation .* sin(2*pi*f_resonance*t); % 加一点随机噪声 signal signal 0.3 * randn(size(t));这里特意把故障特征频率设置为91.5 Hz不取整数就是为了避免跟频谱栅栏对齐模拟真实情况中频率分辨率不足的问题。3.2 4层小波分解与各频带重构wavename db8; N 4; [C, L] wavedec(signal, N, wavename); % 提取最后一层近似系数 cA4 appcoef(C, L, wavename, N); % 提取1-4层细节系数 cD1 detcoef(C, L, 1); cD2 detcoef(C, L, 2); cD3 detcoef(C, L, 3); cD4 detcoef(C, L, 4); % 重构为时域信号 A4 wrcoef(a, C, L, wavename, N); D1 wrcoef(d, C, L, wavename, 1); D2 wrcoef(d, C, L, wavename, 2); D3 wrcoef(d, C, L, wavename, 3); D4 wrcoef(d, C, L, wavename, 4);wrcoef是“重构指定节点”的函数它内部会自动把其他系数补零后再逆变换。所以我得到的D1、D2、D3、D4和A4都是跟原始信号等长度的时域序列。3.3 分析重构分量判断特征频率落在哪一层分别对各个重构分量做FFT看能量集中位置。这一步的核心逻辑是哪一层的频带范围覆盖了你关注的特征频率就在那一层上继续做包络谱分析。% 画各层频谱 figure; for k 1:4 subplot(4,1,k); [Pxx, fvec] myFFT(eval([D num2str(k)]), fs); plot(fvec, Pxx); xlim([0 5000]); title([细节分量 D num2str(k)]); end上面的myFFT是我习惯写的一个局部函数封装了单边幅值谱计算避免重复写。它的核心逻辑见下一节。在D1层2500-5000 Hz会看到2500 Hz共振频率处有一个明显的谱峰说明高频共振成分确实被这一层“抓住”了。3.4 重构任意频带信号后的频谱观察如果把注意力放在D2层和D3层可能发现几乎没有有效峰值说明目标频段不在这两层。这个现象本身就是有效信息它告诉你要重点分析的频段是D1层。提示做小波分解前先对信号做一次带宽分析或直接预览FFT知道能量集中带大概在什么范围能大大减少试错时间。我一般先用pspectrum或FFT快速看一眼再定分解层数。4. 小波包变换把高频段也精细切开特征频率一个都跑不掉小波变换对高频段的处理比较“粗糙”——每层只对低频部分继续分解高频部分一旦拆出来就不再动了。小波包变换对高频部分同样做二分所以你会得到2^N个等宽频带N为层数。这意味着不管你的特征频率落在低频还是高频都可以用一个相对窄的频带去覆盖它。4.1 wpdec和wpfrqord的配合使用MATLAB里做小波包最常用的是wpdec但有个关键坑默认的节点编号顺序不是按频率高低排列的而是按“小波包树”的遍历顺序。直接按节点顺序画频谱会看到频率乱跳这也是很多新手一上来就懵的地方。解决方法是加一行wpfrqord% 3层小波包分解 wpt wpdec(signal, 3, db8); % 按频率顺序重新排列节点 wpt wpfrqord(wpt); nodes leaves(wpt); % 得到按频率顺序排列的叶子节点这一步做完nodes里的顺序才对应从低频到高频的8个频带。每个频带宽度为fs/(2^(N1))代入fs10000、N3就是625 Hz带宽8个频带覆盖0-5000 Hz。4.2 从叶子节点重构窄带信号% 假设我们想重构第5个节点按频率顺序2500-3125 Hz target_idx 5; sig_recon wprcoef(wpt, nodes(target_idx)); % 如果想直接看各节点系数不重构时域信号用 wpcoef coeffs wpcoef(wpt, nodes(target_idx));注意wprcoef和wpcoef的区别wprcoef重构回时域信号长度跟原始信号一致wpcoef给出的是小波包系数序列长度随层数不同会变化。做特征频率提取时大部分情况用wprcoef就够了。4.3 用小波包各频带能量做特征向量小波包除了用来提取特定频带的信号还有一个很常见的用途构造特征向量。做法是对所有叶子节点分别重构信号算每个频带的能量或能量占比拼成特征向量。energies zeros(length(nodes), 1); for k 1:length(nodes) sig_k wprcoef(wpt, nodes(k)); energies(k) sum(sig_k.^2); end energy_ratio energies / sum(energies); % 直接绘制能量占比条形图观察故障特征频带 bar(energy_ratio); set(gca, XTickLabel, arrayfun((x) sprintf(%.0fHz, x), ... (0:length(nodes)-1)*(fs/(2^(31))) fs/(2^(4)), UniformOutput, false));这个条形图在故障诊断里很直观某几个频带的能量明显偏高说明信号能量主要集中在那几个频带。如果你发现2500 Hz附近的频带能量占比异常大那就锁定了共振频带后面在该频带内做包络谱即可。提示能量特征向量很适合作为机器学习或神经网络的输入。你不需要自己设计复杂的滤波器和窗函数小波包分解天然就给了你一套规范化的频带特征。5. 包络谱是怎么做出来的希尔伯特变换提取调制特征的关键操作包络谱这个术语说穿了就是对“包络信号”再做一次FFT。前面已经提到早期故障信号往往表现为“高频共振 低频调幅”。直接FFT只能看到高频载波的大鼓包看不到低频调制成分。把信号经过希尔伯特变换取绝对值后得到的就是这条包络线包络线的频谱中才会露出真正的故障特征频率。5.1 MATLAB里三行搞定包络解调analytic hilbert(sig_recon); % 希尔伯特变换得到解析信号 env abs(analytic); % 取幅值得到包络 env env - mean(env); % 去直流包络中直流分量往往太大很多人不理解为什么每步都要写。解释一下hilbert并不是“算一个希尔伯特谱”它返回一个复数序列虚部是原始信号的希尔伯特变换实部是原始信号本身。对这个复数的模abs()取到的就是信号的瞬时幅值。这一步在数学上等价于AM解调很经典。5.2 包络谱的FFT计算与频率轴对应有了env之后再对它做FFT。这里有一个非常容易出错的点FFT之后频率轴怎么画对、幅值怎么修正。我写成函数方便复用function [P1, f] myFFT(sig, fs) N length(sig); Y fft(sig); P2 abs(Y / N); P1 P2(1:floor(N/2)1); P1(2:end-1) 2 * P1(2:end-1); f fs * (0:floor(N/2)) / N; end这个函数处理了三个关键点对双边谱截取前半部分、对非直流除直流外的频率乘2做幅值修正、构造正确的频率刻度。把上一节重构出来的D1层信号2500-5000 Hz输入这个函数得到包络谱后你会在91.5 Hz处看到明显的谱峰它才是真正要找的故障特征频率。而直接对原始信号做FFT91.5 Hz处通常被低频大能量成分压住或者被噪声淹没得根本看不清。5.3 包络谱里看什么包络谱的X轴通常只画0到几百Hz就够。为什么因为调制频率故障特征频率一般远低于载波频率大多是几Hz到几百Hz。你如果还把包络谱画到几千Hz图上一堆噪声峰反而重点不突出。实用经验是X轴限制到调制频率高值的1.5~2倍就够了比如关注100 Hz的特征频率X轴画到200 Hz。包络谱的谱峰位置对应调制频率也就是轴承故障特征频率BPFO、BPFI、BSF等。实际中还要留意谐波成分如果在f_fault的2倍频、3倍频处也有谱峰可靠性就更高了。6. 一个仿真信号走完全流程三种分析结果逐层对照代码写了这么多关键还是要看它们组合起来的效果。为了让你直观感受到“小波/小波包 包络谱”比“直接FFT”强在哪我用3.1节那个仿真信号走一遍完整流程把每一步的结果拿来做对比。6.1 直接FFT看到了大致轮廓但特征频率不突出直接对signal做FFT画出0-300 Hz范围的频谱25 Hz处有明显谱峰这是转频幅值大很容易看出来。91.5 Hz处有谱峰但被噪声污染旁边一堆干扰峰甚至可能淹没在峰值中。2500 Hz附近有宽泛的能量看不出调制信息。这就是最常见的痛点低频区的故障特征频率会被转频、齿轮啮合频率及其谐波压制高频区只能看到“一坨能量”却不知道这坨能量是被哪个低频频率调制的。6.2 小波分解 重构 包络谱锁定高频频段后特征频率浮现流程如下wavedec做4层分解提取并重构D1层2500-5000 Hz信号能量集中在共振频率附近对D1层重构信号做希尔伯特解调得到包络对包络做myFFTX轴0-200 Hz。结果91.5 Hz的位置出现一个非常干净的谱峰旁边几乎没有什么干扰。25 Hz转频成分因为已经被高频段隔离根本不会进来捣乱。6.3 小波包分解 重构 包络谱频带定位更灵活小波包同样能做这件事3层分解按频率顺序排列找到2500-3125 Hz节点第5个节点重构后用同一套包络谱流程。效果与小波方法基本一致但适用性更广。假如你的共振频带不在D1层覆盖的2500-5000 Hz而是一半落在两个小波细节层的边界上小波包的精细频带可以更好地贴近真实频带边界。6.4 三种方法的对照分析方法能否看到91.5Hz特征边界效应适用频率范围复杂度直接FFT困难易被淹没无全部频段最低小波分解包络谱清晰轻微低频为主中小波包分解包络谱清晰轻微低频高频均可精细定位中一个小细节提醒小波包分解后如果共振频率恰好落在某个节点的边界上能量会被分配到相邻两个节点导致包络谱幅值略低于真实值。解决方法是增加分解层数让频带更窄边界更贴近真实频率。但也要权衡层数增加后短数据段的频带统计稳定性下降所以层数不是越多越好。7. 复盘我踩过的坑和优化技巧文章最后一部分梳理我实际使用小波分析时踩过的几个坑不少坑是文档里不会写的。7.1 边界效应信号两端“翘起来”的假象小波变换本质上是卷积滤波边界处的数据不足会导致重构信号两端出现幅度异常。尤其是数据量比较短的时候这个问题非常明显。解决办法是要么把数据采集长一点丢掉前10%和后10%再分析要么在分解前用wextend做边界扩展例如对信号做‘sym’对称延拓分析完后再截掉延拓部分。signal_ext wextend(1, sym, signal, round(length(signal)*0.2)); % 对signal_ext做分解重构... % 最后截掉前后扩展的部分 sig_result sig_result(ext_len1:end-ext_len);7.2 小波包节点顺序的坑别问我为什么知道这个词要单独拿出来说——wpdec画出来的节点默认顺序跟频率顺序不一致很多人第一次用小波包所有节点频带对不上号最后拿到的特征频率是错的还在那纠结半天。一定要记得wpfrqord这行。如果你用Python的PyWavelets顺序问题也一样存在处理思路是查频带排列再做索引映射。7.3 分解层数和你想要的频带宽度要匹配分解层数不够低频分量混着高频分量层数太多每个频带里的有效点数太少统计特征不稳定。一个实用的判断方法先做一层分解看D1层频谱里最高峰频率是否接近你预期的高频共振频率。如果已经能对上再在该层做包络谱就够了。如果对不上再加一层。7.4 包络谱X轴限制要合适包络谱可以算但画图时X轴范围一定要限制到调制频率附近别一上来全画。很多朋友包络谱画出来一大片全是高频噪声峰反而把特征频率藏掉了。我通常的做法是先按经验值限制在0~200 Hz或0~fs/40有峰值后再缩小范围精确定位。7.5 阈值去噪的正确用法别把小波重构当成万能滤波器小波阈值去噪确实能压噪声但阈值定得太狠会把真实冲击成分也削掉包络谱看起来“干净”了特征频率幅值却也掉了。我的建议是在提取特征频率的场景下先做窄带重构去过掉无关频带再做包络谱除非噪声确实大得离谱否则不要轻易加非线性阈值处理。阈值处理更适合那些需要保留时域波形做进一步分析的场景。7.6 顺手加一个自动化特征频率提取模板很多时候你不需要每次都人肉看图。我常用的一个自动化流程是先小波包分解计算各个叶子节点频带能量锁定能量最大的前两个频带把这两个频带重构后的信号叠加起来做包络谱在包络谱里找出前三个幅值最大的谱峰输出它们对应的频率。这套逻辑可以直接打包成一个函数输入原始信号和fs输出候选特征频率。做批量数据分析时能省很多事。function candidateFreqs autoFaultFreq(signal, fs) wpt wpdec(signal, 3, db8); wpt wpfrqord(wpt); nodes leaves(wpt); energies zeros(size(nodes)); for k 1:length(nodes) sig_k wprcoef(wpt, nodes(k)); energies(k) sum(sig_k.^2); end [~, idx] sort(energies, descend); recon wprcoef(wpt, nodes(idx(1))) wprcoef(wpt, nodes(idx(2))); env abs(hilbert(recon)); env env - mean(env); [P1, f] myFFT(env, fs); [~, pks] findpeaks(P1(1:round(length(f)/10)), SortStr, descend, NPeaks, 3); candidateFreqs f(pks); end这个函数不算长但把“小波包分解 → 频带能量选择 → 重构 → 包络谱 → 候选峰值提取”整条链路串起来了。做批量数据分析时直接调用再去排查候选频率对应的是不是真实故障特征频率。小波变换、小波包变换加包络谱这一套组合不是那种“高大上但落不了地”的学术玩具。它解决的就是工程里最实际问题信号混在一起分不开、噪声太大、特征频率被淹没。拿这个流程去分析你自己手头的数据先看频谱确定大致频带区域再用小波或小波包把目标频带抠出来最后包络谱一锤定音。多跑几组数据你就知道这套东西比肉眼盯FFT好用太多。