張光普，陳日清，付軼帆，吳 劍

(清華大學 深圳研究生院 深圳市無損監(jiān)測與微創(chuàng)醫(yī)學技術(shù)重點實驗室，廣東 深圳 518055)

0 引 言

心電(electrocardiography,ECG)信號是微弱的電生理信號，易受人體自身、外界環(huán)境及測量儀器的影響，產(chǎn)生如工頻噪聲、基線漂移、極化電壓和肌電干擾等類型的干擾[1]。此外多通道(32導(dǎo)聯(lián)及以上)心內(nèi)ECG采集時單通道的采樣率可達2 kHz，增加了數(shù)據(jù)傳輸和存儲的難度。目前，ECG去噪有自適應(yīng)濾波、有限長單位沖激響應(yīng)(finite impulse response,FIR)濾波、模糊邏輯、經(jīng)驗?zāi)Ｊ椒纸夂托〔ㄗ儞Q等方法[2]。其中離散小波變換(discret wavelet transform,DWT)因不需先驗知識，能根據(jù)信號自身選取閾值的特性，近年來被廣泛應(yīng)用于ECG濾波中[3,4]。本文采用基于FIR低通與DWT相結(jié)合的方法，前者在濾波的同時縮短了信號的帶寬，有利于DWT對頻帶進行更為精細的分解和濾波。ECG信號有直接壓縮、變換域壓縮和特征參數(shù)提取壓縮等3種方法[5]。其中變換域壓縮可以去除信號相關(guān)性并減小冗余，因而備受關(guān)注，有Karhunen-Loeve變換(KLT)，傅立葉變換(Fourier transform,FT)，余弦變換(cosine transform,CT)和小波變換(wavelet transform,WT)[6～10]等。其中小波變換具有良好的時頻定位和能量集中特性，本文采用基于DWT的SPIHT算法對ECG信號進行壓縮。SPIHT算法在圖像壓縮領(lǐng)域取得了巨大的成功[11]，經(jīng)改造后對ECG的壓縮也取得了良好的效果[12]，算法原理包括：用子集劃分搜索算法對小波系數(shù)按幅值進行部分排列，利用小波系數(shù)層間的相關(guān)性假設(shè)和順序輸出系數(shù)的比特位[13]。

針對ECG去噪和壓縮存在算法復(fù)雜、處理數(shù)據(jù)量大和計算時間長的問題，提出了一種基于現(xiàn)場可編程門陣列(field programmable gate array,FPGA)的方法來實現(xiàn)多通道ECG信號的實時濾波和壓縮。

1 方 法

基于FPGA對ECG進行濾波和壓縮的實現(xiàn)方法如圖1所示，采用8個8-1多路開關(guān)配合8輸入ADC完成64路ECG的分時采集，減少多通道采集對ADC芯片和FPGA管腳的用量。信號x(n)進入FPGA后，經(jīng)1-64分解器分時傳入64通道的FIR低通濾波、Coif1小波分解和閾值處理模塊，完成64路信號的并行濾波。濾波系數(shù)由64-8多路開關(guān)并行送入8個SPIHT數(shù)據(jù)壓縮模塊進行壓縮。壓縮后由8-1多路開關(guān)經(jīng)通用串行總線(universal serial bus,USB)傳輸模塊送入上位機。

圖1 實現(xiàn)方法框圖

1.1 FPGA中的FIR低通濾波實現(xiàn)

設(shè)計采用31階等波紋FIR低通濾波器去除100 Hz以上噪聲，N階FIR濾波實現(xiàn)的卷積公式為

(1)

式中h(n)為濾波器系數(shù)；x(n)為輸入信號；y(n)為輸出信號；N為濾波器階數(shù)。FIR濾波器在FPGA內(nèi)由分布式算法(distributed algorithm,DA)[14]實現(xiàn)，利用查找表(look-up table,LUT)映射乘積項，代替通用乘法器的使用，廣泛地應(yīng)用在計算乘積和之中，是一項重要的FPGA技術(shù)。2個長度為N的向量c與x的內(nèi)積

(2)

式中y為輸出變量。如果系數(shù)c[n]為已知常數(shù)，x[n]為變量，則B+1位有符號DA系統(tǒng)的x[n]表達式

xb[n]∈[0,1]

(3)

式中xb[n]為x[n]的第b位；x[n]為x的第n次采樣。結(jié)合式(2)和式(3)可得

(4)

