王吉鳴 呂穎瑩 董 晗 吳 韜 包 濤 馮靖祎*

1(浙江大學(xué)醫(yī)學(xué)院附屬第一醫(yī)院醫(yī)學(xué)工程部， 杭州 310003)2(浙江省醫(yī)藥經(jīng)濟(jì)發(fā)展中心， 杭州 310012)

一種濾除高采樣心電工頻干擾的改進(jìn)算法

王吉鳴1呂穎瑩1董 晗2吳 韜1包 濤1馮靖祎1*

1(浙江大學(xué)醫(yī)學(xué)院附屬第一醫(yī)院醫(yī)學(xué)工程部， 杭州 310003)2(浙江省醫(yī)藥經(jīng)濟(jì)發(fā)展中心， 杭州 310012)

為了解決高采樣頻率下50 Hz工頻及其倍頻對(duì)心電波形的影響并保留有用的心電高頻信息，提出一種濾除高采樣心電工頻干擾的改進(jìn)算法。首先，通過(guò)修改零極點(diǎn)得到改進(jìn)的IIR梳狀陷波器，對(duì)帶工頻干擾的ECG信號(hào)進(jìn)行濾波處理；然后，將濾波后的ECG信號(hào)進(jìn)行線性段判別，并對(duì)線性段進(jìn)行平滑濾波處理以消除振鈴現(xiàn)象。通過(guò)R波幅值、信噪比(SNR)和均方誤差(RMSE)對(duì)ECG信號(hào)工頻干擾的消除結(jié)果進(jìn)行評(píng)價(jià)，與其他4種濾波方法(平滑濾波、Notch濾波、自適應(yīng)濾波和Levkov濾波)相比，IIR梳狀陷波器改進(jìn)方法得到的ECG信號(hào)R波幅值與原始未加工頻的ECG信號(hào)最為接近，恢復(fù)率達(dá)到98%以上，且SNR提高約11倍、RMSE降低約30倍。所提出方法不僅能濾除高采樣下的工頻及其倍頻的干擾，且能最大程度保留有用信號(hào)。

高采樣頻率；工頻干擾；IIR梳狀陷波器；線性段；平滑濾波

引言

在心電圖中，大于100 Hz的高頻成分對(duì)許多心臟疾病(如冠心病、心肌病、心肌炎等)都有重要的診斷價(jià)值[1-3]，并已逐漸引起心血管臨床工作者的廣泛關(guān)注。目前，相關(guān)的心電高頻分析技術(shù)如高頻心電圖、心室晚電位等已在臨床上得到應(yīng)用，且正處于逐步深入與擴(kuò)大的研究中。為了獲取心電的高頻信息，需要對(duì)心電信號(hào)進(jìn)行高頻率采樣，但在高采樣頻率下，會(huì)引入更多的噪聲，如肌電、工頻干擾、隨機(jī)噪聲等，從而對(duì)ECG信號(hào)的噪聲濾除提出了更高的要求。特別地，在工頻干擾噪聲的濾除中，由于其倍頻信號(hào)的存在，在濾波器設(shè)計(jì)中必須做針對(duì)性處理。同時(shí)，由于高采樣下的高尖型QRS波(高幅值、短寬度)信號(hào)的頻譜與工頻干擾的頻譜相混疊，在濾波過(guò)程中易產(chǎn)生振鈴現(xiàn)象[4]，從而對(duì)心電圖的診斷造成影響。

為解決工頻干擾的影響，前期研究者已提出多種濾波算法，如平滑濾波法[5]、簡(jiǎn)單整系數(shù)帶阻濾波(Notch濾波)法[6]、自適應(yīng)濾波法[7]、自適應(yīng)相干模板法[8]、小波變換法[9]、Levkov法[10]及其改進(jìn)算法[11-12]等，這些算法均有較好的工頻濾除效果，但在高采樣下進(jìn)行心電濾波時(shí)都存在一定的局限性：一是平滑濾波法算法簡(jiǎn)單，處理速度快，濾波效果較好，但對(duì)心電波形的R波有較大削峰；二是整系數(shù)Notch濾波法具有嚴(yán)格的線性相位且只需較少的計(jì)算量，但延時(shí)較大且在工頻干擾波動(dòng)時(shí)濾波效果下降明顯；三是自適應(yīng)濾波法能夠跟隨干擾頻率和幅度的變化，取得較高的信噪比，但需要一個(gè)與噪聲有關(guān)而與信號(hào)無(wú)關(guān)的參考信號(hào)；四是自適應(yīng)模板法抑制工頻干擾的效果較好，但在低頻處與窄帶干擾信號(hào)處的抑制帶寬必須相等，導(dǎo)致截止頻率高于0.05 Hz的限制；五是小波變換計(jì)算量比較大，不適用于實(shí)時(shí)處理，且小波基和閾值的選擇都需要通過(guò)大量的實(shí)驗(yàn)進(jìn)行研究論證；六是Levkov法計(jì)算量相對(duì)小，易于實(shí)時(shí)處理，但由于其涉及線性段和非線性段的選取，要依據(jù)信噪比來(lái)決定，較難適應(yīng)復(fù)雜情況下的心電信號(hào)處理[13]。

