«上一篇
文章快速检索     高级检索
下一篇»
  哈尔滨工程大学学报  2020, Vol. 41 Issue (9): 1329-1339  DOI: 10.11990/jheu.202003005
0

引用本文  

张淼, 魏国. 心电信号平滑分解阈值去噪方法[J]. 哈尔滨工程大学学报, 2020, 41(9): 1329-1339. DOI: 10.11990/jheu.202003005.
ZHANG Miao, WEI Guo. Smooth decomposition threshold denoising method for ECG[J]. Journal of Harbin Engineering University, 2020, 41(9): 1329-1339. DOI: 10.11990/jheu.202003005.

基金项目

国家自然科学基金项目(51577041)

通信作者

魏国, E-mail:hitweiguo@hit.edu.cn

作者简介

张淼, 男, 博士研究生;
魏国, 男, 教授, 博士生导师

文章历史

收稿日期:2020-03-01
网络出版日期:2020-09-30
心电信号平滑分解阈值去噪方法
张淼 , 魏国     
哈尔滨工业大学 电气工程及自动化学院, 黑龙江 哈尔滨 150001
摘要:小波分解去噪效果取决于小波函数的选择,经验模态分解存在模态混叠,总体经验模态分解需要大量的重复分解计算。为了解决以上问题,本文基于平滑滤波原理提出了平滑分解阈值去噪方法,用于心电信号去噪。首先将含噪心电信号分解为从高频到低频的多个平滑分解分量,建立边界阈值将其分成3组:含高频噪声分量组、有用信号分量组和低频噪声分量组。提出时变阈值和区间阈值,对含高频噪声分量组去噪。由去噪的含高频噪声分量组和有用信号分量组重构心电信号。通过对仿真和MIT-BIH真实心电信号的去噪实验表明:实验范围内的信噪比提升19.656~2.448 dB;相关系数0.885~0.999;能量比接近于1。去噪心电信号的P波、T波和QRS波完整清晰,波形位置、形态和幅值与源信号一致,具有实际临床应用价值。
关键词平滑滤波    平滑分解    经验模态分解    小波分解    时变阈值    区间阈值    心电去噪    阈值去噪    
Smooth decomposition threshold denoising method for ECG
ZHANG Miao , WEI Guo     
School of Electrical Engineering and Automation, Harbin Institute of Technology, Harbin 150001, China
Abstract: The denoising effect of wavelet decomposition depends on the choice of wavelet functions. Empirical mode decomposition has the problem of mode aliasing, and the ensemble empirical mode decomposition requires a lot of repeated decomposition process. In order to solve the above problems, a smoothing decomposition threshold denoising method based on the principle of smoothing filters is proposed for denoising ECG signals. Firstly, the noisy ECG signal is decomposed into several smoothing decomposition components from high frequency to low frequency, and the boundary threshold is established to divide the components into three groups: the high-frequency noise components group, the useful signal components group and the low-frequency noise components group. The time-varying threshold and the interval threshold are proposed to eliminate the high-frequency noise in the first group. The ECG signal is reconstructed by adding the denoised high-frequency noise components and the useful signal components. The denoising experiment results for the simulated signals and the real signals from MIT-BIH database show that the signal-to-noise ratio (SNR) within the experimental range is improved by 19.656 dB~2.448 dB. The correlation coefficient is between 0.885 and 0.999 and the energy ratio is close to 1. P waves, T waves and QRS complexes of the denoised ECGs are complete and clear. The position, shape and amplitude of the waveforms are consistent with the source signals, indicating that the proposed method has practical clinical application value.
Keywords: smooth filtering    smooth decomposition    empirical mode decomposition    wavelet decomposition    time-varying threshold    interval threshold    ECG signal denoising    threshold denoising    

心电图(electrocardiogram, ECG)是人类心脏生理电信号的记录,直接反映心脏的健康状况。对心律失常、室性早搏、心肌缺血和心肌梗死等心脏疾病的诊断具有重要意义。然而,实测心电图总是混有各种噪声:肌电噪声、电极接触噪声、运动伪迹、基线漂移以及电力线干扰等。这些噪声会影响心电图的真实形态,甚至淹没那些重要的心电图细节,导致医生误诊。因此,寻找一个好的心电图降噪方法一直是科学研究者们的奋斗目标。

现有的心电去噪方法主要包括基于傅里叶分析的经典数字滤波、自适应滤波[1-3]、神经网络及支持向量机学习算法[4]、小波去噪算法[5-6]、经验模态分解算法[7-11]等。由于心电信号具有非线性非平稳特性,傅里叶分析不能很好地解决频谱随时间变化的问题。神经网络学习方法则需要大量的实例信号用于训练,而处于不同时间和空间的受试者,其心电信号有很大差异,训练结果的普适性较差。小波分析和经验模态分解为非线性非平稳信号的分析和心电信号去噪提供了有力工具,并成为了目前心电去噪领域的热点。然而,小波函数的选择对小波分析去噪效果有很大影响,目前还没有针对心电信号的小波函数[5]。相对而言,经验模态分解(empirical mode decomposition, EMD)的自适应性,更有利于心电信号的去噪。但是,EMD在分解心电信号时存在模态混叠现象[7, 12],影响信号分析的成功率。总体经验模态分解(ensemble empirical mode decomposition, EEMD)虽然解决了EMD的模态混叠问题,但是总体数量的增加使得重复分解的计算量成倍增长[13]。而辅助噪声的引入,又产生了噪声残留和信号重构误差的问题。