即2N字寬、預(yù)先設(shè)定程序LUT接收N位輸入向量xb=[xb[0],xb[1],…,xb[N-1]]后輸出f(c[n],xB[n])或f(c[n],xb[n])。每個輸出均由相應(yīng)的二次冪加權(quán)并累加，在N次查詢循環(huán)后即完成了內(nèi)積y的計算[14]。由于濾波器系數(shù)h(n)為偶對稱常數(shù)，并可用壓棧來實現(xiàn)x(n-i)中的移位i，因此，卷積可由內(nèi)積實現(xiàn)。Arria V FPGA采用8輸入的LUT，因而31階的濾波器僅需4個LUT即可，F(xiàn)IR低通濾波實現(xiàn)如圖2。

圖2 FIR低通濾波在FPGA中的實現(xiàn)方法

1.2 FPGA中小波濾波實現(xiàn)方法

設(shè)計中，Coif小波基的分析濾波器由Mallat多分辨率實現(xiàn)[15]，既可以有效去除ECG噪聲又有助于實現(xiàn)數(shù)據(jù)壓縮[16]。分解的層數(shù)影響了濾波精度和壓縮性能，設(shè)計采用4層Coif1小波分解，并選取各尺度下相應(yīng)的閾值完成軟閾值去噪，以期獲得更好地濾波效果。方法實現(xiàn)分為2個部分：分解緩存和系數(shù)處理。分解緩存部分使用多相雙信道濾波器組[14]來實現(xiàn)小波分解，其中多相濾波器由簡化加法器圖(reduced adder graph,RAG)實現(xiàn)，不但提高了運行速度同時減少了乘法器使用。Coif1小波的高頻分析濾波器H(z)和低頻分析濾波器G(z)的多相實現(xiàn)如式(5)

(5)

式中H0(z2),H1(z2),G0(z2)和G1(z2)為分解后的多相濾波器。低通濾波后，信號經(jīng)過H(z)輸出的第j層細節(jié)系數(shù)和G(z)輸出的第四層概貌系數(shù)緩存至隨機存取存儲器(random access memory,RAM)中，j=1,2,3,4。當緩存達到預(yù)設(shè)值時，依次從第一至第四層RAM中讀出32個小波系數(shù)，d1(k),k=0,1,…,15;d2(k),k=0,1,…,7;d3(k),k=0,1,2,3;d4(k)和a4(k),k=0,1，如圖3。

圖3 分解緩存在FPGA中的實現(xiàn)方法

讀出的系數(shù)依次傳入系數(shù)處理部分完成去噪，軟閾值去噪如式(6)，當|ω|<λ時，令其為0；否則，令其減去閾值，即

(6)

式中ω與ωλ分別為處理前、后的小波系數(shù)；λ為該尺度下的閾值。

在實現(xiàn)時，小波系數(shù)分2路輸入，一路作為待處理系數(shù)直接進入去噪環(huán)節(jié)，另一路則依次經(jīng)過取絕對值、排序和求閾值環(huán)節(jié)來求取各層的閾值。其中，取絕對值環(huán)節(jié)將系數(shù)由補碼轉(zhuǎn)化為無符號源碼。排序環(huán)節(jié)通過對絕對值|dj(k)|排序獲取各層中位數(shù)，32個系數(shù)輸入后，每層中相鄰兩個系數(shù)構(gòu)成一個偶或奇數(shù)比較對，根據(jù)比較對的比較結(jié)果使能系數(shù)間的交換，根據(jù)總比較結(jié)果決定是否完成三層|dj(k)|的排序。流程如下：

1)對于第j層的系數(shù)，j=1∶3有

a.若|dj(2k+1)|>|dj(2k) |，則|dj(2k+1)|?|dj(2k)|，k=0～24-j-1。即交換相鄰系數(shù)；否則,不變。

b.若|dj(2k)|>|dj(2k-1) |，則|dj(2k)|?|dj(2k-1)|，k=1～24-j-1；否則,不變。

2)如果|dj(0)|≤…≤|dj(25-j-1) |，j=1～3，則排序完成；否則，回到步驟(1)進行下一循環(huán)。如，當|d3(0)|～|d3(3)|為3,1,5和6時，本層的排序流程如圖4，每輪對箭頭所指的奇或偶數(shù)比較對比較交換，經(jīng)過5輪后即可完成由小到大排序。

圖4 第三層排序示意

