一种滤除高采样心电工频干扰的改进算法

王吉鸣 吕颖莹 董 晗 吴 韬 包 涛 冯靖祎*

1(浙江大学医学院附属第一医院医学工程部， 杭州 310003)2(浙江省医药经济发展中心， 杭州 310012)

一种滤除高采样心电工频干扰的改进算法

王吉鸣1吕颖莹1董 晗2吴 韬1包 涛1冯靖祎1*

1(浙江大学医学院附属第一医院医学工程部， 杭州 310003)2(浙江省医药经济发展中心， 杭州 310012)

为了解决高采样频率下50 Hz工频及其倍频对心电波形的影响并保留有用的心电高频信息，提出一种滤除高采样心电工频干扰的改进算法。首先，通过修改零极点得到改进的IIR梳状陷波器，对带工频干扰的ECG信号进行滤波处理；然后，将滤波后的ECG信号进行线性段判别，并对线性段进行平滑滤波处理以消除振铃现象。通过R波幅值、信噪比(SNR)和均方误差(RMSE)对ECG信号工频干扰的消除结果进行评价，与其他4种滤波方法(平滑滤波、Notch滤波、自适应滤波和Levkov滤波)相比，IIR梳状陷波器改进方法得到的ECG信号R波幅值与原始未加工频的ECG信号最为接近，恢复率达到98%以上，且SNR提高约11倍、RMSE降低约30倍。所提出方法不仅能滤除高采样下的工频及其倍频的干扰，且能最大程度保留有用信号。

高采样频率；工频干扰；IIR梳状陷波器；线性段；平滑滤波

引言

在心电图中，大于100 Hz的高频成分对许多心脏疾病(如冠心病、心肌病、心肌炎等)都有重要的诊断价值[1-3]，并已逐渐引起心血管临床工作者的广泛关注。目前，相关的心电高频分析技术如高频心电图、心室晚电位等已在临床上得到应用，且正处于逐步深入与扩大的研究中。为了获取心电的高频信息，需要对心电信号进行高频率采样，但在高采样频率下，会引入更多的噪声，如肌电、工频干扰、随机噪声等，从而对ECG信号的噪声滤除提出了更高的要求。特别地，在工频干扰噪声的滤除中，由于其倍频信号的存在，在滤波器设计中必须做针对性处理。同时，由于高采样下的高尖型QRS波(高幅值、短宽度)信号的频谱与工频干扰的频谱相混叠，在滤波过程中易产生振铃现象[4]，从而对心电图的诊断造成影响。

为解决工频干扰的影响，前期研究者已提出多种滤波算法，如平滑滤波法[5]、简单整系数带阻滤波(Notch滤波)法[6]、自适应滤波法[7]、自适应相干模板法[8]、小波变换法[9]、Levkov法[10]及其改进算法[11-12]等，这些算法均有较好的工频滤除效果，但在高采样下进行心电滤波时都存在一定的局限性：一是平滑滤波法算法简单，处理速度快，滤波效果较好，但对心电波形的R波有较大削峰；二是整系数Notch滤波法具有严格的线性相位且只需较少的计算量，但延时较大且在工频干扰波动时滤波效果下降明显；三是自适应滤波法能够跟随干扰频率和幅度的变化，取得较高的信噪比，但需要一个与噪声有关而与信号无关的参考信号；四是自适应模板法抑制工频干扰的效果较好，但在低频处与窄带干扰信号处的抑制带宽必须相等，导致截止频率高于0.05 Hz的限制；五是小波变换计算量比较大，不适用于实时处理，且小波基和阈值的选择都需要通过大量的实验进行研究论证；六是Levkov法计算量相对小，易于实时处理，但由于其涉及线性段和非线性段的选取，要依据信噪比来决定，较难适应复杂情况下的心电信号处理[13]。

为解决高采样情况下的心电工频干扰滤除方法中存在的问题，本研究基于MATLAB滤波器设计工具箱(FDATool)提出了一种IIR梳状陷波器改进方法，在最大程度保留有效心电高频成分的前提下，能够有效滤除50 Hz工频及其倍频干扰，并通过对ECG波形进行线性段判别和平滑处理，有效地消除了振铃现象的影响。

1 方法

首先对FDATool工具箱设计的IIR梳状陷波器进行了改进，让其仅对50 Hz及其倍频窄带位置陷波，而对0 Hz窄带位置处不进行陷波。然后，针对陷波器在处理具有高尖型QRS波的ECG信号时容易引入振铃现象的问题，提出了一种ECG信号线性段判别方法，对经过IIR梳状陷波器处理的ECG信号的线性段和非线性段分别进行相应处理，进而完成工频干扰的完整滤除。总体流程如图1所示。

