小波包分解MATLAB程序：时频分析与信号降噪的完整实现

简介一套面向信号处理初学者的MATLAB小波包分解程序将小波分解与小波包分解的基本原理同实际代码相结合帮助读者理解非平稳信号的时频分析、小波基选择、阈值去噪以及信号重构流程。压缩包共2个文件均为m脚本大小仅2KB分别针对原始信号与加噪声信号编写代码结构清晰便于直接运行和修改参数。该资源已有3853人下载学习适合用作课程实验、毕业设计或自学入门的配套参考。程序中完整展示了数据加载、预处理、小波与小波包系数分解、分解结果可视化、软硬阈值去噪、信号重构以及原始与重构信号对比分析等环节初学者可借此观察不同小波基和分解层数对效果的影响深入理解小波分析在信号处理、图像分析和故障诊断等场景中的应用。小波包分解在频带细分上更具优势可帮助读者掌握更精细的信号分析手段。1. 小波包分解 MATLAB 程序一套能同时看透时域与频域的信号处理工具箱有一回我拿到一段设备振动信号50 Hz 工频和故障冲击特征几乎叠在一起时域图看不出名堂直接带通滤波又会把冲击前沿滤圆。换成小波分解之后频率成分和时域位置同时拿到了。这套「小波包分解 MATLAB 程序」解决的就是这个核心问题把混合信号按频带切碎还保留每个成分出现的时间信息。它不只是一个演示而是从多层小波分解、小波包分解、阈值降噪、单支重构到能量特征提取的完整可运行程序适合做信号处理的工程师和研究生直接下载、复现、改造。小波分解不是层数越多越好——层数越多频率切得越细但每层系数长度减半时域定位精度跟着下降小波包把高频分支继续拆分但代价是计算量翻倍。怎么选层数、怎么选小波基、怎么处理重构边界这套程序里都给了可调参数。2. 工程结构与运行准备从目录到第一张图的完整链路2.1 目录结构与模块职责拿到程序第一件事是看它的目录。小波分解程序一般不会只有一个.m脚本而是拆成「主入口 各个算法封装 数据生成」几个文件。这套程序常见的组织方式如下├── main_demo.m % 主入口生成仿真信号并串起全流程 ├── data_gen.m % 生成仿真信号正弦叠加 白噪声 ├── wavelet_decomp.m % 封装 wavedec 的多层小波分解 ├── wavelet_denoise.m % 封装 ddencmp wdencmp 的降噪流程 ├── wavelet_packet.m % 封装 wpdec、wpbest、节点重构的小波包流程 └── feature_extract.m % 将小波包能量转成特征向量文件不多但每个文件职责非常明确。主入口负责整个流程的编排和绘图数据生成模块负责造一个「已知答案」的测试信号——你知道它包含 50 Hz、140 Hz 和噪声这样后面看分解结果就能立刻判断算法有没有跑对。真正的算法模块不画图、不打印只负责计算和返回结果。做二次开发时合不合并都无所谓但拆成这种「入口、数据、算法、可视化」四层结构之后改阈值策略不会碰坏读数据部分换小波基也不会让主流程作废。文件核心调用输出用途data_gen.mrandn、sin已知频率成分的仿真信号wavelet_decomp.mwavedec、wrcoef多层小波系数与各层重构信号wavelet_denoise.mddencmp、wdencmp降噪后的干净信号wavelet_packet.mwpdec、wpbest、wprcoef小波包树、节点信号feature_extract.mwenergy各频带能量特征向量如果你拿到的版本只有一个大脚本也建议按这个思路拆开再改不然调一次阈值要把整个脚本跑一遍效率很低而且容易误改到别的模块。2.2 运行环境与被忽略的 Wavelet Toolbox 检查这套程序依赖 MATLAB 的 Wavelet Toolbox。大部分情况下 MATLAB 默认带这个小波工具箱但实验室的授权版本有时会精简掉遇到undefined function wavedec这种报错十有八九不是代码问题而是工具箱没装上。建议在主入口第一行就做一次检测比运行到一半再报错省事得多% 检查 Wavelet Toolbox 是否可用 v ver( wavelet ); if isempty( v ) error( 未检测到 Wavelet Toolbox请先安装小波工具箱再运行本程序 ); end disp( Wavelet Toolbox 检测通过 );ver(wavelet)返回一个结构体数组如果工具箱没安装返回是空的isempty就能拦下来。程序对 MATLAB 版本要求不高R2016 之后的所有版本都能跑因为小波工具包的这些核心函数十几年没有破坏性变化。唯一要注意的是新版 MATLAB 更推荐wdenoise这种新式 API但这套程序沿用的是经典的wavedec/wdencmp接口二者并存不存在兼容问题。2.3 主入口与数据规范先跑 demo 再换自己的数据先看主入口前半段它负责生成一段仿真信号并统一数据格式% main_demo.m 前段生成仿真信号模拟设备振动 噪声 fs 1000; % 采样率 1000 Hz t ( 0:1023 ) / fs; % 1.024 秒的时间轴共 1024 个采样点 x sin( 2*pi*50*t ) ... % 50 Hz 工频分量 0.5 * sin( 2*pi*140*t ) ... % 140 Hz 故障特征分量 0.3 * randn( size( t ) ); % 高斯白噪声标准差 0.3 if size( x, 1 ) 1 x x(:); % 统一转成行向量 end这里有几个参数要解释。fs 1000是采样率它决定了奈奎斯特频率是 500 Hz也就是说这套分析最多能分辨 500 Hz 以下的成分。t从 0 到 1.023 秒步长是 1/fs正好 1024 个点。sin(2*pi*50*t)生成 50 Hz 正弦波0.5是幅值系数140 Hz 分量幅值是 0.5噪声标准差 0.3——信噪比不算高正好用来验证小波降噪的效果。最后一行x(:)很关键。MATLAB 的小波函数严格区分行向量和列向量wavedec两种都能处理但后面做数据拼接、特征组合时行向量和列向量混淆会导致维度错误。统一成行向量是省事且稳妥的习惯。替换成自己的数据时常见做法是% 替换仿真信号从 .mat 或 Excel 里读入你的数据 load( your_signal.mat ); % 确保工作区里有变量 x % x xlsread( your_data.xlsx, Sheet1, A1:A4096 );替换后立刻检查两条第一x必须是一维向量如果你的数据是多通道矩阵先选其中一列第二fs必须改成采集设备的真实采样率这个值直接影响频带划分结果写错的话所有频率分析都是错的。2.4 跑通的标志两张图和一个无报错的工作区运行主入口后程序会依次执行小波分解、小波包分解最后弹出两个 figure 窗口。第一个窗口显示各层细节系数重构信号第二个窗口显示小波包树和能量分布。只要命令窗口没有 Error并且能明显看到 50 Hz 和 140 Hz 成分被分到不同频带就说明环境没问题、流程跑通了。先确认这两张图再深入改代码不要一上来就动阈值。3. 核心模块拆解小波分解从 wavedec 到 wdencmp 的完整实现3.1 多层分解wavedec 的系数排列逻辑小波分解的入口是wavedec程序里对应的封装是这样% wavelet_decomp.m 核心代码 N 4; % 分解层数 wname db4; % Daubechies 4 小波紧支集正交小波 [ C, L ] wavedec( x, N, wname );wavedec的输出只有两个变量C和L新手最容易在这里看晕。C不是矩阵而是一个按顺序拼起来的行向量先放第 N 层的近似系数再依次放第 N 层、第 N-1 层……直到第 1 层的细节系数。L是一个长度 N2 的行向量记录每一段在C中的长度以及原始信号的长度。拆开看就是这样cA4 C( 1 : L(1) ); % 第 4 层近似系数 cD4 C( L(1)1 : L(1)L(2) ); % 第 4 层细节系数 cD1 C( end-L( N1 )1 : end ); % 第 1 层细节系数在 C 的最后L的最后一个元素固定等于原始信号的长度length(x)这个特性在重构和单支提取时非常有用。理解C/L这对变量是整个小波分解模块的地基后面所有阈值、重构、能量计算都建立在它俩之上。wname db4是 Daubechies 4 小波滤波器长度 8。为什么选它而不是更简单的 Haardb1Haar 小波对脉冲特征敏感但重构出来的信号台阶感严重db4更光滑能保留冲击前沿又不至于把信号滤成方块。对振动信号、语音信号来说db4是实践经验里比较稳妥的起点。如果要分析更平滑的缓变信号可以换db8如果要更短的时域支撑可以换sym3。这套程序里wname是一个独立变量后面所有函数都引用它改一处即可全局生效。3.2 单支重构把特定频带从混合信号里抠出来分解完只有系数看不见时域形态。程序里用wrcoef做单支重构把某一层的细节系数映射回原始时间尺度% 提取第 2 层细节信号对应 125~250 Hz 频带 d2 wrcoef( d, C, L, wname, 2 ); % 提取第 4 层近似信号对应 0~31.25 Hz 频带 a4 wrcoef( a, C, L, wname, 4 ); figure; subplot( 2, 1, 1 ); plot( t, d2 ); title( D2 细节信号 125~250 Hz ); subplot( 2, 1, 2 ); plot( t, a4 ); title( A4 近似信号 0~31.25 Hz );wrcoef的第一参数d或a决定提取细节还是近似第二个和第三个必须是wavedec原样返回的C、L最后一个数字是层号。它的好处是不管提取第几层返回的信号长度始终和原始信号一样——因为它内部做了一次子带重构把系数补零后过滤波器。这样plot可以直接画在同一个时间轴上省去自己做上采样的麻烦。对 fs1000、N4 的配置各层频带大致如下层频带范围 (Hz)信号成分D1250 ~ 500噪声主导D2125 ~ 250140 Hz 故障特征D362.5 ~ 125噪声 干扰D431.25 ~ 62.550 Hz 工频A40 ~ 31.25低频趋势这个表不用背但建议贴在手边。跑完程序对照一下140 Hz 特征出现在 D250 Hz 工频出现在 D4说明这套频带划分和预期一致。往后你处理自己的数据时目标频率落在哪一层直接看这个频带表就能推算出来。3.3 ddencmp wdencmp 阈值降噪参数详解降噪模块是小波分解最实用的部分。程序里通过ddencmp自动生成降噪参数再用wdencmp执行阈值处理% wavelet_denoise.m 核心代码 [ thr, sorh, keepapp ] ddencmp( den, wv, x ); xd wdencmp( gbl, C, L, wname, N, thr, sorh, keepapp ); snr_before 10 * log10( sum( x0.^2 ) / sum( (x - x0).^2 ) ); snr_after 10 * log10( sum( x0.^2 ) / sum( (xd - x0).^2 ) ); fprintf( 降噪前 SNR %.2f dB降噪后 SNR %.2f dB\n, snr_before, snr_after );ddencmp(den,wv,x)是自动阈值生成器第一个参数den表示做降噪第二个wv表示按小波而非小波包方式处理输出三个参数——thr是计算出的阈值sorh是阈值方式s软阈值或h硬阈值keepapp是一个标志位1 表示保留近似系数不做任何处理。这三个参数可以直接传给wdencmp。wdencmp的第一个参数gbl表示全局统一阈值即所有层都用同一个thr。如果你希望每一层独立阈值把gbl换成lvdthr改成每层阈值组成的向量。软阈值s把所有系数向零收缩信号更平滑硬阈值h小于阈值的置零大于的保留原值冲击幅值保留得更真实但重构信号可能带振铃。程序默认用软阈值降噪效果更稳定代价是冲击峰值会被压低一点。3.4 用小波包第一次拆高频wdencmp 的 wp 模式除了对wavedec的结果做阈值小波包里还有一个更直接的用法wdencmp的wp模式会先自动做小波包分解再对全部节点统一阈值最后直接返回降噪后的信号。它最适合的场景是噪声不是集中在高频段而是均匀铺满整个频带——这时候普通小波只拆低频高频段的噪声没法单独压。程序里保留了这个选项% 小波包降噪自动分解到第 3 层对全部节点做软阈值 xd_wp wdencmp( wp, x, wname, 3, thr, sorh, keepapp );注意这里传的是原始信号x而不是系数C/L因为小波包路径需要自己重新做分解。第 3 层意味着把 0~500 Hz 频带切成 8 段每个节点带宽 62.5 Hz。相比普通小波的 D1250~500 Hz小波包能把 140 Hz 和 200 Hz 这种相邻不远的频率成分拆到不同节点这是它最值得用的理由。4. 核心模块拆解小波包从熵准则到最优树与节点重构4.1 熵准则wpdec 的第四个参数是如何影响树的小波包分解的核心函数是wpdec它构造的是一棵满二叉树% wavelet_packet.m 核心代码三层满小波包分解 T wpdec( x, 3, db4, shannon );参数依次是信号、分解深度、小波名、熵准则。深度 3 意味着对每个节点都继续向下拆分直到第 3 层总共产生 2^3 8 个叶子节点。熵准则的作用是给每个节点算一个「代价」代价越小说明这个节点的信息越集中、越值得保留代价越大说明这个节点基本上是噪声或冗余信息后面剪枝会优先剪掉它。常用的熵准则有几个熵类型定义思想适合场景shannon对系数能量分布求信息熵持续振动、能量散布广的信号log energy对系数绝对值的对数求和存在明显尖峰脉冲的信号threshold统计超过阈值的系数个数稀疏信号、瞬态冲击norm系数的 p 范数需要自定义压缩程度时程序默认用shannon因为它对「系数是否均匀」最敏感。若信号里有明显冲击特征可以试试log energy它会放大系数小的节点的代价让剪枝更倾向于保留冲击所在的节点。选熵准则没有绝对正确答案但是选错了后面最优树的结构会明显不同这是小波包模块最需要反复试验的地方。4.2 最优小波包wpbest 剪枝与可视化确认满树不一定是最优的。如果某个父节点的代价已经小于等于两个子节点代价之和说明拆分这个节点没有带来信息增益就该剪掉。程序里用wpbest自动完成这个操作% 基于熵代价剪枝得到最优树 T_best wpbest( T ); figure; plot( T_best ); % 可视化最优小波包树plot(T_best)画出来的树不是满二叉树而是部分分支被剪掉的稀疏树。看这张图能直观判断哪些频带被继续细分了哪些频带保留在一个较粗的节点里。在信号处理里这非常有用——噪声频带通常会被剪掉分支而特征频带会保留多层细分。剪枝后画图时要注意wpbest是对整棵树做优化它不会修改原始信号和原始小波配置只是返回一棵新的树对象。你可以同时保留T和T_best做对比看剪枝到底剪掉了哪些节点。4.3 节点系数与单支重构wpcoef 和 wprcoef 的区别小波包树的每个节点都有唯一索引[depth, index]例如[3, 0]表示第 3 层第 0 个节点。提取节点信息有两种方式程序里都有体现% 提取第 3 层第 0 个节点的分解系数 cfs wpcoef( T_best, [ 3, 0 ] ); % 重构第 3 层第 0 个节点的时间序列长度与原始信号相同 s0 wprcoef( T_best, [ 3, 0 ] );wpcoef返回的是该节点的小波包系数长度不是原始信号长度而是该频带内经过二抽样的系数序列适合做特征计算wprcoef返回的是把该节点系数补零后重构出的时域信号长度和原始信号完全一致适合直接画图观察。这两者很容易搞混我早期的理解是「反正都是取节点」结果拿wpcoef的结果去plot横轴长度对不上找半天才发现是函数用错了。要看某个节点对应哪个频带可以用% 读取节点 [3, 2] 的频带信息归一化频率 fb read( T_best, freqband, [ 3, 2 ] );read返回的是归一化频带需要乘以fs/2才能换算成实际 Hz。这一点很容易被忽视直接拿返回值当 Hz 用会差一半。4.4 小波包能量特征从树到分类器输入的最后一公里小波包模块还有一个高频用途把各节点能量归一化成特征向量喂给分类器做故障诊断。程序里用wenergy一步完成% 提取各节点能量百分比 [ E, node ] wenergy( T_best ); feature_vec E(:); % 转成行向量作为一条样本特征 figure; bar( node, E ); xlabel( 节点索引 ); ylabel( 能量百分比 );wenergy对当前树的所有叶子节点计算能量占比E是各节点能量百分比node是对应的节点索引。特征向量feature_vec的维度等于叶子节点数对深度 3 的满树是 8 维剪枝后会少一些。这个特征相比 FFT 频谱的优势在于小波包能量是按频带统计的天然对噪声和幅值波动不那么敏感而且保留了时域局部信息。实际使用中把 8 维能量特征拼接上时域统计量峰值、峭度、均方根再送进 SVM 或随机森林分类效果通常比直接用原始信号好不少。4.5 一次做完的完整小波包处理流程把上面的步骤串起来就得到小波包模块的完整流程先wpdec建满树再wpbest剪枝然后用wprcoef提取目标频带时域信号最后wenergy出特征。顺序不能乱尤其是wenergy一定要在wpbest之后调用否则对满树算能量会把被剪枝节点也统计进去特征维度就和你预期不一致了。5. 常见问题与排查五个最容易把小波分析跑翻车的细节5.1 小波名无效Error using wfilters现象执行wavedec(x, 4, db)或wpdec(x, 3, dbe4, shannon)时命令行报Invalid wavelet name。原因db是整个 Daubechies 小波家族的统称不是具体某个小波分解函数必须指定到具体的阶数dbe4是拼写错误正确写法是db4。常见的还有sym少写阶数、coif和coif4混淆。解决写一个小波名验证确定能加载再进行分解。wname db4; [ Lo_D, Hi_D ] wfilters( wname ); if isempty( Lo_D ) error( 小波名 %s 无效请换成 db/sym/coif 系列并带阶数, wname ); endwfilters返回分解低通和高通滤波器系数如果小波名不可识别返回是空的。养成在分解前跑一次wfilters的习惯能挡掉一大半手误。5.2 分解层数太大Index exceeds array bounds现象wavedec或wpdec时报Index exceeds array bounds或者提示信号长度不足。原因小波分解的每一层都要做一次二抽样如果层数 N 过大到后面某一层信号长度已经小于滤波器长度分解就无法继续。常见于把 100 Hz 采样、长度 256 的信号强行分解 6 层。fs 越小、信号越短最大可分解层数越低。解决在分解前根据信号长度计算最大层数做一个硬约束。N 4; Nmax floor( log2( numel( x ) ) ); if N Nmax error( 分解层数 %d 超过当前信号最大可分解层数 %d, N, Nmax ); endnumel(x)是信号样本数floor(log2(...))是理论上限。这只是粗略估计实际还要考虑滤波器长度但对绝大多数情况已经能提前拦住报错。5.3 手动改 C 后重构失败C 和 L 不匹配现象对C向量里的某些系数直接赋零或改数值后调用waverec报错Vector C and L must ...或者重构出来的信号长度不对。原因wavedec返回的C是一个按段拼接的长向量L记录每段长度。手动修改C时如果增删了元素而没同步更新Lwaverec就无法正确切分每一段自然重构失败。解决不要直接对C/L做手术。要修改某一层用wrcoef先取出该层完整信号改完再与其他层相加或者用wthresh对整个C做阈值处理它会自动保持向量结构不变。需要手动拼接时务必同步维护L让它始终满足sum(L(1:N1)) length(C)且L(end) length(x)。5.4 阈值得当却降噪过头信号变「平头」现象降噪后信号光滑了但原本的冲击特征也被削平幅值明显变小或者反过来阈值太小噪声一点没少。原因ddencmp自动生成的thr是全局统一的它依据的是整个信号的方差估计。如果噪声不是均匀分布在所有层比如高频层噪声大、低频层噪声小全局阈值就会导致高频压得不够、低频压得过头。解决换成每层独立阈值模式lvd并手动为每层算一个通用阈值。常见做法是使用thselect配合各层细节系数thr_vec zeros( 1, N ); for k 1:N dk wrcoef( d, C, L, wname, k ); thr_vec(k) thselect( dk, heursure ); % 启发式 SURE 阈值 end xd wdencmp( lvd, C, L, wname, N, thr_vec, s, 1 );thselect的heursure会综合固定阈值和 SURE 风险估计对每层单独算一个自适应阈值效果比全局gbl稳很多。改完之后把snr_before和snr_after打印出来对比确认降噪有效再继续。5.5 小波包节点选错频带以为取的是 140 Hz 结果全是噪声现象用wpcoef或wprcoef取某个节点画出来的时域信号形态和噪声一样没有预期中的 140 Hz 正弦特征。原因小波包树的节点编号并不是从左到右严格按频率升序排列的而是遵循小波包库的排列规则简单按节点序号拍脑袋取频带大概率选错。解决不要猜用read直接读节点频带信息换算成实际频率后再决定取哪个节点。fb read( T_best, freqband, [ 3, 2 ] ); fb_hz fb * fs / 2; % 归一化频带换算成 Hz fprintf( 节点 [3,2] 实际频带%.1f ~ %.1f Hz\n, fb_hz );如果读出来不是自己想要的频带换节点索引继续读直到命中目标频率范围。这一步虽然多花十几秒但能避免拿到节点后白做一整轮特征提取。6. 二次开发把固定脚本改造成自适应降噪与特征工程管道6.1 每层独立阈值的自适应降噪封装程序默认的降噪是自动阈值加全局处理实际工程里噪声分布往往不均匀。我一般会把降噪模块改成接受一个结构体参数方便切换gbl和lvd% 自适应降噪封装mode 可选 gbl 或 lvd function xd adapt_denoise( x, wname, N, mode ) [ C, L ] wavedec( x, N, wname ); if strcmp( mode, gbl ) [ thr, sorh, keepapp ] ddencmp( den, wv, x ); xd wdencmp( gbl, C, L, wname, N, thr, sorh, keepapp ); else thr_vec zeros( 1, N ); for k 1:N dk wrcoef( d, C, L, wname, k ); thr_vec(k) thselect( dk, heursure ); end xd wdencmp( lvd, C, L, wname, N, thr_vec, s, 1 ); end end这样调用方只需要传一个模式字符串内部自动完成分解、阈值计算和重构。建议把snr计算也塞进去每次降噪后自动打印前后信噪比肉眼判断不如数值判断可靠。6.2 把特征提取拆成可批处理管道做故障诊断时往往是几百段信号要一起处理。我会将wenergy封装成一个函数输入一段信号输出一条特征向量然后用parfor批量跑function feat extract_wpt_feature( x, wname, depth ) T wpdec( x, depth, wname, shannon ); T wpbest( T ); [ E, ~ ] wenergy( T ); feat E(:); end % 批量处理X 是 样本数x信号长度 的矩阵 parfor i 1:size( X, 1 ) features( i, : ) extract_wpt_feature( X( i, : ), db4, 3 ); endparfor要求每个迭代任务独立小波包分解正好满足这个条件。唯一要注意的是不要在每个迭代里画图图窗操作在parfor里会被禁止或拖慢速度。特征矩阵features每一行是一个样本直接接到分类器的fitcsvm或fitcensemble前面就能跑通一条完整的故障诊断链路。6.3 验证改造是否有效的方法每次改完算法我会用一组已知答案的数据做回归验证拿纯信号加固定噪声跑新旧两版降噪代码对比snr_after特征提取则检查同一类样本的特征向量是否聚类比如用gscatter看前两维能量是否有明显的类别分界。如果新版结果不比旧版好立即回滚参数不要恋战。从那以后我每次处理批量信号都强制先把单样本跑通确认小波名、层数、阈值方式没问题再上批量这套检查清单替我省掉了大量返工。希望这套程序的拆解思路能帮到你。本文还有配套的精品资源点击获取