郑敏敏 高小榕 谢海鹤
1(莆田学院机电工程学院,福建 莆田 351100)2(清华大学医学院,北京 100084)
心电信号小波去噪的改进算法研究
郑敏敏1*高小榕2谢海鹤1
1(莆田学院机电工程学院,福建 莆田 351100)2(清华大学医学院,北京 100084)
对心电信号(ECG)中的基线漂移、工频干扰和肌电干扰等噪声进行去除,在波形识别、医疗诊断和治疗等领域具有重要意义。提出用sym5小波函数对心电信号进行8层小波分解。根据有用信号强度在每一层平均分配而噪声强度随分解层数增加而减少的规律,将分解得到的每一层的小波细节系数设置不同的阈值,最后用所提出的新阈值函数进行小波阈值去噪。该阈值函数既能克服硬阈值函数在阈值附近不连续的缺点,又可弥补软阈值函数与原函数之间存在固定差值的不足。以MIT-BIH心电数据库中的101号文件作为原始数据,将整个数据文件进行平均分段,每段有1 200个数据点,对每段数据进行加噪仿真分析,结果表明所提出的去噪算法得到的去噪信号信噪比比硬阈值函数和软阈值函数分别提高2.31%和8.04%,从而证明所提算法的有效性。
心电信号;小波函数;阈值;阈值函数;信噪比
2(BiomedicalEngineeringDepartment,TsinghuaUniversity,Beijing100084,China)
心电信号是人类最早研究并应用于医学临床的生物信号之一,它比其他生物电信号更易检测、规律性更强,可用于对各种心律失常、心室心房肥大、心肌梗死、心律失常、心肌缺血等病症检查。但是,心电信号比较微弱,从人体体表获取的正常心电信号的幅值范围在10 μV~4 mV之间,典型值为1 mV,频率范围在0.05~100 Hz以内,而90%的ECG频谱能量集中在0.25~35 Hz之间。在心电信号采集、数模转换和存储过程中,极易受到环境的影响,如工频干扰、运动伪迹、基线漂移、肌电干扰等[1]。
50/60 Hz的工频干扰可看成由50/60 Hz的正弦信号及其谐波组成,幅值通常与ECG峰峰值相当或更强;运动伪迹幅值通常为几十毫伏;基线漂移的频率范围在0.015~0.300 Hz,基线变化幅度为ECG峰峰值的15%;肌电干扰(EGM)主要能量集中在30~300 Hz范围内,频率范围与心电信号的频率范围重叠。
在去噪过程中,由于心电信号具有非平稳特性且污染噪声分布范围大,限制了传统线性滤波器的使用,所以在过去的几年中,小波分析被广泛地应用于心电信号的去噪中。
郭新蕾等、苏丽等和赵志栋等分别针对阈值的选择策略进行了深入的研究,对阈值选择的局限性以及造成消噪效果不理想的问题提出了解决方法,并取得了较好的结果[2-4]。
国外对心电信号的去噪的研究也非常活跃,先后提出了基于经验模态分解、平移转移算法和非线性滤波器等方法进行去噪[5-9]。
在利用小波变换对心电信号去噪的研究过程中,小波变换阈值法被广泛使用。但是,也存在一些不足:硬阈值函数法在阈值处发生突变,导致去噪重构后的信号出现伪吉布斯现象;软阈值函数法处理后的分解系数与原分解系数存在恒定差值,导致重构信号能量发生衰减。如何尽可能地去除噪声并且保持心电信号波形不失真,是目前尚需进一步探讨的问题。确定最佳小波函数、分解层数、阈值和阈值函数的方法,是心电信号小波阈值去噪算法研究的重点[10]。
本研究针对阈值、小波函数和阈值函数的选择,提出了一种新的小波去噪算法,通过对MIT-BIH心电数据库中的心电信号进行加噪仿真,利用本研究提出的小波去噪算法实现去噪,得到的去噪信号信噪比(SNR)和最小均方误差(MSE)均达到了理想的效果。
1.1 材料
在本研究中,使用的是MIT-BIH心电数据库中的原始心电数据,选取了101号心电数据进行算法的验证。101号心电数据是一组正常人的心电数据,记录时长为30 min左右,采样率为360 Hz,2个导联,记录中掺杂的噪声较少,所以可将101号数据近似看成是无噪声的干净信号,本研究将第一导联中的所有数据进行分段,每段包含1 200个采样点。为了模拟工频干扰、基线漂移和肌电干扰,本研究在每段采样点中分别加入60 Hz正弦波、0.5 Hz正弦波和信噪比为10 dB的高斯白噪声。
1.2 方法
读取数据后,加入3种信号分别仿真工频干扰、基线漂移和肌电干扰。将加噪后的信号用适当的小波函数在合适的层数上进行离散小波分解,然后提取出小波分解的细节系数和近似系数,对每一层细节系数用改进的阈值和阈值函数进行去噪,将去噪后的小波系数重构,即得到去噪后的心电信号。心电信号的去噪效果用信噪比(signal noise ratio,SNR)和均方误差(mean squared error,MSE)进行评价。基于心电信号的小波去噪步骤如图1所示。
图1 改进阈值去噪算法流程Fig.1 Flow chart of improved threshold denoising algorithm
1.2.1 小波函数的选择
在数值分析中,内积的物理含义是两个图形的相似性:若两个图形完全正交,则内积为0;若两个图形完全一样,则系数为1(相对值)。小波系数的本质是小波基函数与原信号相似度度量。小波变换后的系数比较大,就表明了小波和信号的波形相似程度较大;反之,则比较小。
为了使有用信号分解对应的小波系数与噪声分解对应的小波系数差异尽量明显,应该选择与心电信号的波形最相似的小波函数。以去噪后重构信号的信噪比和均方误差作为评价去噪性能的指标,信噪比越大、均方误差越小,说明去噪效果越好,去噪后的信号和原信号越相似。经过大量对比分析,本研究选择sym5小波作为小波分解函数。
信噪比和均方误差分别定义如下:
(1)
(2)
式中,x(n)为原始心电信号,y(n)为小波去噪后的重构信号,N为信号的长度。
1.2.2 分解尺度的选择
在小波分解中,分解尺度i的选择也是非常重要的一步。i取值越大,则噪声和信号表现的不同特性越明显,越有利于二者的分离。但另一方面,分解尺度越大,重构到的信号失真也会越大,在一定程度上又会影响最终去噪的效果。因此,在应用时要格外注意处理好两者之间的矛盾,选择一个合适的分解尺度。
由于离散小波变换(discrete wavelet transform, DWT)都使用二进小波变换,所以通常每一高阶分解的小波系数个数都约是每一低阶的一半,重构时再插值到原始信号的采样点数。所以,可以分析得出,各阶重构时空长度不变,小波系数的个数随阶次增大而以2的幂次减少,每一高阶重构信号的最大频率应是其相邻低阶重构信号最大频率的一半。
本研究中使用的原始心电信号的采样率为360 Hz,根据奈奎斯特采样定理可知,信号的最高频率可为180 Hz。如果进行8层小波分解,那么对应的每层频带范围如图2所示,其中cAi表示近似系数,cDi表示细节系数(i=1,2,…,8)。
图2 小波分解层数与频率关系Fig.2 Relationship between wavelet decomposition level and frequency
从图2可看出,第8层分解得到的cA8近似系数的频带范围正好覆盖了基线漂移的噪声频带范围,所以本研究中选择的分级尺度为8层。
1.2.3 阈值的选择
进行小波去噪时,阈值选择的好坏直接影响到去噪的效果。如果阈值设置过小,噪声不能被完全去除,达不到去噪的目的;如果阈值设置过大,信号中的有用信息也会被当成噪声去除,使去噪后的信号失真。目前,常见的阈值选取规则主要有无偏似然估计、固定阈值估计,启发式阈值估计和极值阈值估计[11]。一般来讲,极值阈值估计和无偏似然估计方法比较保守,当噪声在信号的高频段分布较少时,这两种阈值估计方法去噪效果较好,可以将微弱的信号提取出来。而固定阈值估计法和启发式阈值估计法去噪比较彻底,在去噪时显得更为有效,但是也容易把有用的高频信号误认为噪声而去除掉。
1994年,Donoho在最小最大估计理论的基础上提出了统一阈值去噪方法,阈值T1的取值由下式确定:
(3)
将式(3)得到的阈值T1应用于小波分解的每一层小波系数进行去噪。噪声随着分解层数的增加,强度逐渐降低,用式(3)的统一阈值去噪方法不能反映噪声的这种变化规律,所以本研究中小波分解的每一层的去噪阈值T2j由下式确定:
(4)
式中,j为进行阈值去噪的小波系数所在的层数,σj为第j层小波分解系数中的噪声标准方差,Nj为第j层小波分解的系数个数。
1.2.4 阈值函数的改进
用于心电信号的小波去噪方法最先是由Donoho和Johnstone等人提出的,目前最常用的有硬阈值法和软阈值法。硬阈值法的原理是设定一个固定阈值,大于该阈值的小波系数保留不变,小于该阈值的小波系数直接置零。这种处理方法能够完整地保留大于阈值部分的信息,但是在阈值附近的小波系数在去噪前后由连续变为阶跃,导致硬阈值去噪后重构得到的信号出现伪吉布斯现象。为了避免出现吉布斯现象,软阈值法提出把大于阈值的小波系数进行相应的收缩,小于阈值的小波系数置零,处理后的效果是:在阈值附近的小波系数是连续的,得到的去噪后信号是平滑的。但是,由于进行了收缩,所以软阈值处理后的小波系数与处理前的小波系数之间存在恒定的差异,使重构后的信号失真,同时造成边缘模糊。
为了尽量完整保留原始信号,避免去噪后出现吉布斯现象,本研究提出了一种新的阈值函数,函数值由下式确定:
(5)
假设Tj=2mV时,新阈值函数和软硬阈值函数在同一个坐标系中的对比如图3所示(见下页)。可以看出,所提出的阈值函数在阈值附近具有连续性,同时随着自变量的增大,新阈值函数对应的函数值与原函数值的差异越来越小。所以,该函数综合了软硬阈值函数各自的优点,同时避免了软硬阈值函数的不足。
图3 新阈值函数和软硬阈值函数波形对比Fig.3 Comparison between new threshold function, soft threshold function and hard threshold function
2.1 小波函数
MIT-BIH心电数据库中的101号数据中的一段1 200个采样点在加入仿真噪声前后的波形如图4所示。
图4 加入仿真噪声前后波形。(a)原始信号;(b)加噪信号Fig.4 Waveform before and after adding the simulation noise. (a)Original signal; (b)Noisy signal
对加噪后的信号进行8层小波分解,分解后用本研究提出的阈值分别进行软阈值、硬阈值去噪。表1为用不同小波函数进行分解重构后信号的信噪比和均方误差的分析结果。
表1 基于多种小波函数的软、硬阈值去噪
Tab.1 Soft threshold denoising and hard threshold denoising based on different wavelet functions
小波SNRMSE加噪软阈值硬阈值加噪软阈值硬阈值Haar98514291294000330002000023Db498514121266000330002000024Db698514621337000330001900022Bior2498514201404000330002000021Bior3598513651399000330002100021Bior4498514431378000330002000021Bior5598515051343000330001800022Bior6898514581532000330001900018Sym398515161323000330001800022Sym498513901243000330002100025Sym598515511433000330001700020Sym698514551336000330001900022Sym898514971486000330001800019
从表1可知,经sym5小波函数分解重构得到的信号信噪信噪比最大,均方误差最小,说明sym5小波函数的去噪效果最好。
2.2 阈值
为了进一步说明本研究使用的阈值的去噪效果,根据前面的研究结果,用sym5小波函数对加噪信号进行8层小波分解,分别用不同的阈值进行去噪,得到的去噪效果如表2所示。表中涉及到的相应阈值分别由下式确定:
(6)
(7)
表2 不同阈值的软阈、硬阈值去噪
Tab.2 Soft threshold denoising and hard threshold denoising based on different thresholds
小波SNRMSE加噪软阈值硬阈值加噪软阈值硬阈值T195214531348000340001900022T2j95214891344000340001900022T3j95210421041000340003100020T4j9521164973000340002700034
从表2可知,当软阈值函数和硬阈值函数使用不同的阈值进行去噪时,所得去噪信号的信噪比和均方误差的差异明显,使用本研究中的阈值T2j得到的去噪效果相对较好。
2.3 阈值函数
本研究用式(4)确定每一层小波细节系数的去噪阈值,并将本研究提出的阈值函数和文献[12]中的4种改进阈值函数用于加噪信号的去噪,得到的去噪结果用信噪比(SNR)和均方误差(MSE)进行对比,如表3所示。
表3 本研究和文献[12]改进阈值函数的去噪效果对比
Tab.3 Comparison of denoising effect between proposed threshold function and improved threshold function in literature [12]
阀值函数SNRMSE硬阈值函数14575000019软阈值函数13802600021文献[12]改进114180400020文献[12]改进214847700019文献[12]改进314685100019文献[12]改进414686700019本研究阈值函数14913100019
从表3可以看出,用本研究提出的阈值函数进行去噪得到的信号信噪比和均方误差相比其他阈值函数得到的结果更好。将本研究的阈值函数波形和表3中除改进3、4外的其他阈值函数波形在同一坐标系中进行对比,如图5所示。观察图5并结合表3的去噪结果,可以看出,阈值函数越接近硬阈值函数,去噪后的信号信噪比越高。所以,设计阈值函数,应该在避免出现伪吉布斯现象的同时尽量减少和硬阈值函数的差异。
图5 不同阈值函数波形对比Fig.5 Waveform comparison of different threshold functions
对心电信号进行离散小波变换时,采用的小波函数不同,对后续的去噪、特征提取等操作的影响非常大。小波分解的系数实质上表征了待分解信号和小波函数的相似性,为了将规律性的心电信号和无规律的干扰噪声分开,应该选择与心电信号波形最相似的小波函数进行操作。常见的小波函数中与心电信号较相似的有db系列、sym系列和coif系列小波,通过表1的分析结果可以看出,sym5小波函数用于小波分解得到的结果更有利于后续的去噪处理。
同时,为了综合考虑分析效果和计算量,小波分解的层数选择也是非常重要的。根据小波分解的规律,由表1可以看出,当小波分解到第8层时,基本上可以将基线漂移噪声提取出来,其他噪声(如工频干扰和肌电干扰)的频率分布范围与心电信号的频率范围重叠,只能通过选择合适的阈值进行去除。
根据噪声的强度随着分解层数的增加而减少的规律,本研究用式(4)来确定每一层的阈值,进而用改进的阈值函数进行阈值去噪,使得去噪后的信号信噪比和均方误差都有较好的结果。
对心电信号进行去噪,最大限度地保留有用信号,真实反映患者的心电状态,对参数测量、波形识别、医疗诊断和治疗具有重要意义。笔者针对小波阈值法去噪过程中小波函数、小波分解层数、阈值和阈值函数等关键参数进行研究,提出了用sym5小波函数进行8层小波分解,并根据有用信号和噪声在每一层的分布规律,对每层的小波分解系数设置不同阈值。最后,用本研究提出的阈值函数进行去噪,其去噪效果比其他去噪方法的去噪效果有所提高。
目前,该算法的处理对象为MIT-BIH心电数据库的心电数据,处理方式为离线处理,下一步将重点提高算法的实时性,并将该算法应用到临床心电数据的去噪处理。
[1] 江依法,周青,叶含笑,等.基于扩散模型的心电信号基线漂移去除法[J].中国生物医学工程学报,2013,32(5):631-635.
[2] 郭新蕾,杨开林,郭永鑫.基于阈值自学习小波算法的压力信号去噪方法[J].数据采集与处理,2008,23(3):322-326.
[3] 苏丽,赵国良,张仁彦.基于改进小波阈值法的平移小波不变心电信号去噪[J].哈尔滨工程大学学报,2006,27(6):839-843.
[4] Zhao Zhidong,Chen Yuquan.Research of threshold value selections and applications for generalized threshold functions[J].Chinese Journal of Sensors and Actuators,2007,20(3):601-605.
[5] Blanco-Velasco M, Weng B, Barner KE.ECG signal denoising and baseline wander correction based on the empirical mode decomposition [J].Computers in Biology and Medicine, 2008,38(1):1-13.
[6] Poornachandra S. Wavelet-based denoising using subband dependent threshold for ECG signals[J].Digital Signal Processing, 2008,18(1):49-55.
[7] Yan Jingyu, Lu Yan, Liu Jia. Self-adaptive model-based ECG denoising using features extracted by mean shift algorithm[J].Biomedical Signal Processing and Control, 2010,5(2):103-113.
[8] Mustaffa I, Trenado C, Schwerdtfeger K, et al. Denoising of single-trial matrix representations using 2D nonlinear diffusion filtering[J].Journal of Neuroscience Methods, 2010,185(2):284-292.
[9] Pal S, Mitra M. Empirical mode decomposition based ECG enhancement and QRS dectection[J].Computers in Biology and Medicine, 2012,42(1):83-92.
[10] 张翠娜.心电信号去噪算法的研究[D].桂林: 桂林理工大学, 2012.
[11] 杨建国.小波分析及其工程应用[M].北京:机械工业出版社,2005.
[12] 刘卫东,刘尚合,胡小锋,等.小波阈值去噪函数的改进方法分析[J].高电压技术,2007, 33(10):59-63.
Research on an Improved Algorithm for Wavelet Denoising of ECG
Zheng Minmin1*Gao Xiaorong2Xie Haihe1
1(SchoolofMechanicalandElectricalEngineering,PutianUniversity,Putian351100,Fujian,China)
ECG; wavelet function; threshold; threshold function; signal to noise ratio
10.3969/j.issn.0258-8021. 2017. 01.015
2016-07-07, 录用日期:2016-08-26
福建省自然科学基金(2013J01251);福建省高校现代精密测量与激光无损检测实验室项目(S20160404);莆田学院支持项目(2015032)
R318
A
0258-8021(2017) 01-0114-05
*通信作者(Corresponding author), E-mail: zhengminmin1016@126.com