基於matlab的心電訊號預先處理

標籤：自適應   size   大致   data   ace   axis   未來   plain   影響   

這是前段時間做的一個課程設計，做的比較簡單，沒有考慮到太細，只是初步地達到了想要的效果。這次設計主要是對心電訊號進行預先處理，將其訊號中包含的一些幹擾濾除或者抑制掉。

一、心電訊號

（1）心電訊號的特性

人體心電訊號是非常微弱的生理低頻電訊號，通常最大的幅值不超過5mV，訊號頻率在0．05～100Hz之間。心電訊號是通過安裝在人體皮膚表面的電極來拾取的。由於電極和皮膚組織之間會發生極化現象，會對心電訊號產生嚴重的幹擾。加之人體是一個複雜的生命系統，存在各種各樣的其他生理電訊號對心電訊號產生幹擾。同時由於我們處在一個電磁包圍的環境中，人體就像一根會移動的天線，從而會對心電訊號產生50Hz左右的幹擾訊號。心電訊號具有微弱、低頻、高阻抗等特性，極容易受到幹擾，所以分析幹擾的來源，針對不同幹擾採取相應的濾除措施，是資料擷取重點考慮的一個問題。常見幹擾有如下幾種：

①工頻幹擾②基準漂移③肌電幹擾 

心電訊號具有以下幾個特點：    

  ·訊號極其微弱，一般只有0.05～4mV，典型值為1mV；
  ·頻率範圍較低，頻率範圍為0.1～35Hz，主要集中在5～20Hz；
  ·存在不穩定性。人體內部各器官問的相互影響以及各人的心臟位置、呼吸、年齡、是否經常鍛煉等因素，都會使心電訊號發生相應變化；
  ·幹擾雜訊很強。對心電訊號進行測量時，必然要與外界聯絡，但由於其自身的訊號非常微弱，因此，各種幹擾雜訊非常容易影響測量。
其雜訊可能來自工頻(50Hz)幹擾、電極接觸噪點、運動偽跡、肌電雜訊、呼吸引起的基準漂移和心電幅度變化以及其他電子裝置的機器雜訊等諸多方面。

（2）心電訊號的選擇
本次實驗所採用的心電訊號來自MIT-BIH庫（心律失常資料庫，關於這個庫的介紹可以參考：點擊開啟連結），庫中有48組失常的心電訊號，要在其中找出符合實驗要求的心電訊號（即含有肌電幹擾、工頻幹擾和基準漂移）。
（3）正常心電訊號波形

 

                                           圖1   正常心電訊號波形

圖1是正常心電訊號在一個周期內的波形，由P波、QRS波群和T波組成。

P波是由心房的去極化產生的，其波形比較小，形狀有些圓，幅度約為0.25mV，期間為0.08~0.11s。竇房結去極化發生在心房肌細胞去極化之前，因而在時間上要先於P波，只是竇房結處於心臟內部，其電活動在體表難以採集。

P-R間期是指P波起點和QRS波群起點所跨越的時間，是竇房結產生的興奮，經過右心房、左心房、房室交接區、房室束、左右束支之後，傳到到心室所需要的時間。在正常的體表心電圖中，P-R間期的值為0.12~0.2s，其中大部分時間是興奮在房室交界區內傳導所需要的時間。P-R間期也稱為房室傳導時間。

P-R段是指P波終點和QRS波群起點之間所跨越的時間。在正常的體表心電圖中，P-R段的心電訊號電位值都是接近基準水平的很小點位。在P-R段期間，左右心房同時興奮，因而兩者產生的綜合電場對體表心電圖的影響較小。另外，此時的興奮還處於房室交界區和房室束特殊傳導系統中，沒有到達心室，因而沒有產生較大波動的體表心電圖訊號。

QRS波群是左右心室肌細胞一次發生去極化所產生的膜外負電位在體表的反應。QRS波群的期間為0.06~0.1s。由於心室肌細胞在興奮過程中的綜合電場向量多次發生改變，因而形成了體表心電圖中大小和方向多次發生變化的心電訊號，其中QRS波群中第一個向下的波為Q波，第一個向上的波為R波，R波後面的為S波。

