基于SP061A實(shí)現(xiàn)心電數(shù)據(jù)的FFT與壓縮

在SP061A 單片機(jī)上實(shí)現(xiàn)對(duì)ECG信號(hào)的FFT、濾波和壓縮。合理組織SP061A的硬件資源，并采取數(shù)據(jù)分段長(zhǎng)度可選、避開高頻分量的計(jì)算和簡(jiǎn)易的數(shù)據(jù)壓縮算法，使存儲(chǔ)開銷、運(yùn)算速度和精度滿足實(shí)用要求。

在遠(yuǎn)程心電監(jiān)護(hù)系統(tǒng)中，心電信號(hào)采集器是實(shí)現(xiàn)心電信號(hào)的現(xiàn)場(chǎng)采集、存儲(chǔ)和傳輸?shù)闹匾K端設(shè)備。對(duì)采集器的基本要求之一是：及時(shí)對(duì)采集到的心電信號(hào)進(jìn)行濾波和壓縮等預(yù)處理，以減少存儲(chǔ)器占用量和數(shù)據(jù)遠(yuǎn)程傳輸?shù)筋^端服務(wù)器的開銷。為降低成本，這些任務(wù)一般采用單片機(jī)完成。然而，限于單片機(jī)的資源、運(yùn)算能力和運(yùn)行速度，許多壓縮算法，如周期壓縮法、小波變換壓縮法和神經(jīng)網(wǎng)絡(luò)方法等無法使用，一些缺乏快速算法的頻域變換法也很難達(dá)到實(shí)用的程度[3]。高性價(jià)比的心電信號(hào)采集器的研制一直是一個(gè)熱點(diǎn)問題。

通過研究FFT（快速傅立葉變換）的算法結(jié)構(gòu)和心電信號(hào)的特點(diǎn)發(fā)現(xiàn)，采用分段FFT，保留分析心電波形需要的諧波成分，巧妙地組織單片機(jī)的片內(nèi)RAM資源，可使數(shù)據(jù)運(yùn)算量和RAM開銷大大減少，能實(shí)現(xiàn)數(shù)據(jù)濾波和壓縮，且能達(dá)到實(shí)時(shí)采集與處理所需的運(yùn)算速度。

SP061A 是凌陽科技公司研制的一款16位超低功耗單片機(jī)[1]，片內(nèi)有2K字RAM、10位A/D轉(zhuǎn)換器，CPU時(shí)鐘高達(dá)49.152MHz，且價(jià)格低廉，還特別具有一套精簡(jiǎn)、高效的指令系統(tǒng)和類似于DSP的硬件內(nèi)積運(yùn)算功能。這些特點(diǎn)很適合心電信號(hào)的采集和處理。圖1是作者研發(fā)的心電信號(hào)采集器中有關(guān)硬件的組成框圖：多路ECG模擬信號(hào)送SP061A進(jìn)行A/D轉(zhuǎn)換，轉(zhuǎn)換數(shù)據(jù)送NVRAM DS1265W暫存；待采集完成后，由SP061A進(jìn)行FFT和濾波、壓縮；壓縮結(jié)果送回DS1265W，再適時(shí)通過電話線或計(jì)算機(jī)網(wǎng)絡(luò)送到監(jiān)護(hù)中心處理、診斷。

本文僅討論用SP061A實(shí)現(xiàn)FFT[2]、低通濾波與壓縮。設(shè)對(duì)心電信號(hào)的采樣率為500次／秒，數(shù)據(jù)精度為10位。

1 數(shù)據(jù)分段算法

設(shè)采集到的原始數(shù)據(jù)存于片外RAM中，將這些數(shù)據(jù)分為若干段，逐段讀入片內(nèi)進(jìn)行FFT。各段的變換結(jié)果及時(shí)送回片外RAM中保存。