本文提出平滑分解(smoothing decomposition, SD)方法,用于心电信号去噪,意在避免模态混叠和基本函数的选择,简化分解过程以提高计算效率。SD基于平滑滤波原理,通过改变平滑点数,将含有噪声的心电信号分解为从高频到低频的平滑分解分量(smoothing decomposition component, sdc)。进而建立分界阈值,将sdc分成3组:含高频噪声sdc组、有用信号sdc组和低频噪声sdc组。提出时变阈值(time-varying threshold, TVth)和区间阈值(interval threshold, Ith)概念,对含高频噪声sdc组进行阈值去噪。最后将低频噪声sdc组中的信号作为基线漂移和运动伪迹剔除,用去噪后的含高频噪声sdc组和有用信号sdc组进行信号重构获得去噪心电信号。通过对仿真心电信号和MIT-BIH心律不齐数据库的真实心电信号进行去噪实验,验证了本文方法的有效性,去噪效果良好。

1 平滑分解阈值去噪方法 1.1 平滑分解

根据奈奎斯特定理,采样频率必须大于等于有用信号最高频率的2倍。假设对心电信号的采样满足奈奎斯特采样定理,则实测信号采样频率的1/2为有用心电信号的最高频率。通过三点平滑滤波,可以将频率大于1/2采样频率的信号滤除,将滤除的信号定义为第1阶平滑分解分量(sdc1)。而剩余信号中的最高频率则为前一阶提取频率的下限。继续按剩余信号最高频率的1/2进行新的平滑滤波提取,平滑滤波点数为前一级平滑滤波点数的2倍减1。直到平滑点数大于1/2数据长度且小于数据长度时分解结束。平滑分解方法描述如下。

设实际检测到的心电信号为s(n),n=1,2,…,N为样本序号,N为信号长度。记第k阶平滑分解分量为sdck(n),剩余信号记为rk(n),k=1,2,…,MM=int(lb2(N-1))为最大分解阶数。

1) 当k=1时

$ {{\rm{sd}}{{\rm{c}}_1}(n) = s(n) - \frac{1}{4}\left( {s(n) + \sum\limits_{i = n - 1}^{n + 1} s (i)} \right)} $ (1)
$ {{r_1}(n) = s(n) - {\rm{sd}}{{\rm{c}}_1}(n)} $ (2)

2) 当k≥2时

$ \begin{array}{*{20}{l}} { {\rm{sd}}{{\rm{c}}_k}(n) = {r_{k - 1}}(n) - \frac{1}{{2({2^{k - 1}} + 1)}} \cdot }\\ {\quad \left( {{r_{k - 1}}(n) + \sum\limits_{i = n - {2^{k - 1}}}^{n + {2^{k - 1}}} {{r_{k - 1}}} (i)} \right)} \end{array} $ (3)
$ {r_k}(n) = {r_{k - 1}}(n) - {\rm{sd}}{{\rm{c}}_k}(n) $ (4)

在以上各式中,当i < 1或i>N时,s(i)和rk-1(i)取相应的数据端点值。

3) 信号重构

$ s(n) = \sum\limits_{k = 1}^M { {\rm{sd}}{{\rm{c}}_k}} (n) + {r_M}(n) $ (5)
1.2 sdc分组

借鉴小波分解和经验模态分解去噪方法的经验,需将sdc分成3组:含高频噪声sdc组、有用信号sdc组和低频噪声sdc组。对于含高频噪声sdc组,主要包含随机噪声、肌电噪声和有用信号的高频成分,采取阈值去噪法进行去噪[9-11, 14]。对于有用信号sdc组,噪声很小可以忽略。而低频噪声sdc组的主要成分为基线漂移和运动伪迹,在心电信号重构时被剔除。

1.2.1 含高频噪声sdc组分界

文献[15]指出,高斯白噪声的固有模态函数(intrinsic mode function, imf)仍然服从高斯分布。imf的能量随阶数k的增加近似按1/2k倍的规律减少,而sdc的特性与imf相似,这一点在后面的实验结果分析中可以得到证实。通常当信号的均值趋于零时,方差与平均能量具有相同的变化规律。因此,用各阶sdc的总体方差与相应阶次中包含的噪声估计方差进行依次比较,当两者的比值达到极大值时,表明sdc中噪声能量占比最小。将噪声能量占比最小的sdc阶数定义为含高频噪声sdc组的分界阈值(high-frequency boundary threshold, HBth),阶数小于等于HBth的所有sdc为含高频噪声sdc组。实际应用中,当sdc的阶数较高时,对噪声方差的估计会有较大的偏差。因此,必须对参与噪声能量占比计算的sdc做一定的限制。具体的HBth计算方法如下:

σsdc为sdc的总体方差,$ {\hat \sigma _{{\rm{nsdc}}}} $为sdc中的噪声方差估计,具体为:

$ {\sigma _{{\rm{sdc}}}}(k) = \sqrt {\frac{1}{{N - 1}}\sum\limits_{i = 1}^N {({\rm{sd}}{{\rm{c}}_k}(} i) - \overline {{\rm{sd}}{{\rm{c}}_k}} {)^2}} $ (6)
$ {\hat \sigma _{{\rm{ nsdc }}}}(k) = \frac{{ {\rm{median}} | {\rm{sd}}{{\rm{c}}_k}(i) - \overline { {\rm{sd}}{{\rm{c}}_k}} |}}{{0.674{\kern 1pt} {\kern 1pt} {\kern 1pt} 5}} $ (7)

式中:$ \overline {{\rm{sd}}{{\rm{c}}_k}} = \frac{1}{N}\sum\limits_{i = 1}^N {{\rm{sd}}{{\rm{c}}_k}\left( i \right)} $, 为sdck(n)的均值。

计算整体方差与噪声方差的比值σratio(k)和差值σdiff(k),并考虑到$ {\hat \sigma _{{\rm{nsdc}}}}\left( k \right) \le {\sigma _{{\rm{sdc}}}}\left( k \right) $,则:

$ {\sigma _{{\rm{ ratio }}}}(k) = \left\{ {\begin{array}{*{20}{l}} {{\sigma _{{\rm{sdc}}}}(k)/{{\hat \sigma }_{{\rm{nsdc}}}}(k),}&{{\sigma _{{\rm{sdc}}}}(k) > {{\hat \sigma }_{{\rm{nsdc}}}}(k)}\\ {1,}&{{\sigma _{{\rm{sdc}}}}(k) \le {{\hat \sigma }_{{\rm{nsdc}}}}(k)} \end{array}} \right. $ (8)
$ {\sigma _{{\rm{ diff }}}}(k) = \left\{ {\begin{array}{*{20}{l}} {{\sigma _{{\rm{ sdc }}}}(k) - {{\hat \sigma }_{{\rm{ nsdc }}}}(k),}&{{\sigma _{{\rm{ sdc }}}}(k) > {V_{{\rm{ nsdc }}}}(k)}\\ {0,}&{{\sigma _{{\rm{ sdc }}}}(k) \le {{\hat \sigma }_{{\rm{ nsdc }}}}(k)} \end{array}} \right. $ (9)

首先计算σdiff(k)的第1个极小值,对应的sdc阶数记为Kdiff,然后在[2, Kdiff]区间内计算σratio(k)的第1个极大值,对应的sdc阶数记为Kratio。含高频噪声sdc组的分界阈值HBth由下式确定:

$ {\rm{HBth}} = \left\{ {\begin{array}{*{20}{l}} {{K_{{\rm{ratio}}}},}&{2 \le {K_{{\rm{ratio}}}} \le {K_{{\rm{diff}}}}}\\ {1,}&{{\rm{其他}}} \end{array}} \right. $ (10)
1.2.2 低频噪声sdc组分界

根据实测心电信号的R峰个数和位置,计算出平均心率周期。将各阶sdc的平均周期与心电平均周期比较,小于等于心电平均周期的sdc中,最大的sdc阶数定义为低频噪声sdc组的分界阈值(low-frequency boundary threshold, LBth)。阶数大于LBth的所有sdc确定为低频噪声sdc组。

设心电信号的平均周期为PECG,sdc的平均周期为Psdc,分别由(11)和(12)计算。

$ {{P_{{\rm{ECG}}}} = \frac{{{t_R}({N_R}) - {t_R}(1)}}{{{N_R} - 1}}} $ (11)
$ {{P_{{\rm{sdc}}}}(k) = \frac{{t {\rm{Peak}}{ _k}({N_k}) - t {\rm{Peak}}{ _k}(1)}}{{{N_k} - 1}}} $ (12)

式中:NRR峰个数;tR(1)和tR(NR)分别为第1个和最后一个R峰对应的时间;Nk为sdck的极大值点或极小值点个数;tPeakk(1)和tPeakk(Nk)分别为sdck的第1个极值点和最后一个极值点对应的时间。低频噪声sdc组的分界阈值LBth为:

$ {\rm{LBth}} = \mathop {{\rm{argmax}} [{P_{{\rm{sdc}}}}(k)]}\limits_{{P_{{\rm{sdc}}}}(k) \le {P_{ECG}}} $ (13)
1.3 阈值去噪方法 1.3.1 时变阈值去噪

基于通用阈值去噪方法的思想[16-17],提出时变阈值(time-varying threshold,TVth)的概念,用于适应肌电噪声统计特性的时变性。由于通用阈值去噪会使噪声大的地方产生噪声残留,噪声小的地方会造成心电信号的损失。因此,需要建立自适应变化的时变阈值进行去噪。时变阈值以心电平均周期为时间窗口,在窗口内按通用阈值法计算阈值,当时间窗连续通过信号长度时,得到时变阈值函数曲线,具体方法如下:

时间窗内信号的噪声方差估计为:

$ {\hat \sigma _k}(n) = \frac{{ {\rm{median}} (| {\rm{sd}}{{\rm{c}}_k}(i) - \overline { {\rm{sd}}{{\rm{c}}_k}} (n)|)}}{{0.674{\kern 1pt} {\kern 1pt} {\kern 1pt} 5}} $ (14)
$ {n - \frac{{{P_{{\rm{ECG}}}}}}{2} \le i \le n + \frac{{{P_{{\rm{ECG}}}}}}{2}} $
$ {\overline {{\rm{sd}}{{\rm{c}}_k}} (n) = \frac{1}{{{P_{{\rm{ECG}}}}}}\sum\limits_{i = n - \frac{{{P_{{\rm{ECG}}}}}}{2}}^{n + \frac{{{P_{{\rm{ECG}}}}}}{2}} {\rm{s}} {\rm{d}}{{\rm{c}}_k}(i)} $

式中PECG由(11)式确定。

计算时变阈值曲线:

$ {\rm{TVth}}{ _k}(n) = {\hat \sigma _k}(n)\sqrt {C{\rm{ln}}({P_{{\rm{ECG}}}})} $ (15)

用时变阈值进行去噪:

$ {\rm{sdc}}{ _{TVd,k}}(n) = \left\{ {\begin{array}{*{20}{l}} { {\rm{sdc}}{ _k}(n),}&{| {\rm{sdc}}{ _k}(n)| \ge {\rm{TVth}}{ _k}(n)}\\ {0,}&{| {\rm{sdc}}{ _k}(n)| < {\rm{TVth}}{ _k}(n)} \end{array}} \right. $ (16)

式中:下标k代表第k阶分解分量;系数C由经验确定,本文取为1.5。

图 1为通用阈值和时变阈值的比较。其中水平虚线为通用阈值线,加粗实线为时变阈值线,信号为含有肌电噪声的心电信号。可以看出,时变阈值随肌电噪声的大小而自适应变化。

Download:
图 1 通用阈值和时变阈值的比较 Fig. 1 Comparison of universal and time-varying threshold
1.3.2 区间阈值去噪

心电信号经过时变阈值去噪后,往往还会残留有随机的孤立噪声,这些噪声对信噪比和相关系数影响不大,但是会对局部波形产生干扰。本文提出一种区间阈值(interval threshold, Ith)法,用于消除这些残留噪声。定义由分界sdcHBth去噪后的波形宽度为阈值区间,阈值区间内数据保留,区间外数据剔除。将利用时变阈值去噪后的sdcHBth记为sdcTVd, HBth,并对其进行平滑处理以消除波形不连续点,恢复实际波形宽度。平滑信号记为$ {\widetilde {{\rm{sdc}}}_{{\rm{TVd, HBth}}}} $

$ \widetilde {{\rm{sdc}}}{ _{{\rm{TVd,HBth}}}}(n) = \frac{1}{{2C + 1}}\sum\limits_{i = n - C}^{n + C} {{\rm{ sdc}}{{\rm{ }}_{{\rm{TVd,HBth}}}}} (i) $ (17)

式中参数2C+1为平滑点数,由经验选取,本文取C=3。

区间阈值计算如下:

$ {\rm{Ith}}(n)\left\{ {\begin{array}{*{20}{l}} 1&{{{\widetilde {{\rm{sdc}}}}_{{\rm{TVd,Hth}}}}(n) \ne 0}\\ 0&{{{\widetilde {{\rm{sdc}}}}_{{\rm{TVd,Hth}}}}(n) = 0} \end{array}} \right. $ (18)

区间阈值法剔除残留噪声:

$ \begin{array}{*{20}{c}} { {\rm{sdc}}{ _{Id,k}}(n) = {\rm{sdc}}{ _{{\rm{TVd}},k}}(n) \cdot {\rm{Ith}} (n)}\\ {k = 1,2, \cdots ,{\rm{HBth}} - 1} \end{array} $ (19)
1.4 去噪心电信号的重构

由时变阈值和区间阈值去噪后的含高频噪声sdc组和有用信号sdc组相叠加,重构去噪心电信号:

$ {s_d}(n) = \sum\limits_{k = 1}^{{\rm{HBth}}} { {\rm{sdc}}{ _{Id,k}}} (n) + \sum\limits_{k = {\rm{HBth}} + 1}^{{\rm{LBth}}} { {\rm{sdc}}{ _k}} (n) $ (20)
1.5 SD阈值去噪方法的计算流程

本文提出的SD阈值去噪方法的去噪流程如图 2所示。

Download:
图 2 本文方法去噪流程 Fig. 2 Flow chart of the proposed method for denoising
2 去噪结果评估方法 2.1 信噪比提升

信噪比提升定义为去噪后信号的信噪比与去噪前信号的信噪比之差,表明去噪方法的去噪能力和去噪效果。设SNRimp、SNRd、SNRn分别为信噪比提升、去噪后信噪比和去噪前信噪比:

$ {\rm{SNR}}{ _n} = 10{\rm{lg}}\frac{{\sum\limits_{n = 1}^N {{s^2}} (n)}}{{\sum\limits_{n = 1}^N {({s_n}(} n) - s(n){)^2}}} $ (21)
$ {\rm{SNR}}{ _d} = 10{\rm{lg}}\frac{{\sum\limits_{n = 1}^N {{s^2}} (n)}}{{\sum\limits_{n = 1}^N {({s_d}(} n) - s(n){)^2}}} $ (22)
$ {\rm{SN}}{{\rm{R}}_{{\rm{imp}}}} = {\rm{SN}}{{\rm{R}}_d} - {\rm{SN}}{{\rm{R}}_n} $ (23)

式中:s(n)为不含噪声的源信号;sn(n)为含噪声的测量信号; sd(n)为去噪后的信号。

2.2 相似度指标

CR定义为去噪后信号与无噪声源信号之间的相关系数,用于表达去噪后信号与源信号的相似程度,从波形形态方面评价去噪效果:

$ {C_R} = \frac{{\sum\limits_{n = 1}^N {({s_d}(} n) - \overline {{s_d}(n)} )(s(n) - \overline {s(n)} )}}{{\sqrt {\sum\limits_{n = 1}^N {({s_d}(} n) - \overline {{s_d}(n)} {)^2}} \sqrt {\sum\limits_{n = 1}^N {(s(} n) - \overline {s(n)} {)^2}} }} $ (24)
2.3 幅值度指标

相关系数CR的大小只能反映去噪信号和源信号之间在波形上的相似度。而心电信号的P波、T波、QRS波的大小和位置(间期)则是心脏病理判断的重要指标。为此,用去噪信号与源信号的平均能量之比,作为去噪后心电信号幅值变化的评价指标,用ER表示。当ER等于1时,表明去噪信号的幅值与源信号一致;ER大于1时,说明去噪信号中还有残留噪声;ER小于1时,说明去噪后信号出现了有用信息的损失。能量比ER的计算式为:

$ {\rm{ER}} = \frac{{\frac{1}{N}\sum\limits_{n = 1}^N {s_d^2} (n)}}{{\frac{1}{N}\sum\limits_{n = 1}^N {{s^2}} (n)}} $ (25)
3 实验验证与实验结果分析 3.1 实验数据来源

1) 仿真心电信号。6个成人仿真心电信号Syn00~Syn05,由文献[18]和[19]提供的心电信号工具包(OSET)中的FECGSYN仿真生成。采样频率为360 Hz,信号长度为10 s。如图 3所示。需要指出的是,本文中为了后续计算信噪比等参数,对数据均进行了统一化处理,因此纵坐标没有实际的物理含义。通过图展示的是信号形状,因而不关心振幅的单位。本文后续的图均为此情况。

Download:
图 3 仿真心电信号 Fig. 3 Simulated ECG signals

2) 真实心电信号。真实心电信号来源于MIT-BIH心律不齐数据库(MITDB)[20]。采样频率360 Hz,截取的时间长度为10 s。表 1给出本文使用的数据记录编号、通道名和开始时间位置。

表 1 从MIT-BIH心律不齐数据库选取的真实心电信号 Table 1 Real ECG data chosen from MIT-BIH Arrhythmia Database

3) 实验中的噪声。本文实验涉及了4种噪声,分别为基线漂移(baseline wandering, BW)、肌电噪声(electromyogram, EMG)、白噪声(white noise, WN)和混合噪声(hybrid noise, HN)。其中BW和EMG来源于MIT-BIH的NSTDB数据库[21],采样频率360 Hz,截取长度10 s。白噪声WN由Matlab 2016的随机函数生成。而混合噪声HN则由以上3种噪声混合生成。其中BW经过高频噪声滤除,而EMG则经过基线漂移修正,如图 4所示。

Download:
图 4 本文实验中的4种噪声 Fig. 4 Four categories of noises used in experiments.
3.2 实验结果及分析 3.2.1 平滑分解方法的验证

随机采集了30个单位方差的白噪声信号,采样频率360 Hz,采样时间10 s。分别用本文方法进行平滑分解。针对每个分解分量sdc计算峰值点个数、平均频率以及sdc能量与原始信号能量的百分比。平均结果列于表 2图 5给出了某个白噪声信号的分解结果。结果显示sdc的峰值点数、平均频率和能量百分比都与文献[18]中的EMD分解结果一致。因此平滑分解(SD)方法具有与EMD相似的性能。

表 2 30个随机白噪声平滑分解结果的平均值 Table 2 Mean values of smoothing decomposition results from 30 random white noise
Download:
图 5 白噪声的平滑分解分量(前5 s,1 800采样点) Fig. 5 Smooth decomposition component of white noise (first 5 s, 1 800 sampling points)
3.2.2 时变阈值和区间阈值去噪效果

仿真信号Syn00叠加EMG噪声,改变信噪比从-5~30 dB。用本文的SD方法进行分解后,分别用通用阈值法、时变阈值法、区间阈值法和时变+区间阈值法进行去噪实验。结果列于表 3~5。表中数据显示,时变阈值法去噪效果好于通用阈值法,时变阈值+区间阈值法最好。

表 3 Syn00+EMG的去噪结果(信噪比提升) Table 3 Denoising results of Syn00+EMG(SNRimp)
表 4 Syn00+EMG的去噪结果(相关系数) Table 4 Denoising results of Syn00+EMG(CR)
表 5 Syn00+EMG的去噪结果(能量比) Table 5 Denoising results of Syn00+EMG(ER)

图 6为Syn00+EMG信号在信噪比5 dB时,应用不同阈值法得到的去噪信号波形。结果显示,通用阈值法和时变阈值法,在EMG变化较大的地方,存在噪声残留。区间阈值法则在Q和S波附近有噪声残留。这是因为区间阈值不能去除波群范围内的噪声所致。而时变+区间阈值法综合了两者优势,能够有效弥补其不足,因此去噪效果最好,既没有残留噪声,又能很好地保留信号波形的原有特征,波形失真度最小,适合于肌电噪声的消除。

Download:
图 6 4种阈值去噪方法对Syn00+EMG的去噪结果(5 dB) Fig. 6 Denoising results of Syn00+EMG by four threshold denoising methods (5 dB)
3.2.3 与小波去噪和EMD去噪结果的比较

选择仿真信号Syn00和Syn03分别叠加5 dB的肌电噪声和白噪声进行小波阈值法和EMD阈值法实验。实验结果与本文方法进行比较,验证本文方法的有效性和效果。小波阈值去噪法选择了OriginPro 8 SR4的wtdenoise方法,以避免编程对去噪效果的影响。阈值选择为sqtwolog固定阈值。针对不同的DB小波进行实验,实验结果的局部波形如图 7所示,从中不难看出不同DB小波去噪结果的不同。

Download:
图 7 不同小波函数对小波去噪结果的影响 Fig. 7 Wavelet denoising resuls by different wavelet functions

EMD阈值去噪分2种情况:原始EMD通用阈值去噪和EEMD通用阈值去噪。表 6为EMD、EEMD和本文的SD方法的分解时间和循环次数比较。表中信号采样频率360 Hz,信号长度3 600点。EEMD的总体数量取为50,循环次数可由EMD循环次数乘以50估算。表 6结果显示,本文SD方法的分解时间远远小于EMD和EEMD。

表 6 含噪心电信号分解时间比较 Table 6 Comparison of decomposition time for noisy ECG

图 8为DB4小波阈值去噪、EMD通用阈值去噪、EEMD通用阈值去噪和本文SD通用阈值去噪及SD时变+区间阈值去噪的比较。图中左侧为Syn00叠加肌电噪声的结果,右侧为Syn00叠加白噪声的结果。可以看出,本文方法去噪结果好于其他分解方法,特别是时变+区间阈值方法效果更好。从去噪波形上看,EMD、EEMD和WT更适合于白噪声的去噪,而本文方法对肌电噪声去噪也适合。

Download:
图 8 本文方法SD与WT、EMD、EEMD去噪波形比较(5 dB) Fig. 8 Comparison of denoising waveforms among the proposed SD, WT, EMD and EEMD (5 dB)
3.2.4 仿真心电信号实验结果及分析

通过以上实验,验证了本文方法具有很好的去噪效果,特别是对肌电噪声的滤除具有明显优势。本节将通过仿真心电信号叠加不同信噪比的各种噪声进行去噪实验,进一步验证本文方法的有效性、适应性和去噪效果。

表 7~10分别给出了基线漂移(BW)、肌电噪声(EMG)、白噪声(WN)和混合噪声(HN)与仿真信号Syn00~Syn05叠加后,用本文方法进行去噪,各信号去噪结果的平均值。表中数据表明:在信噪比从0 dB~20 dB范围内,基线漂移滤除效果最好,相关系数大于0.99,能量比在1~1.02,表明没有信号损失,噪声残留较小。肌电噪声去噪效果相对最差,但是信噪比提升大于2.44 dB,相关系数大于0.88,能量比大于0.94小于1,这说明去噪时有一定的信号损失,但损失较小,不影响结果的合理性,仍具有很好的去噪效果。白噪声去噪结果的信噪比提升和相关系数均高于肌电噪声的情况。但是,能量比小于1,并且小于肌电噪声去噪的能量比。这说明白噪声去噪时,产生了较大的幅值损失,但是相关系数大于0.954,具有优良的波形相似性。混合噪声的去噪效果介于肌电噪声和白噪声之间,这是因为多个随机变量的叠加,统计特性趋于高斯分布的原因。

表 7 仿真心电信号基线漂移(BW)滤除结果 Table 7 Results of eliminating BW for simulated ECG
表 8 仿真心电信号肌电噪声(EMG)去噪结果 Table 8 Results of eliminating EMG for simulated ECG
表 9 仿真心电信号白噪声(WN)去噪结果 Table 9 Results of eliminating WN for simulated ECG
表 10 仿真心电信号混合噪声(HN)去噪结果 Table 10 Results of eliminating HN for simulated ECG

图 9给出仿真信号Syn00~Syn05叠加5 dB肌电噪声后,去噪信号与未去噪信号的波形比较。图 10为Syn04分别与BW、HN、EMG和WN 4种噪声叠加的去噪波形。由图 9图 10可以看出,去噪后波形与无噪信号(Syn00~Syn05)波形一致,P波、T波和QRS波群完整清晰。

Download:
图 9 仿真心电信号的肌电噪声去噪结果(5 dB) Fig. 9 Waveforms after eliminating EMG for simulated ECG (5 dB)
Download:
图 10 Syn04分别叠加4种噪声时的去噪结果(5 dB) Fig. 10 Waveforms after eliminating four categories of noise for Syn04 (5 dB)
3.2.5 真实心电信号实验结果及分析

真实心电信号来源于MIT-BIH心律不齐数据库的临床记录,选择其中的10个记录,如表 1所示。本文将其命名为MB100、MB101、MB103、MB112、MB113、MB115、MB117、MB119、MB122、MB123。由于这些临床数据不可避免地含有一定的噪声,首先对其进行高频噪声滤除和基线漂移修正,目的是为了在分析去噪效果时能够定量分析。经过去噪和修正的真实心电信号记为DMB100、DMB101、DMB103、DMB112、DMB113、DMB115、DMB117、DMB119、DMB122、DMB123。以DMB作为干净信号,估计原信号MB的信噪比,结果列于表 11图 11选择给出了3个信号的噪声滤除前后的波形和滤除的噪声曲线。由表 11数据和图 11的波形可以证明,DMB能够反映真实的心电特性。

表 11 MIT-BIH真实心电信号的信噪比估计 Table 11 SNR estimation for MIT-BIH real ECG
Download:
图 11 MB100、MB115和MB122的噪声滤除结果 Fig. 11 Noise filtering results of MB100, MB115 and MB122

以DMB作为无噪真实心电信号,叠加BW、EMG、WN和HN建立含噪信号。然后应用本文方法进行平滑分解和时变+区间阈值去噪。去噪结果按10个DMB信号进行平均,结果给出在表 12~15中。在信噪比0 dB~20 dB范围内,真实信号与仿真信号去噪结果一致。真实心电信号去噪结果中,基线漂移去噪效果最好,信噪比提升大于10 dB,相关系数大于0.99,表明去噪信号与源信号的相似度达到99%,能量比大于1表明去噪信号的能量大于源信号,去噪信号未丢失有用信息。能量比小于1.017表明去噪信号中残留的噪声成分不大于1.7%。肌电噪声去噪效果相对差一些,信噪比提升相对较少,相关系数大于0.89,能量比大于0.93。白噪声去噪结果在信噪比提升和相关系数上优于肌电噪声。但是,能量比大于0.87小于1,比肌电噪声的能量比低,因此有更大的幅值改变。混合噪声的去噪结果和白噪声相当。

表 12 真实心电信号基线漂移(BW)去噪结果 Table 12 Results of eliminating BW for real ECG
表 13 真实心电信号肌电噪声(EMG)去噪结果 Table 13 Results of eliminating EMG for real ECG
表 14 真实心电信号白噪声(WN)去噪结果 Table 14 Results of eliminating WN for real ECG
表 15 真实心电信号混合噪声(HN)去噪结果 Table 15 Results of eliminating HN for real ECG

图 12为真实心电信号叠加5 dB肌电噪声时的去噪信号波形。图 13选择了真实心电信号DMB119与BW、HN、EMG和WN 4种噪声叠加时,本文方法的去噪波形。可以看出,真实心电信号去噪后的P波、T波和QRS波群还原良好。证明本文方法不仅适合于肌电噪声的去噪,而且对其他噪声的去噪效果更好。

Download:
图 12 真实心电信号的肌电噪声去噪结果(5 dB) Fig. 12 Waveforms after eliminating EMG for real ECG (5 dB)
Download:
图 13 DMB119分别叠加4种噪声时的去噪结果(5 dB) Fig. 13 Waveforms after eliminating four categories of noise for DMB119 (5 dB)
4 结论

1) 通过噪声分解分析,证明了本文提出的平滑分解方法与EMD分解具有相似的滤波组特性,并具有更高的分解效率;通过对仿真心电信号和MIT-BIH真实心电信号的去噪实验表明,本文方法的sdc分组、时变阈值和区间阈值,解决了肌电噪声去噪的难题,去噪效果明显好于目前的EMD阈值去噪、EEMD阈值去噪和小波阈值去噪效果;在实验范围内,本文的平滑分解阈值去噪方法具有较高的信噪比提升,相关系数总体上大于0.9以上,去噪后的心电波形失真小,P波、T波和QRS复合波完整清晰,其位置、形态和幅值与源信号一致,达到实际临床应用的水平。

2) 本文方法基于所提出的时变阈值和区间阈值方法,很好地解决了心电信号中肌电噪声问题。平滑滤波分解大大提高了EMD和EEMD的分解效率。避免了小波分解中小波函数选择的不确定性。

