CN114548147A - 一种基于EEMD-Hilbert变换的ECG去噪方法 - Google Patents

Abstract

本发明公开了一种基于EEMD‑Hilbert变换的ECG去噪方法，先采用集成经验模态分解方法对原始ECG_noise进行模态分解，得到由多个固有模态函数构成的一组固有模态函数，然后然后基于得到的一组固有模态函数依次进行Hilbert变换、求解相位以及求解瞬时频率，基于每个固有模态函数对应的一组瞬时频率的中心瞬时频率对每个固有模态函数进行分类，将分类为以基线漂移噪声为主的固有模态函数舍弃，分别去除分类为以信号为主的固有模态函数和以高频噪声为主的固有模态函数中的肌电干扰噪声后进行重构得到最终去噪后的ECG；优点是能够对基线漂移噪声和肌电干扰噪声进行更加有效的降噪，得到的最终去噪后的ECG信号仍然保留着明显的ECG形态特征，去噪效果好。

Description

一种基于EEMD-Hilbert变换的ECG去噪方法

技术领域

本发明涉及一种ECG去噪方法，尤其是涉及一种基于EEMD-Hilbert变换的ECG去噪方法。

背景技术

在ECG信号(Electrocardiogram signal，心电图信号)采集过程中，无法避免的会带有基线漂移噪声以及肌电干扰噪声。而通过可穿戴设备采集得到的ECG信号会受到更加严重的噪声干扰，使得ECG信号无法有效的被利用。

为了得到干扰较小的，能够有效利用的ECG信号，在ECG信号去噪领域内，小波阈值去噪方法被广泛的使用。小波阈值去噪方法的通常具有一下五个步骤：步骤一、使用小波变换将含噪信号(即原始ECG信号)逐级分解为多级不同尺度的细节小波系数；步骤二、分别对所得到的多级细节小波系数进行噪声标准差评估，并根据评估得到的噪声标准差设置小波阈值；步骤三、基于小波阈值提出阈值函数(硬阈值函数、软阈值函数等)；步骤四、利用提出的阈值函数处理多级细节小波系数得到去噪后的细节小波系数；步骤五、将去噪后的细节小波系数进行小波逆变换以此重构去噪信号。

但是，上述小波阈值去噪方法由于阈值评定方法的有效性有限，且肌电干扰噪声的频谱范围与ECG信号的频谱范围存在重叠，在应对肌电干扰噪声强度较强时，采用阈值函数对噪声与信号频谱重叠部分进行去噪，信号部分可能会被当作噪声，使得去噪后的ECG信号不但会产生局部失真的现象，而且ECG信号的形态特征可能会被当成噪声一起被滤除，由此导致最终去噪效果较差。

发明内容

本发明所要解决的技术问题是提供一种去噪效果好的基于EEMD-Hilbert变换的ECG去噪方法。

本发明解决上述技术问题所采用的技术方案为：一种基于EEMD-Hilbert变换的ECG去噪方法，包括以下步骤：

步骤1.采集ECG信号，并将采集到的ECG信号记为ECG_noise，其中，ECG信号的采样频率为f，f等于250Hz，ECG_noise包括n个样本点，即ECG_noise的信号长度为n，n为大于等于1800的整数，任意相邻两个采样时刻之间的间隔为1/f秒，采样时间为n/f秒，将ECG_noise中第L个样本点的采样时刻记为t_L，t_L等于L/f，L＝1，2，…，n；

步骤2.采用集成经验模态分解(Ensemble Empirical Mode Decomposition,EEMD)方法对ECG_noise进行模态分解，得到由多个固有模态函数构成的一组固有模态函数(Intrinsic Mode Function，IMF)，其中采用集成经验模态分解方法对ECG_noise进行模态分解的具体过程为：

S1、先生成100组均值与方差均为0的高斯白噪声序列，且每组高斯白噪声序列所包含的样本点数量与ECG_noise包含的样本点数量相同，将生成的第i组高斯白噪声序列记为w⁽ⁱ⁾，i＝1,2,....,100；

S2、构造x⁽ⁱ⁾＝ECG_noise+βw⁽ⁱ⁾，其中，x⁽ⁱ⁾为在ECG_noise中加入第i组高斯白噪声序列后所构造出的信号，将其称为第i个构造信号，β为信噪比控制系数，β＝0.05；

S3、设定模态分解次数变量为k，对k进行初始化，令k＝0；

S4、进行第k+1次模态分解，具体过程为：

S4-1、将第k次模态分解得到的第i个残差r_k ⁽ⁱ⁾作为第k+1次模态分解的第i个输入信号，计算出第i个输入信号r_k ⁽ⁱ⁾的所有极大值点和极小值点，其中，当k＝0时，r_k ⁽ⁱ⁾＝x⁽ⁱ⁾，将第i个输入信号r_k ⁽ⁱ⁾中每相邻两个极大值点均作为一组极大值点，分别采用三次样条插值法在第i个输入信号r_k ⁽ⁱ⁾的每组极大值点中两个极大值点之间进行插值，然后将第i个输入信号r_k ⁽ⁱ⁾的每个极大值点与此时得到的所有插值依次相连得到第i个输入信号r_k ⁽ⁱ⁾的上包络线e_max-(k+1) ⁽ⁱ⁾，同理，将第i个输入信号r_k ⁽ⁱ⁾的每相邻两个极小值点作为一组极小值点，分别采用三次样条插值法在第i个输入信号r_k ⁽ⁱ⁾的每组极小值点中两个极小值点之间进行插值，然后将每个极小值点与此时得到的所有插值依次相连得到第i个输入信号r_k ⁽ⁱ⁾的下包络线e_min-(k+1) ⁽ⁱ⁾；