為解決高采樣情況下的心電工頻干擾濾除方法中存在的問(wèn)題，本研究基于MATLAB濾波器設(shè)計(jì)工具箱(FDATool)提出了一種IIR梳狀陷波器改進(jìn)方法，在最大程度保留有效心電高頻成分的前提下，能夠有效濾除50 Hz工頻及其倍頻干擾，并通過(guò)對(duì)ECG波形進(jìn)行線性段判別和平滑處理，有效地消除了振鈴現(xiàn)象的影響。

1 方法

首先對(duì)FDATool工具箱設(shè)計(jì)的IIR梳狀陷波器進(jìn)行了改進(jìn)，讓其僅對(duì)50 Hz及其倍頻窄帶位置陷波，而對(duì)0 Hz窄帶位置處不進(jìn)行陷波。然后，針對(duì)陷波器在處理具有高尖型QRS波的ECG信號(hào)時(shí)容易引入振鈴現(xiàn)象的問(wèn)題，提出了一種ECG信號(hào)線性段判別方法，對(duì)經(jīng)過(guò)IIR梳狀陷波器處理的ECG信號(hào)的線性段和非線性段分別進(jìn)行相應(yīng)處理，進(jìn)而完成工頻干擾的完整濾除?？傮w流程如圖1所示。

圖1 去工頻流程Fig.1 The process of elimination of power-line interference

1.1 IIR梳狀陷波器原理及其改進(jìn)

由于數(shù)字濾波器的頻率響應(yīng)就是其單位沖激響應(yīng)在單位圓上的Z變換，在單位圓上相應(yīng)于所需濾波器阻帶位置的頻率處設(shè)置零點(diǎn)可以使濾波器的頻率響應(yīng)在所需阻帶頻率處為零。同時(shí)為了得到非常陡峭的過(guò)渡帶和常數(shù)幅度的通帶特性，必須在Z平面上為每一個(gè)零點(diǎn)再配置一個(gè)相應(yīng)的極點(diǎn)。

基于FDATool工具設(shè)計(jì)IIR梳狀陷波器的傳遞函數(shù)，因?yàn)槊恳粋€(gè)極點(diǎn)都有一個(gè)零點(diǎn)與之相對(duì)應(yīng)，所以能夠用于抑制某一頻率成分及其諧波，有

(1)

式中，N為濾波器階數(shù)，由采樣頻率和工頻決定，參數(shù)a、b的計(jì)算方法為

(2)

(3)

由于IIR梳狀陷波器對(duì)0 Hz窄帶成分也有抑制作用，且抑制帶寬與干擾信號(hào)抑制帶寬相等，而在ECG信號(hào)中這部分信號(hào)成分需要單獨(dú)處理，不能被濾除，因此需要對(duì)該IIR陷波器進(jìn)行改進(jìn)。

對(duì)式(1)進(jìn)行分解，可以得到式(6)，z=1處的零點(diǎn)對(duì)應(yīng)的是0 Hz阻帶，z=a1/N為相應(yīng)極點(diǎn)。去掉0 Hz對(duì)應(yīng)的零極點(diǎn)，可以消除陷波器對(duì)0 Hz處的影響，最后改進(jìn)的IIR梳狀陷波器的傳遞函數(shù)如下：

(6)

(7)

1.2 ECG線性段判別和處理

實(shí)際的信號(hào)處理過(guò)程中振鈴現(xiàn)象的產(chǎn)生具有兩個(gè)條件：從濾波器的角度來(lái)講是因?yàn)槠渚哂袉挝幻}沖響應(yīng)；從信號(hào)的角度來(lái)說(shuō)是因?yàn)檩斎霝楦呒庑盘?hào)[14]，其表現(xiàn)為持續(xù)一段時(shí)間并逐漸衰減的紋波。由于心電波形的QRS波又高又尖，類似于脈沖響應(yīng)，用陷波器對(duì)其處理必然引起振鈴效應(yīng)，從而影響濾波結(jié)果。