每層的中位數(shù)進入求閾值環(huán)節(jié)后，通過移位加實現(xiàn)與噪聲強度估計值的乘積，求取各尺度下的閾值。去噪環(huán)節(jié)利用各尺度下閾值，按照式(6)對系數(shù)進行處理，實現(xiàn)軟閾值的施加。處理后的高頻系數(shù)dj(k)由Di(k)表示,低頻系數(shù)a4(k)用C4(k)表示，其中，i=1,2,3,4。

1.3 SPIHT壓縮

設(shè)計利用FPGA的并行處理和可編程特性，對文獻[13]中的算法進行了改造，大幅降低了原算法的復(fù)雜度并縮短了比較迭代等過程所用時間，提高了處理的實時性。

1.3.1 SPIHT數(shù)據(jù)壓縮實現(xiàn)

處理后的小波系數(shù)從第一層D1(k)(k=0,1…15)至第四層D4(k),C4(k)(k=0,1)依次進入數(shù)據(jù)壓縮模塊。為滿足層間相關(guān)性假設(shè)，系數(shù)先通過由RAM和地址計數(shù)器構(gòu)成的排序環(huán)節(jié)將輸入順序變?yōu)镃4(k)至D1(k)。為便于編碼，轉(zhuǎn)源碼環(huán)節(jié)將排序后的系數(shù)從補碼轉(zhuǎn)化為源碼。之后的取閾值環(huán)節(jié)則將源碼中模最大值的最高二進制數(shù)篩選出來，作為初始閾值和源碼系數(shù)一起傳入編碼環(huán)節(jié)。編碼環(huán)節(jié)由寄存器、選通開關(guān)、計數(shù)器、比較器和邏輯門等組成，構(gòu)成SPIHT算法的實現(xiàn)主體，如圖5，圖中的黑圓點表示兩個子代總比較結(jié)果的邏輯或。

圖5 編碼環(huán)節(jié)結(jié)構(gòu)

圖5中系數(shù)樹共有4層，每個方框代表一個以系數(shù)命名的系數(shù)單元。選通輸出部分由兩級多路開關(guān)組成，用于選通輸出各系數(shù)單元的總比較結(jié)果、自身比較結(jié)果、自身符號位和子代比較結(jié)果。系數(shù)單元結(jié)構(gòu)如圖6，系數(shù)模塊存儲小波系數(shù)；閾值模塊存儲閾值；比較模塊存儲兩者的比較結(jié)果，大于閾值時自身比較結(jié)果為1，稱系數(shù)有意義；否則，為0，無意義。符號模塊是系數(shù)符號位已輸出標志，1表示已輸出，0未輸出；若系數(shù)為正，自身符號位為0；否則，為1。總比較結(jié)果是子代和自身比較結(jié)果的邏輯或，為1時表示自身與其兩個子代中有大于閾值的系數(shù)；0,則沒有。子代比較結(jié)果是其兩個子代總比較結(jié)果的邏輯或。

圖6 系數(shù)單元結(jié)構(gòu)

編碼環(huán)節(jié)的SPIHT算法實現(xiàn)流程為：

1)將初始閾值T和小波系數(shù)C4(k)至D1(k)存入系數(shù)樹中對應(yīng)的系數(shù)單元，并將比較結(jié)果鎖存至各自的比較模塊，同時存儲T到loop_cnt寄存器中作循環(huán)控制，之后執(zhí)行步驟(2)。

2)判斷l(xiāng)oop_cnt是否為零，如果為零，則發(fā)出信號all_zero=1，表示輸入系數(shù)全為零，執(zhí)行步驟(6)；否則，轉(zhuǎn)到步驟(3)。

3)系數(shù)比較輸出。系數(shù)的比較輸出流程偽代碼如下，不同系數(shù)的流程有所區(qū)別。從第四層至第一層，由左向右依次對系數(shù)單元判斷輸出。

a.if(Di(k)T==1){

y=1;

b.if(Di(k)>=T){

y=1;Di(k)=Di(k)-T;

if(Di(k)S==0){

y=Di(k)Sgn;

Di(k)S=1;

}

c.y=Di(k)B;

d. }elsey=0;

e.}elsey=0;

其中，第i層k個系數(shù)用Di(k)表示；Di(k)T代表該系數(shù)的總比較結(jié)果；Di(k)S代表符號位已輸出標志；Di(k)Sgn代表自身符號位；Di(k)B代表子代比較結(jié)果；T為閾值；y代表輸出比特。

①對系數(shù)D4(k)，k=0,1和D1(k)，k=0～15,比較輸出包括了上述步驟(b)和步驟(d)。

②對系數(shù)D4(k)，k=0,1和D3(k)，k=0～3，包括步驟(a)、步驟(b)、步驟(d)和步驟(e)。

