基于MATLAB的心電信號分析心電信號分析(自己做的)帶程序 帶圖片.doc

河北工業(yè)大學信息工程學院 信號與線性系統(tǒng)課程設計 基于MATLAB的心電信號分析摘要： 本課題設計了一個簡單的心電信號分析系統(tǒng)。直接采用Matlab語言編程的靜態(tài)仿真方式、采用Simulink進行動態(tài)建模和仿真的方式，對輸入的原始心電信號，進行線性插值處理，并通過matlab語言編程設計對其進行時域和頻域的波形頻譜分析，根據(jù)具體設計要求完成系統(tǒng)的程序編寫、調(diào)試及功能測試。得出一定的結(jié)論。 關鍵字：matlab、心電信號提取、線性插值、濾波、simulink仿真。一、課題目的及意義心電信號是人類最早研究并應用于醫(yī)學臨床的生物信號之一，它比其它生物電信號更易于檢測，并且具有較直觀的規(guī)律性，因而心電圖分析技術促進了醫(yī)學的發(fā)展。然而，心電圖自動診斷還未廣泛應用于臨床，從國內(nèi)外的心電圖機檢測分析來看，自動分析精度還達不到可以替代醫(yī)生的水平，僅可以為臨床醫(yī)生提供輔助信息。其主要原因是心電波形的識別不準，并且心電圖診斷標準不統(tǒng)一。因此，探索新的方法以提高波形識別的準確率，尋找適合計算機實現(xiàn)又具診斷價值的診斷標準，是改進心電圖自動診斷效果，擴大其應用范圍的根本途徑。如何把心電信號的特征更加精確的提取出來進行自動分析，判斷出其異常的類型成了亟待解決的焦點問題。本課題通過matlab語言編程，對原始心電信號進行一定的分析處理。二、課題任務及要求1、必做部分（1）利用Matlab對MIT-BIH數(shù)據(jù)庫提供的數(shù)字心電信號進行讀取，并還原實際波形。（2）對原始心電信號做線性插值（3）對處理前后的心電信號分別做頻譜分析 利用Matlab軟件對處理前后的心電信號編程顯示其頻譜，分析比較濾 波前后的頻譜，得出結(jié)論。（4） Simulink仿真 根據(jù)前面的設計，進行基于Simulink的動態(tài)仿真設計。實現(xiàn)心電信號 的分析和處理。2、 選作部分（1） 只截取大約2.5s，三個周期左右，大約800個采樣數(shù)據(jù)進行分析。（2）60Hz工頻陷波器設計三、設計技術指標四、設計方案論證1、必做部分2、選作部分五、設計內(nèi)容及結(jié)果分析1、 基于matlab編寫的程序如下：%讀取心電信號并轉(zhuǎn)化成數(shù)組形式function t,Xn=duquexinhao1(w)fid=fopen(w);C=textscan(fid,%8c %f %*f,headerlines,2);%去除前兩行fclose(fid);a=C2;b=C1;k=length(b);for i=1:k c(i)=strread(b(i,:),%*s %f,delimiter,:);endc=c;d=c,a;t=d(:,1); %時間Xn=d(:,2); %幅度%線性插值function t3,Xn3=xianxingchazhi(t,Xn)m=max(t);t3=0:0.001:m;t3=t3;Xn3=interp1(t,Xn,t3);%保存插值前的信號function baocun1(t,Xn)fid = fopen(t.txt,wt);fprintf(fid,%gn,t); fclose(fid);fid = fopen(Xn.txt,wt);fprintf(fid,%gn,Xn); fclose(fid);%保存插值后的信號function baocun2(t1,Xn1)fid = fopen(t1.txt,wt);fprintf(fid,%gn,t1); fclose(fid);fid = fopen(Xn1.txt,wt);fprintf(fid,%gn,Xn1); fclose(fid); %畫初始信號和即插值后信號頻譜function keshehuatu(t,Xn,t1,Xn1)f=1000;T=1/f;m=1:length(Xn);k1=length(Xn1);m1=1:k1;q=f*m/length(Xn);q1=f*m1/k1;subplot(2,2,1)plot(t,Xn)title(初始信號時域波形)subplot(2,2,2)Y=fft(Xn);plot(q,abs(Y)title(初始信號頻譜)subplot(2,2,3)axis(0,1000,0,1000)plot(t1,Xn1)title(插值信號時域波形)Y1=fft(Xn1);subplot(2,2,4)axis(0,1000,0,5000)plot(q1,abs(Y1)title(插值信號頻譜)%低通濾波器function H,f=kesheditonglvboqi(wp,ws,Rp,As,Xn1)T=0.001;f=1/T;N,Wc=buttord(wp,ws,Rp,As,s);b,a=butter(N,Wc,s);f=(0:length(Xn1)-1)*f/length(Xn1);w=f*2*pi;H=freqs(b,a,w);%高通濾波器function H,f=keshegaotonglvboqi(wp,ws,Rp,As,Xn1)T=0.001;fs=1/T;N,Wc=buttord(wp,ws,Rp,As,s);b,a=butter(N,Wc,high,s);f=(0:length(Xn1)-1)*fs/length(Xn1);w=f*2*pi;H=freqs(b,a,w);%帶阻濾波器function H,f=keshedaizulvboqi(wp,ws,p,s,Xn1)T=0.001;f=1/T;N,Wc=buttord(wp,ws,p,s,s);b,a=butter(N,Wc,stop,s);f=(0:length(Xn1)-1)*f/length(Xn1);w=f*2*pi;H=freqs(b,a,w);主函數(shù)如下（1） 、將信號通過低通、高通、帶阻濾波器程序t,Xn=duquexinhao1(117.txt);baocun1(t,Xn) %保存讀取信號t1,Xn1=xianxingchazhi(t,Xn);baocun2(t1,Xn1)%保存插值后信號xy=t1,Xn1; %仿真輸入二維數(shù)組figure(1)keshehuatu(t,Xn,t1,Xn1) %畫原始信號和插值后信號波形和頻譜wp=90*2*pi; %低通濾波器濾波ws=99*2*pi;p=1;s=35;H1,f=kesheditonglvboqi(wp,ws,p,s,Xn1); wp=4*2*pi; %高通濾波器濾波ws=0.25*2*pi;p=1;s=35;H2,f=keshegaotonglvboqi(wp,ws,p,s,Xn1);wp=58,62*2*pi; %帶阻濾波器ws=59.9,60.1*2*pi;H3,f=keshedaizulvboqi(wp,ws,p,s,Xn1);H=abs(H1).*abs(H2).*abs(H3); %低通和高通和帶阻組合的濾波器Y=H.*abs(fft(Xn1); %經(jīng)過濾波后心電信號頻譜y=ifft(Y); %濾波后心電信號時域波形figure(2)subplot(2,2,1)plot(f,abs(H1)axis(0,150,0,1.5)title(低通濾波器)subplot(2,2,2)plot(f,abs(H2)axis(0,50,0,1.5)title(高通濾波器)subplot(2,2,3)plot(f,abs(H3)axis(0,150,0,1.5)title(帶阻濾波器)subplot(2,2,4)plot(f,abs(H)axis(0,100,0,1.5)title(組合后濾波器)figure(3)plot(f,abs(Y)axis(0,100,0,80)title(濾波后心電信號頻譜)figure(4)subplot(2,1,1)plot(t1,Xn1)title(濾波前信號)subplot(2,1,2)plot(t1,y)title(濾波后信號)所出圖形如下結(jié)果分析：（2） 、直接通過帶通濾波器程序t,Xn=duquexinhao1(117.txt);baocun1(t,Xn) %保存讀取信號t1,Xn1=xianxingchazhi(t,Xn);baocun2(t1,Xn1)%保存插值后信號figure(1)keshehuatu(t,Xn,t1,Xn1) %畫原始信號和插值后信號波形和頻譜wp=2,80*2*pi;ws=0.25,99*2*pi;p=1;s=35; H1,f=kesheditonglvboqi(wp,ws,p,s,Xn1);H=abs(H1) ; %帶通Y=H.*abs(fft(Xn1);%經(jīng)過濾波后心電信號頻譜y=ifft(Y); %濾波后心電信號時域波形figure(2)subplot(1,2,1)plot(f,abs(H1)axis(0,200,0,1.5)title(帶通濾波器)subplot(1,2,2)plot(f,abs(Y)axis(0,100,0,80)title(濾波后心電信號頻譜)figure(3)subplot(2,2,1)plot(t1,Xn1)title(濾波前信號)subplot(2,2,2)plot(t1,y)title(濾波后信號)subplot(2,2,3)plot(t1,Xn1)axis(0,1.5,-1.5,1.5)title(濾波前截取一部分信號)subplot(2,2,4)plot(t1,y)axis(0,1.5,-1.5,1.5)title(濾波后截取一部分信號) 所出圖形如下結(jié)果分析：（3） 、將信號通過低通、高通組合成的帶通濾波器程序t,Xn=duquexinhao1(117.txt);baocun1(t,Xn) %保存讀取信號t1,Xn1=xianxingchazhi(t,Xn);baocun2(t1,Xn1)%保存插值后信號figure(1)keshehuatu(t,Xn,t1,Xn1) %畫原始信號和插值后信號波形和頻譜xy=t1,Xn1;wp=0.52*2*pi; %低通濾波器濾波ws=0.62*2*pi;p=1;s=35;H1,f=kesheditonglvboqi(wp,ws,p,s,Xn1);wp=0.10*2*pi; %高通濾波器濾波ws=0.25*2*pi;p=1;s=35;H2,f=keshegaotonglvboqi(wp,ws,p,s,Xn1);H=abs(H1).*abs(H2); %低通和高通組合的帶通Y=H.*abs(fft(Xn1); %經(jīng)過濾波后心電信號頻譜y=ifft(Y); %濾波后心電信號時域波形figure(2)subplot(2,2,1)plot(f,abs(H1)axis(0,1,0,1.5)title(低通濾波器)subplot(2,2,2)plot(f,abs(H2)axis(0,1,0,1.5)title(高通濾波器)subplot(2,2,3)plot(f,abs(H)axis(0,1,0,1.5)title(組合 帶通濾波器)subplot(2,2,4)plot(f,abs(Y)axis(0,1,0,260)title(濾波后心電信號頻譜)figure(3)subplot(2,1,1)plot(t1,Xn1)title(濾波前信號)subplot(2,1,2)plot(t1,y)title(濾波后信號) 所出圖形如下結(jié)果分析：三種方案比較分析：（4）系統(tǒng)零極點分析（在此只以高通濾波器為例）%求高通濾波器的階數(shù)及分子分母系數(shù)wp=0.1*pi;ws=0.25*pi;Rp=1;As=35;T=1;%數(shù)字指標OmegaP=(2/T)*tan(wp/2);%通帶模擬頻率OmegaS=(2/T)*tan(ws/2);%阻帶模擬頻率cs,ds=afd_butt(OmegaP,OmegaS,Rp,As);%歸一化巴特沃斯濾波器原型系統(tǒng)函數(shù)N=ceil(log10(10(Rp/10)-1)/(10(As/10)-1)/(2*log10(wp/ws) OmegaC=wp/(10(Rp/10)-1)(1/(2*N); %求對應于N的3db截止頻率；b,a=u_buttap(N,OmegaC);%去歸一化巴特沃斯濾波器原型系統(tǒng)函數(shù)db,mag,pha,w=freqz_m(b,a);subplot(2,1,1);plot(w/pi,mag);title(digital filter Magnitude Response); axis(0,1,0,0.01)subplot(2,1,2);plot(w/pi,db);title(digital filter Magnitude in DB); axis(0,1,-30,5);%結(jié)果：N = 5% Butterworth Filter Order= 5%OmegaC = 0.3626%b = 0.0063%a = 1.0000 1.1734 0.6884 0.2496 0.0559 0.0063%N=6（5）求上述高通濾波器的系統(tǒng)函數(shù)及其頻譜b=0.0063;a= 1.0000 1.1734 0.6884 0.2496 0.0559 0.0063; h=impz(b,a); %系統(tǒng)的單位取樣響應figure(1);plot(h) %畫出單位取樣響應title(h(n)figure(2)fs=1000; H,f=freqz(b,a,256,fs); %求出系統(tǒng)的頻率響應mag=abs(H); %幅度響應ph=angle(H); %相位響應ph=ph*180/pi;subplot(2,1,1),plot(f,mag);grid %畫出幅度響應xlabel(frequency(Hz);ylabel(magnitude);title(|H(jw)|);subplot(2,1,2);plot(f,ph);grid %畫出相位響應xlabel(frequency(Hz);ylabel(phase);title(相位);figure(3)zr=roots(b) %求出系統(tǒng)的零點pk=roots(a) %求出系統(tǒng)的極點zplane(b,a); %zplane函數(shù)畫出零極點圖%結(jié)果：zr = Empty matrix: 0-by-1 %pk = -0.1120 + 0.3438i % -0.1120 - 0.3438i %-0.3674 %-0.2910 + 0.2156i %-0.2910 - 0.2156i圖形如下系統(tǒng)函數(shù)及級聯(lián)圖：結(jié)果分析：（7）Simulink仿真：（在此只取第一種方案）圖形如下：技術指標及結(jié)果分析：選作部分：1、 只截取大約2.5s，三個周期左右，大約800個采樣數(shù)據(jù)進行分析程序如下：function t,Xn=duqu2(w)fid=fopen(w);C=textscan(fid,%8c %f %*f,headerlines,2);fclose(fid);a=C2;b=C1; k=length(b);for i=1:k c(i)=strread(b(i,:),%*s %f,delimiter,:);endc=c;d=c,a;%截取2.5s的心電信號for i=1:k if c(i)=2.5 e(i,:)=d(i,:); else break; endendt=e(:,1); %時間Xn=e(:,2); %幅度調(diào)用程序：圖形如下：結(jié)果分析：2、60Hz工頻陷波器設計：程序如下：%60Hz工頻陷波器設計wp=58,62*2*pi; ws=59.9,60.1*2*pi;H3,f=keshedaizulvboqi(wp,ws,p,s,Xn1);plot(f,abs(H3)axis(0,150,0,1.5)title(60Hz工頻陷波器設計)圖形如下：分析：課題總結(jié)附錄一、參考文獻1 北京迪陽正泰科技發(fā)展公司.綜合通信實驗系統(tǒng)信號與系統(tǒng)指導書（第二版）. 2006,62 丁玉美.數(shù)字信號處理（第二版）.西安電子科技大學出版社，20013 吳大正. 信號與線性系統(tǒng)分析（第四版）. 高等教育出版社，2005,84 謝嘉奎. 電子線路-線性部分(第四版). 高等教育出版社，2003,25 陳后金. 信號分析與處理實驗. 高等教育出版社，2006，8二、 附錄設計原理1.心電信號的讀取txt格式的數(shù)據(jù)文件內(nèi)容及格式如圖1-1所示（以100.txt為例）。 圖1-1 txt格式的心電數(shù)據(jù)文件其中文件的第一列為采樣時間，第二列是在以MLII這種導聯(lián)方式所得到的采樣數(shù)據(jù)，第三列是以V5這種導聯(lián)方式所得到的采樣數(shù)據(jù)，全文件記錄了約為10s的心電數(shù)據(jù)，3600個采樣數(shù)據(jù)，每一行數(shù)據(jù)之間用Tab符分隔。由于數(shù)據(jù)文件中后兩列數(shù)據(jù)是對同一種心電信號進行不同的導聯(lián)方式所得到的采樣數(shù)據(jù)，所以可以只采用其中的一種采樣數(shù)據(jù)，摒棄另外一種，即可完成對此心電信號的分析。全部的心電文件記錄時間約為10s，共計12個左右周期的心電信號。實際設計心電信號數(shù)據(jù)文件時應注意：2.心電信號的線性插值處理根據(jù)上文中提到的插值公式，以此為原理，設計Matlab程序，對心電信號數(shù)據(jù)做線性插值處理。插值完以后的數(shù)據(jù)應該是時間均勻的、以0.001秒為間隔的。此步驟的實現(xiàn)主要是基于Matlab中的數(shù)組操作函數(shù)來實現(xiàn)，建議一定首先熟悉并掌握Matlab中的所有數(shù)組操作函數(shù)的作用和操作方法。其中一種插值方法的思路是：將第一步中讀取的心電信號數(shù)據(jù)的時間數(shù)據(jù)和幅值數(shù)據(jù)分別存放在一個一維數(shù)組中。然后利用for循環(huán)結(jié)構(gòu)把所有數(shù)據(jù)依次讀取進來。判斷時間數(shù)據(jù)數(shù)組中前后兩個相鄰的數(shù)據(jù)間隔是否為0.001s，如果是則判斷下一對相鄰兩個數(shù)據(jù)；如果間隔大于0.001s則進行一維插值處理。注意對時間數(shù)據(jù)做插值的同時一定不要忘記對幅值數(shù)據(jù)同樣做插值處理，時間數(shù)據(jù)和幅值數(shù)據(jù)一定是相互對應的。3. 濾波器設計3.1 模擬濾波器設計原理（1）模擬巴特沃思濾波器原理 巴特沃斯濾波器具有單調(diào)下降的幅頻特性：在小于截止頻率的范圍內(nèi)，具有最平幅度的響應，而在后，幅頻響應迅速下降。巴特沃思低通濾波器幅度平方函數(shù)為： (3-1)式中N為濾波器階數(shù)，為3dB截止角頻率。將幅度平方函數(shù)寫成s的函數(shù)形式： (3-2) 該幅度平方函數(shù)有2N個等間隔分布在半徑為的圓上的極點， 為了形成穩(wěn)定的濾波器，取左半平面的N個極點構(gòu)成，即： (3-3) 為使設計統(tǒng)一，將頻率歸一化，得到歸一化極點，相應的歸一化系統(tǒng)函數(shù)為： (3-4)多項式形式為： (3-5) （2）模擬切比雪夫濾波器原理 切比雪夫濾波器的幅頻特性具有等波紋特性，有兩種形式，在通帶內(nèi)等波紋、阻帶單調(diào)的是型濾波器，在通帶內(nèi)單調(diào)、在阻帶內(nèi)等波