小波阈值去噪的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_apply：

function 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越接近 1，k越小，阈值越低，保留的细节越多。

这里有个容易被绕进去的概念要澄清：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 后，高频噪声和有效瞬态往往会纠缠在一起，很难干净分离。

几个容易踩的细节。分解层数上限受信号长度限制，N最大不能超过floor(log2(numel(x)))，否则最后一层细节系数长度不够，阈值计算没有统计意义。短信号（几百个点）建议直接Level=3，长信号（几万点）也别一上来就 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); end

garrote的思路是用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。从那以后，“先枚举、再定型、最后上真实数据”成了固定流程。规则数值、小波基、分解层数这三个维度可以分开调，先定规则，再定层数，最后微调小波基，一次只动一个变量，出问题才知道改的是谁。希望帮到你。

本文还有配套的精品资源，点击获取