S4-2、采用公式(1)计算第i个输入信号r_k ⁽ⁱ⁾的平均包络线e_m-(k+1) ⁽ⁱ⁾：

e_m-(k+1) ⁽ⁱ⁾＝(e_max-(k+1) ⁽ⁱ⁾+e_min-(k+1) ⁽ⁱ⁾)/2 (1)

S4-3、设定第k+1代需进一步判断的第i组候选固有模态函数，将其记为d_k+1 ⁽ⁱ⁾，采用公式(2)计算得到d_k+1 ⁽ⁱ⁾：

d_k+1 ⁽ⁱ⁾＝r_k ⁽ⁱ⁾-e_m-(k+1) ⁽ⁱ⁾ (2)

S4-4、判断当前得到的d_k+1 ⁽ⁱ⁾是否为固有模态函数，具体判断过程为：先计算d_k+1 ⁽ⁱ⁾所有极值点数量与过零点数量，然后判断以下两个条件是否同时成立：①d_k+1 ⁽ⁱ⁾的极值点与过零点的数量相同或至多相差一个；②d_k+1 ⁽ⁱ⁾的极大值不存在负数，d_k+1 ⁽ⁱ⁾的极小值不存在正数；如果上述两个条件同时成立，则保存当前d_k+1 ⁽ⁱ⁾，然后通过r_k+1 ⁽ⁱ⁾＝r_k ⁽ⁱ⁾-d_k+1 ⁽ⁱ⁾计算得到第k+1次模态分解得到的第i个残差r_k+1 ⁽ⁱ⁾，并进入步骤S4-5，其中计算时d_k+1 ⁽ⁱ⁾的取值为其当前最新值，如果上述两个条件不能同时成立，则采用当前得到的d_k+1 ⁽ⁱ⁾更新r_k ⁽ⁱ⁾，此时第k次模态分解得到的第i个残差r_k ⁽ⁱ⁾被更新，然后基于更新后的r_k ⁽ⁱ⁾重复步骤S4-1至步骤S4-4，直至计算得到第k+1次模态分解得到的第i个残差r_k+1 ⁽ⁱ⁾；

S4-5、将第k+1次模态分解得到的固有模态函数记为IMF_k+1，采用式(3)计算得到IMF_k+1：

式(3)中，d_k+1 ⁽ⁱ⁾为当前得到的d_k+1 ⁽ⁱ⁾；

S4-6、判断第k+1次模态分解得到的第1个残差r_k+1 ⁽¹⁾至第100个残差r_k+1 ⁽¹⁰⁰⁾是否都满足以下条件：同时具有2个以上极值点以及具有至少1个过零点；如果不是都满足该条件，则模态分解结束，统计当前得到的所有固有模态函数的数量，将其记为y，进入步骤3，如果都满足该条件，则采用k的当前值加1的和更新k的取值，然后返回步骤S4进行下一次模态分解；

步骤3.采用式(4)分别对步骤2得到的y个固有模态函数进行Hilbert变换：

其中，IMF_j(t_L)为第j个固有模态函数IMF_j对应于ECG_noise中第t个样本点采样时刻t_L的幅值，*代表卷积运算，H[·]为Hilbert变换函数，π表示圆周率；

步骤4.使用式(5)所示相位函数计算第j个固有模态函数IMF_j对应于ECG_noise中第t个样本点采样时刻t_L的相位θ_j(t_L)：

此时，得到每个固有模态函数对应于ECG_noise中各样本点采样时刻的相位，其中每个固有模态函数中对应于ECG_noise中所有样本点采样时刻的相位构成该固有模态函数的一组相位；

步骤5.由瞬时频率的定义得到第j个固有模态函数IMF_j对应于ECG_noise中第t个样本点采样时刻t_L的瞬时频率ω_j(t_L)：

此时，得到每个固有模态函数中对应于ECG_noise中各样本点采样时刻的瞬时频率，其中每个固有模态函数中对应于ECG_noise中所有样本点采样时刻的瞬时频率构成该固有模态函数的一组瞬时频率，即得到y组瞬时频率；

步骤6.分别获取每组瞬时频率的中心瞬时频率(即该组瞬时频率的中位数)，并根据中心瞬时频率对每个固有模态函数进行分类：当某组瞬时频率的中心瞬时频率小于1Hz时，该组瞬时频率对应的固有模态函数被判别为基线漂移噪声和肌电干扰噪声为主的固有模态函数，属于第1类，当某组瞬时频率的中心瞬时频率在1Hz～40Hz之间时，该组瞬时频率对应的固有模态函数被判别为信号为主的IMF固有模态函数，属于第2类，当某组瞬时频率的中心瞬时频率大于40Hz时，该组瞬时频率对应的固有模态函数被判别为高频噪声为主的固有模态函数，属于第3类；

步骤7.设定包含n个样本点的信号u，将属于第2类的所有固有模态函数的所有样本点幅值一一对应相加后得到的值一一对应作为信号u的n个样本点的幅值，对信号u采用非局部均值(Non-Local Means,NLM)滤波方法去噪，得到第一个去噪信号

步骤8.设定包含n个样本点的信号g，将属于第3类的所有固有模态函数的所有样本点幅值一一对应相加后得到的值一一对应作为信号g的n个样本点的幅值，对信号g采用阈值方法进行去噪，得到第二个去噪信号

步骤9.将第一个去噪信号

