基于FPGA的多通道心電信號實時濾波與壓縮*

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

(清華大學 深圳研究生院 深圳市無損監測與微創醫學技術重點實驗室，廣東 深圳 518055)

0 引 言

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

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

1 方 法

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

圖1 實現方法框圖

1.1 FPGA中的FIR低通濾波實現

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

(1)

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

(2)

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

xb[n]∈[0,1]

(3)

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

(4)

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

圖2 FIR低通濾波在FPGA中的實現方法

1.2 FPGA中小波濾波實現方法

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

(5)

式中H0(z2),H1(z2),G0(z2)和G1(z2)為分解后的多相濾波器。低通濾波后，信號經過H(z)輸出的第j層細節系數和G(z)輸出的第四層概貌系數緩存至隨機存取存儲器(random access memory,RAM)中，j=1,2,3,4。當緩存達到預設值時，依次從第一至第四層RAM中讀出32個小波系數，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中的實現方法

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

(6)

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

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

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

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

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)進行下一循環。如，當|d3(0)|～|d3(3)|為3,1,5和6時，本層的排序流程如圖4，每輪對箭頭所指的奇或偶數比較對比較交換，經過5輪后即可完成由小到大排序。

圖4 第三層排序示意

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

1.3 SPIHT壓縮

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

1.3.1 SPIHT數據壓縮實現

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

圖5 編碼環節結構

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

圖6 系數單元結構

編碼環節的SPIHT算法實現流程為：

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

2)判斷loop_cnt是否為零，如果為零，則發出信號all_zero=1，表示輸入系數全為零，執行步驟(6)；否則，轉到步驟(3)。

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

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個系數用Di(k)表示；Di(k)T代表該系數的總比較結果；Di(k)S代表符號位已輸出標志；Di(k)Sgn代表自身符號位；Di(k)B代表子代比較結果；T為閾值；y代表輸出比特。

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

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

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

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

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

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

⑥編碼環節結束。

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

表1 系數的判斷轉移環節

2 實驗結果

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

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

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

表2 濾波前后信噪比 dB

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

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

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

3 分析與討論

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

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

4 結 論

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

[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] 貝 斯.數字信號處理的FPGA實現[M].劉 凌,胡永生,譯.北京：清華大學出版社，2006:50-54.

[15] 胡廣書.現代信號處理教程[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.