【滤波专题-第6篇】小波阈值去噪方法看这一篇就明白了~（附MATLAB实现）

小波阈值去噪的算法是近些年比较流行的一种滤波方法，由于其阈值函数有着众多的改进方式和改进空间，改进阈值函数也往往可以作为创新点和亮点写到论文中，所以对于正在搞相关研究的同学们写论文是比较友好的（轻松水论文方式+1）。本篇将用尽量易懂的方式对小波阈值的原理进行讲解，帮大家梳理几个效果还可以的改进阈值函数，并提供一种非常便捷的MATLAB实现方法，供同学们使用。

小波阈值去噪的基础思想还是比较简单的，也就是通过分解+有选择的重构，实现去除噪声成分，留下关键信息的作用。我们从两个角度去理解就可以，谜底就在谜面上，这两个理解角度的关键词就是“小波”和“阈值”。

一、先说“小波”

需要注意的是，这里提到的小波指的是“小波分解”，而不是小波变换、离散小波变换或者小波包分解等等，涉及小波的概念比较多，这里要做好区分。

小波分解在我们之前的文章里也提到过，通过小波分解我们可以得到一串小波分解系数。而小波阈值去噪的理论基础就是[1]：

小波变换具有很强的去数据相关性，他能使信号的能量在小波域集中在一些大的小波系数当中，而噪声的能量却分布于整个小波域内。因此，经小波分解后，信号的小波系数幅值要大于噪声的系数幅值，可以认为，幅值比较大的小波系数一般以信号为主，而幅值比较小的系数在很大程度上是噪声。于是，采用阈值的办法，使大部分噪声系数减小至0，即将小于阈值部分的系数作为干扰噪声去除，然后再进行小波重构，从而实现降噪。

二、然后是“阈值”

在小波阈值去噪中，阈值函数体现了对超过和低于阈值的小波系数的不同处理策略，是阈值去噪中关键的一步。

对于“策略”，主要也包括两部分内容，分别是“阈值”的选择和“阈值函数”的选择。

（一）阈值的选择

所谓阈值，就是设定一个数值，当小波系数的绝对值小于这个阈值的时候，通通赋值为0（注：有些改进方法不是）。

阈值的选择方法常用的有这几种：

'rigsure' — Use the principle of Stein's Unbiased Risk.
'heursure' — Use a heuristic variant of Stein's Unbiased Risk.
'sqtwolog — Use the universal threshold √2ln(length(x)).
'minimaxi' — Use minimax thresholding.
'VisuShrink'— 通用阈值方法

上述前四种是MATLAB内置的阈值选择方法，还有一种比较常用的VisuShrink方法没有收录到MATLAB官方库中，笔者在编程过程中一起把它集成进去了。（代码见下边的文章介绍）

阈值选择的的几种方法这里就不展开说了，其计算公式在论文中都不难找到。

（二）阈值函数的选择

阈值函数的设计通常是使用小波阈值降噪方法的核心关键所在，在写论文时这里可以加以改进，实现更好的滤波效果，妥妥的创新点。

传统的阈值函数有两种，分别是硬阈值和软阈值。

1.硬阈值

所谓硬阈值函数，其实就是“不做处理”。。。换句话说，我们不是已经将小于阈值的小波系数的置0了嘛，大于阈值的小波系数就保持原值（这操作真够硬的）。所以阈值函数画出来就是下图这样的，这样处理就导致了函数的不连续，“不连续”这种事情在信号处理中显然是不够优雅的。

从图上可以看出，转折处就是阈值，在这个例子里也就是2

2.软阈值

而所谓的软阈值函数，其实也没“软”到哪去（核心理念依旧有些生硬），只是将上图中W=2的幅值跳变的位置做了一个平移，他的图像如下图：

阈值依然是2

这样处理带来一个显而易见的问题：小波分解系数值变啦！相当于每个大于阈值的系数都减掉了2（也就是阈值），这毫无疑问会在信号重构的时候引入误差。

