简介:Matlab小波工具箱(Wavelet Toolbox)配套学习资料包,面向从事信号处理、图像分析与数据建模的科研人员和工程师,也适合系统学习小波分析的研究生;内容覆盖离散小波变换、MODWT、多分辨率分析、连续小波变换与小波包分析,涉及去噪、压缩、特征提取、金融时间序列与地震信号处理等场景。压缩包共1152个文件、约17.76MB,以792个.m脚本和227个.mat数据为主体,辅以56张png、35个fig演示图与30张jpg截图,另有少量xml、mexw64等配置与加速文件;脚本对应dwt、modwt、wavedec、cwt等函数的实际调用,mat数据可直接加载验证,fig则保留小波分解界面与波形结果,构成可运行、可对照的示例集合。已有10737人学习下载,可据此梳理工具箱的函数体系,按变换类型对照源码与结果图复现实验,并参考重构误差、边界延拓、尺度选择等问题的实现思路。
1. 小波工具箱不是“一键去噪”:先把信号拆开看
Matlab 小波工具箱(Wavelet Toolbox)在很多人手里就等同于一个wden:一段带噪信号丢进去,出来一条光滑曲线,曲线看着挺美,一放到实测数据上就露馅——阶跃被抹圆、脉冲被削平、两端冒出一串衰减振荡。问题不在工具箱本身,而是小波变换的底层逻辑是“按尺度把信号拆成近似系数和细节系数”,拆得合不合理,取决于小波基、分解层数、边界延拓方式、阈值规则这四个参数,缺一个都会翻车。这份工具箱真正有价值的部分,是它把离散多分辨分析、连续小波变换、小波包、提升小波、小波散射这几条链路都做成了可调用函数,还配了 App 让你先把参数试出来。它适合手里已经有振动、心电、语音、图像、地震这类数据,需要做去噪、压缩、突变点定位、时频分析的人;如果只是想看公式推导,那它帮不上什么忙。下面按“入口怎么选 → 参数怎么定 → 一维二维怎么跑 → 哪一步最容易踩坑”的顺序拆。
2. 三个入口与一张参数表:从 wfilters 到 wdenoise 的调用路径
同一个去噪任务,工具箱里至少有三条路能走通:命令行函数族、Wavelet Analyzer 图形界面、以及 2016b 之后逐步铺开的对象化 API。三条路的底层算法是一致的,差别在于可控性、可批量化程度和可复现性。我一般把 App 当成“探路工具”,把命令行当成“生产代码”,把对象化 API 当成“默认参数更聪明的现代写法”,三者混着用,但交付物里只留后两者。
2.1 命令行函数族:wfilters → wavedec → wthresh → waverec
这条链路是小波工具箱的骨架,所有 App 导出的代码、所有高层封装最终都会落到这几个函数上。它的调用顺序固定:先用wfilters拿到滤波器组,用wavedec做多尺度分解得到系数向量c和长度索引l,用wthresh改系数,最后用waverec重构。中间那一步“改系数”才是你真正做算法的地方,其余都是管道。
% 确认工具箱装好、查看当前边界延拓模式 ver('wavelet') dwtmode('status') % 打印当前模式,默认通常是 sym dwtmode('sym') % 显式设置,避免换机器后结果漂移 % 拿 db4 的四个滤波器:分解低通、分解高通、重构低通、重构高通 [LoD, HiD, LoR, HiR] = wfilters('db4'); fprintf('滤波器长度 = %d\n', numel(LoD)); % 一段 2048 点的合成信号,两个频率成分 L = 2048; t = (0:L-1)' / 1000; x = sin(2*pi*0.013*t*1000) + 0.3*sin(2*pi*0.21*t*1000); % 最大可分解层数由信号长度和滤波器长度共同决定 N = wmaxlev(L, 'db4'); [c, l] = wavedec(x, N, 'db4'); fprintf('实际层数 = %d, 系数总长 = %d, l 长度 = %d\n', N, numel(c), numel(l));l这个向量是新手最容易忽略的东西。它的长度是N+2,第一个元素是第 N 层近似系数的长度,中间依次是第 N 层到第 1 层的细节系数长度,最后一个元素是原始信号长度。waverec只认c和l这一对,你手动改了c却没同步改l,重构必报尺寸不匹配。
提示:
wmaxlev算出来的层数不是“推荐层数”,是“物理上限”。它按floor(log2(L/(Lf-1)))估算,Lf是滤波器长度,db4 是 8。真按上限分解,最深层几乎全是边界延拓出来的假数据。
2.2 Wavelet Analyzer App:试参数可以,交付别用它
命令行敲waveletAnalyzer会拉起图形界面,里面按“一维信号 / 二维图像 / 小波包 / 连续变换”分了几个工作区。它的价值在于可视化:你能一边拖分解层数滑杆,一边看各层细节系数的波形,一眼看出哪一层是噪声、哪一层是故障冲击。这种交互在调参阶段效率很高,尤其是不熟悉某类信号的时候,先看一眼系数分布再决定策略,比盲试wden快得多。
问题出在“导出到命令行”这个按钮。导出的代码通常长这样:小波名、层数、阈值全部硬编码,变量名是x、c、l、s这类通用的,粘回工程里第一件事就是改名,第二件事是把写死的阈值抽成参数。更麻烦的是 App 里某些交互操作(比如手动刷选系数)导出的代码根本还原不了,只能还原成“用某个固定阈值处理全部细节”。所以我的用法是:App 里试出小波基和层数,记下来,回到命令行自己重写一遍。
还有一点,App 对超长信号会降采样显示。你看到的那条漂亮的去噪曲线,可能是抽稀过的版本,实际全分辨率数据上边界振荡更明显。
2.3 小波基与分解层数:四个判据 + 一张可照填的表
小波基的选择在小波分析里几乎是玄学级别的存在,但落到工程上其实就四个判据:支撑长度、消失矩、对称性、正交性。支撑长度决定时域局部化能力,越短越能抓住突变,但频域分辨率差;消失矩阶数越高,对多项式型光滑信号的压制越彻底,压缩场景下表现好;对称性影响相位失真,图像处理里非对称小波会产生方向性偏移;正交性保证能量守恒,重构误差理论上为零,但正交小波基本都不对称(Haar 除外)。
| 小波族 | 支撑长度 | 消失矩 | 对称性 | 典型用途 |
|---|---|---|---|---|
| haar / db1 | 2 | 1 | 对称 | 突变检测、教学演示、快速原型 |
| dbN | 2N | N | 不对称 | 通用去噪,N 取 4~8 最常见 |
| symN | 2N | N | 近似对称 | 对相位敏感的信号,减少波形偏移 |
| coifN | 6N | 2N | 近似对称 | 数据压缩,同等误差下系数更稀疏 |
| biorNr.Nd | 可变 | 可变 | 对称 | 图像压缩,双正交可分别设计分析/重构 |
| morl / mexh | — | — | 对称 | 只用于cwt连续变换,不做重构 |
工程上的默认选择:振动、电流这类一维信号先上sym4或db4;图像去噪先上sym4;图像压缩这类需要对称滤波器的场景上bior4.4或bior3.5;只做时频分析不上重构,直接cwt配morl。层数方面,一维信号的经验值是让最深层细节系数的长度不小于滤波器长度两倍,二维图像一般 3 层够用,再多边缘块会碎得没法看。
2.4 边界延拓模式:sym、sp0、ppd 的实际差别
wavedec每次下采样都会让系数长度缩水,信号两端必须补数据才能继续分解,补的方式就是边界延拓。dwtmode提供了七八种,常用的就三个:sym(对称延拓,默认)、sp0(补零)、per(周期化延拓)。补零会在两端造出一个陡峭的跳变,分解后表现为巨大的边界系数;对称延拓相对温和,但如果信号本身首尾不连续,照样有残留;周期化延拓要求信号长度是2^N的整数倍,长度不对会被截断,好处是首尾连续时完全没有边界伪影。
% 同一段信号,三种延拓模式下的重构边界误差对比 L = 1024; x = sin(2*pi*(0:L-1)'/64); [c, l] = wavedec(x, 5, 'db4'); edgeLen = 64; % 只统计首尾各 64 个点的误差 modes = {'sym', 'sp0', 'per'}; for k = 1:numel(modes) dwtmode(modes{k}, 'nodisp'); xr = waverec(c, l, 'db4'); e = max(abs(x(1:edgeLen) - xr(1:edgeLen))); fprintf('%-4s 模式,边缘最大重构误差 = %.3e\n', modes{k}, e); end dwtmode('sym', 'nodisp'); % 改完记得改回来这段代码的用途不是“选最优模式”,而是让你建立直觉:边缘误差跟你选的模式强相关,而多数人去噪失败就是因为没意识到去噪后的“边缘振荡”其实是延拓方式引起的,不是阈值设错了。dwtmode是全局状态,改一次影响后面所有调用,脚本里如果忘了复位,下次跑出来的结果和上次对不上,这种“改了一行无关代码结果就变了”的情况,八成就是这里。
3. 一维去噪全流程:wdenoise 默认参数为什么偏激进
wdenoise是工具箱里最省事的一维去噪函数,一行调用就能出结果。省事的代价是默认参数为“通用场景”调的,对信噪比低、含冲击成分的信号往往偏激进,表现为细节被磨平、峰值被压低。要判断它到底做了什么,最直接的办法是造一段已知真值的信号,把不同参数组合的 SNR 打出来横向比。
3.1 先造一段带 ground truth 的含噪信号
没有真值就没法量化评价,这是去噪调参的第一原则。合成信号里至少放三类成分:一个低频正弦(考察低频保持)、一个高斯型脉冲(考察冲击保持)、一段高频成分(考察噪声压制)。固定随机种子,保证每次跑的噪声完全一致,否则参数对比毫无意义。
rng(7); % 固定种子,噪声可复现 L = 2048; t = (0:L-1)'; % 三个成分:低频基线、高斯脉冲、后段高频 clean = 0.8*sin(2*pi*0.013*t) ... + 1.0*exp(-((t-900)/25).^2) ... + 0.35*sin(2*pi*0.21*t).*(t > 1200); noise = 0.35*randn(L,1); xn = clean + noise; snr = @(ref, est) 20*log10(norm(ref)/norm(ref - est)); fprintf('加噪前 SNR = %.2f dB\n', snr(clean, xn));这个snr匿名函数后面会反复用。注意它用的是二范数比,不是逐点误差,对冲击成分不敏感——如果你的核心诉求是保住那个脉冲峰值,还得额外看峰值衰减率,单看 SNR 会被低频部分的大能量掩盖。
3.2 wdenoise 默认参数实测与四个可调旋钮
wdenoise的参数不少,真正影响结果的就四个:Wavelet(小波基)、DenoisingMethod(阈值选取策略)、ThresholdRule(阈值施加方式)、NoiseEstimate(噪声估计是否随层变化)。默认小波是sym4,默认层数按log2(L)给,会被滤波器长度自动截断,默认阈值策略是通用阈值加软阈值。通用阈值sqrt(2*log(L))在 L 较大时会明显偏大,软阈值又会把所有系数按比例压缩,两者叠加的结果就是“干净但发闷”。
N = wmaxlev(L, 'sym4'); % 全默认 xd_def = wdenoise(xn); % 显式指定:db4 + 通用阈值 + 软阈值 + 层相关噪声估计 xd_a = wdenoise(xn, N, 'Wavelet', 'db4', ... 'DenoisingMethod', 'UniversalThreshold', ... 'ThresholdRule', 'Soft', ... 'NoiseEstimate', 'LevelDependent'); % sym6 + SURE 阈值 + 硬阈值 + 层无关噪声估计 xd_b = wdenoise(xn, N, 'Wavelet', 'sym6', ... 'DenoisingMethod', 'SURE', ... 'ThresholdRule', 'Hard', ... 'NoiseEstimate', 'LevelIndependent'); fprintf('默认 SNR = %.2f dB\n', snr(clean, xd_def)); fprintf('db4+Univ+Soft SNR = %.2f dB\n', snr(clean, xd_a)); fprintf('sym6+SURE+Hard SNR = %.2f dB\n', snr(clean, xd_b));跑下来通常能看到:默认组合的 SNR 不是最高的,但在低噪段最平滑;SURE配硬阈值在冲击附近保持更好,代价是残留一些“毛刺”;LevelDependent在噪声强度随尺度变化时(比如传感器带宽受限)优势明显,白噪声场景下反而和LevelIndependent差不多。这些不是绝对结论,但能建立一个判断顺序:先看残差里有没有残留的结构性成分,再看脉冲有没有被压,最后才比 SNR 数字。
四个旋钮之外还有Boundary可以切周期/反射延拓,以及可以传入自定义阈值向量做尺度相关阈值。默认走的是周期延拓那一套,首尾不连续的信号在两端会有明显的波纹,这个后面避坑章会展开。
3.3 手写阈值链路:wavedec + wthresh + waverec
封装函数好用,但要做“只对第 3 层以下做阈值”“近似系数完全保留”“按层给不同阈值”这类精细控制,还是得回到手工链路。核心只有一句:噪声标准差用第一层细节系数的 MAD 估计,阈值用通用阈值公式,然后只对细节系数做阈值处理。
wname = 'db4'; N = wmaxlev(L, wname); [c, l] = wavedec(xn, N, wname); % 第一层细节系数的 MAD 估计,除以 0.6745 是无偏化常数 cD1 = detcoef(c, l, 1); sigma = median(abs(cD1)) / 0.6745; % 通用阈值 thr = sigma * sqrt(2*log(L)); fprintf('sigma = %.4f, 阈值 = %.4f\n', sigma, thr); % c 的布局是 [cA_N, cD_N, cD_{N-1}, ..., cD_1] % l(1) 是近似系数长度,后面全是细节系数 cT = c; cT(l(1)+1 : end) = wthresh(c(l(1)+1 : end), 's', thr); xrec = waverec(cT, l, wname); fprintf('手工链路 SNR = %.2f dB\n', snr(clean, xrec));索引l(1)+1 : end是关键。l(1)是第 N 层近似系数的长度,从它之后一直到末尾全是细节系数,这一刀切下去就把近似系数完整保留了下来,等价于wdencmp里的keepapp = 1。如果把近似系数也一起阈值处理,低频基线会被整个削掉,重构出来是一条贴着零轴的曲线,看起来“噪声全没了”,实际上是信号也没了。
wthresh的第二个参数决定阈值方式:'s'软阈值(系数向零收缩,连续但会整体压缩幅值)、'h'硬阈值(小于阈值置零,大于阈值原样保留,保幅值但不连续)、'g'非负绞杀(介于两者之间,小系数压得更狠、大系数基本不动)。选择上没有普适答案:冲击类信号偏'h'或'g',平滑类信号偏's'。
3.4 三种阈值规则横向对比
阈值规则指的是“阈值大小怎么定”,wden把这套逻辑封装成了一个字符串参数,方便快速横向比。四种常用规则里,sqtwolog就是通用阈值,rigrsure是基于 Stein 无偏风险估计的自适应阈值,heursure是前两者的折中,minimaxi按极小极大准则给阈值,对低信噪比更稳健。
rules = {'rigrsure', 'heursure', 'sqtwolog', 'minimaxi'}; names = {'SURE', 'Heuristic', 'Universal', 'Minimax'}; for k = 1:numel(rules) xd = wden(xn, rules{k}, 's', 'sln', N, 'db4'); fprintf('%-10s SNR = %.2f dB, 峰值衰减 = %.1f%%\n', ... names{k}, snr(clean, xd), 100*(1 - max(xd)/max(clean))); endwden的五个参数分别是:输入信号、阈值规则、软硬阈值、噪声尺度处理方式、分解层数、小波名。第三个参数'sln'表示按第一层细节系数估计噪声并做单尺度调整,另外两个选项'one'不做调整、'mln'做多尺度调整。这一段跑完,通常会看到Universal的 SNR 最高但峰值衰减最严重,Minimax峰值保留最好但残留噪声多一些。这种对比做两三次,你对自己手上这类信号的“取舍点”就有数了,比记任何结论都管用。
4. 二维图像与批量任务:wavedec2 系数布局、压缩阈值与 parfor
二维的小波分解和图像压缩是这套工具箱另一条主力链路。很多人卡在第一步:wavedec2返回的C是一维长向量,S是个矩阵,怎么从里面把水平、垂直、对角方向的细节取出来?搞不清这个布局,后面的阈值、压缩、可视化全都无从下手。
4.1 二维分解结构与系数索引
wavedec2(X, N, wname)返回[C, S],C是把所有子带按特定顺序串成的一维向量,S是(N+2)×3的尺寸矩阵。S的第一行是第 N 层近似子带(cA_N)的尺寸,第二到四行依次是第 N 层的水平(cH)、垂直(cV)、对角(cD)子带尺寸,再往下是 N-1 层、N-2 层,最后一行是原图尺寸。
| S 的行 | 内容 | 取系数的方式 |
|---|---|---|
| 1 | 第 N 层近似 cA_N | appcoef2(C, S, wname, N) |
| 2 | 第 N 层水平 cH_N | detcoef2('h', C, S, N) |
| 3 | 第 N 层垂直 cV_N | detcoef2('v', C, S, N) |
| 4 | 第 N 层对角 cD_N | detcoef2('d', C, S, N) |
| 5~(N+1)×3+1 | N-1 层重复上述四项 | 把最后一个参数换成对应层号 |
% 构造一张 256×256 的测试图,避开对图像处理工具箱的依赖 n = 256; [Xg, Yg] = meshgrid(1:n, 1:n); I = 0.5 + 0.3*sin(2*pi*Xg/32) + 0.2*(Yg > 128); In = I + 0.08*randn(size(I)); % 手动加高斯噪声 wname = 'db4'; N2 = 3; [C, S] = wavedec2(In, N2, wname); fprintf('S 的尺寸 = %d×%d,C 的长度 = %d\n', size(S,1), size(S,2), numel(C)); % 按层取子带 cA3 = appcoef2(C, S, wname, 3); [cH3, cV3, cD3] = detcoef2('all', C, S, 3); fprintf('cA3 %s, cH3 %s\n', mat2str(size(cA3)), mat2str(size(cH3))); psnr = @(ref, est) 10*log10(1/mean((ref(:)-est(:)).^2));detcoef2的第一个参数'all'会一次性返回三个方向的细节,比调三次省事。S的尺寸对 N=3 是 5 行,别按N+1去算。这段代码里psnr假设像素值归一化到 0~1,如果你读进来的是uint8,得先转double再除以 255,否则算出来的数字毫无意义。
4.2 wdenoise2 去噪与 wthresh 压缩两条路
二维去噪和二维压缩用的是同一套分解结果,区别只在“怎么改系数”。去噪是按估计的噪声水平做阈值收缩,压缩是按能量保留比例做大阈值截断。wdenoise2把去噪这一路封装好了,压缩这一路得自己写,但也就几行。
% 路线一:封装去噪 Idn = wdenoise2(In, 3, 'Wavelet', 'db4'); fprintf('wdenoise2 后 PSNR = %.2f dB\n', psnr(I, Idn)); % 路线二:手写压缩,按最大系数的比例设全局阈值 thr = 0.35 * max(abs(C)); CT = wthresh(C, 'h', thr); Icomp = waverec2(CT, S, wname); Icomp = min(max(Icomp, 0), 1); % 重构后可能越界,必须裁剪 keepRatio = nnz(CT) / numel(CT); fprintf('保留系数 %.1f%%,PSNR = %.2f dB\n', 100*keepRatio, psnr(I, Icomp));waverec2重构出来的矩阵尺寸可能比原图大几个像素,因为分解时的边界延拓会补长度,重构后多出来的部分需要按原图尺寸裁掉,或者直接靠索引截断。另外重构值会越界(小于 0 或大于 1),不裁剪直接显示就是一片黑或一片白,这个坑几乎每个新手都踩过。
压缩率和 PSNR 的关系不是线性的。阈值设到最大系数的 0.35,通常能保留 5%~15% 的非零系数,PSNR 还能维持在 30 dB 以上;继续加大阈值,系数是少了,但重构图像会出现明显的方块状振铃,尤其在斜边缘附近。判断压缩质量别只看 PSNR,肉眼看一眼对角子带重构的那张图(只保留 cD 层),振铃会非常明显。
4.3 把 App 的操作固化成函数,parfor 批量跑
批量处理是 App 永远做不了的事。一条产线上几百个测点的数据,手工点 App 不现实,得把整个流程封成函数,然后用parfor并行。封装时注意两点:函数的入参必须包含所有可调参数,不能有硬编码的层数和小波名;函数内部不要调用dwtmode这类全局状态函数,否则并行时各 worker 之间会互相干扰。
function [xd, residRatio] = denoiseOne(x, wname, N) % 单通道一维信号去噪,返回去噪结果和残差能量占比 x = double(x(:)); xd = wdenoise(x, N, 'Wavelet', wname, ... 'DenoisingMethod', 'SURE', ... 'ThresholdRule', 'Soft'); residRatio = norm(x - xd) / norm(x); enddataDir = fullfile(pwd, 'raw'); files = dir(fullfile(dataDir, '*.mat')); nFile = numel(files); % 预分配,parfor 里不能动态扩结构体 out = struct('name', cell(1,nFile), 'resid', cell(1,nFile), 'ok', cell(1,nFile)); parfor k = 1:nFile try S = load(fullfile(files(k).folder, files(k).name)); x = S.x; N = wmaxlev(numel(x), 'db4'); [out(k).name, out(k).ok] = deal(files(k).name, true); [~, out(k).resid] = denoiseOne(x, 'db4', N); catch ME out(k).name = files(k).name; out(k).ok = false; out(k).resid = NaN; fprintf('文件 %s 处理失败:%s\n', files(k).name, ME.message); end endtry-catch在批量里是必需品,单个文件格式异常不该让整批任务挂掉。存到out(k).resid里的残差能量占比只是一个粗筛指标,正常应该在 0.05~0.4 之间,明显偏小说明几乎没去噪,偏大说明信号被削掉太多。拿这个指标先粗筛一遍,把异常的那几个挑出来单独看,比逐个波形翻效率高得多。
注意:
parfor里调用的函数如果依赖工具箱的全局状态,结果在并行和串行下可能不一致。上批量之前先用parfor只跑 3 个文件,和串行的结果逐个对一遍,差在小数点后十位以内才算过。
5. 避坑排查:五条血泪记录,每条能省你半天
前面几章的流程跑通了,不代表结果可信。小波变换的失败往往是“看起来正常”的那种失败:曲线光滑、图像干净,但细节早就被抹掉了。下面五条,是我在不同项目里反复遇到的,每条按“现象 → 原因 → 解决”写清楚。
5.1 去噪后两端出现衰减振荡
现象:信号中段去噪效果不错,首尾各几十个采样点出现幅度递减的波纹,越靠边越明显,严重时首尾直接发散。
原因:wavedec的每一层都要在两端补数据,默认的对称延拓在信号本身首尾不连续时,会在边界造出一个阶跃,这个阶跃被逐层放大,重构后就变成振荡。分解层数越深,边界影响向内渗透的范围越大,因为每层下采样都会把边界效应往中间“拉”一截。
解决:三个手段按成本从低到高。先把层数降一到两层,多数情况下振荡就消失了;再把信号长度截断到2^N的整数倍并切到dwtmode('per'),首尾连续时完全没有边界伪影;最后是重构后直接裁掉首尾各一个滤波器长度的采样点(Lf-1个),别把这些点交给下游算法。
5.2 阈值处理之后信号整体漂移
现象:去噪后的曲线形状没问题,但整条线相对原始信号平移了一个常数,或者低频基线被拉平了。
原因:把wthresh直接作用在整个c上,连近似系数一起阈值处理了。近似系数承载的就是低频基线,阈值再小也会削掉一部分,层数越深、近似系数越少,削掉的基线越明显。
解决:只对细节系数做阈值,近似系数原样保留。手工链路里就是cT(l(1)+1:end) = wthresh(c(l(1)+1:end), ...);用wdencmp时把keepapp设成 1。如果确实需要对近似系数做平滑,单独用低通滤波,别混在阈值这一步。
5.3 uint8 图像分解后重构出一片黑
现象:imread读进来的图直接喂给wavedec2,分解正常,waverec2出来之后用imshow显示是一片黑或者满屏噪点。
原因:wavedec2内部会把输入转成双精度做卷积,但不会做数值归一化。uint8的取值是 0~255,小波系数动辄上千,重构后的值会大幅越界;如果你自己又用uint8强转,超出 255 的全被截断成 255,细节全丢。
解决:进分解之前先double(I)/255归一化,重构之后先min(max(x,0),1)裁剪,再按需im2uint8或其他方式转回去。中间所有计算都在双精度下做,别中途转类型。
5.4 waverec 报尺寸不匹配
现象:报错信息类似 “size mismatch” 或 “系数向量长度与索引不一致”,但代码看起来和示例一模一样。
原因:两个高频错误。一是分解用'db4'、重构用'sym4',小波名不一致时滤波器长度不同,l的预期就变了;二是手动截断或拼接了c,却没同步更新l,waverec按旧l去切新c,必然越界。
解决:把小波名抽成变量,分解和重构共用同一个;任何对c的结构性改动都同步维护l,实在不想维护就用appcoef/detcoef逐层取再逐层重构,让工具箱自己管索引。
5.5 cwt 时频图两端出现虚假高能量带
现象:cwt出来的时频图上,最左和最右各有一条能量很高的竖带,看着像信号两端的强突发,实际原始信号两端很平静。
原因:连续小波变换在每个尺度上把母小波平移过整段信号,靠近边界时小波的支撑范围超出信号区间,参与计算的样本不全,归一化后能量被放大。这条边界影响区的宽度由coi(cone of influence,影响锥)刻画,在图上通常用虚线标出。
解决:把coi返回的影响锥之外区域直接掩掉再分析,或者只截取信号中间段做结论。如果必须看边界,把信号两端各补一段数据(对称延拓或者用相邻时段的数据),变换完再裁掉。任何基于cwt系数的量化指标,都应该在掩掉影响锥之后计算,否则边界的假能量会污染统计。
6. 进阶:MODWT 的平移不变性和一个自检习惯
离散小波变换有个硬伤——平移敏感性。同一段信号往后挪一个采样点,各层系数会整体变化,去噪结果也跟着变,这对冲击检测类是致命的,因为脉冲落在不同相位上会被不同程度地压制。modwt(最大重叠离散小波变换)去掉了下采样这一步,各层系数和原信号等长,平移一个点只让系数平移一个点,形态不变。代价是计算量翻几倍,系数冗余度高,不能直接套用 DWT 的阈值。
% MODWT:返回 (N+1)×L 矩阵,第 1 行是尺度系数,2:end 是各层细节 w = modwt(xn, 'db4', 5); fprintf('MODWT 系数矩阵尺寸 = %d×%d\n', size(w,1), size(w,2)); sigma = median(abs(w(2,:))) / 0.6745; thr = sigma * sqrt(2*log(L)); wT = w; wT(2:end, :) = wthresh(w(2:end, :), 's', thr); xr = imodwt(wT, 'db4'); fprintf('MODWT + 阈值 SNR = %.2f dB\n', snr(clean, xr));注意modwt和
本文还有配套的精品资源,点击获取