小波阈值去噪的算法是近些年比较流行的一种滤波方法，由于其阈值函数有着众多的改进方式和改进空间，改进阈值函数也往往可以作为创新点和亮点写到论文中，所以对于正在搞相关研究的同学们写论文是比较友好的（轻松水论文方式+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这两个封装函数的方便之处在于，同学们差不多只需要导入自己要分析的数据，然后直接调用函数就可以了，更多的时间精力可以放在写论文和改进算法上。