按照FFT的要求，段中包含的數(shù)據(jù)個(gè)數(shù)必須為2N，N為FFT變換的層數(shù)?？紤]到SP061A片內(nèi)RAM為2K字，此處取N=9或N=10，即段中數(shù)據(jù)為512或1024，以保證RAM夠用。顯然，段頭和段尾的數(shù)據(jù)大小相等時(shí)，以該段作為一個(gè)周期而無限重復(fù)的波形將無跳躍點(diǎn)。經(jīng)過“FFT變換到頻域”→“丟棄高頻成分”→“IFFT（快速傅立葉反變換，在頭端PC上進(jìn)行）”一系列操作而重建的時(shí)域波形，段與段之間的結(jié)合點(diǎn)將是連續(xù)的。但實(shí)際上，按上述分段幾乎不能做到段頭和段尾的數(shù)據(jù)大小相等。取兩種段長(zhǎng)的目的就是提供兩種可能的選擇——選擇首尾數(shù)據(jù)之差較小的段作FFT。盡管如此，段首尾數(shù)據(jù)之差仍存在，經(jīng)處理、復(fù)原后的波形在段的結(jié)合部位仍將有間斷點(diǎn)。而采用加窗、延拓等辦法在單片機(jī)上又難以實(shí)現(xiàn)。解決問題的策略為：分段時(shí)，各段間的數(shù)據(jù)首、尾各覆蓋10個(gè)數(shù)據(jù)。頭端PC在完成重建后，應(yīng)將首、尾各5個(gè)數(shù)據(jù)丟棄。

2 時(shí)域數(shù)據(jù)的整序與加載

分段后，將該段加載到SP061A的RAM中，以實(shí)施FFT。原始數(shù)據(jù)以采集的時(shí)間先后順序存放，加載時(shí)則應(yīng)“整序”，即改變數(shù)據(jù)的先后順序，以保證變換后的頻域數(shù)據(jù)為正序。

設(shè)RS為指向片外RAM的、待加載的段內(nèi)數(shù)據(jù)的偏移地址，RS=0…2N-1；Rd為指向片內(nèi)RAM的、待寫入數(shù)據(jù)的偏移地址，如圖2。將RS按N位二進(jìn)制逐位高低互換就得到Rd。例如，當(dāng)N=9時(shí)，若RS為011001011B，則Rd為110100110B。為加快計(jì)算速度，將N=9時(shí)Rd的值制表存于FLASH ROM，供整序時(shí)查詢。當(dāng)N=10時(shí)，取RS的B0～B9位查表獲得Rd，再將RS的B10位傳送到Rd的B15位，最后將Rd循環(huán)左移1位。

FFT變換是復(fù)數(shù)運(yùn)算。在將原始數(shù)據(jù)加載到片內(nèi)RAM的同時(shí)，應(yīng)把實(shí)數(shù)轉(zhuǎn)換為復(fù)數(shù)，即令虛部為0。于是，一個(gè)原始數(shù)據(jù)加載到RAM中要占用2個(gè)字。復(fù)數(shù)的存儲(chǔ)格式為：實(shí)部字存于低地址，虛部字存于相鄰的高地址?，F(xiàn)在考察RAM需要量。N=9時(shí)，段長(zhǎng)為512個(gè)數(shù)據(jù)，加載到RAM中要占用 512×2=1024字；N=10時(shí)，段長(zhǎng)為1024個(gè)數(shù)據(jù)，全部加載將占用 1024×2=2048字，超過片內(nèi)RAM的可用容量。此時(shí)，將數(shù)據(jù)分為兩部分，先將第一部分加載到RAM作FFT，得到中間結(jié)果，再將第二部分加載、變換，最后相加合成。

3 FFT變換及低通濾波

FFT將時(shí)域序列{x[i],i∈0…2N}變換為頻域序列{F[i],i∈0…2N}。為了實(shí)現(xiàn)低通濾波，僅須保留{F[i]}中≤75Hz的頻率分量。當(dāng)N=9時(shí)，應(yīng)保留{F[i]}中的前77個(gè)低頻分量；當(dāng)N=10時(shí)，則應(yīng)保留{F[i]}中的前154個(gè)低頻分量。這也同時(shí)減少了計(jì)算量，加快了計(jì)算速度；存放周轉(zhuǎn)量所需的片內(nèi)RAM也能得到保證。

為敘述簡(jiǎn)便，以N=3為例，研究FFT的計(jì)算結(jié)構(gòu)，如圖3所示。

kN=0…2N-1。將W[k]的實(shí)部和虛部都乘214，取整后制成表，存于FLASH ROM中，供程序查表獲得其值；而W[k]與某數(shù)相乘，將32位運(yùn)算結(jié)果右移14位作為積。這就使全部運(yùn)算為整數(shù)運(yùn)算，適應(yīng)SP061A的硬件乘法功能。由圖3知，第一層的計(jì)算僅涉及實(shí)部加減，虛部保持為0，可單獨(dú)進(jìn)行。從第二層開始有復(fù)數(shù)乘，但是，當(dāng)只需計(jì)算{F[i]}中的低頻分量時(shí)，許多中間結(jié)果可不計(jì)算。例如，如果需計(jì)算出F[0]和F[1]（即保留原始信號(hào)的直流分量和1次諧波），則僅需計(jì)算x[0]3、x[4]3和x[1]3、x[5]3。計(jì)算層數(shù)N越多，減少的運(yùn)算也越多。