图1 去工频流程Fig.1 The process of elimination of power-line interference

1.1 IIR梳状陷波器原理及其改进

由于数字滤波器的频率响应就是其单位冲激响应在单位圆上的Z变换，在单位圆上相应于所需滤波器阻带位置的频率处设置零点可以使滤波器的频率响应在所需阻带频率处为零。同时为了得到非常陡峭的过渡带和常数幅度的通带特性，必须在Z平面上为每一个零点再配置一个相应的极点。

基于FDATool工具设计IIR梳状陷波器的传递函数，因为每一个极点都有一个零点与之相对应，所以能够用于抑制某一频率成分及其谐波，有

(1)

式中，N为滤波器阶数，由采样频率和工频决定，参数a、b的计算方法为

(2)

(3)

由于IIR梳状陷波器对0 Hz窄带成分也有抑制作用，且抑制带宽与干扰信号抑制带宽相等，而在ECG信号中这部分信号成分需要单独处理，不能被滤除，因此需要对该IIR陷波器进行改进。

对式(1)进行分解，可以得到式(6)，z=1处的零点对应的是0 Hz阻带，z=a1/N为相应极点。去掉0 Hz对应的零极点，可以消除陷波器对0 Hz处的影响，最后改进的IIR梳状陷波器的传递函数如下：

(6)

(7)

1.2 ECG线性段判别和处理

实际的信号处理过程中振铃现象的产生具有两个条件：从滤波器的角度来讲是因为其具有单位脉冲响应；从信号的角度来说是因为输入为高尖信号[14]，其表现为持续一段时间并逐渐衰减的纹波。由于心电波形的QRS波又高又尖，类似于脉冲响应，用陷波器对其处理必然引起振铃效应，从而影响滤波结果。

为消除振铃现象，本研究提出一种ECG线性段判别方法，并对得到的信号线性段进行平滑滤波处理。具体过程如图2所示，其中阈值T由前一个R波幅值决定，取其1/4。当检测到一个线性段后，判断前面数据段是否为线性段，若是则对连续线性段进行N点平滑滤波，否则仅对该线性段进行N点平滑滤波。当检测到一个非线性段后，则不作任何处理，保持原来的滤波结果。

图2 ECG波形线性段判别标准Fig.2 The linear criterion of ECG waveform

1.3 工频干扰的去噪效果评价

为更加准确地评价工频干扰的去噪效果，本研究采用R波幅值、信噪比(signal to noise ratio，SNR)及标准均方误差(root mean square error，RMSE)等指标进行定量分析。QRS波幅值作为心电图诊断的一项指标，当R波被削峰太多，有些疾病(如心肌梗死)将无法诊断，滤波过程中保留原始R波幅值程度越大，保证相关分析的准确性就越高。

SNR和RMSE是两个最基本的信号处理统计评价指标，其计算公式[15]为

(8)

(9)

式中，xi表示原始含工频干扰的ECG信号，yi表示去工频干扰后的ECG信号，N表示信号长度。

一般地，SNR越高，RMSE越小，表明去噪后的信号就越接近原始信号，去噪的效果也就越好。

2 结果

测试数据均来源于FlukeMPS450标准心电采集设备，其参数设置为采样频率1 000Hz、幅值4mV、心率80。测试数据分两段，第1段是未受工频干扰的原始ECG信号，第2段是在原始信号上叠加50Hz工频干扰的ECG信号。

2.1 IIR梳状陷波器的改进结果

对于采样频率1 000 Hz、工频50 Hz的信号滤波，IIR梳状陷波器的幅频响应和零极点分布如图3所示；改进的IIR梳状陷波器，其幅频响应和零极点分布如图4所示；可以看出，0 Hz窄带阻滞消失。

图3 IIR梳状陷波器。(a)幅频响应；(b)零极点分布Fig.3 IIR notch filter. (a)Amplitude-frequency response； (b)Pole-zero diagram

图4 改进IIR梳状陷波器。(a)幅频响应；(b)零极点分布Fig.4 Improved IIR notch filter. (a)Amplitude-frequency response; (b)Pole-zero diagram

图5 心电波形滤波前后比较。(a) 原始ECG；(b) 工频ECG；(c)改进的IIR梳状陷波滤波和ECG；(d)改进的IIR梳状陷波器处理后再经过线性段平滑后的ECGFig.5 ECG waveform comparison before and after filtering.(a)original ECG；(b) ECG with power-line interference；(c) ECG afterimproved IIR notch filter；(d) ECG after improved IIR notch filter and smoothing filter

