类EMD-小波阈值降噪（以ICEEMDAN-样本熵-小波阈值方法为例）

最后更新于：2024-04-07 21:35:03

一、代码运行环境

MATLAB2018b及更新版本。

二、程序介绍

注：图标代表该m文件为脚本文件，可以直接运行；图标代表函数文件，在没有输入变量的情况下无法直接运行。更详细的解释可以看这里。

1.demoEMDsWavelet文件夹

进行ICEEMDAN-样本熵-小波阈值滤波的脚本文件，函数filEMDsWaveletTh的测试文件，可以直接运行。程序运行完成后，在MATLAB的工作区，双击reSig变量，可以查看滤波后的具体数值。

运行该文件，将绘制以下图片：

另外在命令行窗口将会打印：

2.filEMDsWaveletTh.m文件

ICEEMDAN-样本熵-改进的小波阈值的封装函数，该函数中的ICEEMDAN方法、样本熵、小波阈值的相关参数可以进行修改替换。

function reSig = filEMDsWaveletTh(sig)
% ICEEMDAN-样本熵-改进的小波阈值的封装函数，该函数中的ICEEMDAN方法、样本熵、小波阈值的相关参数可以进行修改替换
% 输入参数:
%   sig - 待滤波信号
%
% 输出参数:
%   reSig - 滤波后的信号

2.filterWaveletTh.m文件

使用小波阈值法实现滤波的主函数。

阈值滤波的不同参数设置全部封装在此函数，只需要输入待滤波数据、小波名称、阈值类型、小波分解水平和阈值选择规则，实现一行代码完成小波阈值滤波。

function s = filterWaveletTh(data,wname,SORH,lev,tptr)

% 使用小波阈值法实现滤波，本文件是为了形式统一方便调用的一层封装，具体实现在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：  阈值类型 可以选择：
%       -'s' 软阈值
%       -'h' 硬阈值
% 
%       SORH ='a1'时，采用改进方法1，改进后的软阈值：
%                     参考论文：孙万麟, 王超. 基于改进的软阈值小波包网络的电力信号消噪[J]. 海军工程大学学报, 2019(4).
%       SORH = 'a2'时，采用改进方法2：
%                     参考论文：刘冲, 马立修, 潘金凤,等. 联合VMD与改进小波阈值的局放信号去噪[J]. 现代电子技术, 2021, 44(21):6.
%       SORH = 'a3'时，采用改进方法4：
%                     参考论文：基于改进小波阈值-CEEMDAN算法的ECG信号去噪研究
%       SORH = 'a4'时，采用改进方法3：
%                     参考论文：基于改进小波阈值函数的语音增强算法研究
%                     采用该方法需要输入两个调节因子，分别是a4_alpha和a4_gamma，取a4_alpha>0，0<a4_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_alpha>0
%          -a4_gamma：改进方法4的gamma值设置 0<a4_gamma<1
%          -a5_alpha：改进方法5的alpha值设置
%          -a5_beta：改进方法5的beta值设置
% 输出：
% s：经过滤波后得到的数据

3.FilterEffectEvaluation.m文件

对滤波效果进行评估。目前包含的评估指标有：信噪比SNR，均方差MSE，波形相似系数NCC。

function [SNR,MSE,NCC] = FilterEffectEvaluation(ori,fil)

% 对滤波效果进行评估
% 目前包含的评估指标有：信噪比SNR，均方差MSE，波形相似系数NCC
% 输入：
% ori：无噪声的原始数据，一维序列
% fil：滤波后的数据，一维序列
% 输出：
% SNR：信噪比
% MSE：均方误差
% NCC：波形相似系数

4.kthselect.m文件

经店长修改后的小波阈值滤波阈值选择函数（原函数为thselect），加入了 VisuShrink 方法，即统一阈值去噪方法。

function thr = kthselect(x,tptr,lev,wname)


% 经 Mr.看海 修改后的小波阈值滤波阈值选择函数（原函数为thselect）
% 加入了 VisuShrink 方法，即通用阈值去噪方法
% 该方法的实现主要参考：https://core.ac.uk/download/pdf/217399263.pdf （Denoising techniques - a comparison）
% 如果要使用VisuShrink方法，需要令TPTR =
% 'visushrink'，注意大小写，并指定lev和wname值。此时输出thr为lev*1维度的数组
% 其中lev：  小波分解水平
%     wname：小波名称

5.kwden.m文件

经店长修改后的小波阈值滤波函数（原函数为kwden），加入了 visushrink 方法，即统一阈值去噪方法。

function [xd,cxd,lxd,thrs] = kwden(in1,in2,in3,in4,in5,in6,in7)

% 经 Mr.看海 修改后的小波阈值滤波函数（原函数为wden）
% 加入了 visushrink 方法，即统一阈值去噪方法
% 加入了改进阈值函数，具体见kwthresh.m文件
% 加入了为满足改进阈值函数额外的调节因子的入口参数

6.kwthresh.m文件

经店长修改过后的阈值函数（原函数为wthresh）。除软阈值、硬阈值外，添加了改进阈值函数类型。

function y = kwthresh(x,sorh,t)

% 经 Mr.看海 修改过后的阈值函数（原函数为wthresh）
% 添加了改进阈值函数类型
% 输入：
% x：阈值函数的输入值，不需修改
% 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
% t：阈值
% 输出：
% y：阈值函数输出值
% 测试阈值函数图线可以用下述代码：
% figure('color','w');plot(-5:0.1:5,kwthresh(-5:0.1:5,'a1',1));
% 测试阈值函数图线与软硬阈值图线对比图可用下述代码：
% figure('color','w');plot(-5:0.1:5,kwthresh(-5:0.1:5,'a1',1));hold on
% plot(-5:0.1:5,kwthresh(-5:0.1:5,'s',1))
% plot(-5:0.1:5,kwthresh(-5:0.1:5,'h',1))
% legend('改进阈值方法','软阈值','硬阈值')

