简介Matlab小波工具箱配套资料面向科研人员、工程师及学习小波分析的学生围绕信号处理、图像分析、去噪、特征提取与模式识别等任务整理了从离散小波变换、多分辨率分析到连续小波变换、MODWT及小波包分析的实践素材。压缩包共1152个文件约17.76MB以792个m脚本/函数、227个mat数据、56个png与35个fig/30个jpg图示为主另有xml、cfg、bin等配置与数据文件便于按主题查找示例并复现实验。已有10737人学习下载说明其在课程练习与工程验证中有较高参考价值。读者可借助其中的函数示例、图形界面文件和测试数据理解DWT、MODWT、wavedec/waverec、cwt及二维小波分解等用法并用于去噪、统计检验和图像处理的对照学习。1. 同一段振动信号两个人跑 wden 差了 20 dB先看清工具箱的五个默认值一段 2 秒、采样率 1 kHz 的振动信号两个人拿同一台设备的采集数据都用 Matlab 小波工具箱里的wden去噪。一个人出来是干净的渐变频另一个人满屏毛刺算出来的信噪比差 20 dB 以上。问题不在数据而在wden的五个入参——阈值规则、软硬阈值、噪声缩放方式、分解层数、小波基每个都有默认值默认值之间又互相耦合随便改一个结果就变一次很多人在这里就觉得小波去噪是门玄学。Matlab 小波工具箱Wavelet Toolbox本质上是一套围绕多分辨率分析搭起来的函数库和交互界面一维/二维离散小波变换、小波包、连续小波变换、去噪、压缩、图像融合、波形特征提取都在里面。它适合手里已经有实测信号、需要快速验证去噪或特征提取方案的从业者也适合刚接触小波变换、想先看到效果再回头补理论的人。这篇按「分解重构怎么算账 → 入口怎么选 → 参数怎么落 → 坑在哪」的顺序拆一遍每一步都给能直接跑的代码。2. 分解重构这笔账滤波器组、层数上限与边界延拓小波变换在工程上落地绕不开三个决策用哪个小波基分几层信号两端怎么补。这三个决策不是各自独立的小波基决定了滤波器的长度滤波器长度又反过来限制了你最多能分几层边界延拓模式则决定了你能不能在重构后原样还原信号长度。下面按这三件事依次拆。2.1 从滤波器组看 dwt 与 wavedec 的差别dwt是单层分解输入一段信号输出近似系数cA和细节系数cD两者长度都约等于输入的一半。wavedec是 N 层分解它把每一层的近似系数继续往下分最后返回一个拼接向量C和一个长度向量L。很多人第一次用wavedec会卡在C和L上因为C不是矩阵是一维拼接得靠appcoef、detcoef按层索引去取。% 构造一段测试信号低频缓变 高频扫频 fs 1000; t (0:1/fs:2-1/fs); % 2 秒2000 个采样点 f1 20 8*sin(2*pi*0.5*t); % 20~28 Hz 缓变分量 f2 120 30*t; % 120~180 Hz 线性扫频 clean 0.7*sin(2*pi*cumsum(f1)/fs) 0.4*sin(2*pi*cumsum(f2)/fs); rng(7); % 固定随机种子保证结果可复现 noisy clean 0.15*randn(size(clean)); % 单层分解看清 cA / cD 的长度关系 [cA1, cD1] dwt(noisy, sym8); fprintf(原信号 %d 点cA1 %d 点cD1 %d 点\n, numel(noisy), numel(cA1), numel(cD1)); % 三层分解[C, L] 是拼接结构不是矩阵 [C, L] wavedec(noisy, 3, sym8); A3 appcoef(C, L, sym8, 3); % 第 3 层近似系数 D3 detcoef(C, L, 3); % 第 3 层细节系数 D1 detcoef(C, L, 1); % 第 1 层细节系数 disp(L); % L 记录了每一段的长度累加即 C 的总长上面这段里L的顺序是从高层往低层排的L(1)是第 3 层近似系数的长度L(2)是第 3 层细节依此类推到L(end)是原信号长度。appcoef的第 4 个参数写层号detcoef的第 3 个参数写层号层号写错不会报错只会取到错的频段这是最隐蔽的一类错误。wrcoef可以从单层系数直接重构回原长度waverec则是从C、L整体重构。两者在去噪里分工很明确如果你只改某一层系数用wrcoef重构该层再看成分如果所有系数都动过用waverec一次性还原。% 只重构第 1 层细节看高频成分长什么样 D1_full wrcoef(d, C, L, sym8, 1); % 整体重构验证无损性没改任何系数时应与 noisy 完全一致 rec waverec(C, L, sym8); err max(abs(rec - noisy)); fprintf(无损重构最大误差%g\n, err);err应该接近机器精度量级1e-12 左右。如果这个值明显不为零说明中间动过C或者延拓模式被改过——这是判断「我的重构流程是不是干净」最便宜的一个自检手段我一般在新环境里第一件事就是跑一遍这个。2.2 分解层数不是 fix(log2(N)) 就完事网上流传最广的一句话是「层数取fix(log2(N))」但这个公式在小波场景里只对 Haar 基近似成立。真实约束是每一层分解后数据长度要大于该小波对应的滤波器长度否则wavedec直接报错。工具箱给了wmaxlev来算这个上限N numel(noisy); for w {haar,db4,sym8,coif5,bior3.5} lev_max wmaxlev(N, w{1}); fprintf(%-8s 允许的最大层数%d\n, w{1}, lev_max); end跑出来你会发现不同小波基允许的层数不一样coif5滤波器长 30能分的层数明显少于haar。这里有个参数经验层数过多时最高层近似系数的长度可能只剩几十个点这一层反映的是信号的趋势而不是细节去噪时把这一层也拿去做阈值处理容易把信号的整体基线削掉一块。一层分解对应的频带是[0, fs/4]和[fs/4, fs/2]第 j 层细节对应[fs/2^(j1), fs/2^j]。反过来推如果你想保留 20 Hz 以上的成分采样率 1 kHz那么第 4 层细节的下边界是 1000/32 ≈ 31 Hz第 5 层是 15.6 Hz所以至少要分到 5 层才能覆盖到 20 Hz。这个反推比死记log2(N)有用得多。层数 j细节频带fs1000 Hz系数长度N2000sym81250 ~ 500 Hz约 10072125 ~ 250 Hz约 510362.5 ~ 125 Hz约 262431.2 ~ 62.5 Hz约 138515.6 ~ 31.2 Hz约 7667.8 ~ 15.6 Hz约 45表里的系数长度是近似值实际长度受延拓模式影响会差一两个点。第 6 层只剩 45 个点这一层做阈值处理意义已经不大除非你的目标就是提取极低频趋势。2.3 边界延拓两端那几十个点为什么总是歪的离散小波变换在每个尺度上都要做下采样信号长度不能被 2 整除时就必须在两端补数据。补法由dwtmode控制这个设置是全局状态一旦改过后面所有dwt、wavedec都跟着变而它不会在函数签名里体现出来属于典型的「黑匣子」式干扰。old_mode dwtmode(status, plt); % 查看当前模式 modes {sym,per,zpd,spd}; for k 1:numel(modes) dwtmode(modes{k}, nodisp); [C, L] wavedec(noisy, 4, sym8); rec waverec(C, L, sym8); n min(numel(rec), numel(noisy)); e max(abs(rec(1:n) - noisy(1:n))); fprintf(%-4s 重构误差 %g重构长度 %d\n, modes{k}, e, numel(rec)); end dwtmode(old_mode, nodisp); % 还原别污染后续脚本四种模式的实际表现差别很大sym对称延拓最常用边界畸变小但重构长度可能和原信号差一个点。per周期延拓重构长度严格等于原长度前提是信号长度是 2 的幂否则要先补零补零那一段会在两端引入跳变。zpd补零实现最简单但对非零均值的信号两端会出现明显台阶。spd一阶平滑延拓对趋势明显的信号比如温度、位移比sym更稳。提示dwtmode是全局状态函数里改了不还原会污染同一 session 里其他脚本。养成「进函数先存、出函数先还原」的习惯比事后排查省事得多。到这里分解重构的三个决策就算有了可执行的判断依据小波基按信号平滑度选sym8、coif5适合连续信号haar适合突变层数按目标频带反推延拓模式按信号两端是否有强趋势选。这三件事定下来后面去噪参数的调整才有意义否则你调阈值调半天改的其实是边界效应。3. GUI 点出来的结果要能复现四个入口、参数映射与批量化工具箱最容易被低估的一点是它有四个完全不同的入口而且这四个入口的参数并不总是一一对应。很多人先在图形界面里试出一个满意结果回头写脚本怎么都复现不出来最后只能截图交差。3.1 函数族、对象族与交互界面各管什么工具箱的入口可以分成四类第一类是过程式函数dwt、wavedec、wden、wdencmp、wpdec、cwt这些输入输出都是数组写脚本、做批处理首选。第二类是较新版本引入的对象式接口比如一维去噪的wdenoise参数用名值对传DenoisingMethod、ThresholdRule、NoiseEstimateMethod、Wavelet、DenoisingLevel。它比wden好读但底层阈值处理和wden不完全等价两者结果对不上是正常的不要拿一个的结果去校验另一个。第三类是交互界面wavemenu较新版本叫waveletAnalyzer一维信号分析、二维图像分析、小波包分析、连续小波分析都在里面。它的价值不在「点着看」而在「点完能导出脚本」。第四类是专用工具箱函数图像去噪、融合、边缘检测那一套接口更封闭适合快速出图不适合细调。选型上我的习惯是先用waveletAnalyzer三十秒看个大概确认小波基和层数的量级然后把界面里的参数抄进wden或wavedec脚本之后所有调整都在脚本里做。界面只用来看不用来交付。3.2 wavemenu 里的参数怎么还原成脚本waveletAnalyzer的一维去噪界面里有一个「File → Generate M-file」的选项能把当前界面上的所有设置导出成一个可运行的函数。这个功能被严重低估了它是解决「界面效果和脚本效果不一致」最直接的办法。导出的文件大概长这样节选函数名和变量名会因版本不同有差异function varargout wdenComp1D(X, varargin) % 由交互界面导出的去噪函数骨架 % X 为输入信号varargin 允许覆盖界面里的设置 % 界面里的设置被固化成了默认值 wname sym8; % 小波基 level 5; % 分解层数 tptr sqtwolog; % 阈值规则 sorh s; % 软阈值 scal sln; % 单层噪声缩放 if nargin 1 % 实际导出代码里是一长串名值对解析这里简化成示例 for k 1:2:numel(varargin) switch lower(varargin{k}) case wavelet, wname varargin{k1}; case level, level varargin{k1}; case tptr, tptr varargin{k1}; end end end xd wden(X, tptr, sorh, scal, level, wname); varargout{1} xd; end这段代码的重点不在实现而在它把界面里那五个我前面说的「隐藏默认值」全部显式写出来了。你可以直接拿它和wden(X,sqtwolog,s,sln,5,sym8)对拍逐项确认是哪一个参数不一样导致结果有差异。这一步做完「界面好用脚本难用」这个问题基本就消掉了。需要注意导出的文件里scal常常是sln单层噪声缩放而很多教程里写的是one不缩放这两者在噪声水平随尺度变化时会给出明显不同的阈值。这是第一类翻车点后面避坑章还会提。3.3 把一次分析固化成批量函数单次调通了下一步就是批量。批量处理的坑不在算法在数据组织文件命名不统一、每条信号长度不一样、采样率写在文件名里、有空文件。下面这个骨架是我常用的形态去掉了业务细节只留结构。function batchDenoise(inDir, outDir, cfg) % 批量一维去噪读目录下所有 .mat 或 .csv逐条去噪后写回 % cfg 字段wname / level / tptr / sorh / scal / fs if ~exist(outDir, dir), mkdir(outDir); end files dir(fullfile(inDir, *.mat)); report struct(name, {}, snr_in, {}, snr_out, {}, level, {}); dwtmode(sym, nodisp); % 统一延拓模式避免逐文件不一致 for i 1:numel(files) S load(fullfile(inDir, files(i).name)); fn fieldnames(S); x S.(fn{1}); % 取第一个变量命名不统一时兜底 x double(x(:)); % 强制列向量防止行向量导致维度错 % 层数按信号长度和配置取小值防止短信号报错 lev min(cfg.level, wmaxlev(numel(x), cfg.wname)); if lev 1 warning(%s 太短跳过, files(i).name); continue; end xd wden(x, cfg.tptr, cfg.sorh, cfg.scal, lev, cfg.wname); % 结构体保存带上处理参数方便回溯 meta struct(wavelet, cfg.wname, level, lev, ... tptr, cfg.tptr, sorh, cfg.sorh, scal, cfg.scal); save(fullfile(outDir, files(i).name), xd, meta); report(end1) struct(name, files(i).name, snr_in, NaN, ... snr_out, NaN, level, lev); %#okAGROW end save(fullfile(outDir, report.mat), report); end几个参数和写法的用意x(:)强制列向量是因为wden对行向量和列向量的处理一致但后续拼接容易出错wmaxlev兜底是因为批量数据里总有几条短信号把meta和处理结果一起存是为了三个月后回头看还能知道当时用的是哪套参数——这是花十行代码买的一份后悔药。调用方式cfg struct(wname,sym8,level,5,tptr,sqtwolog, ... sorh,s,scal,sln,fs,1000); batchDenoise(./raw, ./denoised, cfg);这套骨架替换掉wden那一行就能换成小波包、换成wdenoise对象接口、甚至换成二维图像处理结构不用动。批量化真正的价值不是省时间是让每次处理都留下可复现的参数记录。4. 去噪与压缩的参数落地阈值规则、软硬阈值与两个回归指标参数调整是这套工具箱里最容易陷入「看着差不多就行」的环节。要跳出这个状态得先知道每个参数在数学上到底做了什么再用指标把「差不多」量化。4.1 四类阈值规则分别在什么信号上会翻车wden的tptr参数支持四种阈值选取规则它们对噪声的假设完全不同规则阈值来源适用前提常见翻车场景sqtwolog通用阈值 σ√(2lnN)噪声近似白噪声样本量大信号本身含大量尖峰时尖峰被一起削掉rigrsureStein 无偏风险估计信号稀疏噪声占比不高信号连续且能量分散时阈值偏低噪声残留heursure前两者取折中信噪比未知通用场景N 很小时折中策略失效minimaxi最小最大准则关心最坏情况误差追求平均误差最优时反而偏差大rigrsure和heursure需要估计噪声标准差sqtwolog和minimaxi不需要。这一条决定了它们在短信号上的表现信号只有几百个点时rigrsure的方差估计很不稳结果会跳。sorh只有s和h两个值但影响比想象中大。软阈值把小于阈值的系数直接置零、大于阈值的系数整体往零的方向平移一个阈值量结果是波形更平滑但有系统性收缩幅值会偏小硬阈值只把小于阈值的置零、其余原样保留幅值保住了但会引入振铃。工程上如果后续要做幅值定标比如振动烈度、应变峰值硬阈值更合适如果后续要做频谱分析或者包络软阈值更稳。scal参数是最容易被忽略的一个one不缩放sln按第一层细节系数估的噪声做单层缩放mln做多层独立缩放。噪声在每一层的投影能量本来就不一样mln才是理论上更合理的做法但它对噪声非平稳的信号会过度压制某一层。我在实测数据上一般先用sln打底再对比mln。4.2 一个能直接抄的去噪脚本把前面所有决策串起来下面这段是可以直接换数据跑的完整流程包含参数扫描和指标输出%% 1. 数据准备 fs 1000; t (0:1/fs:2-1/fs); clean 0.7*sin(2*pi*cumsum(208*sin(2*pi*0.5*t))/fs) ... 0.4*sin(2*pi*cumsum(12030*t)/fs); rng(7); noisy clean 0.15*randn(size(clean)); %% 2. 参数网格小波基 × 层数 × 阈值规则 wnames {db4,sym8,coif5}; tptrs {sqtwolog,rigrsure,heursure,minimaxi}; levels 3:6; best struct(snr, -Inf); dwtmode(sym,nodisp); % 统一边界处理保证可比 for wi 1:numel(wnames) for ti 1:numel(tptrs) for li 1:numel(levels) lev min(levels(li), wmaxlev(numel(noisy), wnames{wi})); xd wden(noisy, tptrs{ti}, s, sln, lev, wnames{wi}); % 信噪比按「信号能量 / 残差能量」定义残差 去噪结果 - 真值 snr_out 10*log10(sum(clean.^2) / sum((xd-clean).^2)); if snr_out best.snr best struct(snr, snr_out, wname, wnames{wi}, ... tptr, tptrs{ti}, level, lev); end end end end fprintf(最佳组合%s / %s / %d 层SNR %.2f dB\n, ... best.wname, best.tptr, best.level, best.snr); %% 3. 用最佳参数出图并保存 xd wden(noisy, best.tptr, s, sln, best.level, best.wname); figure(Position,[100 100 900 600]); subplot(3,1,1); plot(t, clean); title(原始信号); ylabel(幅值); subplot(3,1,2); plot(t, noisy); title(含噪信号); ylabel(幅值); subplot(3,1,3); plot(t, xd); title(去噪结果); ylabel(幅值); xlabel(时间 / s);参数说明上wmaxlev的兜底逻辑必须写网格搜索里层数一旦超过上限整个循环就中断snr_out用的是「真值能量 / 残差能量」注意分母是去噪结果减真值不是减含噪信号——后者算出来的只是「改动了多少」不是「改对了多少」这两个指标经常被混用是个不大不小的坑。网格搜索的代价也要说清楚3 个基 × 4 个规则 × 4 个层数 48 次wden调用2 秒 2000 点的信号在普通笔记本上不到两秒跑完可以接受。但如果信号是几万点、基的数量上到十几种这种暴力搜索就不合适了得改成先按频带反推层数再在小范围内搜规则。4.3 用指标把「看着不错」变成可回归的数字去噪效果的评价不能靠眼睛。下面三个指标我每次都会一起看% 指标 1输出信噪比相对真值 snr_out 10*log10(sum(clean.^2)/sum((xd-clean).^2)); % 指标 2均方根误差 rmse sqrt(mean((xd-clean).^2)); % 指标 3残差与原始噪声的相关性越接近 0 说明去噪越「只动噪声」 resid xd - clean; noise_only noisy - clean; corr_resid corrcoef(resid, noise_only); fprintf(SNR %.2f dB, RMSE %.4f, 残差相关 %.3f\n, ... snr_out, rmse, corr_resid(1,2));第三个指标最容易被忽略但它是最能暴露问题的如果残差和噪声高度相关相关性接近 1说明去噪基本没动信号如果相关性接近 0 但 SNR 提升很小说明去噪把信号也削了一块两边抵消了。正常情况应该是 SNR 提升明显、残差与噪声相关性在 0.3 以下。压缩场景用的是另一套指标wdencmp会直接返回两个% 全局阈值压缩THR 是阈值s 软阈值 [C, L] wavedec(noisy, 5, sym8); [xc, perf0, perfl2] wdencmp(gbl, C, L, sym8, 5, 0.3, s); fprintf(保留系数 %.2f%%保留能量 %.2f%%\n, perf0, perfl2);perf0是保留下来的非零系数占比perfl2是保留的系数能量占原能量的比例。压缩比看perf0失真程度看perfl2和重构误差的配合——只看压缩比不看能量保留率是最典型的「压过头」写法。5. 避坑清单五条让我返工的记录从报错到两端畸变下面这五条都是从实际返工里攒出来的现象、原因、解决按顺序写能对上你正在遇到的问题最好。5.1 输入与维度相关的三条第一条wavedec报层数超限。现象是脚本在批量跑的时候绝大多数文件正常个别文件直接抛错退出循环。原因很简单那几条信号的长度不够分到配置的层数wmaxlev的上限被越过了。解决办法是在循环里用lev min(cfg.level, wmaxlev(numel(x), wname))压一层并且对lev 1的短信号直接跳过并打 warning。批量脚本里任何和长度相关的参数都必须做下限保护这一条不止适用于小波。第二条重构长度和原信号差一两个点导致后面按索引对齐的代码全错位。现象是waverec出来的向量比原信号短 1直接用xd - clean算误差时报维度不匹配。原因是延拓模式和信号长度不匹配per模式要求长度是 2 的幂不是的话内部补零补的位置和你想的不一样。解决是统一用dwtmode(sym)并在重构后按min(numel(xd), numel(x))截断再做对比。如果想严格控长度就先x [x; zeros(2^nextpow2(numel(x)) - numel(x), 1)]补到 2 的幂再按原始长度截回来。第三条行向量和列向量混用。现象是单条跑没问题多条拼成矩阵后wden的结果变了形状。原因是工具箱多数函数对一维输入的方向敏感拼接时[a; b]和[a, b]出来的是完全不同的东西。解决是在所有入口统一x double(x(:))输出也强制列向量别让形状沿数据流传播。5.2 重构与评估相关的两条第四条去噪后波形整体幅值偏小。现象是峰值比原始信号低了 10% 以上频谱却没变干净多少。原因是用了软阈值软阈值对超过阈值的系数做了平移所有大幅值系数被系统性地拉小层数越多累积收缩越明显。解决是改sorh h用硬阈值或者在软阈值之后做一次幅值补偿用去噪前后低频近似的能量比做缩放。如果下游要做幅值定标硬阈值基本是唯一选择。第五条两端各十几个点严重畸变中间完全正常。现象是去噪结果在首尾出现明显下沉或上翘越靠近端点越夸张。原因是边界延拓和阈值处理在端点附近互相放大延拓出来的数据本身不是真实信号阈值判断在这些系数上失准小波基支撑越长比如coif5受影响的点数越多。解决办法有三个方向优先用第一个一是换spd或asym延拓模式对带趋势的信号效果明显二是两端各裁掉滤波器长度 - 1个点再分析把端点效应移出关注区三是降低层数层数越少边界影响扩散得越浅。6. 进阶技巧小波包挑频带与 CWT 时频图的两个参数小波包比普通小波多一层自由度普通小波只对近似系数继续分解高频段永远保持粗分辨率小波包对细节系数也继续分解等于把整个频带均匀切细。判断标准很直接——如果你的目标频段落在高频区比如采样率 1 kHz 下想看清 200~400 Hz 的成分普通小波分到第 2 层就停了小波包还能往下切。wpdencmp是最省事的入口参数顺序是信号、软硬阈值、层数、小波基、准则、准则参数、是否保留近似% 小波包去噪3 层threshold 准则下 PAR 就是阈值 [xdp, treed, perf0, perf2] wpdencmp(noisy, s, 3, sym8, threshold, 0.35, 1);第七个参数写1表示保留近似系数不参与阈值处理这是去噪和压缩的关键区别去噪必须保留近似否则信号的整体趋势会被一起削掉压缩可以设0因为趋势本身也在被压缩的范围内。treed返回的是分解树结构可以拿wptree的相关函数去看每个节点被保留还是被裁掉用来确认「频带挑对了没有」。CWT 那边的问题集中在接口版本上。老接口要显式给尺度序列和小波基新接口可以直接传采样率出时频图。老接口的好处是尺度序列完全可控做定量分析时更稳scales 1:128; cw cwt(noisy, scales, sym8); % 复小波系数矩阵 f scal2frq(scales, sym8, 1/fs); % 尺度转实际频率 figure; contour(t, f, abs(cw), 20); % 用等值线看能量脊 set(gca,YScale,log); ylabel(频率 / Hz); xlabel(时间 / s);scal2frq的第三个参数是采样周期的倒数也就是1/fs这里最容易写错成fs写反了频率轴会整体偏掉几个数量级图还是能画出来但读数全错。contour用 20 条等值线是我试下来的经验值太少看不出脊线走向太多整张图糊成一团。如果信号是非平稳的尺度序列不要取等间隔改成2.^(0:0.25:7)这种对数序列低频分辨率会好很多。从那以后我每次碰小波不管是去噪还是时频分析都强制先跑一遍「不改系数直接重构」的误差自检再加一行dwtmode的显式设置最后才动参数。这三步加起来不到十行代码省掉的返工远不止十次。希望帮到你。本文还有配套的精品资源点击获取