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