与第二个去噪信号

重建得到最终的去噪后的ECG信号ECG_d，重建过程表示为：

所述的步骤8中对信号u采用非局部均值(Non-Local Means,NLM)滤波方法去噪，得到第一个去噪信号

的具体过程为：

A、先对信号u前后进行对称填充，得到填充后信号：在信号u初始样本点前填充10个样本点，这10个样本点为信号u的前10个样本点的镜像点，在信号u末尾样本点后填充10个样本点，这10个样本点为信号u的最后10个样本点的镜像点，将填充后信号记为

B、采用式(7)进行去噪：

式(7)中，

为信号u的第v个样本点去噪后的幅值，即填充后信号

的第v+10个样本点去噪后的幅值，v＝1,2,...,n，n为信号u的长度；构造出以填充后信号

的第v+10个样本点为中心再取中心前后各10个样本点的幅值组成一个包含21个样本点的幅值的滤波块，再构造出以填充后信号

的第v+10个样本点的幅值为中心再取中心前后多个样本点的幅值组成一个搜索邻域，将该搜索邻域的范围记为[r₁,r₂]，搜索邻域按照以下规则确定：当v-190≤0且v+210≤n时，r₁取v+10，r₂取v+210；当v-190≤0且v+210>n时，r₁取v+10，r₂取n；当v-190＞0且v+210≤n时，r₁取v-190，r₂取v+210；当v-190＞0且v+210>n时，r₁取v-190，r₂取n；将搜索邻域包含的样本点幅值的总数量记为q，设遍历搜索邻域过程中，第r₁个样本点的幅值至第r₂个样本点的幅值中，除第v+10个样本点的幅值以外的某个样本点的幅值为第p个样本点的幅值，即p＝r₁，r₁+1，...,，r₂，且p≠v+10，分别以第p个样本点的幅值为中心，取中心前后各10个样本点的幅值，在搜索邻域中构造q-1个包含21个样本点的幅值的对比滤波块，Z(v)表示以填充后信号

的第v+10个样本点的幅值为中心所构成滤波块与在搜索邻域内构造出的q-1个对比滤波块的相似滤波权重之和，

为填充后信号

中第p个样本点的幅值，W(v,p)为以填充后信号

的第v+10个样本点的幅值为中心的滤波块与搜索邻域内第p个样本点的幅值为中心的对比滤波块的相似滤波权重，W(v,p)表示为：

其中，exp(·)为指数函数；L_Δ为信号u以第v个样本点的幅值为中心对应的滤波块的大小，取值为21；λ为平滑滤波参数，取值为0.5σ，σ为信号的噪声标准差，采用公式

计算得到，mean(u)表示信号u中所有的样本点的幅值的平均值，u(m)为信号u的第m个样本点的幅值，d²(v,p,ε)为填充后信号

中以第v+10个样本点的幅值为中心再取中心前后各10个样本点的幅值组成一个包含21个样本点的幅值的滤波块和搜索邻域内的第p个样本点的幅值为中心再取中心前后各10个样本点的幅值组成一个包含21个样本点的幅值的对比滤波块的差值的平方，d²(v,p,ε)通过式

计算得到，其中，

为填充后信号

的第v+ε+10样本点的幅值，

为填充后信号

的第p+ε样本点的幅值，ε＝-10,-9,...,9,10；

C.将得到的

至

一一对应作为n个样本点的幅值，采用这n个样本点构成的信号即为第一个去噪信号

所述的步骤9中对信号g采用阈值方法进行去噪，得到第二个去噪信号

的具体过程为：

a.首先计算信号g的噪声标准差γ:

其中，median(·)为取中位数函数；g'为信号g中n个样本点的幅值绝对值的集合。

b.根据噪声标准差计算阈值T：

c.采用阈值函数对信号g进行去噪，阈值函数为：

其中，sgn(·)为符号函数，ln(·)为对数函数，b为阈值函数可调整参数，取0.1，g(a)为信号g的第a个样本点的幅值，||为取绝对值符号，a＝1,2，…，n；

d.将得到的

至

作为n个样本点的幅值，采用这n个样本点构成的信号即为第二个去噪信号

与现有技术相比，本发明的优点在于通过先采用集成经验模态分解方法对原始ECG_noise进行模态分解，得到由多个固有模态函数构成的一组固有模态函数，然后对得到的一组固有模态函数进行Hilbert变换，得到每个固有模态函数的Hilbert变换函数，基于每个固有模态函数的Hilbert变换函数确定每个固有模态函数对应于ECG_noise中每个样本点采样时刻的相位，继而得到每个固有模态函数对应于ECG_noise中每个样本点采样时刻的瞬时频率，采用每个固有模态函数中对应于ECG_noise中所有样本点采样时刻的瞬时频率构成该固有模态函数的一组瞬时频率，即得到每个固有模态函数对应的一组瞬时频率，分别获取每组瞬时频率的中心瞬时频率，并根据中心瞬时频率对每个固有模态函数进行分类，确定以基线漂移噪声为主固有模态函数、以信号为主的固有模态函数以及以高频噪声为主的固有模态函数，然后将以基线漂移噪声为主的固有模态函数舍弃，而由于肌电干扰噪声其频谱范围为20～2000Hz，频谱覆盖较广，被判断为以信号为主的固有模态函数的特点在于以有用的ECG信号为主导的，能够体现出ECG信号的形态特征但仍然包含少部分的与ECG信号频谱范围重叠的肌电干扰噪声，为了去除这少部分与ECG信号频谱范围重叠的肌电干扰噪声，将以信号为主的固有模态函数采用非局部均值滤波方法去噪，在保留ECG信号的形态特征的同时，去除肌电干扰噪声，以高频噪声为主的固有模态函数的特点在于以不与ECG信号频谱范围重叠的肌电干扰噪声为主但包含少部分的ECG信号信息，为了在肌电干扰噪声中提取出这少部分的ECG信号信息，将以高频噪声为主的固有模态函数采用阈值方法去噪，去除肌电干扰噪声，最终将上述两种方法得到的去噪结果重构得到最终去噪后的ECG，由此本发明能够对基线漂移噪声和肌电干扰噪声进行更加有效的降噪，通过对真实采集到的ECG信号进行方法验证，得到的最终去噪后的ECG信号仍然保留着明显的ECG形态特征，去噪效果好。