3.第一种改进的阈值函数

既然软阈值和硬阈值有缺点就好说了，此时我们就可以有针对性地进行改进，而改进的方法无非就是围绕一个核心目的：在消除阈值函数的不连续性的同时，使函数迅速靠近硬阈值函数。

就比如下图这样（蓝线）：

这是一种比较简单且相对有效的阈值函数，表达式是这样的[2]：

式中 λ 就是阈值。

这个改进方法表达式简单明了，阈值函数曲线也比较漂亮。

4.第二种改进的阈值函数

这里再举一个例子，其阈值函数的表达式是这样的[3]：

按照论文里的说法，配合这个阈值函数使用的是VisuShrink阈值计算方法。

这个阈值函数图像是这样的：

这个方法相对于第一种改进方法靠近硬阈值的速度更快一些，不过在阈值位置还是有一个小一些的跳变，这个方法介于硬阈值与软阈值之间，也算是改进方法的一种。

5.其他更多改进方法

还有更多的阈值函数改进方法，不过很多种都要引入一个或多个调节因子，也就是引入了新的变量，这种情况下阈值函数的种类就比较多了，这里不再一一介绍。这里挑几个，贴一下他们的阈值函数图像，并附上参考论文，有需要的同学们可以更进一步查阅。

（1）改进方法3：《基于改进小波阈值-CEEMDAN算法的ECG信号去噪研究》

为了体现效果，这个图像范围画到了10

（2）改进方法4：《基于改进小波阈值函数的语音增强算法研究》，这个方法在小于阈值的位置也做了平滑处理。这个方法引入了两个调节因子alpha和gamma，可以调节阈值函数形状。

（3）改进方法5：《基于VMD与改进小波阈值的地震信号去噪方法研究》，这个方法也引入了两个调节因子，分别是alpha和beta，可以更加灵活地调整图线形状。

这个是alpha和beta都等于2的结果

三、MATLAB案例实现

我把上边讲到的7种阈值方法都进行了MATLAB代码实现。下面采用一个案例来进行演示。

首先我们仿真生成一段信号，并向其中加入白噪声。x为未加入白噪声的纯净信号，sig为加入白噪声后的信号，代码如下：

%% 1.生成仿真信号，x为无噪声信号，sig为添加噪声后的信号
rng(123); %设置随机种子，保证每次生成噪声的一致性
N = 200;
t = linspace(0,1,N);
x = 4*sin(4*pi*t);
x = x - sign(t-0.3) - sign(0.72-t);
sig = x + 0.2*randn(size(t));
figure('Color','w');subplot(2,1,1);plot(t,x);title('未加入白噪声信号')
subplot(2,1,2);plot(t,sig);title('加入白噪声信号')

下边我们尝试做一下软阈值、硬阈值、以及5种不同改进阈值方法的滤波效果对比吧。

（一）MATLAB编程实现方法

在MATLAB中，软、硬阈值的小波阈值滤波方法通过调用wden函数可以实现，但是想要实现改进阈值方法，对于没有MATLAB代码基础的同学们可能就得花费一些周章了，而且MATLAB官方函数里还没有visushrink这个很常用的阈值选择方法，所以笔者改造了wden、thselect和wthresh三个函数文件，并进一步封装成filterWaveletTh函数，延续本专栏中以往代码的风格，实现“一行代码”完成小波阈值去噪的效果。当然啦，这里所说的“一行代码”还需要配合一些参数的设置。

举个例子，如果对上述的sig信号进行软阈值滤波，选取'db3'类型的小波，分解水平为3，阈值选择原则为visushrink，那代码可以这样写：

s = filterWaveletTh(sig,'db3','s',3,'visushrink',[]);  %调用函数进行滤波，filterWaveletTh函数获取方法见文末

这样就可以得到滤波后的数据s。

我们来画图看一下滤波效果：