為消除振鈴現(xiàn)象，本研究提出一種ECG線性段判別方法，并對(duì)得到的信號(hào)線性段進(jìn)行平滑濾波處理。具體過(guò)程如圖2所示，其中閾值T由前一個(gè)R波幅值決定，取其1/4。當(dāng)檢測(cè)到一個(gè)線性段后，判斷前面數(shù)據(jù)段是否為線性段，若是則對(duì)連續(xù)線性段進(jìn)行N點(diǎn)平滑濾波，否則僅對(duì)該線性段進(jìn)行N點(diǎn)平滑濾波。當(dāng)檢測(cè)到一個(gè)非線性段后，則不作任何處理，保持原來(lái)的濾波結(jié)果。

圖2 ECG波形線性段判別標(biāo)準(zhǔn)Fig.2 The linear criterion of ECG waveform

1.3 工頻干擾的去噪效果評(píng)價(jià)

為更加準(zhǔn)確地評(píng)價(jià)工頻干擾的去噪效果，本研究采用R波幅值、信噪比(signal to noise ratio，SNR)及標(biāo)準(zhǔn)均方誤差(root mean square error，RMSE)等指標(biāo)進(jìn)行定量分析。QRS波幅值作為心電圖診斷的一項(xiàng)指標(biāo)，當(dāng)R波被削峰太多，有些疾病(如心肌梗死)將無(wú)法診斷，濾波過(guò)程中保留原始R波幅值程度越大，保證相關(guān)分析的準(zhǔn)確性就越高。

SNR和RMSE是兩個(gè)最基本的信號(hào)處理統(tǒng)計(jì)評(píng)價(jià)指標(biāo)，其計(jì)算公式[15]為

(8)

(9)

式中，xi表示原始含工頻干擾的ECG信號(hào)，yi表示去工頻干擾后的ECG信號(hào)，N表示信號(hào)長(zhǎng)度。

一般地，SNR越高，RMSE越小，表明去噪后的信號(hào)就越接近原始信號(hào)，去噪的效果也就越好。

2 結(jié)果

測(cè)試數(shù)據(jù)均來(lái)源于FlukeMPS450標(biāo)準(zhǔn)心電采集設(shè)備，其參數(shù)設(shè)置為采樣頻率1 000Hz、幅值4mV、心率80。測(cè)試數(shù)據(jù)分兩段，第1段是未受工頻干擾的原始ECG信號(hào)，第2段是在原始信號(hào)上疊加50Hz工頻干擾的ECG信號(hào)。

2.1 IIR梳狀陷波器的改進(jìn)結(jié)果

對(duì)于采樣頻率1 000 Hz、工頻50 Hz的信號(hào)濾波，IIR梳狀陷波器的幅頻響應(yīng)和零極點(diǎn)分布如圖3所示；改進(jìn)的IIR梳狀陷波器，其幅頻響應(yīng)和零極點(diǎn)分布如圖4所示；可以看出，0 Hz窄帶阻滯消失。

圖3 IIR梳狀陷波器。(a)幅頻響應(yīng)；(b)零極點(diǎn)分布Fig.3 IIR notch filter. (a)Amplitude-frequency response； (b)Pole-zero diagram

圖4 改進(jìn)IIR梳狀陷波器。(a)幅頻響應(yīng)；(b)零極點(diǎn)分布Fig.4 Improved IIR notch filter. (a)Amplitude-frequency response; (b)Pole-zero diagram

圖5 心電波形濾波前后比較。(a) 原始ECG；(b) 工頻ECG；(c)改進(jìn)的IIR梳狀陷波濾波和ECG；(d)改進(jìn)的IIR梳狀陷波器處理后再經(jīng)過(guò)線性段平滑后的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 濾波結(jié)果

圖5顯示了本濾波方法的濾波結(jié)果，其中(a)為原始ECG波形，(b)為疊加50 Hz工頻干擾的ECG波形，(c)為經(jīng)過(guò)改進(jìn)的IIR梳狀陷波器后的ECG波形，(d)為經(jīng)過(guò)改進(jìn)的IIR梳狀陷波器處理后再經(jīng)過(guò)線性段平滑后的ECG波形。從(c)可看到，經(jīng)過(guò)改進(jìn)的IIR梳狀陷波器處理后的50 Hz工頻干擾得到了很好的抑制，但是在QRS波后的心電波形出現(xiàn)了振鈴現(xiàn)象。從(d)可看出，經(jīng)過(guò)線性段判別和處理后，振鈴現(xiàn)象基本得到抑制。

