基于matlab的心電信號預(yù)處理

1、基于matlab的心電信號預(yù)處理一、心電信號（1）心電信號的特性人體心電信號是非常微弱的生理低頻電信號，通常最大的幅值不超過5mV，信號頻率在005100Hz之間。心電信號是通過安裝在人體皮膚表面的電極來拾取的。由于電極和皮膚組織之間會發(fā)生極化現(xiàn)象，會對心電信號產(chǎn)生嚴(yán)重的干擾。加之人體是一個復(fù)雜的生命系統(tǒng)，存在各種各樣的其他生理電信號對心電信號產(chǎn)生干擾。同時由于我們處在一個電磁包圍的環(huán)境中，人體就像一根會移動的天線，從而會對心電信號產(chǎn)生50Hz左右的干擾信號。心電信號具有微弱、低頻、高阻抗等特性，極容易受到干擾，所以分析干擾的來源，針對不同干擾采取相應(yīng)的濾除措施，是數(shù)據(jù)采集重點考慮的一個問題。

2、常見干擾有如下幾種：工頻干擾基線漂移肌電干擾 心電信號具有以下幾個特點： ·信號極其微弱，一般只有0.054mV，典型值為1mV； ·頻率范圍較低，頻率范圍為0.135Hz，主要集中在520Hz； ·存在不穩(wěn)定性。人體內(nèi)部各器官問的相互影響以及各人的心臟位置、呼吸、年齡、是否經(jīng)常鍛煉等因素，都會使心電信號發(fā)生相應(yīng)變化； ·干擾噪聲很強。對心電信號進行測量時，必然要與外界聯(lián)系，但由于其自身的信號非常微弱，因此，各種干擾噪聲非常容易影響測量。其噪聲可能來自工頻(50Hz)干擾、電極接觸噪點、運動偽跡、肌電噪聲、呼吸引起的基線漂移和心電幅度變化以及其他電子設(shè)備

3、的機器噪聲等諸多方面。（2）心電信號的選擇本次實驗所采用的心電信號來自MIT-BIH庫，庫中有48組失常的心電信號，要在其中找出符合實驗要求的心電信號（即含有肌電干擾、工頻干擾和基線漂移）。（3）正常心電信號波形圖1是正常心電信號在一個周期內(nèi)的波形，由P波、QRS波群和T波組成。P波是由心房的去極化產(chǎn)生的，其波形比較小，形狀有些圓，幅度約為0.25mV，持續(xù)時間為0.080.11s。竇房結(jié)去極化發(fā)生在心房肌細(xì)胞去極化之前，因而在時間上要先于P波，只是竇房結(jié)處于心臟內(nèi)部，其電活動在體表難以采集。P-R間期是指P波起點和QRS波群起點所跨越的時間，是竇房結(jié)產(chǎn)生的興奮，經(jīng)過右心房、左心房、房室交接區(qū)

4、、房室束、左右束支之后，傳到到心室所需要的時間。在正常的體表心電圖中，P-R間期的值為0.120.2s，其中大部分時間是興奮在房室交界區(qū)內(nèi)傳導(dǎo)所需要的時間。P-R間期也稱為房室傳導(dǎo)時間。P-R段是指P波終點和QRS波群起點之間所跨越的時間。在正常的體表心電圖中，P-R段的心電信號電位值都是接近基線水平的很小點位。在P-R段期間，左右心房同時興奮，因而兩者產(chǎn)生的綜合電場對體表心電圖的影響較小。另外，此時的興奮還處于房室交界區(qū)和房室束特殊傳導(dǎo)系統(tǒng)中，沒有到達(dá)心室，因而沒有產(chǎn)生較大波動的體表心電圖信號。QRS波群是左右心室肌細(xì)胞一次發(fā)生去極化所產(chǎn)生的膜外負(fù)電位在體表的反應(yīng)。QRS波群的持續(xù)時間為0.