③對系數(shù)D2(k)，k=0～7，包括步驟(a)～步驟(e)。

對Di(2k+1),k=0～24-i-1,i=1～4在比較輸出后還有判斷轉(zhuǎn)移環(huán)節(jié)，如表1，表中第一行的數(shù)字代表判斷的順序，第一列為有此環(huán)節(jié)的系數(shù)。每個數(shù)字有2列，第一列為判斷條件，第二列為當判斷條件為真(1)時下個比較輸出的系數(shù)，若為假(0)則檢查下個判斷條件，以此類推，若一行的判斷條件均為假，則執(zhí)行步驟④。例如對D4(1)，若D4(0)T為1則轉(zhuǎn)移至系數(shù)D3(0)；若為0則檢查D4(1)T，若為1則轉(zhuǎn)移至D3(2)；否則，轉(zhuǎn)到步驟④。

④loop_cnt右移一位，并判斷是否為零。若為0則執(zhí)行步驟⑥;否則,轉(zhuǎn)到步驟⑤。

⑤將每個系數(shù)單元中的閾值右移一位，并鎖存新的比較結(jié)果，回到步驟③。

⑥編碼環(huán)節(jié)結(jié)束。

編碼環(huán)節(jié)輸出的位流經(jīng)緩存整合為32 bits，與初始閾值一起由USB模塊送入上位機。上位機則根據(jù)初始閾值和編碼數(shù)據(jù)進行解壓縮與信號重建。

表1 系數(shù)的判斷轉(zhuǎn)移環(huán)節(jié)

2 實驗結(jié)果

設(shè)計使用MIT-BIH T-wave Alternans Challenge和Motion Artifact Contaminated ECG[17]的數(shù)據(jù)來驗證提出的方法，選取的ECG數(shù)據(jù)采樣率為500 Hz，采樣精度為16 bits。FIR低通濾波器系數(shù)由MATLAB的Filter Design Analysis工具生成，采樣率500 Hz，通帶截止頻率100 Hz的31階等波紋濾波器。圖7為加入白噪聲的twa03m樣本濾波前后的對比，圖7(a)為加入白噪聲的信號，圖7(b)為軟件濾波重建的信號，圖7(c)為設(shè)計中硬件濾波重建后的信號。

圖7 twa03m的噪聲信號和重建信號

表2為twa03m樣本加入不同強度白噪聲后濾波前后的信噪比(signal to noise ratio,SNR)。由圖7和表2可知，本文方法對白噪聲具有較強的抑制能力，在較強噪聲環(huán)境下仍能恢復(fù)信號的整體形態(tài)并保留重要特征；相較于軟件濾波，F(xiàn)PGA濾波的效果同樣顯著，兩者去噪能力相當。濾波后SNR平均提高了6.4 dB最高達7.4 dB，達到了有效抑制噪聲的目的。

表2 濾波前后信噪比 dB

表3為通道9的twa03m信號濾波后，分10段壓縮的壓縮比和壓縮時間，其中，1為原始信號，2為幅值加倍后結(jié)果，每段壓縮32個16 bits系數(shù)，輸出為bit，壓縮時間以FPGA的時鐘周期數(shù)來衡量，設(shè)計使用50 MHz時鐘。表4為10個不同通道的twa03m信號壓縮比與壓縮時間，每個通道輸入320個16 bits系數(shù)，輸出為bit。

表3 壓縮比及壓縮時間(f=50 MHz)

從表3和表4可知，壓縮時間隨著壓縮比的降低而增加；壓縮比及壓縮時間和輸入信號幅值有關(guān)，隨著幅值增加壓縮比有所降低，因為需要更多的比特位來表征信號；壓縮時間也會相應(yīng)增加，由于要進行更多輪的系數(shù)比較輸出，但增加的時間有限，仍然能滿足實時壓縮的需求。信號的局部壓縮比起伏較大，整體壓縮比較高且穩(wěn)定。

3 分析與討論

ECG信號經(jīng)過2次數(shù)字濾波，各頻段噪聲得到了有效地抑制，SNR有了較大提升；得益于流水線寄存器的使用，濾波可在100個時鐘周期內(nèi)完成，相對于采集速度完全可以實現(xiàn)實時濾波。由于FPGA實現(xiàn)時濾波器系數(shù)量化為整形，波形的光滑度有所下降，但并不影響信號的整體形態(tài)和特征參數(shù)的獲取。改造后的SPIHT算法復(fù)雜度和壓縮時間顯著降低，并行性和實時性增強，取得了較高的壓縮比，減輕了多通道ECG的傳輸量和存儲量的負擔；壓縮時間通常在幾百個時鐘周期內(nèi)，為實現(xiàn)多通道ECG的實時處理提供了可能。