2.2 滤波结果

图5显示了本滤波方法的滤波结果，其中(a)为原始ECG波形，(b)为叠加50 Hz工频干扰的ECG波形，(c)为经过改进的IIR梳状陷波器后的ECG波形，(d)为经过改进的IIR梳状陷波器处理后再经过线性段平滑后的ECG波形。从(c)可看到，经过改进的IIR梳状陷波器处理后的50 Hz工频干扰得到了很好的抑制，但是在QRS波后的心电波形出现了振铃现象。从(d)可看出，经过线性段判别和处理后，振铃现象基本得到抑制。

2.3 与其他滤波算法结果比较

图6 不同方法处理后的ECG波形。(a)原始ECG；(b)工频ECG；(c)平滑滤波后ECG；(d)Notch滤波后ECG；(e)自适应滤波后ECG；(f)Levkov滤波后ECG；(g)本方法滤波后ECGFig.6 ECG waveform after different methods.(a)Original ECG；(b) ECG with power-line Interference；(c) ECG after smoothing filter；(d) ECG after notch filter；(e) ECG after adaptive filter；(f) ECG after Levkov filter；(g)ECG after the method in this article

图6、7分别为本方法与平滑滤波、Notch滤波、自适应滤波及Levkov滤波对一段和一个周期的ECG信号进行处理的结果比较。从中可看出，平滑滤波和Levkov滤波极大地削弱了R波幅值信息，而Notch滤波和自适应滤波虽然保留了大部分的R波幅值信息，但仍存在一部分噪声干扰。相比其他方法，本滤波方法能够在获得较好的滤波效果的同时，最大程度地保留R波幅值信息。

图7 不同方法处理后的ECG波形(一个周期)。(a)原始ECG；(b)工频ECG；(c)平滑滤波后ECG；(d)Notch滤波后ECG；(e)自适应滤波后ECG；(f)Levkov滤波后ECG；(g)本方法滤波后ECGFig.7 ECG waveform after filtering by different methods(one period). (a)Original ECG；(b) ECG with power-line interference；(c) ECG after smoothing filter；(d) ECG after notch filter；(e) ECG after adaptive filter；(f) ECG after Levkov filter；(g)ECG after the method in this article

2.4 去噪结果评价

表1给出了不同方法处理后的R波幅值、SNR和RMSE等参数对比。可以看出，本方法、Notch滤波法及自适应滤波方法均能最大化保持R波幅值信息，R波幅值恢复率达到98%以上，但是本方法在SNR和RMSE上表现更佳，SNR提高了约11倍而RMSE降低了约30倍，说明有更好的信号去噪能力。

表1 不同方法处理后的ECG信号的R波幅值、SNR和RMSE

Tab.1 Amplitude of R wave SNR and MSE of ECG signals after different treatments

不同处理方法R波幅值/mVSNR/dBRMSE加工频后4262150651平滑滤波后2428780164Notch滤波后39420170051自适应滤波后39020230051Levkov滤波后23272240191本方法滤波后395232140024

3 讨论和结论

在心电分析中，为了提取高频成分，需要更高的采样频率。但采样频率提高后，需要考虑工频干扰的倍频成分对心电波形的影响，使用传统的滤波方法，相对低频的工频陷波的计算量会大大增加。本研究根据滤除50 Hz工频及其倍频并保留心电高频信息的要求，采用了IIR梳状陷波器，并通过改变零极点来改变陷波器结构以消除对0 Hz窄带的影响，从而达到滤波效果，且实现简单。与传统的很多工频干扰滤波方法相比，本方法仅仅滤除了心电信号中的工频及其倍频窄带而保留了其他所有频率信息，以供后续分析。局限是IIR梳状陷波器的阶数是由采样频率和工频决定的，所以要求采样频率必须为工频的整倍数。

本研究对滤波后的振铃现象处理方法参考了Levkov法对线性段进行平均滤波的思想，不同的是先对采集的心电信号进行了窄带宽50 Hz工频陷波，以便提高信号的信噪比并使R波削峰程度最小，再根据前面提到的ECG波形线性段判别标准来确定线性段并进行平滑滤波。这样的处理方式在最大化保持R波幅值信息的同时比Levkov法直接根据信噪比来判定线性段和非线性段的方法稳定很多，进而去振铃的方法也比较稳定。