5、060.1s。由于心室肌細(xì)胞在興奮過程中的綜合電場向量多次發(fā)生改變，因而形成了體表心電圖中大小和方向多次發(fā)生變化的心電信號，其中QRS波群中第一個向下的波為Q波，第一個向上的波為R波，R波后面的為S波。S-T段是指QRS波群終點和T波起點之間所跨越的時間。S-T段期間，左右心室的肌細(xì)胞都處于興奮期間，因而兩者形成的綜合電場向量在體表心電圖中的貢獻(xiàn)非常小，導(dǎo)致S-T段心電信號處于大約基線的水平。T波由心室肌細(xì)胞的復(fù)極化產(chǎn)生，其幅度為0.10.8mV，持續(xù)時間為0.050.25s。由于復(fù)極化差異的存在，T波的方向和QRS波群主波的方向一致。在R波向上的情況下，T波的幅度一般都超過R波幅度的1/10

6、。Q-T間期是指QRS波群起點和T波終點所跨越的時間段，代表心室肌細(xì)胞開始去極化到結(jié)束復(fù)極化所需要的時間，與心率呈負(fù)相關(guān)。二、濾波器的選擇1.肌電干擾的濾除低通濾波器 通常來說，肌電信號的頻率為205000HZ，其主要成分的頻率與肌肉的類型有關(guān)，一般在30300HZ，而心電信號的頻率主要集中在520HZ，所以選擇低通濾波器來濾除肌電干擾。巴特沃斯濾波器的特點是通頻帶內(nèi)的頻率響應(yīng)曲線最為平坦，沒有起伏，而在阻頻帶則逐漸下降為零。巴特沃斯濾波器的振幅對角頻率單調(diào)下降，并且濾波器的階數(shù)越高，在阻頻帶幅度衰減速度越快，其他濾波器高階的振幅對角頻率圖和低階數(shù)的振幅對角頻率有不同的形狀。2.工頻干擾的抑制

7、帶陷濾波器工頻干由于供電網(wǎng)絡(luò)無所不在，因此50Hz的工頻干擾是最普遍的，也是心電信號的主要干擾來源。50HZ陷波器的軟件設(shè)計方法多種多樣，常見方法有小波變換濾波、自適應(yīng)濾波、模板匹配濾波等，但都需要手工計算獲得濾波器的參數(shù)，運算比較復(fù)雜。濾波器設(shè)計中，使用IIR濾波器，可使階數(shù)降低，運算量減少，但破壞了相位特性；使用FIR濾波器既能得到很好的濾波效果，是波形失真達(dá)到最下，而且，F(xiàn)IR濾波器可以做成線性相位特性，這正好是心電信號濾波所需要的。利用MATLAB設(shè)計FIR濾波器的方法有窗函數(shù)法、頻率抽樣法和切比雪夫逼近法等，本次課設(shè)采用窗函數(shù)法設(shè)計50HZ陷波濾波器。窗函數(shù)方法的基本思想是：首先根據(jù)

8、要求選擇一個適當(dāng)?shù)睦硐氲屯V波器，因為其脈沖響應(yīng)是非因果且無限長的，用最優(yōu)化窗結(jié)構(gòu)窗函數(shù)來截取它的脈沖響應(yīng)，從而得到線性相位和因果的FIR濾波器。Kaiser窗是接近最優(yōu)化窗結(jié)構(gòu)的窗函數(shù)，它可以根據(jù)不同的參數(shù)調(diào)整濾波器的各項指標(biāo)，因此采用Kaiser窗函數(shù)進行濾波器設(shè)計擾的抑制帶陷濾波器3.基線漂移的糾正零相移濾波器零相移濾波器是指一個信號序列經(jīng)過該濾波器濾波后相位不發(fā)生變化，即該濾波器系統(tǒng)函數(shù)的相位響應(yīng)為零。顯然，對于因果系統(tǒng)來說是不可能實現(xiàn)零相移的，在事先無法知道信號相位譜的情況下，實現(xiàn)零相移是不可能的。零相移只能是對非因果系統(tǒng)來說的。具體而言，零相移濾波器使用了當(dāng)前信號點前面和后面的信號

9、點所包含的信息，從本質(zhì)上說就是使用了“未來的信息”來消除相位失真。三、程序及結(jié)果1.心電信號讀取 因為對MIT-BIH庫不是很熟悉，在官網(wǎng)上看過之后，還是不懂（全英文，而且是醫(yī)學(xué)方面的。）。所以，此處的心電信號的讀取程序是來自網(wǎng)上的 rddata.m 。如果自己要用的話，在選取好要處理的心電信號后，把路徑更改，并選取合適的樣本數(shù)，就可以了。我選取的是MIT-BIH中的 109，樣本數(shù)為1500，下圖為心電信號讀取后的圖形：從圖2紅色曲線可以看到，波形上存在許多“毛刺”，并且其相位在發(fā)生變化（以波峰為例，各波峰大致不在一條水平線上，即所說的“基線漂移”），部分波形收到的干擾比較嚴(yán)重，比較符合對信