%% 画滤波前后对比图
figure('color','w')
subplot(311);plot(x,'k');title('原始信号（未加入噪声）')
subplot(312);plot(sig,'k');title('原始信号（加入噪声）')
subplot(313);plot(s,'k');title('滤波后信号')

软阈值滤波效果

结合我们之前讲过的计算一下滤波效果评价指标。计算SNR、MSE和NCC这三个指标值，计算方法也同样进行了封装，就像这样调用就行：

ori = x;  %无噪声信号
fil = s;  %滤波后信号
[SNR,MSE,NCC] = FilterEffectEvaluation(ori,fil);
% 对滤波效果进行评估
% 目前包含的评估指标有：信噪比SNR，均方差MSE，波形相似系数NCC
% 输入：
% ori：无噪声的原始数据，一维序列
% fil：滤波后的数据，一维序列

软阈值滤波后的SNR、MSE和NCC值分别为：26.4164、0.021658、0.99887。

这个效果还是可以的，下面我们再计算一下另外几种滤波方法的结果。

（二）硬阈值及五种改进方法的滤波结果

下面我们直接贴上来滤波效果对比图以及评价指标计算结果，具体代码实现都与上边的例子类似，直接给定参数调用就可以。

硬阈值滤波结果：滤波后的SNR、MSE和NCC值分别为：27.1828、0.018154、0.99905

改进方法1：滤波后的SNR、MSE和NCC值分别为：28.4053、0.0137、0.99928

改进方法2：滤波后的SNR、MSE和NCC值分别为：25.1265、0.029149、0.99847

改进方法3：滤波后的SNR、MSE和NCC值分别为：27.6984、0.016122、0.99916

改进方法4：滤波后的SNR、MSE和NCC值分别为：28.4599、0.013529、0.99929

改进方法5：滤波后的SNR、MSE和NCC值分别为：28.3339、0.013927、0.99927

做个表格列一下评价指标，做一个直观的对比：

小波阈值方法	SNR	MSE	NCC
软阈值	26.4164	0.021658	0.99887
硬阈值	27.1828	0.018154	0.99905
改进方法1	28.4053	0.0137	0.99928
改进方法2	25.1265	0.029149	0.99847
改进方法3	27.6984	0.016122	0.99916
改进方法4	28.4599	0.013529	0.99929
改进方法5	28.3339	0.013927	0.99927

这三个指标里，SNR和NCC都是越大越好，MSE是越小越好，所以总体来看改进方法4对于这个信号来说是比较好的小波阈值方法。

这里有三点注意事项需要提醒大家：

第一。改进方法4和改进方法5都引入了调整因子，通过调参有可能得到更好的效果，不过这里仅做演示，就不详细地做参数调试啦。

第二。虽然改进的阈值方法在理论上对阈值函数进行了完善，但是并不意味着改进方法的滤波结果一定优于软/硬阈值方法，就像这个例子里的改进方法2。又或者换了一个分析对象，改进方法4的效果也比不上软/硬阈值也是有可能的，所以在实际应用过程中要灵活分析处理。

第三。SNR、NCC和MSE指标，都是需要知道“纯净信号”才能计算的，对于只知道真实含噪信号，不知道纯净信号的情况，是没办法计算这三个指标的。我在之前的文章里进行过详细的说明：Mr.看海：【滤波专题-第4篇】滤波器滤波效果的评价指标（信噪比SNR、均方误差MSE、波形相似参数NCC）

（三）获取小波阈值的封装函数

上边提到了两个笔者封装的函数，一个是实现小波阈值滤波的filterWaveletTh，还有一个是进行滤波效果评估指标计算的FilterEffectEvaluation，函数的详细说明如下，大家在调用的时候建议详细读一下：