附图说明

图1为本发明的基于EEMD-Hilbert变换的ECG去噪方法的流程图；

图2为本发明实施例的第一组实验结果图；

图3为本发明实施例的第二组实验结果图；

具体实施方式

以下结合附图实施例对本发明作进一步详细描述。

实施例：如图1所示，一种基于EEMD-Hilbert变换的ECG去噪方法，包括以下步骤：

步骤1.采集ECG信号，并将采集到的ECG信号记为ECG_noise，其中，ECG信号的采样频率为f，f等于250Hz，ECG_noise包括n个样本点，即ECG_noise的信号长度为n，n为大于等于1800的整数，任意相邻两个采样时刻之间的间隔为1/f秒，采样时间为n/f秒，将ECG_noise中第L个样本点的采样时刻记为t_L，t_L等于L/f，L＝1，2，…，n；

步骤2.采用集成经验模态分解(Ensemble Empirical Mode Decomposition,EEMD)方法对ECG_noise进行模态分解，得到由多个固有模态函数构成的一组固有模态函数(Intrinsic Mode Function，IMF)，其中采用集成经验模态分解方法对ECG_noise进行模态分解的具体过程为：

S1、先生成100组均值与方差均为0的高斯白噪声序列，且每组高斯白噪声序列所包含的样本点数量与ECG_noise包含的样本点数量相同，将生成的第i组高斯白噪声序列记为w⁽ⁱ⁾，i＝1,2,....,100；

S2、构造x⁽ⁱ⁾＝ECG_noise+βw⁽ⁱ⁾，其中，x⁽ⁱ⁾为在ECG_noise中加入第i组高斯白噪声序列后所构造出的信号，将其称为第i个构造信号，β为信噪比控制系数，β＝0.05；

S3、设定模态分解次数变量为k，对k进行初始化，令k＝0；

S4、进行第k+1次模态分解，具体过程为：

S4-1、将第k次模态分解得到的第i个残差r_k ⁽ⁱ⁾作为第k+1次模态分解的第i个输入信号，计算出第i个输入信号r_k ⁽ⁱ⁾的所有极大值点和极小值点，其中，当k＝0时，r_k ⁽ⁱ⁾＝x⁽ⁱ⁾，将第i个输入信号r_k ⁽ⁱ⁾中每相邻两个极大值点均作为一组极大值点，分别采用三次样条插值法在第i个输入信号r_k ⁽ⁱ⁾的每组极大值点中两个极大值点之间进行插值，然后将第i个输入信号r_k ⁽ⁱ⁾的每个极大值点与此时得到的所有插值依次相连得到第i个输入信号r_k ⁽ⁱ⁾的上包络线e_max-(k+1) ⁽ⁱ⁾，同理，将第i个输入信号r_k ⁽ⁱ⁾的每相邻两个极小值点作为一组极小值点，分别采用三次样条插值法在第i个输入信号r_k ⁽ⁱ⁾的每组极小值点中两个极小值点之间进行插值，然后将每个极小值点与此时得到的所有插值依次相连得到第i个输入信号r_k ⁽ⁱ⁾的下包络线e_min-(k+1) ⁽ⁱ⁾；

S4-2、采用公式(1)计算第i个输入信号r_k ⁽ⁱ⁾的平均包络线e_m-(k+1) ⁽ⁱ⁾：

e_m-(k+1) ⁽ⁱ⁾＝(e_max-(k+1) ⁽ⁱ⁾+e_min-(k+1) ⁽ⁱ⁾)/2 (1)

S4-3、设定第k+1代需进一步判断的第i组候选固有模态函数，将其记为d_k+1 ⁽ⁱ⁾，采用公式(2)计算得到d_k+1 ⁽ⁱ⁾：

d_k+1 ⁽ⁱ⁾＝r_k ⁽ⁱ⁾-e_m-(k+1) ⁽ⁱ⁾ (2)

S4-4、判断当前得到的d_k+1 ⁽ⁱ⁾是否为固有模态函数，具体判断过程为：先计算d_k+1 ⁽ⁱ⁾所有极值点数量与过零点数量，然后判断以下两个条件是否同时成立：①d_k+1 ⁽ⁱ⁾的极值点与过零点的数量相同或至多相差一个；②d_k+1 ⁽ⁱ⁾的极大值不存在负数，d_k+1 ⁽ⁱ⁾的极小值不存在正数；如果上述两个条件同时成立，则保存当前d_k+1 ⁽ⁱ⁾，然后通过r_k+1 ⁽ⁱ⁾＝r_k ⁽ⁱ⁾-d_k+1 ⁽ⁱ⁾计算得到第k+1次模态分解得到的第i个残差r_k+1 ⁽ⁱ⁾，并进入步骤S4-5，其中计算时d_k+1 ⁽ⁱ⁾的取值为其当前最新值，如果上述两个条件不能同时成立，则采用当前得到的d_k+1 ⁽ⁱ⁾更新r_k ⁽ⁱ⁾，此时第k次模态分解得到的第i个残差r_k ⁽ⁱ⁾被更新，然后基于更新后的r_k ⁽ⁱ⁾重复步骤S4-1至步骤S4-4，直至计算得到第k+1次模态分解得到的第i个残差r_k+1 ⁽ⁱ⁾；

