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

一、心電訊號

（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));

figure
freqz(bz,az);title('巴特沃斯低通濾波器幅頻曲線');
      
figure
subplot(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=512
n=0:N-1;
mf=fft(M(:,1),N);               %進行頻譜變換（傅立葉變換）
mag=abs(mf);
f=(0:length(mf)-1)*Fs/length(mf);  %進行頻率變換

figure
subplot(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);
figure
plot(f,P);grid
xlabel('歸一化頻率');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為高通濾波器截止頻率，對應51Hz
wc2=51/Fs*pi     ;%wc2為低通濾波器截止頻率，對應49Hz
h=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);

figure
subplot(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;
   
figure
subplot (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;
  
figure
N=512
subplot(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，階數Me
function 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); 

figure
freqz(b,a);
figure
subplot(211); plot(TIME,m2); 
xlabel('t(s)');ylabel('幅值');title('原始訊號');grid
subplot(212); plot(TIME,result); 
xlabel('t(s)');ylabel('幅值');title('線性濾波後訊號');grid
  
figure
N=512
subplot(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;s
ubplot(2,1,2);plot(abs(fft(result))*2/N);
xlabel('線性濾波後訊號頻譜');ylabel('幅值');grid;

以上程式的結果：

                                   圖10                                                                                  圖11

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