S-T段是指QRS波群終點和T波起點之間所跨越的時間。S-T段期間，左右心室的肌細胞都處於興奮期間，因而兩者形成的綜合電場向量在體表心電圖中的貢獻非常小，導致S-T段心電訊號處於大約基準的水平。

T波由心室肌細胞的複極化產生，其幅度為0.1~0.8mV，期間為0.05~0.25s。由於複極化差異的存在，T波的方向和QRS波群主波的方向一致。在R波向上的情況下，T波的幅度一般都超過R波幅度的1/10。Q-T間期是指QRS波群起點和T波終點所跨越的時間段，代表心室肌細胞開始去極化到結束複極化所需要的時間，與心率呈負相關。

二、濾波器的選擇

1.肌電幹擾的濾除—低通濾波器  

通常來說，肌電訊號的頻率為20~5000HZ，其主要成分的頻率與肌肉的類型有關，一般在30~300HZ，而心電訊號的頻率主要集中在5~20HZ，所以選擇低通濾波器來濾除肌電幹擾。

巴特沃斯濾波器的特點是通頻帶內的頻率響應曲線最為平坦，沒有起伏，而在阻頻帶則逐漸下降為零。巴特沃斯濾波器的振幅對角頻率單調下降，並且濾波器的階數越高，在阻頻帶幅度衰減速度越快，其他濾波器高階的振幅對角頻率圖和低階數的振幅對角頻率有不同的形狀。

2.工頻幹擾的抑制—帶陷濾波器

工頻幹由於供電網路無所不在，因此50Hz的工頻幹擾是最普遍的，也是心電訊號的主要幹擾來源。50HZ陷波器的軟體設計方法多種多樣，常見方法有小波變換濾波、自適應濾波、模板匹配濾波等，但都需要手工計算獲得濾波器的參數，運算比較複雜。
濾波器設計中，使用IIR濾波器，可使階數降低，運算量減少，但破壞了相位特性；使用FIR濾波器既能得到很好的濾波效果，是波形失真達到最下，而且，FIR濾波器可以做成線性相位特性，這正好是心電訊號濾波所需要的。
利用MATLAB設計FIR濾波器的方法有窗函數法、頻率抽樣法和切比雪夫逼近法等，本次課設採用窗函數法設計50HZ陷波濾波器。窗函數方法的基本思想是：首先根據要求選擇一個適當的理想低通濾波器，因為其脈衝響應是非因果且無限長的，用最佳化窗結構窗函數來截取它的脈衝響應，從而得到線性相位和因果的FIR濾波器。Kaiser窗是接近最佳化窗結構的窗函數，它可以根據不同的參數調整濾波器的各項指標，因此採用Kaiser窗函數進行濾波器設計擾的抑制—帶陷濾波器

3.基準漂移的糾正—零相移濾波器

零相移濾波器是指一個訊號序列經過該濾波器濾波後相位不發生變化，即該濾波器系統函數的相位響應為零。顯然，對於因果系統來說是不可能實現零相移的，在事先無法知道訊號相位譜的情況下，實現零相移是不可能的。零相移只能是對非因果系統來說的。具體而言，零相移濾波器使用了當前訊號點前面和後面的訊號點所包含的資訊，從本質上說就是使用了“未來的資訊”來消除相位失真。

三、程式及結果

1.心電訊號讀取

        因為對MIT-BIH庫不是很熟悉，在官網上看過之後，還是不懂（全英文，而且是醫學方面的。。。）。所以，此處的心電訊號的讀取程式是來自網上的  rddata.m ，詳細地可以參考點擊開啟連結。如果自己要用的話，在選取好要處理的心電訊號後，把路徑更改，並選取合適的樣本數，就可以了。我選取的是MIT-BIH中的 109，樣本數為1500，為心電訊號讀取後的圖形：

                                        

                                                                圖2  MIT-BIH庫109心電訊號圖（雙通道）

        從圖2紅色曲線可以看到，波形上存在許多“毛刺”，並且其相位在發生變化（以波峰為例，各波峰大致不在一條水平線上，即所說的“基準漂移”），部分波形收到的幹擾比較嚴重，比較符合對訊號處理的要求。

