信號與系統(tǒng)課程設(shè)計報告心電信號分析系統(tǒng)的設(shè)計與仿真

1、信號與系統(tǒng)課程設(shè)計報告課題二 心電信號分析系統(tǒng)的設(shè)計與仿真班級：姓名：學(xué)號：組號及同組人：指導(dǎo)教師：日期：2015年1月5日心電信號分析系統(tǒng)的設(shè)計與仿真摘要：本文利用matlab對mit-bih數(shù)據(jù)庫中的心電信號進行分析，利用matlab軟件、simulink平臺、gui圖形用戶界面、labview軟件對心電信號進行讀取、插值、高通低通濾波等處理。并畫出時域、頻域波形進行比較分析。同時將濾波器的系統(tǒng)函數(shù)進行讀取，分析，畫出濾波的信號流程圖，并畫出系統(tǒng)的沖擊響應(yīng)、幅頻響應(yīng)、相位響應(yīng)和零極點圖來判斷系統(tǒng)的穩(wěn)定性。關(guān)鍵詞：matlab，simulink，心電信號，數(shù)字濾波器，gui，labviewa

2、bstract:this article makes use of matlab to analyze ecg signal of mit-bih  ecg database .to ecg signal .we collect it first.then we make linear interpolation.finally we carry a variable of  filter including  lowpass and high pass.we will compare differences after p

3、ainting the time domain and frequency domain waveform .at the same time we read and analyze the system function of filter with painting its the flow chart of the signal.fanally we paint system shock response along with amplitude-frequency response and phase response.we judge system stability by zero

4、 pole figure.key words:matlab, simulink, ecg signal, digital filter, gui, labview1、 課程設(shè)計目的、意義本設(shè)計課題主要研究數(shù)字心電信號的初步分析方法及濾波器的設(shè)計與應(yīng)用。通過完成本課題的任務(wù)，擬主要達到以下幾個目的：1了解matlab軟件的特點和使用方法，熟悉基于simulink的動態(tài)建模和仿真的步驟和過程；2. 了解labview虛擬儀器軟件的特點和使用方法，熟悉采用labview進行信號分析、系統(tǒng)設(shè)計及仿真的方法。3了解人體心電信號的時域特征和頻譜特征；4通過設(shè)計具體的濾波器進一步加深對濾波器性能的

5、理解；5掌握數(shù)字心電信號的分析方法，學(xué)會系統(tǒng)設(shè)計與軟件仿真方法；6通過本課題的訓(xùn)練，培養(yǎng)學(xué)生運用所學(xué)知識分析和解決實際問題的能力。二、課程設(shè)計任務(wù)及要求（一）基于matlab的簡單心電信號分析系統(tǒng)設(shè)計1.對原始數(shù)字心電信號進行讀取，由數(shù)字信號數(shù)據(jù)繪制出其時域波形并加以分析。2.對數(shù)字信號數(shù)據(jù)做一次線性插值，使其成為均勻數(shù)字信號，以便后面的信號分析。3.根據(jù)心電信號的頻域特征（自己查閱相關(guān)資料），設(shè)計相應(yīng)的濾波器去除噪聲。4.繪制進行信號處理前后的頻譜，做頻譜分析，得出相關(guān)結(jié)論。5.使用gui進行系統(tǒng)的圖形用戶界面設(shè)計，（包含以上功能）。（二）基于labview虛擬儀器的簡單心電信號分析系統(tǒng)設(shè)計