S4-5、将第k+1次模态分解得到的固有模态函数记为IMF_k+1，采用式(3)计算得到IMF_k+1：

式(3)中，d_k+1 ⁽ⁱ⁾为当前得到的d_k+1 ⁽ⁱ⁾；

S4-6、判断第k+1次模态分解得到的第1个残差r_k+1 ⁽¹⁾至第100个残差r_k+1 ⁽¹⁰⁰⁾是否都满足以下条件：同时具有2个以上极值点以及具有至少1个过零点；如果不是都满足该条件，则模态分解结束，统计当前得到的所有固有模态函数的数量，将其记为y，进入步骤3，如果都满足该条件，则采用k的当前值加1的和更新k的取值，然后返回步骤S4进行下一次模态分解；

步骤3.采用式(4)分别对步骤2得到的y个固有模态函数进行Hilbert变换：

其中，IMF_j(t_L)为第j个固有模态函数IMF_j对应于ECG_noise中第t个样本点采样时刻t_L的幅值，*代表卷积运算，H[·]为Hilbert变换函数，π表示圆周率；

步骤4.使用式(5)所示相位函数计算第j个固有模态函数IMF_j对应于ECG_noise中第t个样本点采样时刻t_L的相位θ_j(t_L)：

此时，得到每个固有模态函数对应于ECG_noise中各样本点采样时刻的相位，其中每个固有模态函数中对应于ECG_noise中所有样本点采样时刻的相位构成该固有模态函数的一组相位；

步骤5.由瞬时频率的定义得到第j个固有模态函数IMF_j对应于ECG_noise中第t个样本点采样时刻t_L的瞬时频率ω_j(t_L)：

此时，得到每个固有模态函数中对应于ECG_noise中各样本点采样时刻的瞬时频率，其中每个固有模态函数中对应于ECG_noise中所有样本点采样时刻的瞬时频率构成该固有模态函数的一组瞬时频率，即得到y组瞬时频率；

步骤6.分别获取每组瞬时频率的中心瞬时频率(即该组瞬时频率的中位数)，并根据中心瞬时频率对每个固有模态函数进行分类：当某组瞬时频率的中心瞬时频率小于1Hz时，该组瞬时频率对应的固有模态函数被判别为基线漂移噪声为主的固有模态函数，属于第1类，当某组瞬时频率的中心瞬时频率在1Hz～40Hz之间时，该组瞬时频率对应的固有模态函数被判别为信号为主的IMF固有模态函数，属于第2类，当某组瞬时频率的中心瞬时频率大于40Hz时，该组瞬时频率对应的固有模态函数被判别为高频噪声为主的固有模态函数，属于第3类；

步骤7.设定包含n个样本点的信号u，将属于第2类的所有固有模态函数的所有样本点幅值一一对应相加后得到的值一一对应作为信号u的n个样本点的幅值，对信号u采用非局部均值(Non-Local Means,NLM)滤波方法去噪，得到第一个去噪信号

步骤8.设定包含n个样本点的信号g，将属于第3类的所有固有模态函数的所有样本点幅值一一对应相加后得到的值一一对应作为信号g的n个样本点的幅值，对信号g采用阈值方法进行去噪，得到第二个去噪信号

步骤9.将第一个去噪信号

与第二个去噪信号

重建得到最终的去噪后的ECG信号ECG_d，重建过程表示为：

本实施例中，步骤8中对信号u采用非局部均值(Non-Local Means,NLM)滤波方法去噪，得到第一个去噪信号

的具体过程为：

A、先对信号u前后进行对称填充，得到填充后信号：在信号u初始样本点前填充10个样本点，这10个样本点为信号u的前10个样本点的镜像点，在信号u末尾样本点后填充10个样本点，这10个样本点为信号u的最后10个样本点的镜像点，将填充后信号记为

B、采用式(7)进行去噪：

式(7)中，

为信号u的第v个样本点去噪后的幅值，即填充后信号

的第v+10个样本点去噪后的幅值，v＝1,2,...,n，n为信号u的长度；构造出以填充后信号

的第v+10个样本点为中心再取中心前后各10个样本点的幅值组成一个包含21个样本点的幅值的滤波块，再构造出以填充后信号

的第v+10个样本点的幅值为中心再取中心前后多个样本点的幅值组成一个搜索邻域，将该搜索邻域的范围记为[r₁,r₂]，搜索邻域按照以下规则确定：当v-190≤0且v+210≤n时，r₁取v+10，r₂取v+210；当v-190≤0且v+210>n时，r₁取v+10，r₂取n；当v-190＞0且v+210≤n时，r₁取v-190，r₂取v+210；当v-190＞0且v+210>n时，r₁取v-190，r₂取n；将搜索邻域包含的样本点幅值的总数量记为q，设遍历搜索邻域过程中，第r₁个样本点的幅值至第r₂个样本点的幅值中，除第v+10个样本点的幅值以外的某个样本点的幅值为第p个样本点的幅值，即p＝r₁，r₁+1，...,，r₂，且p≠v+10，分别以第p个样本点的幅值为中心，取中心前后各10个样本点的幅值，在搜索邻域中构造q-1个包含21个样本点的幅值的对比滤波块，Z(v)表示以填充后信号

