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

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

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

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

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

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

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

7、速度越快，其他濾波器高階的振幅對(duì)角頻率圖和低階數(shù)的振幅對(duì)角頻率有不同的形狀。2. 工頻干擾的抑制帶陷濾波器工頻干由于供電網(wǎng)絡(luò)無(wú)所不在，因此50Hz 的工頻干擾是最普遍的，也是心電信號(hào)的主要干擾來(lái)源。 50HZ 陷波器的軟件設(shè)計(jì)方法多種多樣，常見(jiàn)方法有小波變換濾波、自適應(yīng)濾波、模板匹配濾波等，但都需要手工計(jì)算獲得濾波器的參數(shù)，運(yùn)算比較復(fù)雜。濾波器設(shè)計(jì)中，使用IIR 濾波器，可使階數(shù)降低，運(yùn)算量減少，但破壞了相位特性；使用FIR 濾波器既能得到很好的濾波效果，是波形失真達(dá)到最下，而且，F(xiàn)IR 濾波器可以做成線性相位特性，這正好是心電信號(hào)濾波所需要的。利用 MATLAB設(shè)計(jì) FIR 濾波器的方法有窗

8、函數(shù)法、頻率抽樣法和切比雪夫逼近法等，本次課設(shè)采用窗函數(shù)法設(shè)計(jì)50HZ 陷波濾波器。窗函數(shù)方法的基本思想是：首先根據(jù)要求選擇一個(gè)適當(dāng)?shù)睦硐氲屯V波器， 因?yàn)槠涿}沖響應(yīng)是非因果且無(wú)限長(zhǎng)的，用最優(yōu)化窗結(jié)構(gòu)窗函數(shù)來(lái)截取它的脈沖響應(yīng)，從而得到線性相位和因果的FIR 濾波器。 Kaiser窗是接近最優(yōu)化窗結(jié)構(gòu)的窗函數(shù)，它可以根據(jù)不同的參數(shù)調(diào)整濾波器的各項(xiàng)指標(biāo)，因此采用Kaiser 窗函數(shù)進(jìn)行濾波器設(shè)計(jì)擾的抑制帶陷濾波器3. 基線漂移的糾正零相移濾波器零相移濾波器是指一個(gè)信號(hào)序列經(jīng)過(guò)該濾波器濾波后相位不發(fā)生變化，即該濾波器系統(tǒng)函數(shù)的相位響應(yīng)為零。 顯然， 對(duì)于因果系統(tǒng)來(lái)說(shuō)是不可能實(shí)現(xiàn)零相移的，在事先無(wú)法知

9、道信號(hào)相位譜的情況下， 實(shí)現(xiàn)零相移是不可能的。 零相移只能是對(duì)非因果系統(tǒng)來(lái)說(shuō)的。具體而言，零相移濾波器使用了當(dāng)前信號(hào)點(diǎn)前面和后面的信號(hào)點(diǎn)所包含的信息，從本質(zhì)上說(shuō)就是使用了“未來(lái)的信息”來(lái)消除相位失真。三、程序及結(jié)果1. 心電信號(hào)讀取因?yàn)閷?duì) MIT-BIH 庫(kù)不是很熟悉， 在官網(wǎng)上看過(guò)之后，還是不懂 （全英文， 而且是醫(yī)學(xué)方面的。）。所以，此處的心電信號(hào)的讀取程序是來(lái)自網(wǎng)上的 rddata.m 。如果自己要用的話，在選取好要處理的心電信號(hào)后，把路徑更改，并選取合適的樣本數(shù)，就可以了。我選取的是MIT-BIH 中的 109 ，樣本數(shù)為1500，下圖為心電信號(hào)讀取后的圖形：從圖 2 紅色曲線可以看到

10、，波形上存在許多“毛刺” ，并且其相位在發(fā)生變化 （以波峰為例，各波峰大致不在一條水平線上， 即所說(shuō)的 “基線漂移” ），部分波形收到的干擾比較嚴(yán)重，比較符合對(duì)信號(hào)處理的要求。2. 心電信號(hào)的預(yù)處理（1）肌電信號(hào)的濾除plain view plain copy在 CODE上查看代碼片派生到我的代碼片clc;%-低 通 濾 波 器 濾 除 肌 電 信 號(hào)-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是確定巴特沃斯模

