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

1、基于matlab的心電信號(hào)預(yù)處理一、心電信號(hào)（1）心電信號(hào)的特性人體心電信號(hào)是非常微弱的生理低頻電信號(hào)，通常最大的幅值不超過5mV，信號(hào)頻率在005100Hz之間。心電信號(hào)是通過安裝在人體皮膚表面的電極來(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)考慮的一個(gè)問題。

2、常見干擾有如下幾種：工頻干擾基線漂移肌電干擾 心電信號(hào)具有以下幾個(gè)特點(diǎn)： ·信號(hào)極其微弱，一般只有0.054mV，典型值為1mV； ·頻率范圍較低，頻率范圍為0.135Hz，主要集中在520Hz； ·存在不穩(wěn)定性。人體內(nèi)部各器官問的相互影響以及各人的心臟位置、呼吸、年齡、是否經(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è)備

3、的機(jī)器噪聲等諸多方面。（2）心電信號(hào)的選擇本次實(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)過右心房、左心房、房室交接區(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)中，沒有到達(dá)心室，因而沒有產(chǎn)生較大波動(dòng)的體表心電圖信號(hào)。QRS波群是左右心室肌細(xì)胞一次發(fā)生去極化所產(chǎn)生的膜外負(fù)電位在體表的反應(yīng)。QRS波群的持續(xù)時(shí)間為0.

5、060.1s。由于心室肌細(xì)胞在興奮過程中的綜合電場(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波群主波的方向一致。在R波向上的情況下，T波的幅度一般都超過R波幅度的1/10

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

7、帶陷濾波器工頻干由于供電網(wǎng)絡(luò)無(wú)所不在，因此50Hz的工頻干擾是最普遍的，也是心電信號(hào)的主要干擾來(lái)源。50HZ陷波器的軟件設(shè)計(jì)方法多種多樣，常見方法有小波變換濾波、自適應(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濾波器的方法有窗函數(shù)法、頻率抽樣法和切比雪夫逼近法等，本次課設(shè)采用窗函數(shù)法設(shè)計(jì)50HZ陷波濾波器。窗函數(shù)方法的基本思想是：首先根據(jù)

8、要求選擇一個(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)過該濾波器濾波后相位不發(fā)生變化，即該濾波器系統(tǒng)函數(shù)的相位響應(yīng)為零。顯然，對(duì)于因果系統(tǒng)來(lái)說是不可能實(shí)現(xiàn)零相移的，在事先無(wú)法知道信號(hào)相位譜的情況下，實(shí)現(xiàn)零相移是不可能的。零相移只能是對(duì)非因果系統(tǒng)來(lái)說的。具體而言，零相移濾波器使用了當(dāng)前信號(hào)點(diǎn)前面和后面的信號(hào)

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

10、號(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'是確定巴特沃斯模擬濾波器階次和3dB 截止模擬頻率 z,P,k=buttap(n); %設(shè)計(jì)歸一化巴特沃斯模擬低通濾波器，z為極點(diǎn)，p為零點(diǎn)

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); %對(duì)模擬濾波器雙線性變換 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('原始心電信號(hào)波形');grid; subplot(2,1,2); plot(TIME,m); xlabel('t(s)');ylabel('mv');title('低通濾波后的時(shí)域圖形');grid; N=512 n=0:N-1; mf=fft(M(:,1),N); %進(jìn)行頻譜變換（傅里葉變換） mag=abs(mf); f=(0:length(mf)-1)*Fs/len

13、gth(mf); %進(jìn)行頻率變換 figure subplot(2,1,1) plot(f,mag);axis(0,1500,1,50);grid; %畫出頻譜圖 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; %畫出頻

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

15、形圖，圖5是在頻域?yàn)V波前后心電信號(hào)的頻譜圖，圖6是心電信號(hào)的功率譜圖（2）工頻干擾的抑制plain view plain copy 在CODE上查看代碼片派生到我的代碼片%-帶陷濾波器抑制工頻干擾- %50Hz陷波器：由一個(gè)低通濾波器加上一個(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)51Hz wc2=51/Fs*pi ;%wc2為低通濾波器截止頻率，對(duì)應(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('原始信號(hào)');grid; subplot(2,1,2);plot(TIME,m2); xlabel('t(s)');ylabel('幅值');title

18、('帶阻濾波后信號(hào)');grid; figure N=512 subplot(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; 其中，

19、ideal_lp（）函數(shù)在另一個(gè)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是在時(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(N,Rp,Rs,Wn,'high'); %求橢圓濾波器的系數(shù) hw,w=freqz