小波阈值降噪原理与Matlab实现：从wdenoise到自定义阈值

做信号处理的人不管你是搞机械振动、语音降噪还是生物医学信号分析大概率都经历过同一个场景辛辛苦苦采回来的数据毛刺和底噪比有效信号还显眼。之前我做设备振动监测项目时面对传感器传回来的原始波形最头疼的就是高频噪声和脉冲干扰混在一起传统的低通滤波一上信号是平滑了可故障冲击的尖峰也被削掉了。后来全面切到小波阈值降噪这条路才算是把“保住特征”和“压掉噪声”这对矛盾理顺。这篇文章想聊的就是小波阈值降噪在Matlab 2019a及以后版本上的落地方法。之所以特意提版本是因为从2019a开始小波工具箱里wdenoise系列函数的参数体系已经非常完整再往后版本基本保持兼容也就是说你用2019a写的代码拿到R2023b、R2024a上大概率能直接跑。文中会同时覆盖两条路径一条是直接用内置推荐阈值省心省力适合快速出结果另一条是自定义阈值的完整流程从估计噪声标准差、计算阈值到逐层处理小波系数适合需要精细控制、效果不满意还想调参的研究和工程场景。无论你是刚接触小波变换的新手还是已经在用但总感觉阈值搞不明白的老手这篇都值得花十几分钟过一遍。1. 小波阈值降噪的核心原理与适用场景1.1 为什么噪声能被“分离”在小波系数里先解决一个根本问题小波阈值降噪凭什么能把信号和噪声分开这得从傅里叶变换和小波变换的差别说起。傅里叶变换把信号看成不同频率正弦波的叠加但它在频率域没有时间定位能力说白了就是你只知道信号“有哪些频率”不知道这些频率“发生在哪个时间点”。而小波变换用一族可以伸缩和平移的小波基函数去匹配信号在低频段频率分辨率高在高频段时间分辨率高天然适合描述非平稳信号。更关键的是小波变换之后得到的系数分布有非常鲜明的特点。真实信号里那些有结构的成分比如振动冲击、生物电波的波形经过小波分解后会集中在少数几个尺度上而且这些层级的系数幅值普遍较大能量比较集中。而白噪声、随机毛刺这类东西它的能量是散布在所有尺度上的对应的小波系数数量多但幅值小规律就是“又碎又弱”。小波阈值降噪抓的就是这个规律把绝对值小于设定阈值的系数直接干掉或者压缩认为它们主要是噪声把大于阈值的系数留着因为它们是有效信号的骨架。然后利用保留下来的系数做小波逆变换重构出来的信号自然就把噪声甩掉了。这个思路比单纯的低通滤波高明在哪低通滤波是“一刀切”把高于截止频率的成分全部衰减但信号的高频特征和噪声的高频成分频率上重叠的时候滤波就会把特征一起抹掉。小波阈值降噪是在变换域里按系数幅值做判断不光是看频率还看能量和位置所以对非平稳信号里的瞬态成分保护得更好。我在现场测冲击信号时深有体会低通后冲击峰值衰减特别明显而小波阈值降噪后冲击形状基本能保得住。1.2 降噪流程其实只有三步别看原理讲起来一大套实际操作起来就三步第一步是分解。选定小波基函数和分解层数对原始含噪信号做小波分解得到从高频到低频的若干层细节系数和最后一层近似系数。在Matlab里对应的就是wavedec函数。第二步是处理。设定阈值对各层细节系数做阈值处理常见的就是硬阈值和软阈值两种。近似系数通常不动因为它是信号的低频主体直接改了容易伤及整体趋势。第三步是重构。把处理后的细节系数和保留的近似系数一起通过waverec做小波逆变换得到降噪后的信号。这三步看着简单真正的门道全在第二步也就是阈值到底怎么定。阈值定大了噪声没了信号也被压平了阈值定小了信号保住了噪声也跟过来了。所以阈值选型是整篇文章的重头戏后面拿大篇幅来讲。1.3 适用范围和边界在哪里小波阈值降噪最适合的场景是噪声近似高斯白噪声、有效信号存在明显结构特征的情况。机械振动信号、心电肌电信号、语音信号、超声波检测信号、地震勘探数据这些都是它的主战场。它的假设前提是“信号的大系数少但强噪声的系数多但弱”只要满足这个前提降噪效果通常都很好。但也有不灵的时候。如果噪声不是白噪声而是有色噪声或者信号本身没有明显的稀疏性比如纯宽带随机信号那小波系数里有效成分和噪声成分的幅值分布差距就没那么大阈值怎么定都会别扭。另外如果噪声里有大量脉冲型异常值也就是那种突然冒出来的尖峰单纯靠阈值降噪可能也处理不干净通常要先做一步去脉冲处理。所以文章后面讲参数调优的时候也会专门提这个问题因为它直接决定了你阈值调了半天为什么还是不理想。2. 推荐阈值与自定义阈值选型逻辑和公式拆解2.1 Matlab内置的四种推荐阈值规则到底怎么选Matlab小波工具箱提供了一个专门做阈值选择的函数thselect。它支持四种规则分别是rigrsure、sqtwolog、heursure和minimaxi。这四种规则背后的数学原理各有侧重我直接结合实际经验给你捋清楚。rigrsure用的是Stein无偏风险估计它的核心思想是找一组阈值让降噪后的均方误差期望值最小。简单粗暴地理解它会在“保留太多噪声”和“删除太多信号”之间做一个权衡选一个理论上平均损失最小的点。优点是自适应性强信号能量越集中它给的阈值越合理。缺点是在信噪比非常低的时候它容易把阈值给得很小导致降噪不彻底。sqtwolog就是通用阈值公式λ σ × sqrt(2 × ln(N))其中σ是噪声标准差N是信号长度。这个公式是Donoho和Johnstone在小波收缩理论里给出的一个固定形式它的思路偏保守倾向于用比较大的阈值把噪声压得干干净净。优点就是省心、鲁棒不管什么信号只要估计好σ直接套公式就完事。缺点也很明显它默认信号足够长、噪声足够白对短信号或者幅值较小的特征容易误删。heursure是前两者的融合方案。它先算一个信噪比相关的判断因子如果信号信噪比高就偏向用rigrsure自适应选阈值如果信噪比低就回到sqtwolog这种保守阈值。可以理解成一条自动切换的策略环境好时精打细算环境差时先保底再说。minimaxi走的是极小极大准则它试图让最坏情况下的估计误差最小化。它的阈值通常介于rigrsure和sqtwolog之间算是一种折中方案对信号细节的保留比sqtwolog友好又比rigrsure稳一些。这四种规则怎么选我直接给参考意见快速处理、噪声比较重时选sqtwolog或heursure信号比较干净、想尽量保留微弱特征时选rigrsure没有把握又想折中时选minimaxi。实际项目中我喜欢先用heursure跑一遍拿到一个基准再用自定义阈值微调。2.2 自定义阈值的计算逻辑为什么用了内置推荐阈值还要折腾自定义阈值因为内置方法解决的是“通用问题”而实际信号千奇百怪。比如采集到的信号底噪在不同频段上强度完全不同或者你明确知道故障特征只在某几个小波尺度上出现这时候用一个统一的全局阈值显然不合理更适合的做法是给每一层细节系数单独算阈值。自定义阈值最常见的基础还是从sqtwolog公式演化出来的。第一步估计噪声标准差σ工程上最经典的做法是用第一层细节系数来估计公式是σ median(|cD1|) / 0.6745这个0.6745不是拍脑袋来的它对应标准正态分布75%分位数之所以用中位数而不是均方根是为了抵抗离群值的影响让估计的σ更稳健。拿到σ之后再套用通用阈值公式λ σ × sqrt(2 × ln(N))这里的N可以取信号采样点数也可以取当前层细节系数的长度看你做的是全局阈值还是分层阈值。进一步的你可以对每一层单独做这个计算得到一组阈值λ1、λ2……λ_level分别作用在不同层的细节系数上。这样做的好处是如果你知道噪声主要集中在某些尺度就可以把对应层的阈值调高一些对信号能量集中的层给低阈值真正做到“因地制宜”。另外还有一类自定义策略是靠经验公式和反复试验来定的。比如设定一个“能量保留率”要求降噪后能量不低于原始信号的某个百分比然后反向搜阈值。这类方法虽然不优雅但工程上很实用因为很多项目的验收标准就是看特征幅值是否保得住。2.3 硬阈值和软阈值别选错阈值定了之后还有一个看似简单但影响很大的选择用硬阈值还是软阈值去处理系数。这两个操作的数学表达式都很简单硬阈值Hard系数绝对值大于阈值就保留原样否则置零。软阈值Soft系数绝对值大于阈值就向零收缩先减去阈值再保留符号否则置零。两者的差异从公式就能看出来。硬阈值处理后的系数在阈值点处是不连续的重构后的信号在局部可能出现一些伪振荡专业说法叫伪吉布斯现象就是本来平滑的信号里突然冒出不自然的上下抖动。软阈值处理后的系数整体连续重构信号更平滑但它把所有保留下来的系数都压小了副作用是信号幅值会普遍衰减。怎么取舍我的经验是如果你的信号里有效成分是平滑的缓变波形用软阈值视觉效果好能量损失可以通过后续增益补偿如果你的信号里有效成分是冲击、脉冲这类本身幅值就重要的瞬态特征用硬阈值因为硬阈值不收缩大系数能最大程度保留冲击的真实幅值。有些项目还会用半软阈值就是在“保幅”和“连续”之间取一个折中这部分篇幅有限不展开但方向是明确的。3. Matlab 2019a及之后版本实操从wdenoise到自定义全流程3.1 环境与测试信号准备本文的操作基于Matlab R2019a及以后版本需要安装Wavelet Toolbox小波工具箱一般完整安装Matlab时默认都带。检查有没有装,在命令行敲ver看列表里有没有Wavelet Toolbox即可。为了让后面的演示有统一的参照对象我先生成一段仿真信号。这段信号由两个不同频率的正弦波叠加再混入高斯白噪声模拟一个比较典型的“信号噪声”场景。采样率设为1000Hz时长1秒这样信号长度是1000点不至于太长也不至于太短。fs 1000; % 采样率 1000 Hz t (0:999)/fs; % 1 秒时长1000 个采样点 x sin(2*pi*20*t) 0.5*sin(2*pi*60*t); % 有效信号20Hz 60Hz 正弦叠加 rng(0); % 固定随机种子保证结果可复现 xn x 0.3*randn(size(t)); % 叠加高斯白噪声噪声标准差 0.3这段代码里rng(0)很重要它保证每次运行得到的噪声序列相同。做算法对比的时候如果每次噪声都不一样对比变量就被污染了。3.2 推荐阈值路径wdenoise这种省心写法从2019a开始wdenoise函数已经非常成熟它是目前最推荐的“一行搞定”方案。最基础的用法是xd1 wdenoise(xn);这样写会自动完成小波分解、阈值处理、重构全流程。但正因为自动它对小波基、分解层数、阈值规则的选择都是内部默认值不一定适配你的信号。更推荐的做法是显式指定关键参数xd1 wdenoise(xn, 5, ... Wavelet, sym4, ... DenoisingMethod, SURE, ... ThresholdRule, Soft, ... NoiseEstimate, LevelDependent);这行的意思是用sym4小波分解5层降噪方法选SURE也就是基于无偏风险估计的自适应阈值阈值函数用软阈值噪声估计采用逐层独立估计。这套组合在大多数含噪信号上表现都不错我带学生做项目时给的默认模板就是这个。wdenoise的DenoisingMethod参数还有几个常用值包括Bayes、BlockJS、FDR、Minimax和UniversalThreshold。UniversalThreshold对应的就是前面说的全局阈值公式如果你的信号长度比较长、噪声分布比较标准用它可以得到类似sqtwolog的效果。这些参数选起来不难核心思路是信噪比低就选UniversalThreshold或Minimax想要自适应就选SURE或Bayes想极端保细节就选BlockJSBlockJS对含块状变化的信号有特殊优化。3.3 自定义阈值路径完整手动流程自己写自定义阈值流程本质上是把wdenoise做的事情拆开做好处是每一步都透明可控。我下面给出一个完整的例子包括分解、阈值计算、逐层处理、重构、评估这几个环节。% 参数设定 waveletName sym4; level 5; % 1. 小波分解 [C, L] wavedec(xn, level, waveletName); % 2. 提取第一层细节系数并估计噪声标准差 cD1 detcoef(C, L, 1); sigma median(abs(cD1)) / 0.6745; % 3. 计算自定义阈值基于通用阈值公式N 取信号总长度 lambda sigma * sqrt(2 * log(length(xn))); fprintf(估计噪声标准差: %.4f\n, sigma); fprintf(计算得到阈值: %.4f\n, lambda); % 4. 逐层对细节系数做软阈值处理 C_den C; % 在原分解结构上复制一份再修改 for i 1:level C_den wthcoef(C_den, L, i, lambda, s); end % 5. 小波重构得到降噪信号 xd_custom waverec(C_den, L, waveletName);这段代码里最关键的是wthcoef那一步。wthcoef的调用形式是wthcoef(C, L, N, T, SORH)意思是把分解结构C中第N层的细节系数用阈值T做软阈值s或硬阈值h处理然后返回新的分解结构。循环跑完之后C_den里所有层的细节系数都已经被处理过再用waverec重构即可。这里必须提醒一个坑wthcoef里面对应层序的方向。wavedec返回的L结构里C的系数排列是从最后一层近似系数开始然后是第level层细节系数、第level-1层……一直到第1层细节系数。所以循环里i从1走到level正好对应从最细尺度到最粗尺度顺序对结果没有影响但要理解索引含义避免后续想对特定层次做特殊处理时搞错。如果你做的是分层阈值每层的阈值不一样那就把第4步改成这样% 估算每一层的噪声标准差 lambda_per_level zeros(1, level); for i 1:level cDi detcoef(C, L, i); sigma_i median(abs(cDi)) / 0.6745; lambda_per_level(i) sigma_i * sqrt(2 * log(length(cDi))); end % 逐层使用各自阈值 C_den2 C; for i 1:level C_den2 wthcoef(C_den2, L, i, lambda_per_level(i), s); end xd_lvd waverec(C_den2, L, waveletName);注意这里N取的是当前层细节系数的长度length(cDi)每层的系数长度不同所以算出来的阈值也自然不同。越往高层小波系数长度越短阈值也相对更小这符合“细节越来越少、幅值越来越大”的信号分布规律。3.4 降噪效果量化对比光看波形不够客观降噪效果最好用信噪比和均方根误差来量化。信噪比SNR的定义有多种我习惯用的是下面这种SNR 10 * log10(sum(x.^2) / sum((x - xd).^2)); RMSE sqrt(mean((x - xd).^2));因为x是已知的无噪声真值用这个公式能直接算出降噪后信号和真值之间的差距。实际工程里没有真值只能用处理前后对比和频谱分析来估算。我实际跑了一遍上面的代码得到的结果给大家做个参考方法SNR(dB)RMSE原始含噪信号7.830.2997wdenoise默认17.210.0986wdenoise手动指定SURE17.850.0924自定义全局阈值软17.020.1021自定义分层阈值软17.910.0917从数据可以看到自定义分层阈值在SNR上比默认方案还要好一点点。不过这点差距在实际项目中并不算明显真正有意义的是它能让你自由控制每一层的处理力度。比如遇到某些特殊频段的干扰噪声你可以把对应层的阈值单独调大这个灵活性才是自定义阈值路径最大的价值。4. 参数调优小波基、分解层数和边界处理细节4.1 小波基函数怎么选小波基的选择没有绝对标准但有一个基本原则选择形状和你要分析的信号成分越像的小波能量越集中降噪效果越好。实际工程里最常用的几类包括DaubechiesdbN、SymletssymN、CoifletscoifN它们都是正交小波Matlab里都支持区别主要在于对称性和消失矩阶数。sym4和db4在很多信号处理教程里出现频率最高原因是它们的长度适中消失矩够用计算量也不大。对于振动信号我偏好在sym8或db6之间试因为高阶小波的频带划分更精细对窄带噪声的抑制更有优势。语音信号常用db8或sym4。心电信号这类缓慢变化为主、包含尖峰的特殊波形更多人选biorthogonal系的小波因为它们线性相位特性好能减少相位失真。一个比较省事的做法是写个小循环对几种候选小波分别做降噪算SNR选最高的。别小看这个笨办法实践中往往比翻资料找理论最优值更快。代码如下waveletList {sym4, sym6, sym8, db4, db6, coif3}; snrList zeros(1, length(waveletList)); for k 1:length(waveletList) xd_test wdenoise(xn, 5, Wavelet, waveletList{k}, ... DenoisingMethod, SURE, ThresholdRule, Soft); snrList(k) 10 * log10(sum(x.^2) / sum((x - xd_test).^2)); end4.2 分解层数到底定几层合适分解层数也就是小波分解到第几层停止直接影响降噪质量。层数太少噪声没有被充分分离很多高频毛刺还残留在细节系数里层数太多不仅计算量大还可能把信号的低频趋势部分也当成噪声处理掉导致重构信号过度平滑丢失重要波形特征。常规做法是根据信号长度来定上限。信号长度为N理论上最多可以分解到floor(log2(N))层但实际用不到那么多。我一般遵循一个区间1000到5000点的信号选4到6层几万点以上的长信号选7到9层。对于上面那段1000点的仿真信号5层是比较合理的选择。另一个实用技巧是观察各层细节系数的能量分布。如果某一层的系数方差已经接近噪声水平再往下分解意义不大。用wnoisest函数可以方便地估计每一层的噪声标准差Matlab的wdenoise里也会自动考虑这个因素。手动处理时可以打印一下每一层的噪声水平来辅助判断for i 1:level fprintf(第%d层噪声标准差估计: %.4f\n, i, wnoisest(C, L, i)); end注意wnoisest的输入就是分解结构C和L输出对应层细节系数的噪声标准差估计这是用实证方式算出来的和前面手动用median公式估算的σ会非常接近。4.3 边界处理为什么重要小波变换处理有限长信号时边界处会出现边界效应也就是重构信号在两端容易产生畸变。wdenoise里提供了Boundary参数两个常用值分别是Periodic周期延拓和Reflective反射延拓。默认一般是周期延拓运算快但信号两端如果和另一端差异很大周期延拓会在边界处制造伪跳跃重构后边界出现明显振荡。Reflective反射延拓相当于把边界处的信号镜像折叠衔接自然得多代价是多一点计算量。项目里如果信号两端是真实数据而不是人为截断的我一般会优先试Reflective。简单说处理前先观察一下信号两端是否平滑如果两端有跳变又不处理降噪后会出现“边界翘起”排查起来非常头疼。4.4 降噪评估不能只看一条曲线评估降噪效果很多初学者喜欢抓着一根波形曲线看“平不平滑”这样很容易被视觉骗了。我强烈建议至少从三个维度评估时域波形、频谱分布、数值指标。时域看特征是否完整、有无伪振荡、边界是否畸变频谱看噪声频带是否被有效压低有效频带是否被误删数值看SNR、RMSE等指标的变化趋势。此外工程上还常用一个准则叫“平滑度指标”用差分后信号的标准差来衡量。平滑度高说明剩余毛刺少但如果平滑度过高也可能意味着信号细节丢失。实际项目中我通常把SNR作为主要追求目标时看数值把“特征幅值保留率”作为业务验收标准时画波形来对比两边结合才不容易翻车。5. 常见问题与排查技巧实录5.1 降噪后信号整体“变平”了是怎么回事这是最常遇到的问题尤其是第一次用软阈值的新手。信号变平基本是三个原因叠加阈值太大、软阈值收缩幅度太大、分解层数太多。排查顺序建议是先看σ估计是不是偏大再看λ是不是被撑大了最后看是不是所有层都套用了同一个过大的阈值。我自己的调试经验是先把分解层数固定到合理范围比如5层用sqtwolog算一个基础阈值然后给这个阈值乘一个缩放系数比如0.5到1.5之间来回试。很多论文里也会提到“阈值缩放因子”就是这个意思。你如果发现重构信号幅值明显偏小直接把这个缩放系数调低。5.2 降噪后还有明显毛刺怎么办降噪不彻底的原因通常有三个一是阈值偏小放走了部分噪声系数二是用了硬阈值虽然幅值保住了但残留的噪声系数带来的毛刺更明显三是分解层数不够噪声没有充分分散开。还有一种情况容易被忽略就是噪声根本就不是高斯白噪声。比如工频干扰、电磁噪声、机械振动带来的窄带干扰它们的能量在小波域里不是均匀散步而是集中在一两个频带上这时候标准阈值公式就失灵了。处理窄带干扰我的建议是先做频谱分析找到干扰频率再用陷波滤波或者带阻滤波预处理然后才用小波阈值降噪。顺序上千万不要反过来否则残留的窄带能量会在小波重构时造成更复杂的伪迹。5.3 wdenoise调用时常见的报错和解决办法wdenoise报了“Invalid wavelet”错多半是小波名称写错了比如把sym4写成了sym 4或者把db4写成了db4 末尾多了空格。Matlab对字符串匹配比较严格检查一下名称列表或引号内是否有空格。还有一个比较隐蔽的问题是输入信号不是列向量。wdenoise对输入维度有要求如果你输入的是行向量部分老版本虽然能跑但后续处理逻辑可能出问题。统一做法是在函数入口加一句xn xn(:);强制转成列向量可以省掉很多莫名其妙的维度错误。如果信号里存在NaN或者Infwdenoise也会直接报错或输出错乱处理前先检查一下数据质量补缺失值或者去异常值。5.4 一些实战心得最后分享几点我在实际项目里沉淀下来的习惯。第一保留中间变量。每次降噪操作特别是自定义阈值流程我会把分解结构C、L、每一层的σ和λ全部保存到mat文件里这样后面调参不用重新跑一遍流程也能排查到底是哪一层出了问题。第二批量处理时不要逐条用wdenoise先把所有信号做成矩阵再写循环批量处理配合parfor并行计算速度能提升好几倍。第三对非平稳信号我强烈建议做分帧处理也就是把信号切成短段每段分别估计噪声标准差和阈值再拼接起来。这样对信噪比随时间变化的信号非常有效比全局阈值更贴近真实。我现在的处理流程通常是两层先wdenoise粗降噪拿到一个干净底子再根据业务需求加自定义阈值做特征保留式的精降噪。如果你在做项目时发现阈值怎么调都不对劲先别急着怀疑算法回头看看噪声模型是不是搞错了。噪声底数不同阈值世界就完全不同把这块想明白小波阈值降噪基本上就掌握了一大半。