的第v+10个样本点的幅值为中心所构成滤波块与在搜索邻域内构造出的q-1个对比滤波块的相似滤波权重之和，

为填充后信号

中第p个样本点的幅值，W(v,p)为以填充后信号

的第v+10个样本点的幅值为中心的滤波块与搜索邻域内第p个样本点的幅值为中心的对比滤波块的相似滤波权重，W(v,p)表示为：

其中exp(·)为指数函数；L_Δ为信号u以第v个样本点的幅值为中心对应的滤波块的大小，取值为21；λ为平滑滤波参数，取值为0.5σ，σ为信号的噪声标准差，采用公式

计算得到，mean(u)表示信号u中所有的样本点的幅值的平均值，u(m)为信号u的第m个样本点的幅值，d²(v,p,ε)为填充后信号

中以第v+10个样本点的幅值为中心再取中心前后各10个样本点的幅值组成一个包含21个样本点的幅值的滤波块和搜索邻域内的第p个样本点的幅值为中心再取中心前后各10个样本点的幅值组成一个包含21个样本点的幅值的对比滤波块的差值的平方，d²(v,p,ε)通过式

计算得到，其中，

为填充后信号

的第v+ε+10样本点的幅值，

为填充后信号

的第p+ε样本点的幅值，ε＝-10,-9,...,9,10；

C.将得到的

至

一一对应作为n个样本点的幅值，采用这n个样本点构成的信号即为第一个去噪信号

本实施例中，步骤9中对信号g采用阈值方法进行去噪，得到第二个去噪信号

的具体过程为：

a.首先计算信号g的噪声标准差γ:

其中，median(·)为取中位数函数；g'为信号g中n个样本点的幅值绝对值的集合。

b.根据噪声标准差计算阈值T：

c.采用阈值函数对信号g进行去噪，阈值函数为：

其中，sgn(·)为符号函数，ln(·)为对数函数，b为阈值函数可调整参数，取0.1，g(a)为信号g的第a个样本点的幅值，||为取绝对值符号，a＝1,2，…，n；

d.将得到的

至

作为n个样本点的幅值，采用这n个样本点构成的信号即为第二个去噪信号

为了证明本发明的基于EEMD-Hilbert变换的ECG去噪方法的优异性能，以下设置两组实验来进行验证。第一组实验中通过几乎无噪声的ECG信号上添加混合噪声来形成ECG_noise，然后分别采用本发明方法和现有的小波阈值去噪方法进行去噪，实验对比结果如图2所示。图2中，子图(a)为几乎无噪声的ECG信号的示意图；子图(b)为几乎无噪声的ECG信号上添加混合噪声来形成ECG_noise的示意图；子图(c)是采用现有的小波阈值去噪方法对ECG_noise进行去噪得到的去噪效果图，子图(d)为采用本发明的方法对ECG_noise进行去噪得到的去噪效果图。分析图2可知本发明的去噪方法相对于小波阈值去噪方法，具有较好的去噪效果，去噪的ECG信号中ECG的形态特征得到完整的保留。第二组实验中通过基于AD8232心电信号传感器与STM32实现的心电信号采集平台采集得到的ECG信号作为ECG_noise，然后分别采用本发明方法和现有的小波阈值去噪方法进行去噪，实验对比结果如图3所示。图3中，子图(a)为采集的ECG信号的示意图，子图(b)为采用现有的小波阈值去噪方法对ECG_noise进行去噪得到的去噪效果图，子图(c)为采用本发明的方法对ECG_noise进行去噪得到的去噪效果图。分析图3可知，对于真实采集的ECG信号，本发明相对于现有的相对于小波阈值去噪方法，仍然具有较好的去噪效果，去噪的ECG信号中ECG的形态特征得到完整的保留。

Claims (3)

1.一种基于EEMD-Hilbert变换的ECG去噪方法，其特征在于包括以下步骤：

步骤1.采集ECG信号，并将采集到的ECG信号记为ECG_noise，其中，ECG信号的采样频率为f，f等于250Hz，ECG_noise包括n个样本点，即ECG_noise的信号长度为n，n为大于等于1800的整数，任意相邻两个采样时刻之间的间隔为1/f秒，采样时间为n/f秒，将ECG_noise中第L个样本点的采样时刻记为t_L，t_L等于L/f，L＝1，2，…，n；

步骤2.采用集成经验模态分解(Ensemble Empirical Mode Decomposition,EEMD)方法对ECG_noise进行模态分解，得到由多个固有模态函数构成的一组固有模态函数(IntrinsicMode Function，IMF)，其中采用集成经验模态分解方法对ECG_noise进行模态分解的具体过程为：

S1、先生成100组均值与方差均为0的高斯白噪声序列，且每组高斯白噪声序列所包含的样本点数量与ECG_noise包含的样本点数量相同，将生成的第i组高斯白噪声序列记为w⁽ⁱ⁾，i＝1,2,....,100；

S2、构造x⁽ⁱ⁾＝ECG_noise+βw⁽ⁱ⁾，其中，x⁽ⁱ⁾为在ECG_noise中加入第i组高斯白噪声序列后所构造出的信号，将其称为第i个构造信号，β为信噪比控制系数，β＝0.05；

S3、设定模态分解次数变量为k，对k进行初始化，令k＝0；