復(fù)數(shù)乘可利用SP061A的內(nèi)積功能實(shí)現(xiàn)。例如，要計(jì)算x[i]×W[j]，設(shè)x[i]×W[j]=(a+jb)×(c+jd)=ac+(-bd)+j(bc+ad)。顯然，結(jié)果的實(shí)部和虛部均為內(nèi)積形式，只是設(shè)置操作數(shù)時(shí)須注意符號(hào)和排列順序。

上述方法使計(jì)算量顯著減少。以512點(diǎn)FFT為例，計(jì)算出全部頻率分量需要512×log2512=4608次運(yùn)算，其中含有2048次復(fù)數(shù)乘。若計(jì)算77個(gè)低頻分量，則只有3611次運(yùn)算，其中含有1767次復(fù)數(shù)乘。

當(dāng)N=10時(shí)，計(jì)算點(diǎn)數(shù)達(dá)1024，片內(nèi)RAM不夠用。此時(shí)，應(yīng)按1024點(diǎn)的整序次序取數(shù)，先對(duì)x[0]1～x[511]1進(jìn)行FFT，算出F1[0]～F1[153]，暫存于片內(nèi)RAM中的一個(gè)緩沖區(qū)；再對(duì)x[512]1～x[1023]1進(jìn)行FFT，算出F2[0]～F2[153]；則最終結(jié)果為:F[i]=F1[i]+F2[i]，i=0…153。

為避免計(jì)算中產(chǎn)生數(shù)據(jù)溢出，從第三層開始，對(duì)x[i]4～x[i]9都算術(shù)右移1位。操作的累積結(jié)果使F[i]縮小了64倍，故在重建時(shí)應(yīng)擴(kuò)大64倍。如此操作實(shí)際上降低了運(yùn)算精度，但實(shí)驗(yàn)表明，重建的波形完全滿足醫(yī)學(xué)觀察要求。
4 數(shù)據(jù)壓縮

采取如下簡(jiǎn)易格式實(shí)現(xiàn)數(shù)據(jù)壓縮：

對(duì)于F[0]，因虛部為0，僅用一個(gè)字存放實(shí)部，重建時(shí)默認(rèn)虛部為0；

對(duì)于F[i],i>0，若實(shí)部在-64～63范圍內(nèi)且虛部在-128～127范圍內(nèi)，則用2個(gè)字節(jié)存放，格式如下：

兩種格式由第1字節(jié)的最高位區(qū)分。

5 實(shí)驗(yàn)結(jié)果與分析

用自行研發(fā)的心電信號(hào)采集器進(jìn)行實(shí)驗(yàn)，對(duì)采集到的4個(gè)樣本進(jìn)行處理，實(shí)驗(yàn)結(jié)果如表1。表1中，PRD為均方根誤差，CC為相關(guān)系數(shù)，計(jì)算公式為：

數(shù)據(jù)系列及其平均值，平均值。

處理時(shí)間為SP061A完成FFT與壓縮花費(fèi)的時(shí)間, CPU時(shí)鐘設(shè)置為49.152MHz。

實(shí)驗(yàn)表明，本方法用價(jià)格低廉的單片機(jī)實(shí)現(xiàn)了復(fù)雜的FFT與數(shù)據(jù)壓縮，計(jì)算耗時(shí)少，所得結(jié)果滿足實(shí)用要求。由圖4可見，重建后的波形在段間結(jié)合點(diǎn)無畸變。噪聲較弱時(shí)PRD和CC參數(shù)較為理想；而當(dāng)噪聲很強(qiáng)時(shí)，如圖4（a）、4（b），因?yàn)V除了高頻噪聲而使得重建波形與原始波形差距較大，PRD和CC參數(shù)已不能說明問題。壓縮算法簡(jiǎn)便，CR約為4。順便說明，本方法未實(shí)現(xiàn)50Hz干擾濾波、肌電干擾濾波和基線漂移，這些處理可在頭端PC上進(jìn)行。