表4 twa03m的10個通道壓縮比及壓縮時間(f=50 MHz)

4 結(jié) 論

通過濾波有效地抑制了來自人體自身和外界環(huán)境的噪聲干擾，提高了SNR；SPIHT壓縮則減少了傳輸?shù)臄?shù)據(jù)量，保證了傳輸速率，提高了實時處理的能力。為進一步增強濾波與壓縮性能，可以選取高階的Coiflet小波基[18]或增加分解層數(shù)，若結(jié)合FPGA的硬核處理系統(tǒng)(hard-core processing system,HPS)部分，則可構(gòu)成一個功能更為完善的ECG采集與處理系統(tǒng)。

[1] Bai Yingwen,Chu Wenyang,Chen Chienyu,et al. Adjustable 60 Hz noise reduction by a notch filter for ECG signal[C]∥International and Measurement Technology Conference,Italy:IEEE,2004:1091-5281.

[2] Sarang L Joshi,Rambabu A Vatti,Rupali V Tornekar.A survey on ECG signal denoising techniques[C]∥2013 International Confe-rence on Communication Systems and Network Technologies,Washington DC:IEEE Computer Society,2013:60-64.

[3] Chmelka L,Kozumplik J.Wavelet-based Wiener filter for electrocardiogram signal denoising[C]∥Computers in Cardiology,2005:771-774.

[4] Mithun P,Prem C Pandey,Sebastian T,et al.A wavelet-based technique for suppression of EMG noise and motion artifact in ambulatory ECG[C]∥2011 Annual International Conference of the IEEE Engineering in Medicine and Biology Society,Boston:IEEE,2011:7087-7090.

[5] Andrews C A,Davies J M,Schwarz G R. Adaptive data compression[J].Proceedings of the IEEE,1967,55(3):267-277.

[6] Ahmed S M,Al-Shrouf A,Abo-Zzahhad M.ECG data compression using optimal non-orthogonal wavelet transform[J].Medical Engineering Physics,2000,22(1):39-46.

[7] Miaou S G,Lin C L.A quality-on-demand algorithm for wavelet-based compression of electrocardiogram signals[J].IEEE Trans on Biomed Eng,2002,49(3):233-239.

[8] Al-Shrouf A,Abo-Zahhad M,Ahmed Sabah M.A novel compression algorithm for electrocardiogram signals based on the linear prediction of the wavelet coefficients[J].Digital Signal Process,2003,13(4):604-622.

[9] Ahmed S M.ECG data compression algorithm based on the combination of singular value decomposition and discrete wavelet transform[J].J Eng Sci,2005,33(6):2267-2280.

[10] Sabah M A.ECG signal compression using combined modified discrete-cosine and discrete-wavelet transforms[J].J Eng Sci,2006,34(1):215-226.

[11] Said A,Pearlman W A.A new,fast and efficient image CDEC based on set partitioning in hierarchical trees[J].IEEE Trans on Circuits and Systems for Video Technology,1996,6(3):243-250.

[12] Lu Zhitao,Pearlman W A.An efficient,low-complexity audio coder delivering multiple levels of quality for interactive applications[C]∥Proceedings of 1998 IEEE the Second Workshop on Multimedia Signal Processing,1998:529-534.

[13] Lu Zhitao,Kim Dongyoun,Pearlman W A.Wavelet compression of ECG signals by the set partitioning in hierarchical trees(SPIHT)algorithm[J].IEEE Transactions on Biomedical Engineering,2000,47(7):849-856.

[14] 貝 斯.數(shù)字信號處理的FPGA實現(xiàn)[M].劉 凌,胡永生,譯.北京：清華大學出版社，2006:50-54.

[15] 胡廣書.現(xiàn)代信號處理教程[M].北京：清華大學出版社，2004：246-375.

[16] Stephane Mallat.A wavelet tour of signal processing[M].3 ed.Manhattan:Academic Press,2008:191-444.

[17] PhysioNet.PhysioBank[DB/OL].(2016—10—17),http:∥www.physionet.org.

[18] Daubechies I,Lagarias J C.Tow-scale difference equations:IL.Local regularity,infinite products of matrices and fractals[J].SIAM J of Math Anal,1992,23(4):1031-1079.