摘 要：在SP061A 單片機上實現(xiàn)對ECG信號的FFT,、濾波和壓縮,。合理組織SP061A的硬件資源，并采取數(shù)據(jù)分段長度可選,、避開高頻分量的計算和簡易的數(shù)據(jù)壓縮算法,，使存儲開銷、運算速度和精度滿足實用要求,。
關鍵詞：ECG數(shù)據(jù) SP061A FFT 濾波 壓縮
在遠程心電監(jiān)護系統(tǒng)中,，心電信號采集器是實現(xiàn)心電信號的現(xiàn)場采集、存儲和傳輸?shù)闹匾K端設備,。對采集器的基本要求之一是：及時對采集到的心電信號進行濾波和壓縮等預處理,，以減少存儲器占用量和數(shù)據(jù)遠程傳輸?shù)筋^端服務器的開銷。為降低成本,，這些任務一般采用單片機完成,。然而,，限于單片機的資源、運算能力和運行速度,，許多壓縮算法,，如周期壓縮法、小波變換壓縮法和神經(jīng)網(wǎng)絡方法等無法使用,，一些缺乏快速算法的頻域變換法也很難達到實用的程度[3],。高性價比的心電信號采集器的研制一直是一個熱點問題。
通過研究FFT(快速傅立葉變換)的算法結構和心電信號的特點發(fā)現(xiàn),，采用分段FFT,，保留分析心電波形需要的諧波成分，巧妙地組織單片機的片內RAM資源,，可使數(shù)據(jù)運算量和RAM開銷大大減少,，能實現(xiàn)數(shù)據(jù)濾波和壓縮，且能達到實時采集與處理所需的運算速度,。
SP061A 是凌陽科技公司研制的一款16位超低功耗單片機[1],，片內有2K字RAM、10位A/D轉換器,，CPU時鐘高達49.152MHz,，且價格低廉，還特別具有一套精簡,、高效的指令系統(tǒng)和類似于DSP的硬件內積運算功能,。這些特點很適合心電信號的采集和處理。圖1是作者研發(fā)的心電信號采集器中有關硬件的組成框圖：多路ECG模擬信號送SP061A進行A/D轉換,，轉換數(shù)據(jù)送NVRAM DS1265W暫存,；待采集完成后，由SP061A進行FFT和濾波,、壓縮,；壓縮結果送回DS1265W，再適時通過電話線或計算機網(wǎng)絡送到監(jiān)護中心處理,、診斷,。
圖1 硬件組成框圖
本文僅討論用SP061A實現(xiàn)FFT[2]、低通濾波與壓縮,。設對心電信號的采樣率為500次/秒,，數(shù)據(jù)精度為10位。
1 數(shù)據(jù)分段算法
設采集到的原始數(shù)據(jù)存于片外RAM中,，將這些數(shù)據(jù)分為若干段,，逐段讀入片內進行FFT。各段的變換結果及時送回片外RAM中保存,。
按照FFT的要求,，段中包含的數(shù)據(jù)個數(shù)必須為2N,，N為FFT變換的層數(shù)?？紤]到SP061A片內RAM為2K字,，此處取N=9或N=10，即段中數(shù)據(jù)為512或1024,，以保證RAM夠用,。顯然，段頭和段尾的數(shù)據(jù)大小相等時,，以該段作為一個周期而無限重復的波形將無跳躍點,。經(jīng)過“FFT變換到頻域”→“丟棄高頻成分”→“IFFT(快速傅立葉反變換，在頭端PC上進行)”一系列操作而重建的時域波形,，段與段之間的結合點將是連續(xù)的,。但實際上，按上述分段幾乎不能做到段頭和段尾的數(shù)據(jù)大小相等,。取兩種段長的目的就是提供兩種可能的選擇——選擇首尾數(shù)據(jù)之差較小的段作FFT,。盡管如此，段首尾數(shù)據(jù)之差仍存在,，經(jīng)處理、復原后的波形在段的結合部位仍將有間斷點,。而采用加窗,、延拓等辦法在單片機上又難以實現(xiàn)。解決問題的策略為：分段時,，各段間的數(shù)據(jù)首,、尾各覆蓋10個數(shù)據(jù)。頭端PC在完成重建后,，應將首,、尾各5個數(shù)據(jù)丟棄。
2 時域數(shù)據(jù)的整序與加載
分段后,，將該段加載到SP061A的RAM中,，以實施FFT。原始數(shù)據(jù)以采集的時間先后順序存放,，加載時則應“整序”,，即改變數(shù)據(jù)的先后順序，以保證變換后的頻域數(shù)據(jù)為正序,。
設RS為指向片外RAM的,、待加載的段內數(shù)據(jù)的偏移地址，RS=0…2N-1,；Rd為指向片內RAM的,、待寫入數(shù)據(jù)的偏移地址,，如圖2。將RS按N位二進制逐位高低互換就得到Rd,。例如,，當N=9時，若RS為011001011B,，則Rd為110100110B,。為加快計算速度，將N=9時Rd的值制表存于FLASH ROM,，供整序時查詢,。當N=10時，取RS的B0～B9位查表獲得Rd,，再將RS的B10位傳送到Rd的B15位,，最后將Rd循環(huán)左移1位。
圖2 整序與數(shù)據(jù)加載操作
FFT變換是復數(shù)運算,。在將原始數(shù)據(jù)加載到片內RAM的同時,，應把實數(shù)轉換為復數(shù)，即令虛部為0,。于是,，一個原始數(shù)據(jù)加載到RAM中要占用2個字。復數(shù)的存儲格式為：實部字存于低地址,，虛部字存于相鄰的高地址?，F(xiàn)在考察RAM需要量。N=9時,，段長為512個數(shù)據(jù),，加載到RAM中要占用 512×2=1024字；N=10時,，段長為1024個數(shù)據(jù),，全部加載將占用 1024×2=2048字，超過片內RAM的可用容量,。此時,，將數(shù)據(jù)分為兩部分，先將第一部分加載到RAM作FFT,，得到中間結果,，再將第二部分加載、變換,，最后相加合成,。
3 FFT變換及低通濾波
FFT將時域序列{x[i]，i∈0…2N}變換為頻域序列{F[i],，i∈0…2N},。為了實現(xiàn)低通濾波,，僅須保留{F[i]}中≤75Hz的頻率分量。當N=9時,，應保留{F[i]}中的前77個低頻分量,；當N=10時，則應保留{F[i]}中的前154個低頻分量,。這也同時減少了計算量,，加快了計算速度；存放周轉量所需的片內RAM也能得到保證,。為敘述簡便,，以N=3為例，研究FFT的計算結構,，如圖3所示,。
圖3 8點FFT的蝶形計算結構
圖3中，W[k]是復因子,，k=0…2N-1,。將W[k]的實部和虛部都乘214，取整后制成表,，存于FLASH ROM中,，供程序查表獲得其值；而W[k]與某數(shù)相乘,，將32位運算結果右移14位作為積,。這就使全部運算為整數(shù)運算，適應SP061A的硬件乘法功能,。由圖3知,，第一層的計算僅涉及實部加減,，虛部保持為0,，可單獨進行。從第二層開始有復數(shù)乘,，但是,，當只需計算{F[i]}中的低頻分量時，許多中間結果可不計算,。例如,，如果需計算出F[0]和F[1](即保留原始信號的直流分量和1次諧波)，則僅需計算x[0]3,、x[4]3和x[1]3,、x[5]3。計算層數(shù)N越多,，減少的運算也越多,。
復數(shù)乘可利用SP061A的內積功能實現(xiàn),。例如，要計算x[i]×W[j],，設x[i]×W[j]=(a+jb)×(c+jd)=ac+(-bd)+j(bc+ad),。顯然，結果的實部和虛部均為內積形式,，只是設置操作數(shù)時須注意符號和排列順序,。
上述方法使計算量顯著減少。以512點FFT為例,，計算出全部頻率分量需要512×log2512=4608次運算,，其中含有2048次復數(shù)乘。若計算77個低頻分量,，則只有3611次運算,，其中含有1767次復數(shù)乘。
當N=10時,，計算點數(shù)達1024,，片內RAM不夠用。此時,，應按1024點的整序次序取數(shù),，先對x[0]1～x[511]1進行FFT，算出F1[0]～F1[153],，暫存于片內RAM中的一個緩沖區(qū),；再對x[512]1～x[1023]1進行FFT，算出F2[0]～F2[153],；則最終結果為：F[i]=F1[i]+F2[i],，i=0…153。
為避免計算中產生數(shù)據(jù)溢出,，從第三層開始,，對x[i]4～x[i]9都算術右移1位。操作的累積結果使F[i]縮小了64倍,，故在重建時應擴大64倍,。如此操作實際上降低了運算精度，但實驗表明,，重建的波形完全滿足醫(yī)學觀察要求,。
4 數(shù)據(jù)壓縮
采取如下簡易格式實現(xiàn)數(shù)據(jù)壓縮：
對于F[0]，因虛部為0,，僅用一個字存放實部,，重建時默認虛部為0；
對于F[i]，i>0,，若實部在-64～63范圍內且虛部在-128～127范圍內,，則用2個字節(jié)存放，格式如下：
否則,，用3個字節(jié)存放,，格式如下：
兩種格式由第1字節(jié)的最高位區(qū)分。
5 實驗結果與分析
用自行研發(fā)的心電信號采集器進行實驗,，對采集到的4個樣本進行處理,，實驗結果如表1。表1中,，PRD為均方根誤差,，CC為相關系數(shù)，計算公式為：
表1 實驗結果數(shù)據(jù)
式中m為樣本系列的數(shù)據(jù)個數(shù),，x[n],、分別為原始數(shù)據(jù)系列及其平均值，y[n],、分別為重建后的系列及其平均值,。CR為壓縮比，,。
處理時間為SP061A完成FFT與壓縮花費的時間,， CPU時鐘設置為49.152MHz。
實驗表明,，本方法用價格低廉的單片機實現(xiàn)了復雜的FFT與數(shù)據(jù)壓縮,，計算耗時少，所得結果滿足實用要求,。由圖4可見,，重建后的波形在段間結合點無畸變。噪聲較弱時PRD和CC參數(shù)較為理想,；而當噪聲很強時,，如圖4(a)、4(b),，因濾除了高頻噪聲而使得重建波形與原始波形差距較大,，PRD和CC參數(shù)已不能說明問題,。壓縮算法簡便,，CR約為4。順便說明,，本方法未實現(xiàn)50Hz干擾濾波,、肌電干擾濾波和基線漂移，這些處理可在頭端PC上進行,。
圖4 實驗結果