2.3 與其他濾波算法結(jié)果比較

圖6 不同方法處理后的ECG波形。(a)原始ECG；(b)工頻ECG；(c)平滑濾波后ECG；(d)Notch濾波后ECG；(e)自適應(yīng)濾波后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濾波、自適應(yīng)濾波及Levkov濾波對(duì)一段和一個(gè)周期的ECG信號(hào)進(jìn)行處理的結(jié)果比較。從中可看出，平滑濾波和Levkov濾波極大地削弱了R波幅值信息，而Notch濾波和自適應(yīng)濾波雖然保留了大部分的R波幅值信息，但仍存在一部分噪聲干擾。相比其他方法，本濾波方法能夠在獲得較好的濾波效果的同時(shí)，最大程度地保留R波幅值信息。

圖7 不同方法處理后的ECG波形(一個(gè)周期)。(a)原始ECG；(b)工頻ECG；(c)平滑濾波后ECG；(d)Notch濾波后ECG；(e)自適應(yīng)濾波后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 去噪結(jié)果評(píng)價(jià)

表1給出了不同方法處理后的R波幅值、SNR和RMSE等參數(shù)對(duì)比?？梢钥闯?，本方法、Notch濾波法及自適應(yīng)濾波方法均能最大化保持R波幅值信息，R波幅值恢復(fù)率達(dá)到98%以上，但是本方法在SNR和RMSE上表現(xiàn)更佳，SNR提高了約11倍而RMSE降低了約30倍，說(shuō)明有更好的信號(hào)去噪能力。

表1 不同方法處理后的ECG信號(hào)的R波幅值、SNR和RMSE

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

不同處理方法R波幅值/mVSNR/dBRMSE加工頻后4262150651平滑濾波后2428780164Notch濾波后39420170051自適應(yīng)濾波后39020230051Levkov濾波后23272240191本方法濾波后395232140024

3 討論和結(jié)論

在心電分析中，為了提取高頻成分，需要更高的采樣頻率。但采樣頻率提高后，需要考慮工頻干擾的倍頻成分對(duì)心電波形的影響，使用傳統(tǒng)的濾波方法，相對(duì)低頻的工頻陷波的計(jì)算量會(huì)大大增加。本研究根據(jù)濾除50 Hz工頻及其倍頻并保留心電高頻信息的要求，采用了IIR梳狀陷波器，并通過(guò)改變零極點(diǎn)來(lái)改變陷波器結(jié)構(gòu)以消除對(duì)0 Hz窄帶的影響，從而達(dá)到濾波效果，且實(shí)現(xiàn)簡(jiǎn)單。與傳統(tǒng)的很多工頻干擾濾波方法相比，本方法僅僅濾除了心電信號(hào)中的工頻及其倍頻窄帶而保留了其他所有頻率信息，以供后續(xù)分析。局限是IIR梳狀陷波器的階數(shù)是由采樣頻率和工頻決定的，所以要求采樣頻率必須為工頻的整倍數(shù)。

本研究對(duì)濾波后的振鈴現(xiàn)象處理方法參考了Levkov法對(duì)線性段進(jìn)行平均濾波的思想，不同的是先對(duì)采集的心電信號(hào)進(jìn)行了窄帶寬50 Hz工頻陷波，以便提高信號(hào)的信噪比并使R波削峰程度最小，再根據(jù)前面提到的ECG波形線性段判別標(biāo)準(zhǔn)來(lái)確定線性段并進(jìn)行平滑濾波。這樣的處理方式在最大化保持R波幅值信息的同時(shí)比Levkov法直接根據(jù)信噪比來(lái)判定線性段和非線性段的方法穩(wěn)定很多，進(jìn)而去振鈴的方法也比較穩(wěn)定。

另外，在心電波形的實(shí)際采集過(guò)程中工頻頻率并不是固定不變的，工頻波動(dòng)會(huì)影響濾波結(jié)果。本研究提出的線性段平滑濾波方法不僅能消除高尖ECG信號(hào)濾波后的振鈴現(xiàn)象，相同原理對(duì)由于工頻干擾波動(dòng)引起的濾波后出現(xiàn)的振鈴現(xiàn)象也有很好的抑制作用。如圖8所示，A、B、C所指分別為工頻波動(dòng)處、陷波后振蕩及線性段平滑濾波對(duì)振蕩的抑制，可以看出，C相對(duì)于B，振蕩小了很多。

