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