做故障诊断这些年我一开始也是拿FFT硬扛。直到有一次处理现场采集的振动信号故障特征频率完全被淹没在宽频噪声里频域图上除了几个工频分量什么都看不出来才老老实实回来研究小波这一套东西。后来花了不少时间把MATLAB里的小波变换、小波包变换、信号分解重构、特征频率提取再到频谱图和包络谱图绘制整条链路跑通才觉得自己真正入了时频分析的门。这篇文章就是把这条链路完整写出来从原理讲到能直接跑的代码再讲到我踩过的坑。如果你正准备做振动信号分析、设备故障诊断或者刚接触小波变换被各种术语卡住这篇文章应该能帮你省下不少摸索时间。1. 先搞明白小波和小波包在信号处理里到底解决什么问题1.1 傅里叶变换的短板时频无法兼顾傅里叶变换的核心作用是告诉你信号里“有哪些频率”但它有一个天然缺陷丢掉了一切时间信息。一个信号在第1秒出现一个高频冲击另一段信号在末尾出现同样的冲击两者的FFT频谱几乎一样你根本分辨不出来。对很多实际信号来说频率成分是什么很重要但出现在哪个时刻、持续多久同样重要设备故障、语音、心电、电力暂态信号都有这个特点。短时傅里叶变换STFT想了个办法把信号切成小段分别做FFT。可窗口长度是固定的这就陷入一个两难——窗口短了时间分辨率高但频率分辨率差窗口长了频率分辨率好时间定位又模糊。生活里好比用固定大小的篮子去装不同大小的水果你只能选一个篮子尺寸总有水果装不合适。小波变换的出现就是为了解决这个“固定篮子”的问题。小波变换的思路是让分析用的波形可以伸缩和平移用一个长度有限、均值为零的小波基函数通过尺度参数调节它的宽窄通过平移参数扫描时间位置。频率越高就用越窄的小波去匹配频率越低就用越宽的小波去匹配于是低频段获得高频率分辨率高频段获得高时间分辨率两者不再互相牵制。这套“可变窗口”的想法在工程上非常实用。拿设备振动信号来说正常状态下的信号相对平稳用FFT看得很清楚一旦出现局部缺陷比如轴承滚道上的剥落会产生周期性的瞬态冲击这种信号既需要定位冲击发生的时刻又需要分析冲击激励起的共振频率小波变换恰好两头都能兼顾。1.2 小波变换的多分辨率思路小波变换在工程里通常分为连续小波变换CWT和离散小波变换DWT。连续小波变换的尺度是连续变化的得到的是时间-尺度平面上的二维图谱视觉上很直观但计算量大、数据冗余多。实际程序里更多用的是离散小波变换它通过滤波器组实现计算效率高并且支持完美的信号重构这也是MATLAB里wavedec等函数的基础。离散小波变换的每一次分解本质上是把信号同时通过一个低通滤波器和一个高通滤波器低通输出得到近似系数对应低频概况高通输出得到细节系数对应高频细节然后对低频部分继续做同样操作。所以它的多分辨率特性是“只对低频一层层细分”每往下一层频率分辨率翻倍但时间分辨率减半。这种结构特别适合低频占主导、高频只是叠加成分的信号。理解这一层逻辑很重要因为后面选小波、定层数、画频谱都要回到这个图景里。很多人拿到代码就开跑结果分解出来一堆系数不知道每个频带对应什么根源就是对这条“滤波器组-频带划分”的主线没有概念。1.3 场景判断什么信号该用小波包而不是小波离散小波变换有个隐含倾向它每次只细分低频部分高频部分分解一次就不再往下分。如果你的目标频率本来就在较低频段这没有问题但很多故障信号的高频部分恰恰是信息最丰富的区域比如滚动轴承早期故障冲击会激励起轴承结构的高频共振故障特征频率藏在高频带里。这时候离散小波变换的高频细节带宽太宽根本分不清里面的成分。小波包变换就解决了这个问题它不只对低频做细分高频部分同样继续分解。每一层都会把上一层的每个频带再一分为二经过n层分解得到2的n次方个等宽频带。相当于把整个分析频带切得又细又均匀你想看高频细节还是低频趋势都有一块对应的“窄带”给你。我自己的判断习惯是先想清楚目标频率大概落在哪个频段。如果是齿轮箱啮合频率、轴转频这类偏中低频的成分离散小波变换足够如果是轴承早期故障、松动冲击这类容易激起高频共振的情况直接用小波包更省心。至于更精细的细节取舍到了程序里再看实际效果调整。2. 动手前的三个关键参数采样率、小波基与分解层数2.1 采样率和分析频带先对齐写代码之前第一件事是确认信号的采样率Fs。根据奈奎斯特采样定理采样率Fs的信号可分析的频带范围是0到Fs/2Fs/2以上就是混叠区程序跑得再漂亮也没用。常见现场采集到的振动加速度信号采样率多为5120Hz、10240Hz、25600Hz对应分析带宽分别是2560Hz、5120Hz、12800Hz。这个参数直接影响后面所有计算但很多人习惯随手写一个数字或者直接从文件读了一个Fs都没核对结果画出来的频谱图频率轴是错的。我的习惯是先做一个最基本的自检生成一段已知频率的正弦信号比如100Hz采样率5120长度1秒跑一遍FFT确认峰值落在100Hz再开始分析真实数据。这个自检步骤只花一分钟能省掉后面大量排查时间。选采样率还有一个实际考虑不是越高越好。采样率太高同样时长下数据点数暴涨小波分解尤其是小波包分解的计算量会大得让你怀疑人生。后面我会专门讲大数据量怎么处理但源头上的建议是根据目标最高分析频率反推所需采样率够用就行。2.2 小波基怎么选db4与sym8的取舍小波基函数不是只有一个而是一个大家族。MATLAB里常见的有db系列、sym系列、coif系列、bior系列。选哪个直接决定分解结果的质量但没有哪个小波是“全局最优”的都是匹配场景选出来的。我做振动信号分析时最常用的是db4和sym8。db4属于Daubechies小波长度短时间定位好对瞬态冲击比较敏感计算量也小适合故障诊断里常见的突变信号。sym8是Symlets小波近似对称性更好频带的局部化能力更强波形畸变更小适合你关心频带能量分布而不是冲击位置的场景。如果用生活类比小波基相当于一把尺子尺子刻度密的适合量精细局部刻度疏的适合看整体轮廓。db4的时域支撑短跟短促冲击更“对得上”sym8的频域特性更好分析窄带成分更利索。一般我建议从一个常用组合起步比如db4三层分解先看结果如果频带划分不理想再换sym8对比不要让选型变成纠结。2.3 分解层数用公式算不靠拍脑袋分解层数n决定频带切到多细。离散小波变换第j层细节分量对应的频带是从Fs/2^(j1)到Fs/2^j近似分量对应0到Fs/2^j。所以你想让目标频率落在某一层就可以反推层数。举个例子采样率Fs5120Hz目标特征频率是100Hz左右。如果选n5第5层细节频带是80到160Hz目标正好落在里面如果选n4第4层细节频带是160到320Hz100Hz就落到近似分量里去了。所以在这个例子里n5是合理起点。目标频率越低需要的层数越多。但层数不是越多越好。层数每加一计算量翻倍边界效应也更明显重构出来的信号两端更容易出现异常摆动。我一般建议n取3到5除非目标频率非常低否则很少超过7层。小波包分解同理频带宽度是Fs/2^(n1)用同一套逻辑推算即可。3. MATLAB分解与重构实操从wavedec到waverec3.1 三步完成离散小波分解MATLAB做离散小波分解很简单核心函数就一个wavedec。先假设你有一段信号x采样率Fs5120Hz代码如下Fs 5120; n 5; % 分解层数 wname db4; % 小波基 [C, L] wavedec(x, n, wname);C是分解后所有系数拼接成的一个长向量L是每一层系数长度的记录数组。理解C和L的组织方式很重要因为后面wrcoef、waverec都要靠它知道哪些系数属于哪个频带。C的排列顺序是第1层细节系数、第2层细节系数……第n层细节系数最后是第n层近似系数。L则依次记录这些段的长度。这个结构看起来枯燥但踩过一次坑就知道它的价值。我第一次用的时候直接对C做处理结果出来的信号长度和频带全乱了后来才明白必须通过L来定位每一段系数的范围。如果你只是想要“分解-重构”这条链路可以不用手动拆C直接用wrcoef按层提取。3.2 wrcoef与waverec的区别重构谁、怎么重构分解之后最常见的需求是根据某一层系数重建出对应频带的时域信号看看这一频带里到底长什么样。wrcoef就是干这个的a5 wrcoef(a, C, L, wname, 5); % 近似分量对应低频 d1 wrcoef(d, C, L, wname, 1); % 第1层细节对应最高频段 d5 wrcoef(d, C, L, wname, 5); % 第5层细节越低层对应频率越低注意wrcoef返回的长度和原始信号x完全一致这是因为它内部做了系数补零和重构滤波。所以你可以把这几个分量加起来检验完整性。xiaorong的经验是把a5和d1到d5全部相加结果应该和x基本一样误差通常在10的负12次方量级这叫做完美重构。如果误差很大优先检查小波基是否一致或者分解层数是否改动过。waverec则不一样它不让你选某一层而是直接把全部系数重新合成为原始信号调用方式是x_recon waverec(C, L, wname);它的用处主要是验证整个分解链路有没有问题或者配合阈值处理后的系数实现滤波去噪。理解两者的区别就能少犯糊涂waverec是“原样还原”wrcoef是“挑某个频带单独看”。3.3 结合软阈值实现信号去噪分解重构最常见的工程应用就是去噪。小波去噪的基本逻辑是噪声能量分散在所有小波系数上而有用信号的系数幅值通常较大所以把绝对值较小的系数压掉再重构就能滤掉一部分噪声。MATLAB里最简单的方式是直接调用wdenoisex_den wdenoise(x, n, Wavelet, wname, DenoisingMethod, SoftThreshold);如果想更精细地控制阈值可以自己选层、选阈值规则thr thselect(d1, heursure); % 基于第1层细节估计阈值 d1_den wthresh(d1, s, thr); % s为软阈值h为硬阈值软阈值会把系数向零收缩结果更平滑硬阈值保留幅值但截断容易留下毛刺。做故障诊断的时候要特别注意阈值设太狠有用的冲击特征会被一起削掉包络谱上看不出故障频率。我自己的建议是去噪时保守一点阈值宁低勿高或者干脆在选定频带内做带通滤波不轻易用全局阈值。4. 小波包变换把高频细节再切细4.1 小波包如何细化频带小波包变换在MATLAB里的核心是wpdec。它的思路和离散小波的差别一句话就能说清离散小波只把低频继续分小波包把低频和高频都继续分。经过n层小波包分解频带被均匀切成2的n次方份每份带宽是Fs/2^(n1)。这意味着如果你关心的高频信号在某个窄带里小波包能把这块窄带单独拎出来重构而离散小波做不到这种粒度。拿Fs5120Hz、n3举例小波包把0到2560Hz分成8个带宽为320Hz的频带每个节点对应其中一段。你就能针对性地分析某个320Hz带宽内的信号做包络谱时干扰成分会少很多。这里有一个工程直觉要建立小波包的“均匀切分”不是免费的。频带越细每个频带内的数据点越少时间分辨率越差而且节点越多处理每个节点的计算量线性上升。所以小波包层数不是越多越好够用就行通常n3到5。4.2 wpdec与wprcoef的完整调用示例小波包分解的代码非常简洁wp wpdec(x, 3, db4); % 3层小波包分解wp是小波包树对象分解完成后用wprcoef提取某个节点的重构信号。节点用两个数表示第一是层数第二是该层内的节点索引rec wprcoef(wp, [3 1]); % 第3层第1个节点的重构信号如果想一次性取出第3层所有8个节点写个循环rec_cell cell(1, 8); for i 0:7 rec_cell{i1} wprcoef(wp, [3 i]); end每个节点重构信号的时长和原始信号一样但是对应的频带不同。拿到这些重构信号后最常做的事情是计算每个频带的能量作为特征向量用于后续分类或识别energy zeros(1, 8); for i 0:7 node_sig rec_cell{i1}; energy(i1) sum(node_sig.^2); end energy energy / sum(energy); % 归一化能量特征这套能量特征在故障诊断和机器学习特征工程里非常常用。小波包能量特征的本质是“信号的频率成分分布”不同状态下能量在不同频带的分布会发生明显变化用这个特征向量训练分类器效果往往比直接用原始波形好。4.3 小波包节点频率排序的坑这里必须提醒一个非常容易翻车的地方小波包树的节点索引不是按频率从低到高排的。第3层的节点索引0到7看起来从左到右应该频率递增但实际频率顺序是乱的这跟滤波器组的频带交错有关。MATLAB的常见节点索引对应的频率顺序在不同小波下会有差异。最稳妥的办法是用wpfrqord函数做频带排序或者绘制节点的频谱图确认。我自己的习惯是提取某个节点信号后顺手对它做一次FFT看它的中心频率落在哪里再决定这个节点是不是我要找的频带。实践中第3层的频率排列常出现类似0、1、3、2、7、6、4、5的非单调顺序所以你只靠索引猜频带十有八九会看错。这个坑如果没踩过很容易出现“明明提取了某个节点包络谱上却看不到想要的故障频率”的困惑然后怀疑算法写错了。其实只是频带没对上。遇到这种问题先确认节点实际频段再继续往下走。5. 频谱图与包络谱图提取特征频率的关键一步5.1 一张规范的FFT频谱图怎么画画频谱图看起来人人都会但频率轴算错的人非常多。我先给一个标准的单边频谱绘制代码N length(x); Xf fft(x); f (0:N-1) * Fs / N; % 频率轴单位Hz Amp 2 * abs(Xf) / N; % 单边幅值 Amp(1) Amp(1) / 2; % 直流分量不乘2 plot(f(1:floor(N/2)), Amp(1:floor(N/2))); xlabel(频率 / Hz); ylabel(幅值);这里面有几个关键细节。第一FFT结果是对称的只看前一半即可第二幅值要乘以2再除以N否则纵轴数值完全不对第三直流分量也就是第1个点不需要乘2。频率分辨率是Fs/N比如Fs5120HzN5120时频率分辨率刚好1Hz能区分相差1Hz的两个峰如果分辨率不足两个峰叠在一起根本分不开。画频谱图之前建议先做去均值处理否则直流分量会特别大纵轴一压缩其他频率成分全看不见。常用做法是x x - mean(x)简单有效。还可以加窗减小频谱泄漏但加窗后要做幅值校正细节处理比较复杂初学时先用不加窗或矩形窗跑通别一上来就堆一堆窗函数。5.2 包络谱分析Hilbert解调流程频谱图很多时候不够用。比如滚动轴承的故障特征本质上是周期性冲击对高频共振信号的调制直接看原始信号的FFT高频段是一团宽峰故障特征频率藏在调制边频里肉眼根本分不清。这时候就要做包络谱分析也就是Hilbert解调。标准的MATLAB流程分成三步。第一步选择目标频带通常用小波包重构出故障冲击所在的高频带信号第二步做Hilbert变换求包络第三步对包络信号做FFT得到包络谱。代码是这样x_band wprcoef(wp, [3 2]); % 选一个包含故障冲击的频带 env abs(hilbert(x_band)); % 希尔伯特变换取包络 env_spec fft(env); N_env length(env); f_env (0:N_env-1) * Fs / N_env; env_amp 2 * abs(env_spec) / N_env; env_amp(1) env_amp(1) / 2; plot(f_env(1:floor(N_env/2)), env_amp(1:floor(N_env/2))); xlabel(频率 / Hz); ylabel(包络幅值);包络的本质是提取信号的瞬时幅值变化。调制信号经过Hilbert变换解调后高频载波被剥掉剩下的包络信号里就只剩下调制频率成分故障特征频率在包络谱上会以清晰谱峰的形式出现。用生活类比相当于把一个调频广播信号先解调成音频你才能听到里面的话否则听到的只是载波的嗡嗡声。5.3 如何从包络谱中确认特征频率拿到包络谱之后下一步是找出谱峰对应的频率并跟理论特征频率对照。可以用findpeaks自动挑峰[pks, locs] findpeaks(env_amp, MinPeakHeight, 0.01, MinPeakDistance, 3); peak_freqs (locs - 1) * Fs / N_env;MinPeakHeight要按实际幅值水平调整MinPeakDistance至少是几个频率分辨率避免把同一个峰的旁瓣当成多个峰。得到峰频率列表后把前几个最大的峰值频率和处理论值对比允许的偏差一般不超过转频的1%到2%或者不超过一个频率分辨率。如果是滚动轴承外圈故障频率、内圈故障频率等都有理论公式需要预先根据轴承参数和转速算出来。这里特别要注意特征频率不是只有一个基频往往还有2倍、3倍频确认时优先看基频再看谐波是否成系列。另外转频、轴的转频及其整数倍也是常见干扰别把转频当成故障频率。我踩过的坑是看到一个大谱峰就激动地当作故障频率结果一算只是转频白高兴一场。所以最后对照理论值时一定把所有候选峰值都列出来跟理论值表一起比对差距大就要回头检查前面的频带选择是否合理。6. 常见问题排查与技术要点速查6.1 重构信号端点异常与边界效应用wrcoef重构某个频带后经常发现信号两端出现明显震荡幅值比中间大很多这是小波变换的边界效应。根本原因是滤波器在处理信号边界时缺少足够的数据只能靠延拓补齐延拓方式带来的误差会在端点集中体现。处理办法有几个一是分析时舍去两端各一小段数据比如各丢弃小波滤波器长度一半的样本点二是换用处理边界能力更强的小波比如sym系列三是在分解之前对信号做对称延拓分析完再截掉延拓部分。我在实际项目中用前两种比较多因为简单直接。如果你做的是实时在线分析端点的短暂异常不会影响整体特征频率判断不必过度纠结。6.2 频率轴标定错误频谱图和包络谱上峰值位置与理论预期偏差大第一反应先查频率轴相关的几个变量采样率是不是读对了、FFT点数是不是用了段长而不是整个信号长度、幅值变换是不是忘了除以N。最让人头疼的是那种“看起来一切正常但频率差了一倍”的情况多数是单边谱处理没做对或者采样率填成了原来的两倍。排查这类问题我有一个固定套路生成一个5Hz、50Hz、500Hz三个正弦叠加的测试信号采样率比如5120Hz跑完整套流程看看频谱图三个峰是不是落在正确位置。测试信号能通过说明程序链路没问题问题出在真实信号的参数设置上测试信号不过说明算法流程本身有bug直接修流程。6.3 计算资源与大数据处理小波包变换的计算量跟节点数成正比层数到5层就是32个节点6层就是64个节点。如果信号长度达到百万点级别全节点重构速度会很慢内存也容易爆。程序卡死、风扇狂转、MATLAB直接报内存不足这些情况我都遇到过。应对措施按优先级排序第一控制层数不要盲目追求细分第二只提取需要的几个节点别全节点都重构第三对数据做降采样前提是降采样之后目标频率仍在分析带宽内第四分块处理每块单独做小波包分析和特征提取再合并结果。最后一条适用于很多工业场景因为故障特征本来就具有周期性截取几段足够长的数据分别分析即可。6.4 小波基与层数的试错方法调试时如果发现特征频率不明显我建议不要在一个配置上反复调整参数而是系统地做一轮比较。先把层数从3试到5小波基在db4、db8、sym8之间切换输出每个配置下的包络谱对比哪个谱峰最清晰。这个对比过程建议写一个小脚本自动化把参数和结果组织成表格或者并排图直观程度远好于改一行跑一次。判断“哪个结果更好”的标准我一般看三个维度目标频率处的谱峰幅值是否足够突出、干扰谱峰是否更少、重构信号是否还保留明显的冲击特征。这个取舍没有绝对的对错但有了标准就不会七改八改找不到方向。另外每次切换参数后务必重新确认分解层数和目标频带是否匹配两者是联动的只改小波基不改层数也未必达到最好效果。最后再分享一个实用的封装思路。我会把“加载信号-设置参数-小波包分解-提取特征频带-包络谱分析”这一整套流程写成一个函数输入只有信号、采样率和目标频率范围输出直接是特征频率列表和谱图。这样换一组数据时只需要改一行参数不用每次从头复制粘贴脚本也不容易漏改某个关键设置。这算是整个调试过程中让我最受益的一步比单纯堆代码函数更省心也让整套方法更容易复用到新的数据上。
