简介：Matlab小波工具箱（Wavelet toolbox）面向科研人员、工程师及数据分析学习者，用于信号处理、图像分析、噪声去除与模式识别等任务，提供从离散小波变换、多分辨率分析到连续小波变换、模极大重叠离散小波变换的完整函数与交互示例。压缩包共1152个文件，约17.76MB，以792个m函数脚本、227个mat数据、56个png与35个fig图形界面为主，另含jpg、xml、mexw64等配套文件，便于查阅调用与实验复现。已有10737人学习下载，适合作为小波分析入门与工程实践的参考材料。内容覆盖DWT、MODWT、MRA、CWT及小波包、二维图像小波分解与系数统计检验，读者可据此理解不同小波基的选择、分解重构流程和典型应用场景，并借助示例数据与图形界面快速验证算法效果。
1. 小波工具箱的落点：为什么 FFT 看不出的冲击，它能定位到毫秒
一段电机电流信号，前两秒只有干净的 50 Hz 基波，第 2.7 秒突然叠进持续 20 ms 的高频冲击。FFT 上看，这只是一片模糊的边带，你分不清它是噪声还是故障；滤波又容易把冲击本身削平。这类非平稳、短时、突变的信号，就是 Matlab 小波工具箱最典型的落点。它把信号按尺度拆成「近似 + 细节」，低频趋势留在近似系数里，突变和瞬态被压进某几层细节系数，时间位置基本不丢。做旋转机械故障诊断、心电肌电预处理、地震波分析、图像压缩的人，都会碰到它。不过先说明一句：这个工具箱不是「一键出图」的按钮，小波基、分解层数、边界模式任何一个设错，图上照样能画出好看的色块，结论却是错的。
2. 小波基与分解层数：选错基，后面的阈值调得再细也白搭
小波分析里最容易被跳过的一步，就是选基。很多人拿到信号直接wavedec(x, 5, 'db4')，跑通了就以为完事，后面去噪效果不好就反复调阈值，调到最后发现是基和层数的问题。这一章把选型的几个硬指标拆开说，再落到一条能跑通的分解-重构闭环上。
2.1 小波基的四个性质：支撑长度、消失矩、对称性、正交性
小波基不是「名字好听就选哪个」，它由四个性质共同决定，而这四个性质在信号处理里是互相拉扯的。
支撑长度指小波滤波器的长度。dbN的滤波器长度是2N，db4就是 8 个抽头。支撑越短，时域定位越准，越适合抓瞬态冲击；但频域选择性差，频率分不开。支撑越长，频域越窄，但边界效应会往信号内部扩散得更远，而且计算量随之上升。
消失矩决定小波对多项式趋势的压制能力。dbN有 N 阶消失矩，意味着它能把 N-1 次多项式完全压到近似系数里，细节系数上不残留。做趋势去除、基线漂移抑制时，这个指标比什么都重要。心电信号里的基线漂移近似是分段线性，选db4以上就比haar干净得多。
对称性影响相位失真。Daubechies 系列除了db1（也就是haar）之外都不对称，重构信号会有相位偏移，对图像处理不友好。图像压缩和边缘检测常见做法是选双正交小波biorNr.Nd，它用两组不同的滤波器做分解和重构，换来严格的线性相位。
正交性保证能量守恒，重构误差在数值上可以压到机器精度。dbN、symN、coifN、dmey都是正交的；bior系列是双正交，不满足能量守恒，但换来了对称性。
| 家族 | 写法 | 正交 | 对称 | 典型落点 |
|---|---|---|---|---|
| Haar | haar/db1 | 是 | 是 | 阶跃、教学演示 |
| Daubechies | dbN | 是 | 否 | 通用去噪、故障冲击 |
| Symlets | symN | 是 | 近似 | 相位要求不高的重构 |
| Coiflets | coifN | 是 | 近似 | 信号压缩 |
| 双正交 | biorNr.Nd | 否 | 是 | 图像、边缘 |
| 离散 Meyer | dmey | 是 | 是 | 高分辨率分析 |
| 连续小波 | morl/mexh/cmor | 不适用 | 不适用 | cwt时频图 |
这张表里最要紧的一列是最后一列：连续小波只能进cwt，离散小波才能进wavedec。后面避坑章节会专门说这个翻车现场。
2.2 用 wavemngr、waveinfo、wfilters 把家族摸清楚
记不住名字和性质的时候，不要靠猜。工具箱自带三个查询入口，先把它们跑一遍，比翻文档快。
% 列出工具箱里所有可用小波的短名，确认拼写 wavemngr('read') % 打印 db 家族的完整性质表：消失矩、支撑长度、正交性 waveinfo('db') % 取 db4 的四个滤波器：分解低通、分解高通、重构低通、重构高通 [LoD, HiD, LoR, HiR] = wfilters('db4'); fprintf('db4 分解低通长度 = %d\n', numel(LoD)); fprintf('db4 滤波器系数：\n'); disp(LoD(:)');wavemngr('read')输出的是短名清单，db后面跟数字、bior后面跟两个数字，拼写错一个字符就会在后面的调用里报Invalid wavelet name。waveinfo('db')会直接打印一张表，里面有消失矩阶数、支撑长度、是否正交、是否双正交，选型前扫一眼这张表，比凭印象可靠。wfilters返回的四个向量长度都是2N，它们就是离散小波变换真正在卷积的东西；理解这一点之后，边界延拓为什么会产生额外样本、为什么重构长度可能比原信号长，就都能对上号了。
我一般会在项目开头把这几个查询跑一次，把选定的基名和滤波器长度记在脚本注释里。半年后回来改参数，看见注释就知道当时为什么选它，省掉一轮重新试错。
2.3 分解层数的两条硬约束：长度上限与物理分辨率
层数不是拍脑袋定的，它同时被两个东西卡住。
第一条是长度约束。工具箱提供wmaxlev(N, wname)算理论最大层数，N 是信号长度。它保证在不超过这个层数时，分解不会因为系数长度不够而出问题。常见做法是先跑wmaxlev，再往下取一到两层，给边界延拓留余量。
第二条是物理约束，这条更关键，也更常被忽略。第 j 层细节系数对应的频带是[fs/2^(j+1), fs/2^j]，这是二进抽取的必然结果。如果关心的冲击主频在 400 Hz、采样率 2000 Hz，那主战场是第 2 层（250～500 Hz），不是第 4 层也不是第 5 层。层数越高，细节系数的时域长度越短，时间分辨率越低；一个 20 ms 的冲击如果被摊到一层跨度 80 ms 的系数上，峰值就被平均掉了。
N = 8000; wname = 'db4'; levMax = wmaxlev(N, wname); fprintf('理论最大层数 = %d\n', levMax); fs = 2000; for j = 1:min(levMax, 8) fprintf('D%d 对应频带 [%6.1f, %6.1f] Hz\n', j, fs/2^(j+1), fs/2^j); end这段代码先算出长度允许的上限，再把每一层的频带打印出来。参数说明：wmaxlev的第一个入参是样本数，不是时间长度；第二个入参必须和后面wavedec用的基一致，换基就要重算。循环里用fs/2^(j+1)和fs/2^j是因为每层做一次二抽取，奈奎斯特带宽逐层减半。打印出来之后，对着自己的关注频段挑层数，比默认写 5 层靠谱得多。
2.4 一次可复现的 wavedec 到 waverec 闭环
选完基和层数，先做一次「分解再重构」，确认这条链路本身是通的，再往上叠加去噪和特征提取。这是我一直坚持的顺序：先把无损闭环跑通，再谈有损处理。
rng(7); % 固定随机种子，保证每次结果一致 fs = 2000; N = 8000; t = (0:N-1)'/fs; x = sin(2*pi*50*t) + 0.2*randn(N,1); idx = 5400:5440; % 第 2.7 秒附近 x(idx) = x(idx) + 1.2*sin(2*pi*400*t(idx)); wname = 'db4'; lev = 4; [c, l] = wavedec(x, lev, wname); xr = waverec(c, l, wname); xr = xr(1:N); % 截断到原始长度，抵消边界延拓 fprintf('重构最大误差 = %.3e\n', max(abs(x - xr))); E = zeros(lev,1); for j = 1:lev d = detcoef(c, l, j); E(j) = sum(d.^2); fprintf('D%d 长度 %5d 能量占比 %.2f%%\n', j, numel(d), 100*E(j)/sum(E)); end这段代码的逻辑要拆开说。wavedec返回的c是一个长向量，排列顺序是[A_lev, D_lev, D_{lev-1}, ..., D_1]；l是配套的长度向量，l(1)是近似系数的长度，后面依次是各层细节的长度。waverec用c和l还原，但边界延拓会让输出比原信号多几个样本，所以用xr(1:N)截断。重构最大误差正常应该在1e-12量级，如果到了1e-3以上，先查dwtmode有没有被别的脚本改过。
detcoef(c, l, j)的第 3 个参数是层号，j=1对应最细的那层（原信号最高频段），j=lev对应最粗的细节。能量占比这一列用来判断冲击落在哪层：如果 D2 的能量占比明显高于其他层，说明 400 Hz 的冲击确实被它抓住了，后面做特征提取就从这一层下手。
这套闭环跑完，再去碰wden、wdencmp这些高层函数，心里就有底了。
3. 连续小波变换：把「尺度」换成频率，时频图才读得懂
连续小波变换的输出是一张时频图，横轴时间、纵轴频率、颜色是幅值。这张图看起来直观，读错的人却不少，最常见的错误是把纵轴当成「尺度」直接标上去，或者忘了传采样率，结果频率轴整体偏移一个系数。这一章把cwt的接口语义、尺度到频率的换算、以及两路信号的相干分析讲清楚。
3.1 cwt 的接口换代与参数语义
cwt在近几个版本里换过一次接口，这是很多老脚本迁移时踩坑的地方。旧写法是cwt(x, scales, 'wname')，由你给出尺度向量；新写法是cwt(x, fs)或cwt(x)，工具箱自己选尺度，直接返回频率向量。新接口省事，但代价是默认参数藏在函数里，不显式设就可能和预期不一致。
fs = 1000; t = (0:1/fs:4-1/fs)'; x = chirp(t, 10, 4, 200, 'logarithmic') + 0.3*randn(size(t)); [cfs, f] = cwt(x, fs, 'morl'); figure; imagesc(t, f, abs(cfs)); axis xy; % 让低频在下、高频在上 xlabel('时间 (s)'); ylabel('频率 (Hz)'); title('Morlet 连续小波变换幅值'); colorbar;这段代码里，cfs是复数矩阵，行对应频率、列对应时间，abs(cfs)取幅值；f是和行一一对应的频率向量。imagesc默认把第一行画在顶部，必须加axis xy翻转，否则低频跑到上面去，读图时容易误判。'morl'是 Morlet 小波，时频聚焦性好，适合看振荡成分；如果信号里是尖锐脉冲，'mexh'（墨西哥帽）的时域支撑更短，定位更锐。
两个常用参数靠 name-value 对传进去：VoicesPerOctave控制每个倍频程里的频率点数，默认 10，调到 16 或 24 图会更细腻，但计算量按比例涨；FrequencyLimits限定分析频带，长信号上能省掉大量无用计算。我一般先扫一遍全频带看整体，再限频带出细节图。
提示：
cwt的默认小波在新版本里已经不是'morl'，显式写出来更稳妥，别依赖默认值。
3.2 尺度到频率的换算与频率轴校准
旧接口只给你尺度，要自己换算成频率。公式是f = fc / (a * dt)，其中a是尺度，dt = 1/fs是采样间隔，fc是小波的中心频率。Morlet 在单位采样率下的中心频率约为 0.8125，所以尺度 1、采样率 1000 Hz 对应约 812.5 Hz。
scales = 2.^(0:1/12:7); % 每倍频程 12 个尺度 freq = scal2frq(scales, 'morl', 1/fs); disp(table(scales(:), freq(:), 'VariableNames', {'尺度','频率Hz'}));scal2frq就是干这个换算的，第三个参数是采样间隔1/fs，不是fs，这里写错频率轴会整体差一个数量级。打印出来的表可以直接对：尺度翻倍，频率减半，这是二进关系。如果自己手写公式，fc一定要查清楚，morl、mexh、cmor的中心频率各不相同，混用会让整张图偏移。
新接口下这一步被封装了，但理解换算关系仍然必要。比如你看到时频图上某个亮斑落在 120 Hz，想回去查它在哪一层离散细节里，就得靠这个关系把频率映射回层号。
3.3 用 wcoherence 看两路信号的相位关系
单路信号的时频图只能看「有什么频率、什么时候出现」，看不出两路信号之间的关系。想知道振动和电流在故障频段上是否同步，用wcoherence。
y = x + 0.5*randn(size(x)); % 构造一路相关信号 [cxy, f] = wcoherence(x, y, fs); figure; imagesc(t, f, abs(cxy)); axis xy; colorbar; xlabel('时间 (s)'); ylabel('频率 (Hz)'); title('小波相干幅值');cxy是复数，幅值在 0 到 1 之间，越接近 1 表示两路信号在该时刻该频率上越相干；相位可以用angle(cxy)取出来，用来看超前滞后。参数上，采样率同样必须给对，否则频率轴又偏。相干分析对噪声敏感，输入信号最好先去趋势、去均值，否则低频段会出现一片假的高相干区，那是直流分量在做功。
3.4 长信号的分段处理与内存控制
cwt的复杂度大致是O(N * 频率点数)，N 到几十万点时内存会明显吃紧，再叠上VoicesPerOctave调高，很容易把机器拖死。工程上的常见做法是分段。
segLen = 2*fs; % 每段 2 秒 ovl = round(0.5*segLen); % 50% 重叠，避免段边界丢事件 starts = 1:ovl:(numel(x)-segLen); for k = 1:numel(starts) seg = x(starts(k):starts(k)+segLen-1); seg = seg - mean(seg); % 去均值，压掉边界低频大值 [cfs, f] = cwt(seg, fs, 'morl'); % 这里按需拼接、取峰值或统计各频带能量 end分段长度按关注事件的持续时间定：事件 20 ms，段长 1～2 秒足够；段太长反而浪费。重叠率取 50% 是常见折中，能保证跨段边界的事件至少被一段完整覆盖。每段先减均值这一步别省，不去均值的话，段首段尾会出现一大片低频高幅值，看起来像故障特征，其实是边界效应。拼接时只取每段的中间部分，边缘各丢ovl/2个点，图会更干净。
4. 小波去噪：阈值规则、软硬阈值、噪声估计的三参数联动
去噪是小波工具箱里用得最多的功能，也是最容易「跑出结果但结果不对」的地方。三个参数互相牵制：阈值规则决定阈值多大，软硬阈值决定怎么削，噪声估计决定 sigma 从哪来。三个都调对了，去噪才站得住。
4.1 三参数模型：规则、方式、噪声估计
阈值规则有四条常用路径。sqtwolog用sqrt(2*log(n))乘噪声标准差，是通用阈值；信号越长，这个阈值越大，所以长信号上容易过平滑。rigrsure基于 Stein 无偏风险估计，逐点算风险再取最小，对稀疏信号更友好。heursure是前两者的启发式组合，信噪比低时偏向sqtwolog，高时偏向rigrsure。minimaxi用极大极小准则，适合信噪比本身就不高、又不想过度削峰的场景。
阈值作用方式分软硬两种。软阈值's'把所有系数往零的方向收缩，输出连续、光滑，但幅值会被系统性压低；硬阈值'h'只把小系数置零，大系数原样保留，保幅值但输出可能有跳变。工程上默认软阈值居多，需要保冲击幅值时改用硬阈值或折中的惩罚阈值。
噪声估计分'sln'和'mln'。'sln'只用第一层细节系数估 sigma，假设噪声主要落在最高频段；'mln'用多层细节联合估计，噪声分布较宽时更稳。
sqtwolog的阈值和信号长度直接相关，这一点在实测数据上经常被忽略。同样噪声水平，N 从 4096 涨到 65536，阈值会涨一截，原本干净的细节被削掉，看起来像「去噪太狠」。遇到这种情况，换heursure或rigrsure，比死磕sqtwolog有效。
4.2 wden 与 wdencmp 两条路
工具箱提供两条去噪路径：wden是一键式，内部完成分解、阈值、重构；wdencmp把阈值和方式交给你，可控性更强。
rng(11); fs = 1000; N = 4096; t = (0:N-1)'/fs; xc = sin(2*pi*30*t) + 0.5*sin(2*pi*120*t); % 干净信号 xn = xc + 0.6*randn(N,1); % 含噪信号 lev = 4; wname = 'db4'; xd1 = wden(xn, 'sqtwolog', 's', 'sln', lev, wname); xd2 = wden(xn, 'heursure', 's', 'sln', lev, wname); xd3 = wden(xn, 'minimaxi', 'h', 'sln', lev, wname);wden的入参顺序是：信号、阈值规则、软硬方式、噪声估计、层数、小波名。它返回的就是重构好的去噪信号，中间系数拿不到。想分析各层细节被削了多少，就得换wdencmp。
[thr, sorh, keepapp] = ddencmp('den', 'wv', xn); xd4 = wdencmp('gbl', xn, wname, lev, thr, sorh, keepapp); fprintf('ddencmp 给出阈值 %.4f，方式 %s，保留近似 %d\n', thr, sorh, keepapp);ddencmp按信号长度和噪声水平自动给一组默认参数，'den'表示去噪、'wv'表示小波。wdencmp的第一个参数'gbl'表示全局统一阈值，换成'lvd'就是逐层不同阈值，后者对噪声在各层分布不均的信号更合适。keepapp为 1 表示近似系数不参与阈值处理，保留低频趋势；设成 0 会把趋势一起削，除非你本来就想滤掉低频。
我自己做实测数据时更倾向wdencmp加'lvd'，因为现场噪声很少是严格白噪声，各层细节的能量分布不均匀，统一阈值要么削不干净要么削过头。
4.3 用 SNR 和 RMSE 做客观评估
去噪效果不能只看图。仿真数据有干净信号作参考，直接算信噪比和均方根误差，比肉眼判断可靠得多。
snr = @(a,b) 10*log10(sum(a.^2)/sum((a-b).^2)); rmse = @(a,b) sqrt(mean((a-b).^2)); fprintf('含噪信号 SNR = %6.2f dB, RMSE = %.4f\n', snr(xc,xn), rmse(xc,xn)); fprintf('sqtwolog SNR = %6.2f dB, RMSE = %.4f\n', snr(xc,xd1), rmse(xc,xd1)); fprintf('heursure SNR = %6.2f dB, RMSE = %.4f\n', snr(xc,xd2), rmse(xc,xd2)); fprintf('minimaxi SNR = %6.2f dB, RMSE = %.4f\n', snr(xc,xd3), rmse(xc,xd3));匿名函数snr的分子是干净信号能量，分母是误差能量，单位是 dB；rmse直接反映误差绝对量级。评估时两个指标一起看：SNR 涨了但 RMSE 也涨，说明信号被整体缩放，软阈值的系统性压低在作怪；SNR 降了，说明阈值过大，把信号成分当噪声削了。
实测数据没有干净参考，这套指标用不了。替代做法是看细节层能量的变化：去噪前后各层细节能量的比值，如果在某一层突然掉到 0.1 以下，那一层大概率被削过头了。
4.4 小波包：当高频段也值得继续拆
普通 DWT 只对低频近似系数继续分解，高频细节不再动。如果故障特征恰好落在高频段、又需要更细的频率分辨，就该上小波包。
T = wpdec(xn, 3, 'db4'); % 三层小波包分解 T = wdencmp('gbl', T, 'db4', 3, thr, 's', 1); % 对整个包树做阈值 xw = wprec(T); % 重构 fprintf('小波包 SNR = %6. <p> <a href="https://download.csdn.net/download/zhuan_shu/10979411" style="color:#ec7500;font-size:14px;"> 本文还有配套的精品资源，点击获取 </a> <img alt="menu-r.4af5f7ec.gif" src="https://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif" style="width:16px;margin-left:4px;vertical-align:text-bottom;cursor:text;"> </p>