10、號處理的要求。2.心電信號的預(yù)處理（1）肌電信號的濾除plain view plain copy 在CODE上查看代碼片派生到我的代碼片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); %設(shè)計歸一化巴特沃斯模擬低通濾波器，z為極點，p為零點

11、和k為增益 bp,ap=zp2tf(z,P,k) %轉(zhuǎn)換為Ha(p),bp為分子系數(shù)，ap為分母系數(shù) bs,as=lp2lp(bp,ap,wp) %Ha(p)轉(zhuǎn)換為低通Ha(s)并去歸一化，bs為分子系數(shù)，as為分母系數(shù) hs,ws=freqs(bs,as); %模擬濾波器的幅頻響應(yīng) bz,az=bilinear(bs,as,Fs); %對模擬濾波器雙線性變換 h1,w1=freqz(bz,az); %數(shù)字濾波器的幅頻響應(yīng) m=filter(bz,az,M(:,1); figure freqz(bz,az);title('巴特沃斯低通濾波器幅頻曲線'); figure subp

12、lot(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/len

13、gth(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; %畫出頻

14、譜圖 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('心電信號的功率譜'); 以上程序的結(jié)果如下：圖3是所設(shè)計的巴特沃斯數(shù)字低通濾波器的幅頻響應(yīng)曲線，圖3是在時域濾波前后心電信號的波

15、形圖，圖5是在頻域濾波前后心電信號的頻譜圖，圖6是心電信號的功率譜圖（2）工頻干擾的抑制plain view plain copy 在CODE上查看代碼片派生到我的代碼片%-帶陷濾波器抑制工頻干擾- %50Hz陷波器：由一個低通濾波器加上一個高通濾波器組成 %而高通濾波器由一個全通濾波器減去一個低通濾波器構(gòu)成 Me=100; %濾波器階數(shù) L=100; %窗口長度 beta=100; %衰減系數(shù) Fs=1500; wc1=49/Fs*pi; %wc1為高通濾波器截止頻率，對應(yīng)51Hz wc2=51/Fs*pi ;%wc2為低通濾波器截止頻率，對應(yīng)49Hz h=ideal_lp(0.132*pi

16、,Me)-ideal_lp(wc1,Me)+ideal_lp(wc2,Me); %h為陷波器 沖擊響應(yīng) w=ser(L,beta); y=h.*rot90(w); %y為50Hz陷波器沖擊響應(yīng)序列 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)

17、-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

18、('帶阻濾波后信號');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; 其中，

19、ideal_lp（）函數(shù)在另一個M文件中，具體如下： %理想低通濾波器 %截止角頻率wc，階數(shù)Me function hd=ideal_lp(wc,Me) alpha=(Me-1)/2; n=0:Me-1; p=n-alpha+eps; %eps為很小的數(shù)，避免被0除 hd=sin(wc*p)./(pi*p); %用Sin函數(shù)產(chǎn)生沖擊響應(yīng) 以上程序的結(jié)果如下：圖7是帶陷濾波器的幅度譜和功率譜，從圖中可以看到在50Hz處，濾波器的幅度很大，而且功率在-150以下，說明帶陷性能較好。圖8是在時域濾波前后的心電信號圖，可以看出，濾波后波形有了略微的改善。圖19是在頻域濾波前后的心電信號頻譜圖。（3）基線漂移的糾正plain view plain copy 在CODE上查看代碼片派生到我的代碼片%-IIR零相移數(shù)字濾波器糾正基線漂移- Wp=1.4*2/Fs; %通帶截止頻率 Ws=0.6*2/Fs; %阻帶截止頻率 devel=0.005; %通帶紋波 Rp=20*log10(1+devel)/(1-devel); %通帶紋波系數(shù) Rs=20; %阻帶衰減 N Wn=ellipord(Wp,Ws,Rp,Rs,'s'); %求橢圓濾波器的階次 b a=ellip(N,Rp,Rs,Wn,'high'); %求橢圓濾波器的系數(shù) hw,w=freqz