2.心電訊號的預先處理

（1）肌電訊號的濾除

clc;%------------------------------低通濾波器濾除肌電訊號------------------------------Fs=1500;                        %採樣頻率fp=80;fs=100;                    %通帶截止頻率，阻帶截止頻率rp=1.4;rs=1.6;                    %通帶、阻帶衰減wp=2*pi*fp;ws=2*pi*fs;   [n,wn]=buttord(wp,ws,rp,rs,‘s‘);     %‘s‘是確定巴特沃斯類比濾波器階次和3dB                               截止類比頻率[z,P,k]=buttap(n);   %設計歸一化巴特沃斯類比低通濾波器，z為極點，p為零點和k為增益[bp,ap]=zp2tf(z,P,k)  %轉換為Ha(p),bp為分子係數，ap為分母係數[bs,as]=lp2lp(bp,ap,wp) %Ha(p)轉換為低通Ha(s)並去歸一化，bs為分子係數，as為分母係數[hs,ws]=freqs(bs,as);         %類比濾波器的幅頻響應[bz,az]=bilinear(bs,as,Fs);     %對類比濾波器雙線性變換[h1,w1]=freqz(bz,az);         %數字濾波器的幅頻響應m=filter(bz,az,M(:,1));figurefreqz(bz,az);title(‘巴特沃斯低通濾波器幅頻曲線‘);      figuresubplot(2,1,1);plot(TIME,M(:,1));xlabel(‘t(s)‘);ylabel(‘mv‘);title(‘原始心電訊號波形‘);grid;subplot(2,1,2);plot(TIME,m);xlabel(‘t(s)‘);ylabel(‘mv‘);title(‘低通濾波後的時域圖形‘);grid;   N=512n=0:N-1;mf=fft(M(:,1),N);               %進行頻譜變換（傅裡葉變換）mag=abs(mf);f=(0:length(mf)-1)*Fs/length(mf);  %進行頻率變換figuresubplot(2,1,1)plot(f,mag);axis([0,1500,1,50]);grid;      %畫出頻譜圖xlabel(‘頻率(HZ)‘);ylabel(‘幅值‘);title(‘心電訊號頻譜圖‘);mfa=fft(m,N);                    %進行頻譜變換（傅裡葉變換）maga=abs(mfa);fa=(0:length(mfa)-1)*Fs/length(mfa);  %進行頻率變換subplot(2,1,2)plot(fa,maga);axis([0,1500,1,50]);grid;  %畫出頻譜圖xlabel(‘頻率(HZ)‘);ylabel(‘幅值‘);title(‘低通濾波後心電訊號頻譜圖‘);    wn=M(:,1);P=10*log10(abs(fft(wn).^2)/N);f=(0:length(P)-1)/length(P);figureplot(f,P);gridxlabel(‘歸一化頻率‘);ylabel(‘功率(dB)‘);title(‘心電訊號的功率譜‘);

以上程式的結果如下：

                            

                                                          圖3                                                                               圖4

                             

                                                  圖5                                                                                            圖6

圖3是所設計的巴特沃斯數字低通濾波器的幅頻響應曲線，圖3是在時域濾波前後心電訊號的波形圖，圖5是在頻域濾波前後心電訊號的頻譜圖，圖6是心電訊號的功率譜圖

（2）工頻幹擾的抑制