3) 本文方法没有给出心电周期的理论确定方法,尚需观察给出。对于包含在有用信号sdc中的中低频接触噪声这一目前各种去噪方法均未涉及的难题,没有进行研究,有望于进一步研究解决。

参考文献
[1]
CUOMO S, DE PIETRO G, FARINA R, et al. A novel O(n) numerical scheme for ECG signal denoising[J]. Procedia computer science, 2015, 51: 775-784. DOI:10.1016/j.procs.2015.05.198 (0)
[2]
KRISHNAN S R, SEELAMANTULA C S. On the selection of optimum savitzky-golay filters[J]. IEEE transactions on signal processing, 2013, 61(2): 380-391. (0)
[3]
VULLINGS R, DE VRIES B, BERGMANS J W M. An adaptive kalman filter for ECG signal enhancement[J]. IEEE transactions on biomedical engineering, 2011, 58(4): 1094-1103. DOI:10.1109/TBME.2010.2099229 (0)
[4]
MOEIN S. An MLP neural network for ECG noise removal based on Kalman filter[M]//ARABNIA H R. Advances in Computational Biology. New York: Springer, 2010: 109-116. (0)
[5]
朱荣亮, 陶晋宜. 基于改进小波阈值去噪算法的心电信号处理及仿真[J]. 数学的实践与认识, 2019, 49(5): 143-150.
ZHU Rongliang, TAO Jinyi. ECG signal processing and simulation based on improved wavelet threshold denoising algorithm[J]. Mathematics in practice and theory, 2019, 49(5): 143-150. (0)
[6]
PHINYOMARK A, LIMSAKUL C, PHUKPATTARANONT P. EMG denoising estimation based on adaptive wavelet thresholding for multifunction myoelectric control[C]//Proceedings of 2009 Innovative Technologies in Intelligent Systems and Industrial Applications. Monash, Malaysia, 2009: 171-176. (0)
[7]
HUANG N E, SHEN Zheng, LONG S R, et al. The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proceedings of the royal society a, 1998, 454(1971): 903-995. DOI:10.1098/rspa.1998.0193 (0)
[8]
尹丽, 陈富民, 张琦, 等. 采用集合经验模态分解和改进阈值函数的心电自适应去噪方法[J]. 西安交通大学学报, 2020, 54(1): 101-107.
YIN Li, CHEN Fumin, ZHANG Qi, et al. ECG adaptive denoising method based on EEMD and improved threshold function[J]. Journal of Xi'an JiaoTong University, 2020, 54(1): 101-107. (0)
[9]
KOPSINIS Y, MCLAUGHLIN S. Development of EMD-based denoising methods inspired by wavelet thresholding[J]. IEEE transactions on signal processing, 2009, 57(4): 1351-1362. DOI:10.1109/TSP.2009.2013885 (0)
[10]
NGUYEN P, KIM J M. Adaptive ECG denoising using genetic algorithm-based thresholding and ensemble empirical mode decomposition[J]. Information sciences, 2016, 373: 499-511. DOI:10.1016/j.ins.2016.09.033 (0)
[11]
MOHGUEN W, BEKKA R E. EMD-based denoising by customized thresholding[C]//Proceedings of 2017 International Conference on Control, Automation and Diagnosis. Hammamet, Tunisia, 2017: 19-23. (0)
[12]
TORRES M E, COLOMINAS M A, SCHLOTTHAUER G, et al. A complete ensemble empirical mode decomposition with adaptive noise[C]//Proceedings of 2011 IEEE International Conference on Acoustics, Speech and Signal Processing. Prague, Czech Republic, 2011: 4144-4147. (0)
[13]
WU Zhaohua, HUANG N E. Ensemble empirical mode decomposition:a noise-assisted data analysis method[J]. Advances in adaptive data analysis, 2009, 1(1): 1-41. (0)
[14]
PAL S, MITRA M. Empirical mode decomposition based ECG enhancement and QRS detection[J]. Computers in biology and medicine, 2012, 42(1): 83-92. DOI:10.1016/j.compbiomed.2011.10.012 (0)
[15]
WU Zhaohua, HUANG N E. A study of the characteristics of white noise using the empirical mode decomposition method[J]. Proceedings of the royal society A, 2004, 460(2046): 1597-1611. DOI:10.1098/rspa.2003.1221 (0)
[16]
DONOHO D L, JOHNSTONE I M. Ideal spatial adaptation by wavelet shrinkage[J]. Biometrika, 1994, 81(3): 425-455. DOI:10.1093/biomet/81.3.425 (0)
[17]
AWAL A, MOSTAFA S S, AHMAD M, et al. An adaptive level dependent wavelet thresholding for ECG denoising[J]. Biocybernetics and biomedical engineering, 2014, 34(4): 238-249. DOI:10.1016/j.bbe.2014.03.002 (0)
[18]
MCSHARRY P E, CLIFFORD G D, TARASSENKO L, et al. A dynamical model for generating synthetic electrocardiogram signals[J]. IEEE transactions on biomedical engineering, 2003, 50(3): 289-294. DOI:10.1109/TBME.2003.808805 (0)
[19]
SAMENI R, CLIFFORD G D, JUTTEN C, et al. Multichannel ECG and noise modeling:application to maternal and fetal ECG signals[J]. EURASIP journal on advances in signal processing, 2007, 2007: 043407. DOI:10.1155/2007/43407 (0)
[20]
MOODY G B, MARK R G. The MIT-BIH arrhythmia database on CD-ROM and software for use with it[C]//Proceedings Computers in Cardiology. Chicago, USA, 1990: 185-188. (0)
[21]
MOODY G B, MULDROW W E, MARK R G. A noise stress test for arrhythmia detectors[J]. Computers in cardiology, 1984, 11: 381-384. (0)