6、1.進行心電信號的頻譜分析，根據(jù)心電信號的頻域特征（自己查閱相關(guān)資料），設(shè)計相應(yīng)的濾波器去除噪聲。要求給出系統(tǒng)的前面板和框圖，并記錄仿真結(jié)果。 2.根據(jù)心電信號的特征，針對系統(tǒng)進行功能拓展，記錄仿真結(jié)果，并進行相應(yīng)的分析。三、設(shè)計方案過程及論證（一）matlab部分 1.設(shè)計流程：開始采集原始心電信號線性插值帶通濾波器低通濾波器帶阻濾波器高通濾波器帶阻濾波器繪圖結(jié)束2.程序（1）m文件%讀取心電信號并轉(zhuǎn)化為數(shù)組形式function t,xn=duqushuju(w)fid=fopen(w)c=textscan(fid,'%8c %f %*f','headerlines&

7、#39;,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 h,f=ditonglvboqi(wp,

8、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=gaotonglvboqi(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

9、=(0:length(xn1)-1)*fs/length(xn1);w=f*2*pi;h=freqs(b,a,w);%帶通濾波器function h,f=ditonglvboqi(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=daizulvboqi(wp,ws,p,s,xn1)t=0.001;f

10、=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);%濾波器的幅值響應(yīng)、相位響應(yīng)及群延遲響應(yīng)function db,mag,pha,w=freqz_m(b,a);h,w=freqz(b,a,1000,'whole'); %在0-2*pi之間選取n個點計算頻率響應(yīng)h=(h(1:501)' %頻率響應(yīng) w=(w(1:501)' %頻率m

11、ag=abs(h); %響應(yīng)幅度db=20*log10(mag+eps)/max(mag); %增益pha=angle(h); %變直接形式為級聯(lián)形式 function b0,b,a=dir2cas(b,a)b0=b(1);b=b/b0;a0=a(1);a=a/a0;b0=b0/a0; %以上步驟求出系數(shù)b0m=length(b); n=length(a);if n>m b=b zeros(1,n-m);elseif m>n a=a zeros(1,m-n);else nm=0;endk=floor(n/2); b=zeros(k,3); a=zeros(k,3);if k*2=n

12、 b=b 0; a=a 0;end broots=cplxpair(roots(b); %以下程序?qū)⒚績蓚€極點和兩個零點組合成二階因子aroots=cplxpair(roots(a); % roots：求多項式的根for i=1:2:2*k brow=broots(i:1:i+1,:); brow=real(poly(brow); b(fix(i+1)/2,:)=brow; arow=aroots(i:1:i+1,:); arow=real(poly(arow); a(fix(i+1)/2,:)=arow;end%讀取2.5s心電信號并轉(zhuǎn)化為數(shù)組形式function t,xn=duqushuj

13、u1(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;for i=1:k if c(i)<=2.5; %讀取2.5秒 e(i,:)=d(i,:); else break; endendt=e(:,1); %時間xn=e(

14、:,2); %幅度（2）（3） 主程序%主函數(shù) %將信號通過低通、高通、帶阻濾波器程序t,xn=duqushuju('117.txt'); %讀取原心電信號fid=fopen('xn.txt','wt'); %保存原信號fprintf(fid,'%gn',xn);fclose(fid);t1,xn1=xianxingchazhi(t,xn); %線性插值fid=fopen('xn1.txt','wt'); %保存插值后信號 fprintf(fid,'%gn',xn1);fclose(

15、fid);shuru=t1,xn1 figure(1)subplot(2,2,1)plot(t,xn)title('初始信號時域波形') %原心電信號時域波形subplot(2,2,2)y=fft(xn);plot(abs(y)title('初始信號頻譜') %原時域信號頻譜subplot(2,2,3)plot(t1,xn1)title('插值信號時域波形')y1=fft(xn1);subplot(2,2,4)plot(abs(y1)title('插值信號頻譜')wp=90*2*pi; %低通濾波器濾波ws=99*2*pi;p=1

16、;s=30;h1,f=ditonglvboqi(wp,ws,p,s,xn1);wp=1*2*pi; %高通濾波器濾波ws=0.25*2*pi;p=1;s=80;h2,f=gaotonglvboqi(wp,ws,p,s,xn1);wp=48,53*2*pi; %帶阻濾波器ws=49,51*2*pi;p=1;s=35;h3,f=daizulvboqi(wp,ws,p,s,xn1);h=abs(h1).*abs(h2).*abs(h3); %低通和高通和帶阻組合的濾波網(wǎng)絡(luò)y=h'.*abs(fft(xn1); f1=ifft(y); %濾波后心電信號時域波形figure(2)subplot(

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

18、e('濾波后心電信號頻譜')figure(4)subplot(2,1,1)plot(t1,xn1)title('濾波前信號')subplot(2,1,2)plot(t1,abs(f1)title('濾波后信號')axis(0,10,0,0.5)圖 1圖 2 圖 3圖 4圖 5%直接通過帶通濾波器程序后通過50hz工頻陷波器t,xn=duqushuju('117.txt');fid=fopen('xn.txt','wt');fprintf(fid,'%gn',xn);fclose(fi

19、d);t1,xn1=xianxingchazhi(t,xn);fid=fopen('xn1.txt','wt');fprintf(fid,'%gn',xn1);fclose(fid); wp=1,86*2*pi;ws=0.25,99*2*pi;p=1; %帶通s=30;h1,f=daitonglvboqi(wp,ws,p,s,xn1); wp=48,53*2*pi;ws=49,51*2*pi;h2,f=daizulvboqi(wp,ws,p,s,xn1); %50hz工頻陷波器設(shè)計h=h1.*h2 y=h1'.*h2'.*abs(

20、fft(xn1); %經(jīng)過濾波后心電信號頻譜y=ifft(y); %濾波后心電信號時域波形figure(1)subplot(1,2,1)plot(f,abs(h1)axis(0,200,0,1.5)title('帶通濾波器')subplot(1,2,2)plot(f,abs(h2)axis(0,150,0,1.5)title('50hz工頻陷波器設(shè)計')figure(2)plot(f,abs(y)axis(0,100,0,80)title('濾波后心電信號頻譜')figure(3)subplot(2,1,1)plot(t1,xn1)title(&

21、#39;濾波前信號')subplot(2,1,2) %6階帶阻濾波器plot(t1,y)title('濾波后信號') axis(0,10,0,0.5) 圖 6 圖 7圖 8 圖 9%直接通過帶通濾波器濾波程序t,xn=duqushuju('117.txt');fid=fopen('xn.txt','wt');fprintf(fid,'%gn',xn); fclose(fid);t1,xn1=xianxingchazhi(t,xn);fid=fopen('xn1.txt','wt

22、9;);fprintf(fid,'%gn',xn1);fclose(fid);shuru=t1,xn1 ;figure(1)subplot(2,2,1)plot(t,xn)title('初始信號時域波形')subplot(2,2,2)y=fft(xn);plot(abs(y)title('初始信號頻譜')subplot(2,2,3)plot(t1,xn1)title('插值信號時域波形')y1=fft(xn1);subplot(2,2,4)plot(abs(y1)title('插值信號頻譜')wp=1,86*2*p

23、i;ws=0.25,99*2*pi;p=1;s=30;h1,f=daitonglvboqi(wp,ws,p,s,xn1);h=abs(h1); %29階帶通濾波器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)subplo

24、t(2,2,1)plot(t1,xn1)title('濾波前信號')subplot(2,2,2)plot(t1,abs(y)axis(0,10,0,0.5)title('濾波后信號')subplot(2,2,3)plot(t1,xn1)axis(0,1.5,-1.5,1.5)title('濾波前截取的一部分信號')subplot(2,2,4)plot(t1,abs(y)axis(0,1.5,-1.5,1.5)title('濾波后截取一部分信號') 圖 10圖 11%將信號通過低通、高通組合成的帶通濾波器程序t,xn=duqushu

25、ju('117.txt');fid=fopen('xn.txt','wt');fprintf(fid,'%gn',xn);fclose(fid);t1,xn1=xianxingchazhi(t,xn);fid=fopen('xn1.txt','wt');fprintf(fid,'%gn',xn1);fclose(fid); %畫原始信號和插值后信號波形和頻譜 xy=t1,xn1; wp=90*2*pi; %低通濾波器濾波 ws=99*2*pi; p=1; s=30; h1,f=dit

26、onglvboqi(wp,ws,p,s,xn1); wp=1*2*pi; %高通濾波器濾波 ws=0.25*2*pi; p=1; s=80; h2,f=gaotonglvboqi(wp,ws,p,s,xn1); h=abs(h1).*abs(h2); %低通和高通組合的帶通 y=h'.*abs(fft(xn1); %經(jīng)過濾波后心電信號頻譜 y=ifft(y); %濾波后心電信號時域波形 figure(1) subplot(2,2,1) plot(f,abs(h1) axis(0,100,0,1.5) title('低通濾波器') subplot(2,2,2) plot(

27、f,abs(h2) axis(0,100,0,1.5) title('高通濾波器') subplot(2,2,3) plot(f,abs(h) axis(0,100,0,3) title('組合 帶通濾波器') subplot(2,2,4) plot(f,abs(y) axis(0,100,0,260) title('濾波后心電信號頻譜') figure(2) subplot(2,1,1) plot(t1,xn1) title('濾波前信號') subplot(2,1,2) plot(t1,y) title('濾波后信號&

28、#39;) axis(0,10,0,0.5)圖 12圖 13圖 14wp=1,90*2*pi;ws=0.25,99*2*pi;p=1; %帶通s=30;n,wc=ellipord(wp./1000,ws./1000,p,s);b,a=ellip(n,p,s,wc) h1=impz(b,a); figure(1)subplot(1,2,1)plot(h1)axis(0,10,-0.3,0.5)title('帶通h(n)')subplot(1,2,2)dstep(b,a,200)title('帶通階躍響應(yīng))')figure(2)db,mag,pha,w=freqz_

29、m(b,a);subplot(1,2,1)plot(mag)axis(0,500,0,2)title('帶通幅頻特性')subplot(1,2,2)b1=roots(b) %求出系統(tǒng)的零點a1=roots(a) %求出系統(tǒng)的極點zplane(b,a) %zplane函數(shù)畫出零極點圖title('帶通幅零極點圖') b0,b,a=dir2cas(b,a) p=roots(a1) %求h(z)的極點 pm=abs(p); %求h(z)的極點的模 if max(pm)<1 disp('系統(tǒng)因果穩(wěn)定'), else disp('系統(tǒng)不因果穩(wěn)

30、定'), end b0 = 0.1672b = 1.0000 1.0560 1.0000 1.0000 0.6125 1.0000 1.0000 -1.9997 1.0000 1.0000 -1.9998 1.0000 1.0000 0.0000 -1.0000a = 1.0000 0.3973 0.9370 1.0000 0.0433 0.6182 1.0000 -1.9883 0.9889 1.0000 -1.9983 0.9987 1.0000 -1.2649 0.2960 運行結(jié)果：系統(tǒng)不因果穩(wěn)定 圖 15圖 16（4） 擴展題%主函數(shù) %將信號通過低通、高通、帶阻濾波器程序t

31、,xn=duqushuju1('117.txt'); %讀取2.5s原心電信號fid=fopen('xn.txt','wt'); %保存原信號fprintf(fid,'%gn',xn);fclose(fid);t1,xn1=xianxingchazhi(t,xn); %線性插值fid=fopen('xn1.txt','wt'); %保存插值后信號 fprintf(fid,'%gn',xn1);fclose(fid);shuru=t1,xn1 wp=90*2*pi; %低通濾波器濾波w

32、s=99*2*pi;p=1;s=30;h1,f=ditonglvboqi(wp,ws,p,s,xn1);wp=1*2*pi; %高通濾波器濾波ws=0.25*2*pi;p=1;s=80;h2,f=gaotonglvboqi(wp,ws,p,s,xn1);wp=48,53*2*pi; %帶阻濾波器ws=49,51*2*pi;p=1;s=35;h3,f=daizulvboqi(wp,ws,p,s,xn1);h=abs(h1).*abs(h2).*abs(h3); %低通和高通和帶阻組合的濾波網(wǎng)絡(luò)y=h'.*abs(fft(xn1); f1=ifft(y); %濾波后心電信號時域波形figu

33、re(1)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(2)plot(f,abs(y)axis(0,

34、100,0,80)title('濾波后心電信號頻譜')figure(3)subplot(2,1,1)plot(t1,xn1)title('濾波前信號')subplot(2,1,2)plot(t1,abs(f1)title('濾波后信號') 圖17 圖18圖 19（4）simulink部分圖 20圖 21圖 22圖 23（5）gui部分圖 24圖 25（2） labview部分1.設(shè)計流程 設(shè)計一個簡單的心電信號分析系統(tǒng)。其基本功能包括：輸入原始心電信號，對其做一定的數(shù)字信號處理，進行時域顯示、分析及頻譜分析。其具體流程如下：讀取原始心電信號對原始

35、信號進行線性插值開始對插值后的信號進行濾波觀察及分析各階段信號時域波形及頻譜波形結(jié)束 2.整體設(shè)計框圖及整體前面板展示圖 26前面板：圖 273.各部分設(shè)計原理及程序框圖，濾波器參數(shù)設(shè)計 （一）信號讀取 a.設(shè)計原理 麻省理工學(xué)院提供的mit-bih數(shù)據(jù)庫，采用txt格式的數(shù)據(jù)文件作為我們的原心電信號數(shù)據(jù)。利用labview提供的文件i/o函數(shù)，讀取txt數(shù)據(jù)文件中的信號，并且還原實際波形。此數(shù)據(jù)庫的數(shù)據(jù)文件的前兩行為解釋說明文字，不是真正的信號數(shù)據(jù)，讀取信號程序需要自動忽略前兩行文字，只讀取真正的數(shù)字信號數(shù)據(jù)。另外labview默認(rèn)的從文本文件中讀取的數(shù)據(jù)都是字符串，因此在使用心電信號數(shù)據(jù)前需要將其轉(zhuǎn)換為數(shù)值才可以。第一列時間數(shù)據(jù)均為0:00.001這種格式，因此需要將字符串0:00.001先轉(zhuǎn)化為字符串0.001，即去除字符串中冒號以前的部分，然后再將其轉(zhuǎn)為數(shù)值。最后利用已經(jīng)轉(zhuǎn)為數(shù)值的分別代表心電信號時間和幅值的兩個一維數(shù)組，圖形化還原原始心電信號波形，圖形化時可使用控件xy圖。 b.程序框圖圖 28