基于相对熵优化VMD的非局部均值去噪方法

前言

本博客对论文 "优化VMD与NLM结合的信号去噪"中的算法进行介绍与仿真。该论文为解决人体心电信号易出现噪声干扰的情况，提出一种自选取参数的变分模态分解(VMD)与非局部均值(NLM) ECG信号去噪方法。针对人为选取变分模态分解参数不准确的问题，运用一种自适应的参数确定方法，根据不同的原始信号，确定不同的参数并对信号进行分解，得到一组信号分量，计算各个分量的样本熵，根据样本熵的值，选取出噪声主导分量和有效分量，对噪声主导信号进行NLM去噪，将去噪后的信号分量与剩余的有效信号分量进行重构得到去噪信号。

论文中所用算法及其MATLAB代码见以下链接：

VMD算法

相对熵优化变分模态分解(KL-VMD)

排列熵、模糊熵、近似熵、样本熵的原理及MATLAB实现之样本熵

一、非局部均值(NLM)算法

1 算法背景

非局部均值算法是近年来图像处理领域的新兴算法之一，是 Buades 等为了解决传统消噪方法在图像信号处理出现的无法很好的保留图像细节及边缘的问题而提出的一种新的消噪方法。NLM 算法的原理是利用图像中存在的众多相似结构的特性，以像素点所在领域为单元，在图像中搜索与该单元相似的成分，并对这些相似成分加权平均操作，以此消除图像中的异常噪声。Buades 等提出NLM 算法的同时证明了其一致性，即NLM 算法并非只对特定的像素点有效，而是适用于图像中的任一像素点，即在自然图像中，任一结构均存在与之类似的结构样本。这些相似特性不仅表现在二维图像中，在一维信号中也具有广泛体现。TRACEY 等和胡新海等应用 NLM 算法分别处理ECG信号和地震信号，去噪后的信号均具有较大的信噪比。

2 在一维信号中的应用步骤

（1）假设1个一维信号x(t)被1个高斯白噪声信号n(t)污染，则含噪声信号为：

（2）NLM算法是通过计算y(t)中全部相似块的加权平均 X(t) 对原始信号 x(t) 进行估计进而实现去噪，计算过程如下：

式中：w(s,t) 为以 s 和 t 为中心的2个搜索块之间的相似度；N(s)为搜索域中全部点的集合；Z(s)为归一化常数，表示全部搜索块相似度之和 。

（3）Z(s) 为归一化因子，其计算公式为：

（4）w(s,t) 需要满足如下2个条件：

（5）w(t,s)的计算过程如下：

式中：λ为滤波器带宽参数；Δ为以 s 为中心的搜索块；LΔ为以 t为中心的搜索块；为以 s 和 t为中心的2个搜索块之间欧氏距离的平方和，其数值越小，表示相似度越大； 控制权重 w(s,t) 的衰减速度，即直接决定信号的最终滤波程度。

二、非局部均值算法对一维信号去噪

1 MATLAB程序

% 论文复现：优化VMD和NLM结合的信号去噪
% 利用K-L散度(相对熵)确定VMD分解信号的K值和惩罚因子alpha，得到一组信号分量;
% 计算各个分量的样本熵，根据样本熵的值，选取出噪声主导分量和有效分量;
% 对噪声主导信号进行非局部均值(NLM)去噪;
% 将去噪后的信号分量与剩余的有效信号分量进行重构得到去噪信号.
clc;
clear;
close all;
tic;

%% 仿真信号
load data;
x=data(1,1:2000)';
N=length(x);
y=awgn(x,20,'measured');  % 添加信噪比为20dB的高斯白噪声
%% 利用K-L散度(相对熵)确定VMD分解信号的最优K值和最优惩罚因子alpha
alpha = 1000;  % 惩罚因子固定
[Best_K,Best_alpha] = KL_VMD(y,alpha);

%% VMD分解
% y：为待分解的时域信号
% u：分解模式的集合
% omega：估计模式中心频率
% Best_K：最优模态数
% Best_alpha：最优惩罚因子

% VMD分解的输入参数
tau = 0;            % 噪声容忍度
DC = 0;             % 无直流分量
init = 1;           % 初始化中心频率为均匀分布
tol = 1e-7;         % 收敛准则容忍度

[u, ~, omega] = VMD(y, Best_alpha, tau, Best_K, DC, init, tol);

[~, sortIndex] = sort(omega(end,:));
u=u(sortIndex,:);
X=u';
%% 样本熵
r0=0.2;   % r为相似容限
dim=2;    % dim为嵌入维数
sampEn=zeros(1,Best_K);
CC=zeros(1,Best_K);
figure;
for i=1:Best_K
    r=r0*std(u(i,:));
    % 计算样本熵值
    sampEn(i) = SampleEntropy( dim, r, u(i,:));
    % 计算相关系数
    CC(i)=corr(X(:,i),y,'type','Pearson');   % 相关系数
    % 各个IMF分量图
    subplot(Best_K,1,i);plot(X(:,i));
    xlabel('t/s');ylabel(['IMF',num2str(i)]);
end
%% 非局部均值(NLM)去噪
PatchHW=10;  % patch half-width：size of smallest feature, in samples
P = 1000;    % neighborhood search width：wider=more computation but more channces to find a similar patch
lambda = 0.6*.02; % set bandwidth for NLM - here set by an "eyeball estimate"
denoise_data=zeros(4,N);
for i=6:7
    [denoise_data(i-5,:),debugEcg]= NLM( u(i,:),lambda,P,PatchHW );
end
%% 信号重构
Y=sum(u(1:5,:))+sum(denoise_data);

%% 画图
% 时域波形对比图
figure;
subplot(3,1,1);plot(x);xlabel('t/s');ylabel('幅值');title('原始信号x(t)');
subplot(3,1,2);plot(y);xlabel('t/s');ylabel('幅值');title('含噪信号y(t)');
subplot(3,1,3);plot(Y);xlabel('t/s');ylabel('幅值');title('去噪后信号Y(t)');

toc;

非局部均值算法的MATLAB代码链接：非局部均值(NLM)信号去噪方法

2 去噪效果

 图1所示是利用NLM算法对原始心电信号进行去噪。

图1 NLM算法对原始心电信号进行去噪

对原始心电信号添加20dB高斯白噪声，再利用相对熵优化VMD的非局部均值去噪方法对含噪声心电信号进行去噪，去噪效果如图2所示。

图2 去噪效果对比图

结束语

大家的点赞和关注是博主最大的动力，博主所有博文中的代码文件都可分享给您（除了少量付费资源），如果您想要获取博文中的完整代码文件，可通过C币或积分下载，没有C币或积分的朋友可在关注、点赞和评论博文后，私信发送您的邮箱，我会在第一时间发送给您。博主后面会有更多的分享，敬请关注哦！

参考文献

[1] 尹佳璠, 陈小奇, 李世林. 优化VMD与NLM结合的信号去噪[J]. 计算机工程与设计, 2021, 42(4):1135-1142.

[2] Tracey Brian H, Miller Eric L. Nonlocal means denoising of ECG signals[J]. IEEE Transactions on Biomedical Engineering, 2012,59(9):2383-2386.

[3] 唐晓红, 胡俊锋, 熊国良, 等. 自适应非局部均值及在轴承故障检测中的应用[J]. 振动、测试与诊断, 2019, 39(1):61-67.

[4] 万书亭, 彭勃. 基于非局部均值去噪和快速谱相关的滚动轴承早期故障诊断方法[J]. 中南大学学报(自然科学版), 2020, 51(1):76-85.