硬阈值、软阈值、改进阈值方法1、改进阈值方法2的阈值函数分别如下图：

软阈值

硬阈值

改进阈值方法1

改进阈值方法2

改进阈值方法3

改进阈值方法4

改进阈值方法5

7.testGenFeaEn.m文件

该脚本文件原本为熵特征提取算法的主测试脚本，放在此处是为了方便大家在替换熵特征时作为参考。

特征提取函数genFeatureEn的测试文件，可以直接运行。程序运行完成后，在MATLAB的工作区，双击fea变量，可以查看求得的具体数值。

8.pICEEMDAN.m

封装好的ICEEMDAN画图程序。

函数参数说明：

function imf = pICEEMDAN(data,FsOrT,Nstd,NE,MaxIter)
% 画信号ICEEMDAN分解图
% 输入：
% data为待分解信号
% FsOrT为采样频率或采样时间向量，如果为采样频率，该变量输入单个值；如果为时间向量，该变量为与y相同长度的一维向量。如果未知采样频率，可设置为1
% Nstd为附加噪声标准差与Y标准差之比
% NE为对信号的平均次数
% MaxIter：最大迭代次数
% 输出：
% imf为经ICEEMDAN分解后的各imf分量值
% 例1：（FsOrT为采样频率）
% fs = 100;
% t = 1/fs:1/fs:1;
% data = sin(2*pi*5*t)+2*sin(2*pi*20*t);
% imf = pICEEMDAN(data,fs,0.2,100);
% 例2：（FsOrT为时间向量，需要注意此时FsOrT的长度要与y相同）
% t = 0:0.01:1;
% data = sin(2*pi*5*t)+2*sin(2*pi*20*t);
% imf = pICEEMDAN(data,t,0.2,100);20

画图效果如下（参考）：

9.kCEEMDAN.m

调整后的CEEMDAN函数，整合了iceemdan函数，并针对工具箱进行了统一性的调整。

函数参数说明：

function imf = kICEEMDAN(data,Nstd,NE,MaxIter)
% 统一封装后的ICEEMDAN函数，原始代码来源：http://bioingenieria.edu.ar/grupos/ldnlys/metorres/re_inter.htm
% 
% 输入：
% data：待分解的数据
% Nstd：为附加噪声标准差与data标准差之比
% NE：对信号的平均次数
% MaxIter：最大迭代次数
% 
% 输出：
% imf：内涵模态分量，统一为n*m格式，其中n为模态数，m为数据点数。例如 imf(1,:)即IMF1，imf(end,:)即为残差

% 注意：在使用该代码之前，请务必安装工具箱：https://khsci.com/docs/index.php/2020/04/09/1/

三、快速开始

1.运行测试脚本

先在MATLAB里打开下载好的文件夹，然后运行demoEMDsWavelet.m程序，程序运行完毕后如果没报错，且正常画出滤波对比图像，则说明运行环境正常，程序正确。

2.修改仿真数据/导入数据

复制一个demoEMDsWavelet.m的文件副本，在这个副本里做如下修改：

（1）第一种情况，你可能想要对你自己要研究的仿真数据进行滤波测试，此时你需要对 demoEMDsWavelet.m 脚本文件中的第1小节内容进行修改替换即可。需要注意的是，请最好保持变量名的一致，即将无噪声信号命名为ps，添加噪声后的信号命名为x，x和ps都应该是一维数据。此时数据替换完成。

（2）第二种情况，你可能是想对一段真实采集的数据进行滤波，此时需要根据你的文件类型的不同（excel,txt,csv等），将数据导入MATLAB的方法有所不同。同学们可以看博主针对常用文件的导入方法的这个教程，教程上未包含的数据类型，大家可以再参考这个文档。导入完成后请讲这个待滤波信号命名为x。此时你可能没有纯净信号（在大多数实际应用中是没有纯净信号的），这不影响你是用小波阈值滤波进行信号去噪，但是你将无法计算滤波评价指标。

3.替换模态分解方法、熵特征类型、小波阈值参数等（可选）

双击打开filEMDsWaveletTh.m函数文件。

该函数文件的第1部分，目前为ICEEMDAN算法，你可以选择替换成其他的模态分解算法，店铺中有相关代码，大家可以在下边链接中找到。

该函数文件的第2部分，可以替换为其他熵特征。大家可以参考testGenFeaEn.m文件中的调用方法。

该函数文件的第3部分，大家可以根据实际情况修改阈值。

该函数文件的第4部分，可以根据需要选取小波名称、阈值函数、分解层数、阈值选择规则等参数。

上述替换过程大家有不清楚的，我录制了视频供参考（需完整版）。

4.运行程序

此时运行程序即可。

（需要注意，对于使用真实采集的数据进行滤波，且无纯净信号的情况下，需要删除代码中的第3节滤波评价指标计算，并将第4节滤波对比图像按照注释进行修改）

四、关于完整版与公开版代码

五、获取公开版程序（需使用电脑浏览器打开）

注：公开版代码需使用MATLAB2022a及以上版本。

六、获取完整版程序（使用电脑浏览器或者手机浏览器打开）

获取通道一（淘宝）：点击此处获取完整版程序

获取通道二（本页面）：点击下面“立即支付”按钮，付款后获取完整版代码下载链接和售后联系方式~本通道处于测试阶段，使用该通道可以额外优惠（仅需118元）。付款完成后刷新一下本页面即可看到下载链接。

（注意支付跳转失败的话，请使用浏览器打开本页面）

七、完整版代码重要更新

20240315 完成初版代码

八、常见问题

无。