S4、进行第k+1次模态分解，具体过程为：

S4-1、将第k次模态分解得到的第i个残差r_k ⁽ⁱ⁾作为第k+1次模态分解的第i个输入信号，计算出第i个输入信号r_k ⁽ⁱ⁾的所有极大值点和极小值点，其中，当k＝0时，r_k ⁽ⁱ⁾＝x⁽ⁱ⁾，将第i个输入信号r_k ⁽ⁱ⁾中每相邻两个极大值点均作为一组极大值点，分别采用三次样条插值法在第i个输入信号r_k ⁽ⁱ⁾的每组极大值点中两个极大值点之间进行插值，然后将第i个输入信号r_k ⁽ⁱ⁾的每个极大值点与此时得到的所有插值依次相连得到第i个输入信号r_k ⁽ⁱ⁾的上包络线e_max-(k+1) ⁽ⁱ⁾，同理，将第i个输入信号r_k ⁽ⁱ⁾的每相邻两个极小值点作为一组极小值点，分别采用三次样条插值法在第i个输入信号r_k ⁽ⁱ⁾的每组极小值点中两个极小值点之间进行插值，然后将每个极小值点与此时得到的所有插值依次相连得到第i个输入信号r_k ⁽ⁱ⁾的下包络线e_min-(k+1) ⁽ⁱ⁾；

S4-2、采用公式(1)计算第i个输入信号r_k ⁽ⁱ⁾的平均包络线e_m-(k+1) ⁽ⁱ⁾：

e_m-(k+1) ⁽ⁱ⁾＝(e_max-(k+1) ⁽ⁱ⁾+e_min-(k+1) ⁽ⁱ⁾)/2 (1)

S4-3、设定第k+1代需进一步判断的第i组候选固有模态函数，将其记为d_k+1 ⁽ⁱ⁾，采用公式(2)计算得到d_k+1 ⁽ⁱ⁾：

d_k+1 ⁽ⁱ⁾＝r_k ⁽ⁱ⁾-e_m-(k+1) ⁽ⁱ⁾ (2)

S4-4、判断当前得到的d_k+1 ⁽ⁱ⁾是否为固有模态函数，具体判断过程为：先计算d_k+1 ⁽ⁱ⁾所有极值点数量与过零点数量，然后判断以下两个条件是否同时成立：①d_k+1 ⁽ⁱ⁾的极值点与过零点的数量相同或至多相差一个；②d_k+1 ⁽ⁱ⁾的极大值不存在负数，d_k+1 ⁽ⁱ⁾的极小值不存在正数；如果上述两个条件同时成立，则保存当前d_k+1 ⁽ⁱ⁾，然后通过r_k+1 ⁽ⁱ⁾＝r_k ⁽ⁱ⁾-d_k+1 ⁽ⁱ⁾计算得到第k+1次模态分解得到的第i个残差r_k+1 ⁽ⁱ⁾，并进入步骤S4-5，其中计算时d_k+1 ⁽ⁱ⁾的取值为其当前最新值，如果上述两个条件不能同时成立，则采用当前得到的d_k+1 ⁽ⁱ⁾更新r_k ⁽ⁱ⁾，此时第k次模态分解得到的第i个残差r_k ⁽ⁱ⁾被更新，然后基于更新后的r_k ⁽ⁱ⁾重复步骤S4-1至步骤S4-4，直至计算得到第k+1次模态分解得到的第i个残差r_k+1 ⁽ⁱ⁾；

S4-5、将第k+1次模态分解得到的固有模态函数记为IMF_k+1，采用式(3)计算得到IMF_k+1：

式(3)中，d_k+1 ⁽ⁱ⁾为当前得到的d_k+1 ⁽ⁱ⁾；

S4-6、判断第k+1次模态分解得到的第1个残差r_k+1 ⁽¹⁾至第100个残差r_k+1 ⁽¹⁰⁰⁾是否都满足以下条件：同时具有2个以上极值点以及具有至少1个过零点；如果不是都满足该条件，则模态分解结束，统计当前得到的所有固有模态函数的数量，将其记为y，进入步骤3，如果都满足该条件，则采用k的当前值加1的和更新k的取值，然后返回步骤S4进行下一次模态分解；

步骤3.采用式(4)分别对步骤2得到的y个固有模态函数进行Hilbert变换：

其中，IMF_j(t_L)为第j个固有模态函数IMF_j对应于ECG_noise中第t个样本点采样时刻t_L的幅值，*代表卷积运算，H[·]为Hilbert变换函数，π表示圆周率；

步骤4.使用式(5)所示相位函数计算第j个固有模态函数IMF_j对应于ECG_noise中第t个样本点采样时刻t_L的相位θ_j(t_L)：

此时，得到每个固有模态函数对应于ECG_noise中各样本点采样时刻的相位，其中每个固有模态函数中对应于ECG_noise中所有样本点采样时刻的相位构成该固有模态函数的一组相位；

步骤5.由瞬时频率的定义得到第j个固有模态函数IMF_j对应于ECG_noise中第t个样本点采样时刻t_L的瞬时频率ω_j(t_L)：

式(6)中，

表示θ_j(t_L)对t_L求导数，此时，得到每个固有模态函数中对应于ECG_noise中各样本点采样时刻的瞬时频率，其中每个固有模态函数中对应于ECG_noise中所有样本点采样时刻的瞬时频率构成该固有模态函数的一组瞬时频率，即得到y组瞬时频率；