11、擬濾波器階次和3dB截止模擬頻率z,P,k=buttap(n); %設(shè)計(jì)歸一化巴特沃斯模擬低通濾波器，z 為極點(diǎn)， p 為零點(diǎn)和 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); %對(duì)模擬濾波器雙線性變換h1,w1=freqz(bz,az);%數(shù)字濾波器的幅頻響應(yīng)m=filter(bz,az,M(:,1

12、);figurefreqz(bz,az);title(巴特沃斯低通濾波器幅頻曲線 );figuresubplot(2,1,1);plot(TIME,M(:,1);xlabel(t(s);ylabel(mv);title(原始心電信號(hào)波形 );grid;subplot(2,1,2);plot(TIME,m);xlabel(t(s);ylabel(mv);title(低通濾波后的時(shí)域圖形 );grid;N=512n=0:N-1;mf=fft(M(:,1),N);%進(jìn)行頻譜變換（傅里葉變換）mag=abs(mf);f=(0:length(mf)-1)*Fs/length(mf);%進(jìn)行頻率變換fig

13、uresubplot(2,1,1)plot(f,mag);axis(0,1500,1,50);grid;%畫(huà)出頻譜圖xlabel( 頻率 (HZ);ylabel(幅值 );title(心電信號(hào)頻譜圖 );mfa=fft(m,N);%進(jìn)行頻譜變換（傅里葉變換）maga=abs(mfa);fa=(0:length(mfa)-1)*Fs/length(mfa);%進(jìn)行頻率變換subplot(2,1,2)plot(fa,maga);axis(0,1500,1,50);grid; %畫(huà)出頻譜圖xlabel( 頻率 (HZ);ylabel(幅值 );title(低通濾波后心電信號(hào)頻譜圖 );wn=M(:,

14、1);P=10*log10(abs(fft(wn).2)/N);f=(0:length(P)-1)/length(P);figureplot(f,P);gridxlabel(歸一化頻率 );ylabel(功率 (dB);title(心電信號(hào)的功率譜);以上程序的結(jié)果如下：圖 3 是所設(shè)計(jì)的巴特沃斯數(shù)字低通濾波器的幅頻響應(yīng)曲線，圖 3 是在時(shí)域?yàn)V波前后心電信號(hào)的波形圖，圖5 是在頻域?yàn)V波前后心電信號(hào)的頻譜圖，圖6 是心電信號(hào)的功率譜圖（2）工頻干擾的抑制plain view plain copy在 CODE上查看代碼片派生到我的代碼片%-帶陷濾波器抑制工頻干擾-%50Hz陷波器：由一個(gè)低通濾波器

15、加上一個(gè)高通濾波器組成%而高通濾波器由一個(gè)全通濾波器減去一個(gè)低通濾波器構(gòu)成Me=100;%濾波器階數(shù)L=100;%窗口長(zhǎng)度beta=100;%衰減系數(shù)Fs=1500;wc1=49/Fs*pi;%wc1為高通濾波器截止頻率，對(duì)應(yīng)51Hzwc2=51/Fs*pi;%wc2為低通濾波器截止頻率，對(duì)應(yīng)49Hzh=ideal_lp(0.132*pi,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

16、,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( 原始

17、信號(hào) );grid; subplot(2,1,2);plot(TIME,m2);xlabel(t(s);ylabel(幅值 );title(帶阻濾波后信號(hào));grid;figureN=512subplot(2,1,1);plot(abs(fft(m)*2/N);axis(0 100 0 1);xlabel(t(s);ylabel(幅值 );title(原始信號(hào)頻譜);grid;subplot(2,1,2);plot(abs(fft(m2)*2/N);axis(0 100 0 1);xlabel(t(s);ylabel(幅值 );title(帶阻濾波后信號(hào)頻譜);grid;其中， ideal_l

18、p（）函數(shù)在另一個(gè)M文件中，具體如下：%理想低通濾波器%截止角頻率wc，階數(shù) Mefunction 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 以下，說(shuō)明帶陷性能較好。圖8 是在時(shí)域?yàn)V波前后的心電信號(hào)圖，可以看出，濾波后波形有了略微的改善。圖19 是在頻域?yàn)V波前后的心電信號(hào)頻譜圖。（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