function s = filterWaveletTh(data,wname,SORH,lev,tptr,options)
% 使用小波阈值法实现滤波，本文件是为了形式统一方便调用的一层封装，具体实现在kwden.m、kwdencmp.m、kwthresh.m、kthselect.m文件中
% 这4个文件均为经本人调整过后的函数文件，原函数分别为wden,wdencmp,wthresh,thselect，所以关于这几个函数更多介绍可以查看对应的官方帮助文档
% 输入：
% data：待滤波数据
% wname：小波名称，可选范围参考这里：https://ww2.mathworks.cn/help/wavelet/ref/wfilters.html?searchHighlight=wname&s_tid=srchtitle_wname_2#d123e130597
% SORH：阈值函数类型
%       sorh ='a1'时，采用改进方法1，改进后的软阈值：
%                     参考论文：孙万麟, 王超. 基于改进的软阈值小波包网络的电力信号消噪[J]. 海军工程大学学报, 2019(4).
%       sorh = 'a2'时，采用改进方法2：
%                     参考论文：刘冲, 马立修, 潘金凤,等. 联合VMD与改进小波阈值的局放信号去噪[J]. 现代电子技术, 2021, 44(21):6.
%       sorh = 'a3'时，采用改进方法4：
%                     参考论文：基于改进小波阈值-CEEMDAN算法的ECG信号去噪研究
%       sorh = 'a4'时，采用改进方法3：
%                     参考论文：基于改进小波阈值函数的语音增强算法研究
%                     采用该方法需要输入两个调节因子，分别是alpha和gamma，取alpha>0，0<gamma<1
%       sorh = 'a5'时，采用改进方法3：
%                     参考论文：基于VMD与改进小波阈值的地震信号去噪方法研究
%                     采用该方法需要输入两个调节因子，分别是alpha和beta% lev：  小波分解水平
% tptr： 阈值选择规则，有可选类型：
%       '' — Adaptive threshold selection using the principle of Stein's Unbiased Risk Estimate (SURE).
%       'sqtwolog' — Fixed-form threshold is sqrt(2*log(length(X))).
%       'heursure' — Heuristic variant of 'rigrsure' and 'sqtwolog'.
%       'minimaxi' — Minimax thresholding.
%       'visushrink' - 通用阈值去噪方法
% options：部分阈值函数需要补充的参数设置，即某些调节因子，需要用结构体方式调用幅值，如options.a4_alpha=2。如果不需要设置，可以赋值为options=[]。具体如下：
%          -a4_alpha：改进方法4的alpha值设置
%          -a4_gamma：改进方法4的gamma值设置
%          -a5_alpha：改进方法5的alpha值设置
%          -a5_beta：改进方法5的beta值设置
% 输出：
% s：经过滤波后得到的数据
% 理论讲解见：https://zhuanlan.zhihu.com/p/579187348/

function [SNR,MSE,NCC] = FilterEffectEvaluation(ori,fil)
% 对滤波效果进行评估
% 目前包含的评估指标有：信噪比SNR，均方差MSE，波形相似系数NCC
% 输入：
% ori：无噪声的原始数据，一维序列
% fil：滤波后的数据，一维序列
% 输出：
% SNR：信噪比
% MSE：均方误差
% NCC：波形相似系数
% 理论讲解见：https://zhuanlan.zhihu.com/p/558808890

这两个封装函数的方便之处在于，同学们差不多只需要导入自己要分析的数据，然后直接调用函数就可以了，更多的时间精力可以放在写论文和改进算法上。

需要上边这个函数文件以及测试代码的同学，可以下边链接获取：

小波阈值滤波算法（软阈值、硬阈值、改进的阈值） | 工具箱文档

end.关于滤波专题

目前计划要讲到的包括：

-1.FIR有限冲激响应、IIR无限冲激响应数字滤波算法的理论讲解、滤波器设计方法和MATLAB代码实现
-2.滤波器滤波效果的评价指标
-3.“类EMD”方法（即包括EMD、EEMD、CEEMD、VMD等一系列方法）的滤波算法讲解与实现
-4.“类EMD”方法与ICA结合的滤波算法讲解与实现
-5.小波阈值滤波方法讲解与实现
-6.卡尔曼滤波方法讲解与实现
...（其他方法随时补充）