步骤6.分别获取每组瞬时频率的中心瞬时频率(即该组瞬时频率的中位数)，并根据中心瞬时频率对每个固有模态函数进行分类：当某组瞬时频率的中心瞬时频率小于1Hz时，该组瞬时频率对应的固有模态函数被判别为基线漂移噪声为主的固有模态函数，属于第1类，当某组瞬时频率的中心瞬时频率在1Hz～40Hz之间时，该组瞬时频率对应的固有模态函数被判别为信号为主的IMF固有模态函数，属于第2类，当某组瞬时频率的中心瞬时频率大于40Hz时，该组瞬时频率对应的固有模态函数被判别为高频噪声为主的固有模态函数，属于第3类；

步骤7.设定包含n个样本点的信号u，将属于第2类的所有固有模态函数的所有样本点幅值一一对应相加后得到的值一一对应作为信号u的n个样本点的幅值，对信号u采用非局部均值(Non-Local Means,NLM)滤波方法去噪，得到第一个去噪信号

步骤8.设定包含n个样本点的信号g，将属于第3类的所有固有模态函数的所有样本点幅值一一对应相加后得到的值一一对应作为信号g的n个样本点的幅值，对信号g采用阈值方法进行去噪，得到第二个去噪信号

步骤9.将第一个去噪信号

与第二个去噪信号

重建得到最终的去噪后的ECG信号ECG_d，重建过程表示为：

2.根据权利要求1所述的一种基于EEMD-Hilbert变换的ECG去噪方法，其特征在于所述的步骤8中对信号u采用非局部均值(Non-Local Means,NLM)滤波方法去噪，得到第一个去噪信号

的具体过程为：

A、先对信号u前后进行对称填充，得到填充后信号：在信号u初始样本点前填充10个样本点，这10个样本点为信号u的前10个样本点的镜像点，在信号u末尾样本点后填充10个样本点，这10个样本点为信号u的最后10个样本点的镜像点，将填充后信号记为

B、采用式(7)进行去噪：

式(7)中，

为信号u的第v个样本点去噪后的幅值，即填充后信号

的第v+10个样本点去噪后的幅值，v＝1,2,...,n，n为信号u的长度；构造出以填充后信号

的第v+10个样本点为中心再取中心前后各10个样本点的幅值组成一个包含21个样本点的幅值的滤波块，再构造出以填充后信号

的第v+10个样本点的幅值为中心再取中心前后多个样本点的幅值组成一个搜索邻域，将该搜索邻域的范围记为[r₁,r₂]，搜索邻域按照以下规则确定：当v-190≤0且v+210≤n时，r₁取v+10，r₂取v+210；当v-190≤0且v+210>n时，r₁取v+10，r₂取n；当v-190＞0且v+210≤n时，r₁取v-190，r₂取v+210；当v-190＞0且v+210>n时，r₁取v-190，r₂取n；将搜索邻域包含的样本点幅值的总数量记为q，设遍历搜索邻域过程中，第r₁个样本点的幅值至第r₂个样本点的幅值中，除第v+10个样本点的幅值以外的某个样本点的幅值为第p个样本点的幅值，即p＝r₁，r₁+1，...,，r₂，且p≠v+10，分别以第p个样本点的幅值为中心，取中心前后各10个样本点的幅值，在搜索邻域中构造q-1个包含21个样本点的幅值的对比滤波块，Z(v)表示以填充后信号

的第v+10个样本点的幅值为中心所构成滤波块与在搜索邻域内构造出的q-1个对比滤波块的相似滤波权重之和，

为填充后信号

中第p个样本点的幅值，W(v,p)为以填充后信号

的第v+10个样本点的幅值为中心的滤波块与搜索邻域内第p个样本点的幅值为中心的对比滤波块的相似滤波权重，W(v,p)表示为：

其中exp(·)为指数函数；L_Δ为信号u以第v个样本点的幅值为中心对应的滤波块的大小，取值为21；λ为平滑滤波参数，取值为0.5σ，σ为信号的噪声标准差，采用公式

计算得到，mean(u)表示信号u中所有的样本点的幅值的平均值，u(m)为信号u的第m个样本点的幅值，d²(v,p,ε)为填充后信号

中以第v+10个样本点的幅值为中心再取中心前后各10个样本点的幅值组成一个包含21个样本点的幅值的滤波块和搜索邻域内的第p个样本点的幅值为中心再取中心前后各10个样本点的幅值组成一个包含21个样本点的幅值的对比滤波块的差值的平方，d²(v,p,ε)通过式

计算得到，其中，

为填充后信号

的第v+ε+10样本点的幅值，

为填充后信号

的第p+ε样本点的幅值，ε＝-10,-9,...,9,10；

C.将得到的

至

一一对应作为n个样本点的幅值，采用这n个样本点构成的信号即为第一个去噪信号

3.根据权利要求1所述的一种基于EEMD-Hilbert变换的ECG去噪方法，其特征在于所述的步骤9中对信号g采用阈值方法进行去噪，得到第二个去噪信号

的具体过程为：

a.首先计算信号g的噪声标准差γ:

其中，median(·)为取中位数函数；g'为信号g中n个样本点的幅值绝对值的集合。

b.根据噪声标准差计算阈值T：

c.采用阈值函数对信号g进行去噪，阈值函数为：

其中，sgn(·)为符号函数，ln(·)为对数函数，b为阈值函数可调整参数，取0.1，g(a)为信号g的第a个样本点的幅值，| |为取绝对值符号，a＝1,2，…，n；

d.将得到的

至

作为n个样本点的幅值，采用这n个样本点构成的信号即为第二个去噪信号

Images