圖8 工頻干擾波動(dòng)(A—工頻波動(dòng)；B—陷波后振蕩；C—線性段平滑濾波對(duì)振蕩的抑制)Fig.8 Power-line interference fluctuation (A: Power frequency fluctuation, B: Oscillation after notch filter, C: Inhibition of oscillation by smoothing the linear segment)

心電高頻成分能夠提供對(duì)疾病有診斷價(jià)值的信息，但在對(duì)心電高采樣的同時(shí)會(huì)引入更多的噪聲，從而對(duì)ECG信號(hào)的噪聲濾除方法有更高的要求。作為其中的一個(gè)重要干擾因素，工頻干擾在高采樣頻率下的濾除仍存在問(wèn)題。本研究提出一種改進(jìn)的工頻干擾濾除算法，通過(guò)結(jié)合改進(jìn)的IIR梳狀陷波器以及線性段判別處理，可在最大化保持R波幅值信息的前提下，有效地濾除ECG信號(hào)中的工頻及其倍頻成分，并可消除由于陷波器造成的振鈴現(xiàn)象，從而為后續(xù)心電的準(zhǔn)確分析提供了保證。

[1] 張開(kāi)滋，肖傳實(shí)，郭繼紅，等.中國(guó)心電信息學(xué)圖解集成：下冊(cè)[M]. 長(zhǎng)沙：湖南科技出版社，2010：1656.

[2] 盛海民．心電高頻信息的臨床應(yīng)用[J]. 醫(yī)學(xué)綜述，1996，2(6)：269-270.

[3] 劉恒亮，單新中．心電高頻信息研究的現(xiàn)狀與展望[J]. 心血管病學(xué)進(jìn)展，2000，21(1)：52-55.

[4] 蘭瑞芬，胡廣書(shū)．高采樣率下簡(jiǎn)單整系數(shù)工頻陷波器的設(shè)計(jì)[J]. 航天醫(yī)學(xué)與醫(yī)學(xué)工程，2008(2)：152-156.

[5] 王鎮(zhèn)，蔡萍．用于去除心電信號(hào)中工頻干擾數(shù)字濾波技術(shù)[J]. 電子測(cè)量技術(shù)，2000(2)：43-47.

[6] 周靜．心電信號(hào)中工頻干擾的消除[J]. 生物醫(yī)學(xué)工程研究，2003，22(4)：61-64.

[7] 林紹杰，賴麗娟，吳效明．基于阻抗檢測(cè)的自適應(yīng)消除心電運(yùn)動(dòng)偽跡方法[J]. 生物醫(yī)學(xué)工程學(xué)雜志，2010，27(3)：529-532.

[8] 于超，李剛，林凌．基于簡(jiǎn)化FFT可跟蹤工頻干擾的自適應(yīng)相干模板法的實(shí)現(xiàn)[J]. 生物醫(yī)學(xué)工程學(xué)雜志，2007，24(4)：780-784.

[9] 張雪，王海燕，李保軍．臨床心電信號(hào)工頻干擾小波去噪方法對(duì)比分析[J]. 計(jì)算機(jī)測(cè)量與控制，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] 孫京霞，白延強(qiáng)，楊玉星. 一種抑制心電信號(hào)50 Hz 工頻干擾的改進(jìn)Levkov方法[J].航天醫(yī)學(xué)與醫(yī)學(xué)工程.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] 劉勝洋，張根選，曹陽(yáng)．心電信號(hào)中的工頻干擾濾除方法[J].生物醫(yī)學(xué)工程學(xué)雜志，2014(3)：577-582.

[14] 羅旻暉．心電圖機(jī)自動(dòng)分析算法的研究與實(shí)現(xiàn)[D]. 哈爾濱：哈爾濱工業(yè)大學(xué)，2012.

[15] 張喜紅，王玉香．心電信號(hào)工頻干擾陷波器的設(shè)計(jì)與實(shí)現(xiàn)[J]. 湖南工業(yè)大學(xué)學(xué)報(bào)，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

# 中國(guó)生物醫(yī)學(xué)工程學(xué)會(huì)高級(jí)會(huì)員(Senior member, Chinese Society of Biomedical Engineering)

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