%-----------------帶陷濾波器抑制工頻幹擾-------------------%50Hz陷波器：由一個低通濾波器加上一個高通濾波器組成%而高通濾波器由一個全通濾波器減去一個低通濾波器構成Me=100;               %濾波器階數L=100;                %視窗長度beta=100;             %衰減係數Fs=1500;wc1=49/Fs*pi;     %wc1為高通濾波器截止頻率，對應51Hzwc2=51/Fs*pi     ;%wc2為低通濾波器截止頻率，對應49Hzh=ideal_lp(0.132*pi,Me)-ideal_lp(wc1,Me)+ideal_lp(wc2,Me); %h為陷波器                                                            衝擊響應w=kaiser(L,beta);y=h.*rot90(w);         %y為50Hz陷波器衝擊響應序列m2=filter(y,1,m);figuresubplot(2,1,1);plot(abs(h));axis([0 100 0 0.2]);xlabel(‘頻率(Hz)‘);ylabel(‘幅度(mv)‘);title(‘陷波器幅度譜‘);grid;N=512;P=10*log10(abs(fft(y).^2)/N);f=(0:length(P)-1);subplot(2,1,2);plot(f,P);xlabel(‘頻率(Hz)‘);ylabel(‘功率(dB)‘);title(‘陷波器功率譜‘);grid;   figuresubplot (2,1,1); plot(TIME,m);xlabel(‘t(s)‘);ylabel(‘幅值‘);title(‘原始訊號‘);grid;subplot(2,1,2);plot(TIME,m2);xlabel(‘t(s)‘);ylabel(‘幅值‘);title(‘帶阻濾波後訊號‘);grid;  figureN=512subplot(2,1,1);plot(abs(fft(m))*2/N);axis([0 100 0 1]);xlabel(‘t(s)‘);ylabel(‘幅值‘);title(‘原始訊號頻譜‘);grid;subplot(2,1,2);plot(abs(fft(m2))*2/N);axis([0 100 0 1]);xlabel(‘t(s)‘);ylabel(‘幅值‘);title(‘帶阻濾波後訊號頻譜‘);grid;  其中，ideal_lp（）函數在另一個M檔案中，具體如下：%理想低通濾波器%截止角頻率wc，階數Mefunction hd=ideal_lp(wc,Me)alpha=(Me-1)/2;n=[0:Me-1];p=n-alpha+eps;              %eps為很小的數，避免被0除hd=sin(wc*p)./(pi*p);       %用Sin函數產生衝擊響應

以上程式的結果如下：

                         

                                                        圖7                                                                                                  圖8

                                                                                      

                                                                                                                   圖9

圖7是帶陷濾波器的幅度譜和功率譜，可以看到在50Hz處，濾波器的幅度很大，而且功率在-150以下，說明帶陷效能較好。圖8是在時域濾波前後的心電訊號圖，可以看出，濾波後波形有了略微的改善。圖19是在頻域濾波前後的心電訊號頻譜圖。

（3）基準漂移的糾正

%------------------IIR零相移數字濾波器糾正基準漂移-------------------Wp=1.4*2/Fs;     %通帶截止頻率 Ws=0.6*2/Fs;     %阻帶截止頻率 devel=0.005;    %通帶紋波 Rp=20*log10((1+devel)/(1-devel));   %通帶紋波係數  Rs=20;                          %阻帶衰減 [N Wn]=ellipord(Wp,Ws,Rp,Rs,‘s‘);   %求橢圓濾波器的階次 [b a]=ellip(N,Rp,Rs,Wn,‘high‘);       %求橢圓濾波器的係數 [hw,w]=freqz(b,a,512);   result =filter(b,a,m2); figurefreqz(b,a);figuresubplot(211); plot(TIME,m2); xlabel(‘t(s)‘);ylabel(‘幅值‘);title(‘原始訊號‘);gridsubplot(212); plot(TIME,result); xlabel(‘t(s)‘);ylabel(‘幅值‘);title(‘線性濾波後訊號‘);grid  figureN=512subplot(2,1,1);plot(abs(fft(m2))*2/N);xlabel(‘頻率(Hz)‘);ylabel(‘幅值‘);title(‘原始訊號頻譜‘);grid;subplot(2,1,2);plot(abs(fft(result))*2/N);xlabel(‘頻率(Hz)‘);ylabel(‘幅值‘);title(‘線性濾波後‘);grid;subplot(2,1,2);plot(abs(fft(result))*2/N);xlabel(‘線性濾波後訊號頻譜‘);ylabel(‘幅值‘);grid;

以上程式的結果：

                                   圖10                                                                                  圖11

圖10是在時域濾波前後的心電訊號圖，可以看出，濾波後基準漂移得到了改善，圖11是在頻域濾波前後的心電訊號頻譜圖。

基於matlab的心電訊號預先處理