小波阈值去噪的MATLAB实现:软/硬阈值与5种改进规则
简介面向需要做信号去噪的 MATLAB 用户这份小波阈值函数包将 wden、thselect、wthresh 等原生函数改造后统一封装为 filterWaveletTh一行代码即可切换软阈值、硬阈值及 5 种改进阈值规则适用于振动分析、故障诊断、语音/传感器信号预处理等场景。资源共 7 个文件由 5 个 p 函数文件核心算法与效果评估、1 个 m 演示脚本和 1 份 docx 说明文档组成整体仅 18KB轻量易集成。目前已有 1320 人浏览学习。拿到手后可先运行示例脚本感受不同阈值方法的滤波差异再结合说明文档明确各参数含义并在自己数据上快速替换调用省去从零编写阈值去噪代码的时间借助效果评估函数可对比去噪前后指标。1. 小波阈值MATLAB代码到底解决什么问题从一次信号去噪的翻车现场说起做信号去噪的工程师迟早会撞上小波阈值这条技术路线。某一次我在处理某采集系统的振动数据时噪声把有效脉冲淹得只剩下轮廓先试了带通滤波结果噪声没去干净信号里的陡峭特征也被削平了。后来把思路换成小波阈值去噪才意识到这类问题的关键不在“滤掉哪个频带”而在“把信号和噪声按尺度拆开再对噪声主导的细节系数做阈值处理”。这篇笔记要交付的东西很具体一套封装好的MATLAB函数能直接跑软阈值、硬阈值以及5种改进阈值规则。适合正在做信号预处理、故障诊断、传感器数据清洗的工程师也适合刚接触小波分析、想找一份能抄的代码的研究生。先说一个反直觉的结论硬阈值实现最简单翻车率却不低标题里的“改进阈值”也不是越多越好选错规则甚至比不处理更糟。下面从原理、封装、代码到坑一层层拆开。2. 小波阈值去噪的原理与封装设计为什么一个小函数能顶掉一套滤波流程2.1 小波阈值去噪的最小闭环分解、阈值、重构三步走小波阈值去噪和非平稳信号是天生一对。普通滤波器在频域里做全局操作信号一旦和噪声在频带上有重叠滤波就会连同信号一起削掉。小波变换把信号按尺度拆成近似系数和细节系数噪声能量主要落在细节系数里而信号的瞬态特征会集中在少数几个大系数上。于是思路变成分解出各层细节系数把绝对值小于阈值的系数压掉或置零再重构回去。噪声被压掉大系数代表的信号特征留下来。这段逻辑对应的最小闭环代码是这样% 小波阈值去噪最小闭环分解 - 阈值处理 - 重构 [C, L] wavedec(x, N, wname); % 小波分解C 为级联系数L 记录各层长度 d detcoef(C, L, 1:N); % 取出第 1 到第 N 层细节系数 d_new cell(1, N); for j 1:N d_new{j} d{j} .* 0; % 演示用把细节系数全部置零 end C_new appcoef(C, L, wname, N); % 保留最后一层近似系数 for j N:-1:1 C_new [C_new, d_new{j}]; % 按 C 原有顺序拼回细节系数 end x_den waverec(C_new, L, wname); % 重构信号这段代码有两点必须说明。第一点是C的结构wavedec返回的系数顺序是最后一层近似系数、最后一层细节系数、倒数第二层细节系数……一直到第一层细节系数L记录每一段的长度。拼回C_new时顺序必须严格一致从第 N 层往第 1 层倒着拼接一旦顺序错位重构出来就是一段错乱的波形。第二点是waverec重构时会沿用分解时的边界延拓方式所以小波基和延拓模式全程要保持一致。这里的“置零”是阈值处理的极端情况实际去噪时我们会用阈值函数决定哪些系数保留、哪些压缩下面第 3 章细说。2.2 封装函数的接口设计一个函数覆盖全部阈值模式有了最小闭环下一步就是把它封装成可复用的函数。封装的目的很直接调用者不需要关心C、L的结构不需要记得细节系数拼接顺序只需要给信号、选规则、定参数。我把阈值处理拆成两个维度一个是阈值函数决定系数怎么被处理也就是软阈值和硬阈值另一个是阈值规则决定阈值本身怎么算对应 5 种改进方法。两个维度互相组合就是标题里“软阈值、硬阈值以及 5 种改进阈值方法”的由来。主函数框架如下子函数的完整实现留到第 3 章function [x_den, info] wthresh_denoise(x, varargin) % 小波阈值去噪封装函数 % Rule 可选: % soft / hard 软阈值、硬阈值 % sqtwolog / rigrsure / heursure / % minimaxi / layered 5 种改进阈值规则 % 组合写法: hard-rigrsure 表示用 rigrsure 规则取阈值、硬阈值处理 p inputParser; addRequired(p, x, (v) isnumeric(v) isvector(v)); addParameter(p, Wavelet, sym4, ischar); % 小波基默认 sym4 addParameter(p, Level, [], isnumeric); % 分解层数默认自动 addParameter(p, Rule, heursure, ischar); % 阈值规则 parse(p, x, varargin{:}); x x(:); if isempty(p.Results.Level) N min(5, floor(log2(numel(x)))); % 默认层数受信号长度限制 else N p.Results.Level; end [threshFunc, threshRule] splitRule(p.Results.Rule); [C, L] wavedec(x, N, p.Results.Wavelet); d detcoef(C, L, 1:N); sigma zeros(N, 1); lambda zeros(N, 1); d_new cell(1, N); for j 1:N sigma(j) median(abs(d{j})) / 0.6745; % MAD 噪声估计抗尖峰 lambda(j) thrule_estimate(d{j}, threshRule, sigma(j)); d_new{j} thresh_apply(d{j}, threshFunc, lambda(j)); end C_new appcoef(C, L, p.Results.Wavelet, N); for j N:-1:1 C_new [C_new, d_new{j}]; % 按 wavedec 顺序拼回 end x_den waverec(C_new, L, p.Results.Wavelet); info struct(sigma, sigma, lambda, lambda, ... Level, N, Wavelet, p.Results.Wavelet); end function [threshFunc, threshRule] splitRule(rule) % 解析 hard-rigrsure 这类组合写法 if contains(rule, -) parts strsplit(rule, -); threshFunc parts{1}; threshRule parts{2}; elseif strcmp(rule, soft) || strcmp(rule, hard) threshFunc rule; threshRule sqtwolog; % 纯软/硬阈值时阈值固定用通用规则 else threshFunc soft; threshRule rule; end end几个参数的设计理由值得展开。Wavelet默认用sym4因为 sym 系列小波近似对称在信号端点处的重构畸变比 db 系列小工程上省去很多边界处理的麻烦。Level默认取min(5, floor(log2(numel(x))))这是经验值分解超过 5 层后细节系数长度变得太短阈值估计不稳定而且计算量成倍上涨。Rule的默认值我选heursure因为它会根据信噪比自动在固定阈值和风险估计之间切换对未知信号的适应性最好。噪声估计用median(abs(d{j})) / 0.6745这是 MAD 估计比std稳健得多——噪声里一旦出现尖峰std会被拉高导致阈值偏大有效信号被误伤。info结构体把每层的 sigma 和 lambda 输出这是排查问题的“后悔药”第一次跑完先去info里看阈值是否合理再决定调参方向。3. 软阈值、硬阈值与5种改进阈值每一条规则都能直接抄的MATLAB实现3.1 软阈值与硬阈值连续性与稀疏性的取舍软阈值和硬阈值是两种最基础的阈值函数也是 5 种改进规则的落点。硬阈值的逻辑是绝对值大于阈值的系数原样保留小于等于阈值的置零。软阈值则更温和大于阈值的系数先减去阈值小于负阈值的加上阈值中间区域置零。换句话说软阈值做了整体收缩硬阈值只做选择。两种方式各有代价。硬阈值在阈值点处不连续重构出的信号容易出现局部抖动看起来像波形表面粘了一层细毛刺软阈值连续重构波形更平滑但所有保留的系数都被削掉了一个阈值量幅值会系统性偏小。如果业务场景对幅值保真要求高比如测振幅、测峰值硬阈值更合适如果对波形平滑度敏感比如做频谱分析前的预处理软阈值更稳。我一般把软阈值当默认硬阈值作为备选。这段逻辑落到代码里就是统一入口thresh_applyfunction d_out thresh_apply(d, func, lambda) % 统一的阈值函数入口func 支持 soft / hard switch func case soft % 软阈值向零收缩连续但幅值有偏差 d_out sign(d) .* max(abs(d) - lambda, 0); case hard % 硬阈值直接置零保留大系数原值但阈值点不连续 d_out d; d_out(abs(d) lambda) 0; otherwise error(未知阈值函数: %s, func); end end代码逻辑不复杂但有几个实现细节值得注意。软阈值用sign(d) .* max(abs(d) - lambda, 0)比逐元素if-else快得多MATLAB 向量化一次搞定硬阈值先复制d再置零避免直接修改原数组防止后续需要原始系数做对比时数据被污染。至于lambda怎么来就是下面 5 种改进规则要做的事两者是解耦的——这也是为什么要封装的原因想换规则只改thrule_estimate想换处理方式只改thresh_apply。3.2 5种改进阈值规则的实现固定阈值、无偏风险估计到分层自适应5 种改进规则分别是sqtwolog、rigrsure、heursure、minimaxi和layered。其中前四种是业界验证过的经典规则第五种layered是给那些“信号能量在不同层分布差异极大”的场景准备的自适应方案。说句实话MATLAB 自带工具箱里也有阈值规则函数但只覆盖前四种且内部实现是黑匣子没法插自己的改进逻辑。自己写的好处就是每个 case 对应一条公式想调参数、想加规则改动范围被限制在一个分支里。整段实现如下function lambda thrule_estimate(d, rule, sigma) % 5 种改进阈值规则 固定规则返回本层阈值 n numel(d); switch rule case sqtwolog % 固定阈值渐近最优适合平稳信号叠加白噪声 lambda sigma * sqrt(2 * log(n)); case rigrsure % 无偏风险估计用小波系数的 SURE 准则选阈值 lambda sure_threshold(d, sigma); case heursure % 启发式信噪比低时兜底到固定阈值高时交给 SURE eta (sum(d.^2) - n) / n; crit (log2(n))^1.5 / sqrt(n); if eta crit lambda sigma * sqrt(2 * log(n)); else lambda sure_threshold(d, sigma); end case minimaxi % 极大极小准则的工程近似适合弱信号保护 lambda sigma * (0.3936 0.1829 * log2(n)); case layered % 分层自适应按本层疑似信号占比动态缩放阈值 lambda sigma * sqrt(2 * log(n)) * layered_factor(d, sigma); otherwise error(未知阈值规则: %s, rule); end end function lambda sure_threshold(d, sigma) % SURE 无偏风险估计对归一化系数遍历候选阈值取风险最小者 d_abs sort(abs(d(:)) / sigma); n numel(d_abs); risk zeros(n, 1); for k 1:n t d_abs(k); risk(k) n sum(min(d_abs, t).^2) - 2 * sum(d_abs t); end [~, idx] min(risk); lambda sigma * d_abs(idx); end function k layered_factor(d, sigma) % 本层显著系数占比越高说明信号成分越多阈值越要放松 sparsity mean(abs(d) 2 * sigma); k 2.5 - 1.5 * sparsity; k min(max(k, 0.6), 2.5); % 限幅防止阈值过小或过大 end逐个规则说适用场景。sqtwolog的公式是sigma*sqrt(2*log(n))它的理论依据是高斯噪声条件下渐近最优简单粗暴适合信噪比稳定、噪声分布均匀的信号比如平稳振动信号的预处理。rigrsure通过遍历候选阈值、最小化 SURE 风险来选阈值对信号中的细节成分保护最好代价是计算量随系数长度近似平方增长信号很长时建议只在高频层使用或者先降采样再做。heursure本质是个开关先算一个信噪比估计量eta低于阈值就走sqtwolog高于阈值就走rigrsure它适合你对信号先验一无所知的时候做默认规则。minimaxi的公式是极大极小准则的工程近似它把“最坏情况下的均方误差”压到最小对弱信号保护明显代价是噪声残留比sqtwolog多一点。layered是我在项目里常用的自研变体核心思想是如果某一层细节系数里大系数占比高说明这层以信号为主阈值要放松如果大系数寥寥无几说明这层基本是噪声阈值要收紧。layered_factor里sparsity越接近 1k越小阈值越低保留的细节越多。这里有个容易被绕进去的概念要澄清5 种改进规则解决的是“阈值取多少”的问题软硬阈值解决的是“拿到阈值后系数怎么处理”的问题。标题里把它们并列是因为实际使用时要两两组合。Rule参数写成soft或hard时内部会默认搭配sqtwolog规则写成hard-rigrsure这样带前缀的组合就能明确指定阈值函数和阈值规则。这种设计让 7 种命名模式覆盖了 10 种组合调用时不用记两套参数。4. 用封装函数快速去噪最小调用示例与参数搭配表4.1 最小调用示例三行代码从带噪信号到去噪结果封装函数最终要能“一把梭”。最常见的用法是信号进来选一个小波基、定一个层数、给一个规则输出去噪结果和每层的参数信息。下面这个例子可以在 MATLAB 里直接跑通只需要 Wavelet Toolbox 支持wavedec和waverec这是绕不开的依赖动手前先确认许可证。% 生成测试信号50Hz 与 180Hz 正弦叠加高斯白噪声 fs 2000; t (0:1999) / fs; x sin(2*pi*50*t) 0.5*sin(2*pi*180*t); x_noisy x 0.3*randn(size(x)); % 调用封装函数软阈值 启发式规则 [x_den, info] wthresh_denoise(x_noisy, ... Wavelet, sym6, Level, 4, ... Rule, heursure); % 三行绘图对比 subplot(3,1,1); plot(t, x); title(纯净信号); subplot(3,1,2); plot(t, x_noisy); title(带噪信号); subplot(3,1,3); plot(t, x_den); title(去噪结果);跑完先别急着看波形先看控制台里info的内容。info.sigma是每层噪声标准差估计info.lambda是每层实际用的阈值。如果第一层的lambda比细节系数的中位数还大好几倍说明阈值定高了信号可能被削如果lambda小到只有 sigma 的一半说明去噪力度不够。这个习惯能帮你从“调参靠猜”变成“看数调参”。初次上手建议把Rule从heursure换成hard再跑一次对比软阈值和硬阈值在波形毛刺上的差异你会直观感受到第 3 章讲的连续性问题的分量。4.2 参数怎么搭配小波基、分解层数与阈值规则的推荐表参数搭配没有万能解但有一套我摸出来的初始经验值可以让你少走一半弯路。小波基的选择看信号形态信号波形平滑、变化缓慢选 db 系列信号有陡峭突变、需要保边缘选 sym 系列。分解层数不是越多越好4 到 6 层覆盖了绝大多数工程场景层数超过 7 后高频噪声和有效瞬态往往会纠缠在一起很难干净分离。信号场景推荐小波基分解层数推荐阈值规则备注平稳振动信号 白噪声db45~7sqtwolog噪声分布均匀固定阈值最稳冲击 / 突变信号sym44~5rigrsure保护瞬态特征优先心电 / 脑电等生理信号sym66~8rigrsure细节波形容错低语音 / 音频降噪db84~5minimaxi避免过平滑导致听感发闷传感器混合噪声sym43~5layered信噪比波动大自适应更稳未知信号首次处理sym45heursure先摸底再按指标细调几个容易踩的细节。分解层数上限受信号长度限制N最大不能超过floor(log2(numel(x)))否则最后一层细节系数长度不够阈值计算没有统计意义。短信号几百个点建议直接Level3长信号几万点也别一上来就 7 层先 5 层起步看info.lambda再决定要不要加深。Rule的组合前缀很好用比如hard-minimaxi适合既要幅值保真、又怕过度平滑的弱信号场景但前提是你已经确认硬阈值的毛刺不影响后续分析。参数表只是起点真正的定参依据是第 6 章的量化对比。5. 小波阈值去噪避坑指南5个高频翻车现场与排查方法5.1 去噪后端点发散波形两端翘起现象重构信号中间平滑但两头出现大幅度的假波像是信号被硬生生掰弯了。原因大部分情况下是重构时小波基和分解时不一致或者边界延拓方式在分解、重构之间被隐式换了。wavedec用 sym4 分解waverec如果写成了 db4两者边界延拓特性不同重构结果在两端会严重失真。解决小波基名称全程使用同一个变量不要二次手敲如果用了自定义延拓参数分解和重构必须保持一致。这个坑我见过不止一次通常是“重构代码复制粘贴时改了名字”。排查方法很简单把info.Wavelet打印出来和waverec第 4 个参数比对。5.2 阈值算得太大有效信号被拦腰削减现象去噪后波形确实平滑了但峰值明显变矮原本陡峭的上升沿被磨成了缓坡。原因噪声方差估计用了std而噪声里一旦有少量尖峰标准差会被大幅拉高阈值随之膨胀把有效系数一起压掉了。解决改用 MAD 估计取第一层细节系数的中位数再除以 0.6745。这里是整个封装里我最坚持的一个设计因为 MAD 对尖峰不敏感工程数据里的野值很难骗到它。如果换了 MAD 后峰值还是偏小把Rule换成rigrsure它选阈值时不直接依赖 sigma对幅度保护更好。5.3 硬阈值去噪后波形表面布满毛刺现象用hard模式去噪后噪声大块是没了但波形表面残留细密抖动像锯齿一样。原因硬阈值函数在正负阈值处不连续系数被生硬置零后重构会在阈值点附近激起振荡分量。解决换soft如果幅值保真要求高我常用一个介于软硬之间的 garrote 函数兜底——大系数保留幅值小系数平滑收缩。代码很短function y th_garrote(x, lambda) % garrote 收缩介于软硬阈值之间兼顾幅值与连续性 y zeros(size(x)); mask abs(x) lambda; y(mask) (x(mask).^2 - lambda^2) ./ x(mask); endgarrote的思路是用x - lambda^2/x代替x - lambda幅值损失比软阈值小连续性比硬阈值好。它不属于标题里的 5 种改进规则但作为后备方案值得留在工具库里。使用时注意mask之外的系数要置零不能沿用原值。5.4 分解层数一味加深波形出现“台阶化”现象去噪结果变得过分平坦像楼梯一样一段一段的局部细节全部消失。原因层数太多时细节系数被大量清零只剩下近似系数在撑场面重构信号自然丢失高频微变。不少新手以为层数越多去噪越干净实际上 7 层以上效果大概率倒退。解决从 4 层起步观察info.lambda和波形如果高频毛刺还在加一层如果波形开始发僵减一层。短信号尤其别贪几百个点用 3 层就够。5.5 阈值规则选错越“改进”效果越差现象同一段信号sqtwolog出来很干净rigrsure出来还有明显噪声换个场景结果反过来。原因每种规则的信噪比假设不同没有银弹。sqtwolog假设噪声均匀rigrsure适合高信噪比minimaxi偏保守。解决别在真实数据上试错先在模拟信号上跑一遍枚举对比。构造一段和真实信号形态相似的干净信号加上相同强度的噪声把 7 种Rule全跑一遍用 SNR 排序选出最优再套到真实数据上。这套流程第 6 章给完整代码。6. 进阶用法用SNR和MSE把阈值规则选择从玄学变成指标参数表给的是经验起点真正决定用哪条规则要靠量化指标说话。第一步先写一个评估函数function [snr, mse] eval_denoise(x_clean, x_den) % 参考信号存在时用 SNR 和 MSE 量化去噪效果 e x_clean - x_den; mse mean(e.^2); snr 10 * log10(sum(x_clean.^2) / sum(e.^2)); end用法很直接构造一段和真实业务信号形态接近的模拟干净信号叠加同等级噪声然后对Rule的所有取值跑一遍循环记录每组snr。哪个规则得分高就用哪个上真实数据。这个方法远比“逐个规则用眼睛看波形”靠谱因为它把参数选择从感觉变成了指标。真实数据没有参考x_clean时可以退一步用残差统计对去噪后的信号做一次高阶差分看差分能量是否接近噪声水平能给出一个无参考的定性判断。我个人的教训是某一次处理一批传感器采集数据嫌麻烦直接用了默认的heursure结果有效脉冲的峰值被削掉了一截后来在模拟信号上枚举了一次规则才发现那个场景下rigrsure比heursure高出近 3 个 dB。从那以后“先枚举、再定型、最后上真实数据”成了固定流程。规则数值、小波基、分解层数这三个维度可以分开调先定规则再定层数最后微调小波基一次只动一个变量出问题才知道改的是谁。希望帮到你。本文还有配套的精品资源点击获取