另外，在心电波形的实际采集过程中工频频率并不是固定不变的，工频波动会影响滤波结果。本研究提出的线性段平滑滤波方法不仅能消除高尖ECG信号滤波后的振铃现象，相同原理对由于工频干扰波动引起的滤波后出现的振铃现象也有很好的抑制作用。如图8所示，A、B、C所指分别为工频波动处、陷波后振荡及线性段平滑滤波对振荡的抑制，可以看出，C相对于B，振荡小了很多。

图8 工频干扰波动(A—工频波动；B—陷波后振荡；C—线性段平滑滤波对振荡的抑制)Fig.8 Power-line interference fluctuation (A: Power frequency fluctuation, B: Oscillation after notch filter, C: Inhibition of oscillation by smoothing the linear segment)

心电高频成分能够提供对疾病有诊断价值的信息，但在对心电高采样的同时会引入更多的噪声，从而对ECG信号的噪声滤除方法有更高的要求。作为其中的一个重要干扰因素，工频干扰在高采样频率下的滤除仍存在问题。本研究提出一种改进的工频干扰滤除算法，通过结合改进的IIR梳状陷波器以及线性段判别处理，可在最大化保持R波幅值信息的前提下，有效地滤除ECG信号中的工频及其倍频成分，并可消除由于陷波器造成的振铃现象，从而为后续心电的准确分析提供了保证。

[1] 张开滋，肖传实，郭继红，等.中国心电信息学图解集成：下册[M]. 长沙：湖南科技出版社，2010：1656.

[2] 盛海民．心电高频信息的临床应用[J]. 医学综述，1996，2(6)：269-270.

[3] 刘恒亮，单新中．心电高频信息研究的现状与展望[J]. 心血管病学进展，2000，21(1)：52-55.

[4] 兰瑞芬，胡广书．高采样率下简单整系数工频陷波器的设计[J]. 航天医学与医学工程，2008(2)：152-156.

[5] 王镇，蔡萍．用于去除心电信号中工频干扰数字滤波技术[J]. 电子测量技术，2000(2)：43-47.

[6] 周静．心电信号中工频干扰的消除[J]. 生物医学工程研究，2003，22(4)：61-64.

[7] 林绍杰，赖丽娟，吴效明．基于阻抗检测的自适应消除心电运动伪迹方法[J]. 生物医学工程学杂志，2010，27(3)：529-532.

[8] 于超，李刚，林凌．基于简化FFT可跟踪工频干扰的自适应相干模板法的实现[J]. 生物医学工程学杂志，2007，24(4)：780-784.

[9] 张雪，王海燕，李保军．临床心电信号工频干扰小波去噪方法对比分析[J]. 计算机测量与控制，2010，18(4)：902-905.

[10] Levkov C, Michov G, Ivanov R, et al. Subtraction of 50Hz interference from the electrocardiogram [J]. Med Bio Eng Compute, 1984, 22(4):371-373.

[11] 孙京霞，白延强，杨玉星. 一种抑制心电信号50 Hz 工频干扰的改进Levkov方法[J].航天医学与医学工程.2000，13(3)：196-199.

[12] Levkov C, Michov G, Ivanov R, et al. Removal of power-line interference from the ECG: a review of the subtraction produre[J]. BioMedical Engineering OnLine 2005, 4(1): 18-36.

[13] 刘胜洋，张根选，曹阳．心电信号中的工频干扰滤除方法[J].生物医学工程学杂志，2014(3)：577-582.

[14] 罗旻晖．心电图机自动分析算法的研究与实现[D]. 哈尔滨：哈尔滨工业大学，2012.

[15] 张喜红，王玉香．心电信号工频干扰陷波器的设计与实现[J]. 湖南工业大学学报，2014(2)：72-76.

An Improved Algorithm for Removing the Power-Line Interference from ECG Signals in High Sampling Rate

Wang Jiming1Lv Yingying1Dong Han2Wu Tao1Bao Tao1Feng Jingyi1*

1(ClinicalEngineeringDepartment，TheFirstAffiliatedHospital，SchoolofMedicine，ZhejiangUniversity，Hangzhou310003，China)2(ZhejiangPharmaceuticalEconomicDevelopmentCenter，Hangzhou310012，China)

high sampling rate; power-line interference; IIR comb notch filter; linear segment; smoothing filter

10.3969/j.issn.0258-8021. 2016. 06.015

2015-12-20， 录用日期:2016-08-21

R318

D

0258-8021(2016) 06-0744-05

# 中国生物医学工程学会高级会员(Senior member, Chinese Society of Biomedical Engineering)

*通信作者(Corresponding author)， E-mail: casper_feng@163.com