MATLAB小波分析实战：从去噪到故障特征提取

简介《MATLAB小波分析超级学习手册》的配套程序代码面向需要借助MATLAB开展信号处理、图像处理、噪声过滤与数据压缩等工作的科研人员和工程师。代码覆盖小波基选择与可视化、一维信号分解重构、二维/三维小波变换、小波包分析、连续小波变换及阈值去噪等常见场景并提供可直接运行的m脚本与示例数据。资源共119个文件以113个m脚本为主辅以4个wav语音样本和2个fig图形文件压缩包仅504KB便于快速下载和对照学习。目前已有499人学习/下载。读者可借助这些代码梳理小波分析从原理到落地的完整流程如通过wavedec/waverec理解多尺度分解利用阈值处理实现信号去噪结合cwt观察时频分布进而迁移到故障诊断、图像压缩等实际任务。1. 小波分析这门手艺为什么 MATLAB 里它比傅里叶好用做信号处理的人迟早会遇到一个让傅里叶失效的场景振动信号里混着一段瞬时冲击生理信号里叠着高频干扰裂缝信号刚到一半就断了。傅里叶变换能告诉你信号里有哪些频率却说不清这些频率出现在哪一刻因为它的基函数是全局的正弦波。小波分析不一样它用一簇可以伸缩和平移的波去匹配信号高频时窗口自动收窄、低频时窗口自适应放宽既能看频率又能定位时间。MATLAB 自带的小波工具箱是这个方向投入产出比最高的入口函数都已经封装好难点不在理论公式而在命令怎么拼、阈值怎么取、分解层数怎么定。这篇文章就按这条动手路线把核心概念、可直接复现的去噪与特征提取代码、以及新手必踩的坑一次讲透。2. 破开黑匣子小波变换的核心概念与 MATLAB 工具箱选型2.1 傅里叶的短板与小波的设计逻辑时频分辨率为什么不能兼得傅里叶变换把信号拆成无限长的正弦基函数得到的是整个时间段内的平均频谱。对平稳信号这个描述是完整的但真实工程信号几乎没有平稳的——轴承磨损初期只在某个转速区间激起高频共振肌电信号只在肌肉发力瞬间出现特定节律雷达回波里目标一出现频率就跳变。这些场景下傅里叶频谱只能告诉你“存在这些频率”不能告诉你“这个频率是何时出现的”也就是失去了时间定位能力。短时傅里叶变换试图补上这个缺陷做法是在信号上滑一个固定长度的窗口对每个窗口内的片段做傅里叶变换。问题也随之而来窗口短时间分辨率高但频率分辨率差窗口长频率分得清但时间定位变模糊。这个矛盾是固定窗函数的固有属性调整窗口长度只是在两头之间换位置并不能同时兼顾。小波变换解决这个问题的思路是换掉基函数不再用无限伸展的正弦波而是用一个有限长度、会衰减的母小波对它做缩放得到不同频率的子波做平移得到不同时刻的子波。缩放因子越大子波在频域越窄、时域越宽缩放因子越小子波在时域越窄、频域越宽。于是高频细节用窄波看低频趋势用宽波看时频分辨率随频率自动调整这就是多分辨率分析的基本逻辑。MATLAB 里对应着两套计算体系。连续小波变换cwt输出一张时频图适合肉眼观察信号里某一频段的成分随时间的变化也适合做时频峰值提取离散小波变换dwt及其多尺度版本wavedec则是把信号按 2 的幂次逐层分解成近似分量和细节分量适合后续做定量计算比如去噪、压缩、特征提取。两者不冲突先用 cwt 看全局心里有数再用 wavedec 落定量特征。这里有个常见误解以为小波变换需要自己实现滤波器组。实际 MATLAB 的 wavedec 底层就是 Mallat 算法信号每分解一层通过一对低通和高通滤波器后做隔点抽样得到长度减半的近似系数和细节系数。这个“每层都减半”的性质直接决定了分解层数的上限层数不能超过 log2(信号长度)否则系数长度会变成 0程序直接报错。日常选型可以按这个标准判断信号里关心的是瞬态突变、冲击、尖峰优先用短支撑、低消失矩的小波比如 db2、haar它们的时域支撑短能紧贴突变点信号是平滑曲线加噪声优先用消失矩高一些的小波比如 sym8、coif4它们能用更少的系数表达光滑部分信号带有对称特征或需要保证重构后不产生相位畸变可以关注 bior 双正交族。理解这一层后面调参时才不会盲目试。2.2 MATLAB 小波工具箱函数地图六个高频函数与选型建议刚接触这个工具箱的人最容易困惑的不是理论而是面对几十个函数不知道从哪个开始。按我平时带项目的经验先掌握六个函数就够覆盖八成以上的工作。函数作用典型使用场景wavedec多层一维离散小波分解输出系数数组 C 和长度记录 L拿到多尺度系数做去噪、特征提取waverec对分解系数做重构还原出信号和 wavedec 配对使用验证处理后信号wthresh对系数做软/硬阈值处理手动阈值去噪的核心步骤wdenoise一站式去噪分解、估算噪声、阈值化、重构一步完成快速验证去噪效果不深究内部细节wpdec小波包分解返回小波包树对象需要对细节分量继续细分时使用cwt连续小波变换返回时频系数观察信号时频分布定位异常时刻做选型时有一个习惯值得养成先用 wavemngr 查看当前 MATLAB 支持哪些小波族再用 wavefun 画出小波函数和尺度函数的形状确认它的振荡次数和支撑长度。这两个命令不产生业务结果但对理解手头的小波基非常有帮助。% 查看 sym 族里可用的所有小波名称 wavemngr(read, sym); % 画出 sym8 的尺度函数 phi 和小波函数 psi精度取 8 层迭代 [phi, psi, xval] wavefun(sym8, 8); plot(xval, psi, LineWidth, 1.5); title(sym8 小波函数形状);wavefun 的第二个参数是迭代次数表示用尺度方程迭代多少轮来生成函数曲线工程上看 6 到 8 次就足够稳定。画出来之后重点看两点一是小波的振荡次数它对应消失矩阶数振荡越多表达平滑信号的能力越强二是支撑长度也就是函数非零部分的宽度它决定了对突变点的定位精度。这两者总是矛盾的支撑越短越能定位突变但越没能力描述光滑曲线所以选基本质上是在“定位精度”和“表达平滑度”之间做取舍。对应到落地经验做齿轮箱故障诊断时冲击特征明显我一般先用 db2 或 db3做心电信号去噪时基线漂移和工频干扰混在一起sym8 比 db 系列更稳做表面粗糙度这类缓变信号分析coif4 往往比短支撑小波少了很多不必要的细节波动。当然这只是起点第 6 章会给一个十分钟跑完的定量选基脚本用重构误差说话而不是靠感觉猜。3. 一维信号小波去噪落地从 wavedec 到阈值参数的手把手代码3.1 最小可运行去噪wdenoise 一行跑通如果只想要一个能立刻看到效果的最小代码wdenoise 是最短路径。它把多层分解、噪声估计、阈值化、重构全封装在内部适合在拿到新信号时先跑一遍观察去噪结果是否符合预期再决定要不要深入做手动阈值。% 假设 x 是一维含噪信号采样率 Fs 1000 Hz % 第 1 个参数是原始信号第 2 个参数是分解层数 [xden, ~, info] wdenoise(x, 5, ... Wavelet, sym8, ... % 使用 sym8 小波基 DenoisingMethod, Bayes, ... % 基于贝叶斯准则估计阈值 ThresholdRule, Soft); % 使用软阈值 % info 结构里可以直接读到噪声标准差估计值和各层阈值 disp(info);这段代码里层数取 5适用于长度几千到几万点的信号Wavelet 指定小波族换成 db4、coif4 效果会不同DenoisingMethod 控制阈值如何生成Bayes 是我个人比较常用的一个选项它假设细节系数服从某种分布然后估计最优阈值对信噪比不高的信号更宽容。ThresholdRule 的 Soft 即软阈值系数超过阈值后向零收缩Hard 即硬阈值超过阈值的系数原样保留。去噪后的信号可以立即做一次快速验证subplot(2, 1, 1); plot(x); title(原始含噪信号); subplot(2, 1, 2); plot(xden); title(wdenoise 去噪结果);运行这一段如果发现端点出现异常波动或者噪声没完全压住不用急着调 wdenoise 内部参数先切换到下一小节的“手动阈值”流程把每一步拆开来控制。这里有一个容易误解的地方wdenoise 虽然是“一行函数”但它内部的噪声估计依赖的是第一层细节系数的统计量如果原始信号本身包含强冲击成分这个估计会被冲击污染导致阈值偏高、把有用的瞬态成分一起削弱。所以 wdenoise 适合做快速原型验证不适合直接当作生产链路里不可控的一环。落地时我通常用它做 baseline基线拿到一个可对比的去噪效果后再用手动流程固定参数。3.2 手动阈值全链路wavedec 拆解、阈值计算、系数处理与重构把去噪过程拆开的意义在于能准确知道每一步对结果的影响方便针对不同信号调参数。完整流程分为四步——分解、估计噪声标准差、计算并施加阈值、重构。第一步和第四步分别由 wavedec 和 waverec 完成中间的阈值部分可以自己控制。% 1. 多层小波分解5 层sym8 小波基 [C, L] wavedec(x, 5, sym8); % 2. 用第一层细节系数估计噪声标准差 % C 的排列顺序是 [第 5 层近似, 第 5 层细节, ..., 第 1 层细节] % L(end-1) 是第一层细节系数的长度所以取 C 末尾那段就是 cD1 d1 C(end - L(end-1) 1 : end); sigma median(abs(d1)) / 0.6745; % 3. 计算固定阈值并施加软阈值 thr wthrmngr(sqtwolog, C, L, sigma); % 统一阈值规则 Cden C; % 细节系数部分做软阈值近似系数保持不变 Cden(L(1) 1 : end) wthresh(C(L(1) 1 : end), s, thr); % 4. 重构去噪信号 xden waverec(Cden, L, sym8);第 2 步的噪声标准差估计是整段代码的关键。0.6745 这个常数来自标准正态分布的四分位距与标准差的换算关系用中位数绝对偏差反推标准差比直接用 std 稳健得多——即使细节系数里还残留个别大幅值冲击中位数也不会被整体抬高。这也是为什么小波去噪在处理含瞬态冲击的信号时优于很多常规滤波器的原因之一它先做一个稳健的噪声水平估计再决定阈值。第 3 步里 wthrmngr 的作用是计算阈值。函数名是“阈值管理器”的缩写它支持多种规则sqtwolog 是其中最常用的一种固定阈值形式公式是 sqrt(2 * log(n)) 乘上 sigman 是细节系数的长度。把 sigma 作为第四个参数传入阈值就被自动缩放到和系数同一个量级。这一步如果漏传 sigma 或误传一个固定值阈值就会偏大或偏小去噪效果直接翻车。细节系数处理完毕后第 4 步的 waverec 会把系数还原成信号注意它必须和分解时使用同一个小波基否则重构结果会南辕北辙。阈值规则特点适用场景sqtwolog固定阈值公式简单阈值偏大噪声均匀、信号较长的场景rigrsure基于 Stein 无偏风险估计阈值偏保守弱信号、信噪比低怕削掉有效成分heursure混合规则自动在固定与无偏之间选择不确定信号特征时的默认选择minimaxi极小极大准则类似固定阈值但更温和细节系数中有少量强突变时选阈值规则时不必迷信默认可以先跑一版 sqtwolog看残差里是否还有明显噪声如果还有换 rigrsure 再跑一遍对比残差的白噪声程度。残差序列应当尽量接近纯噪声这一段如果呈现出周期性或局部相关说明阈值没压干净或阈值过度作用导致信号成分被截断。3.3 分解层数与阈值的三条经验规则和一个验证手段分解层数是另一个直接影响结果的参数。第一条规则层数不能超过 floor(log2(N))N 是信号点数这是算法边界超过直接报错。第二条规则不是层数越多越好每一层会把信号再次隔点抽样信息量减半如果层数过深最低频的近似系数只剩几个点表达趋势的能力大大下降。第三条规则层数和阈值要一起考虑因为每层细节系数的噪声水平相同但信号成分的能量分布不同浅层细节里往往是高频噪声占大头深层细节里可能藏着有效信息盲目对每层用同一个阈值会误伤深层细节。验证去噪效果的一个常用手段是计算信噪比提升量。如果手里有干净的参考信号 xclean可以用下面的代码量化snr_before 10 * log10(sum(xclean.^2) / sum((xclean - x).^2)); snr_after 10 * log10(sum(xclean.^2) / sum((xclean - xden).^2)); fprintf(去噪前 SNR %.2f dB去噪后 SNR %.2f dB\n, snr_before, snr_after);如果没有干净参考退而求其次看残差residual x - xden 应当近似白噪声用 autocorr 画前几十个滞后除了滞后 0 以外不应该有明显规律。这个验证手段比肉眼看图可靠因为人眼很容易被平滑的曲线迷惑忽略了信号细节是否真的被保留。4. 用小波包做故障特征提取节点排序与能量占比的代码实现4.1 为什么小波包比小波分解多一层“细化”空间小波分解有个固有限制每一层只对近似分量继续分解细节分量不再处理。这让它在去噪场景很合理因为噪声主要落在细节系数里压掉它们就是去噪逻辑但在特征提取场景这个限制就变成短板——故障信号的有用特征往往就藏在高频细节的某个窄带里不继续分解就提取不到。小波包分解的思路很直接细节分量同样继续做低通和高通拆分每一层都得到完整的二分频谱。代价自然是系数的数量翻倍但换来的是对频段的精细刻画能力。MATLAB 里小波包分解的入口是 wpdec它返回一个小波包树对象节点用“层数 节点索引”来定位例如 [3 2] 表示第 3 层第 2 个节点。读取某个节点的系数用 wpcoef单独重构某个节点的时域信号用 wprcoef。下面的代码建立一棵 3 层小波包树并读出第 3 层全部节点的系数% 建立 3 层小波包树shannon 熵准则控制最优分解 wpt wpdec(x, 3, db4, shannon); % 遍历第 3 层 0 到 7 号节点逐个读取系数 coefs cell(1, 8); for k 0 : 7 coefs{k 1} wpcoef(wpt, [3 k]); end % 绘制第 3 层第 4 个节点系数的时域波形观察其振荡特征 plot(coefs{5});wpdec 的第三个参数是小波基第四个参数是熵准则。熵准则决定了小波包分解时“最优树”的选择方式常用的是 shannon 熵也可以换成 log energy 熵大多数场景下用 shannon 就足够稳定。这里注意 wpdec 返回的是自适应选出的分解树不代表所有叶节点都来自完全展开如果后续要固定结构做批量特征提取可以在 wpdec 之后统一扩展到指定深度以保证不同信号之间特征维度对齐。4.2 节点频率顺序反直觉探路法与能量占比特征向量小波包的经典坑在于节点的编号顺序和频率由低到高的顺序并不一致。因为每一层隔点抽样会让高频分量的频谱发生翻转节点的编号遵循的是滤波器组的二分树顺序不是频率轴的递增顺序。第一次用这个功能的人经常会遇到一种现象按节点序号 0、1、2、3 去画系数波形发现振荡频率并没有单调变化第三、第四个节点的波形看起来比第二个节点更“快”。这就是节点顺序错位了。处理方式有两种。第一种是直接用已知结论做节点重排三层小波包在频率由低到高时对应的节点顺序是 [0 4 6 2 3 7 5 1]四层是这个规律继续向外扩展。第二种方式是探路法不想背结论时可以构造一个由若干不同频率正弦波叠加的合成信号分解后统计每个节点的系数的过零率过零率高频率就高排序结果一目了然。探路法代码见下前面那个结论不够用时用这个办法自己确定顺序。% 探路信号三个频率成分叠加频率逐渐升高 t 0 : 0.001 : 1; probe sin(2 * pi * 20 * t) sin(2 * pi * 80 * t) sin(2 * pi * 200 * t); wptProbe wpdec(probe, 3, db4, shannon); % 统计第 3 层每个节点的过零率粗略判断频率高低 for k 0 : 7 c wpcoef(wptProbe, [3 k]); crossings sum(abs(diff(sign(c))) 0); zeroRate(k 1) crossings / numel(c); end [~, freqOrder] sort(zeroRate, ascend); freqOrder freqOrder - 1; % 转成 0 起始的节点编号探路法的输出 freqOrder 就是从低频到高频的节点顺序。把它保存下来之后处理同类型信号时直接复用。做完节点重排后下一步就是构造能量占比特征向量这是故障诊断里最常用的特征之一计算方式简洁但有区分度% 计算第 3 层每个节点的能量 E zeros(1, 8); for k 0 : 7 c wpcoef(wpt, [3 k]); E(k 1) sum(abs(c).^2); end % 按探路法得到的频率顺序重排能量向量 E_sorted E(freqOrder 1); % 归一化为百分比作为特征向量 E_feature E_sorted / sum(E_sorted) * 100;能量占比特征的本质是把信号在不同窄带内的能量分布转成一组百分比归一化之后与绝对振幅无关只看相对分布因此对不同采集条件下的同类信号有更强的可比性。使用时需要注意系数向量的长度可能不同但能量计算本身不受长度影响如果担心长度差异带来的偏差可以对每个节点用 numel(c) 做平均能量归一。特征向量构造出来后可以用 PCA 或 t-SNE 压缩可视化也可以直接送入分类器配合标注样本做故障类型识别。需要补充的是小波包树上的节点能量并不具有正交基之间的绝对独立性相邻节点之间会因滤波器非理想矩形特性产生交叠所以能量占比特征是“工程特征”而非“严格物理量”但它到分类模型的落地上通常够用。追求更高精度时可以再追加每个节点的幅度峰值、均方根值或谱峭度特征维度从 8 扩到 16 至 24。5. 小波分析避坑手册边界效应、小波基误选与阈值翻车的四条记录5.1 去噪后端点冒出假毛刺现象去噪后的信号中间段很干净但首尾几十个采样点出现无规则的明显抖动幅度甚至超过原始信号峰值看起来像人为拼接了一段假数据。原因离散小波分解默认对信号做边界延拓MATLAB 里常见的延拓方式是对称延拓。信号端点往往不是周期连续的延拓产生的边界不连续会让边界附近的细节系数出现异常大值阈值处理后又不能完全消除最终在重构时把边界区域的伪波动带回信号。解决先在分解前对信号做一次显式的边界延拓处理完再截回原长度。常见做法是用 wextend 把信号左右各延长原长度的 10% 到 20%延拓模式根据信号长周期或短周期特性选 sym 或 ppd。处理完成后用 wkeep 截取中间段。另一种补救是换用支撑更短的小波比如从 sym8 换到 db3边界效应会减轻但对平滑信号的表达力也会下降所以先延拓再处理是更稳妥的顺序。5.2 选了 db1 做平滑信号去噪趋势直接变形现象用 db1即 Haar 小波对一个缓慢变化的温度曲线去噪得到的信号变成阶梯状折线原始趋势的圆滑形状消失高频细节里也没有提取出有效信息只剩一堆矩形脉冲。原因db1 只有一阶消失矩也就是说它只能精确表达常数信号对一次以上的多项式趋势无能为力。平滑趋势在小波基看来是高频成分会被当成噪声处理掉。这是小波基的消失矩与信号平滑度不匹配导致的典型问题。解决对平滑信号至少选择消失矩为 4 或 6 的小波基sym4、sym6、coif4 在这个场景都比 db1 合适。判断方法可以先用 wavefun 看小波函数的振荡次数再对原始信号做一次分解重构算重构误差。如果平滑信号用 db1 分解再重构误差会明显大于 sym8这就是一种不需要理论计算就能完成的选基降错方式。5.3 手动阈值没按层自适应把有效信号削没了现象手动去噪时设了个固定阈值 0.2处理完发现噪声确实小了但信号里的尖峰脉冲也一起变平了能量和峰值都下降明显看起来像把信号“削”了一遍。原因不同分解层的细节系数尺度不同浅层细节系数幅值小且密集深层细节系数幅值大但稀疏。用一个全局固定阈值去处理所有层浅层可能阈值偏大、抑制过度深层又可能阈值偏小、漏掉噪声。阈值必须和每一层系数的统计尺度匹配。解决改为按层独立计算阈值。对每层细节系数单独估计 sigma再计算该层的阈值然后只对那一层做阈值化。实现上只需把 3.2 节里的单次阈值过程放进一个循环逐层处理细节系数近似系数保持不动。这样每层用了与自身尺度匹配的阈值信号恢复度明显更好。5.4 特征提取时按节点编号直接用换了信号特征就乱现象按小波包节点 0 到 7 顺序提取能量特征画出来频率分布看起来乱序换一个信号之后乱序规律好像又不同送入分类模型后准确率始终不高。原因前面说过小波包节点编号不是按频率递增排列而是按滤波器的二分树顺序排列。直接按编号顺序整理特征等于把一个固定顺序的物理量序列随机打乱了再交给模型。模型虽然理论上可能自己学出映射关系但对样本量小的训练集这种顺序错位会显著增加学习难度。解决先用探路法或者三层固定结论 [0 4 6 2 3 7 5 1] 把节点顺序修正到频率递增再构造特征向量。修正后的特征随频率单调分布可解释性也强很多故障频段对应的特征索引变得有物理意义。如果是更高层数的小波包需要重排建议分两个阶段验证先画能量柱状图确认峰值确实出现在你预期的频段再进入批量特征提取。6. 小波基与分解层数怎么快速敲定一个十分钟跑完的评估脚本6.1 以重构误差为判据的选基脚本小波基的选择方法很多但工程上最直接有效的判据是重构误差。思路是对信号本身做一次分解再完整重构如果小波基能很好地表达这个信号重构回来的信号应该和原始信号几乎一样如果表达不充分重构误差就会偏大。注意这一步要在原始信号上做不要先加噪声再去噪那样会把噪声影响混进选基判断里。wnames {db2, db4, sym5, sym8, coif4}; levels 5; % 统一使用 5 层分解 mseVal zeros(size(wnames)); for i 1 : numel(wnames) [C, L] wavedec(x, levels, wnames{i}); xr waverec(C, L, wnames{i}); mseVal(i) mean((x - xr).^2); end [bestVal, bestIdx] min(mseVal); fprintf(最优小波基%s重构MSE %.3e\n, wnames{bestIdx}, bestVal);跑完之后顺手把第二好、第三好的小波基也记录下来因为最优解和次优解在工程效果上往往差别很小而次优解可能支撑更短、计算更快更适合长序列批量处理。层数的敲定用另一个简单指标从第 1 层开始逐层增加算每一层近似系数与原始信号的相关系数。相关系数降到 0.9 以下说明该层丢失了太多趋势信息层数不应再加深。for lev 1 : 6 [C, L] wavedec(x, lev, sym8); approx C(1 : L(1)); % 第 lev 层近似系数 approxUp upcoef(a, approx, sym8, lev); % 重构到原长度 corrVal(lev) corr(x(:), approxUp(:)); end stem(corrVal); ylabel(近似系数与原始信号的相关系数);这个脚本的妙处在于它把两个最玄学的参数变成了可量化的折线图。我自己的习惯是选基看重构 MSE定层看相关系数两个脚本加起来不到十分钟却省掉了后期无数次“为什么这个特征提取效果不如预期”的返工。有一次做某设备的振动数据分析我直接沿用了之前的 db4结果某特征在正常样本和故障样本之间完全没有区分度排查到最后发现 db4 对该信号的支持长度不足换用 sym8 后特征分离度立刻改善。所以从那以后任何新数据集到手我都先跑选基脚本再谈参数优化。希望这一套流程也能帮到你让你少走一段凭感觉调参数的弯路。本文还有配套的精品资源点击获取