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

1、信號(hào)與系統(tǒng)課程設(shè)計(jì)報(bào)告課題二 心電信號(hào)分析系統(tǒng)的設(shè)計(jì)與仿真班級(jí)：姓名：學(xué)號(hào)：組號(hào)及同組人：指導(dǎo)教師：日期：2015年1月5日心電信號(hào)分析系統(tǒng)的設(shè)計(jì)與仿真摘要：本文利用matlab對(duì)mit-bih數(shù)據(jù)庫(kù)中的心電信號(hào)進(jìn)行分析，利用matlab軟件、simulink平臺(tái)、gui圖形用戶(hù)界面、labview軟件對(duì)心電信號(hào)進(jìn)行讀取、插值、高通低通濾波等處理。并畫(huà)出時(shí)域、頻域波形進(jìn)行比較分析。同時(shí)將濾波器的系統(tǒng)函數(shù)進(jìn)行讀取，分析，畫(huà)出濾波的信號(hào)流程圖，并畫(huà)出系統(tǒng)的沖擊響應(yīng)、幅頻響應(yīng)、相位響應(yīng)和零極點(diǎn)圖來(lái)判斷系統(tǒng)的穩(wěn)定性。關(guān)鍵詞：matlab，simulink，心電信號(hào)，數(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è)計(jì)目的、意義本設(shè)計(jì)課題主要研究數(shù)字心電信號(hào)的初步分析方法及濾波器的設(shè)計(jì)與應(yīng)用。通過(guò)完成本課題的任務(wù)，擬主要達(dá)到以下幾個(gè)目的：1了解matlab軟件的特點(diǎn)和使用方法，熟悉基于simulink的動(dòng)態(tài)建模和仿真的步驟和過(guò)程；2. 了解labview虛擬儀器軟件的特點(diǎn)和使用方法，熟悉采用labview進(jìn)行信號(hào)分析、系統(tǒng)設(shè)計(jì)及仿真的方法。3了解人體心電信號(hào)的時(shí)域特征和頻譜特征；4通過(guò)設(shè)計(jì)具體的濾波器進(jìn)一步加深對(duì)濾波器性能的

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

6、1.進(jìn)行心電信號(hào)的頻譜分析，根據(jù)心電信號(hào)的頻域特征（自己查閱相關(guān)資料），設(shè)計(jì)相應(yīng)的濾波器去除噪聲。要求給出系統(tǒng)的前面板和框圖，并記錄仿真結(jié)果。 2.根據(jù)心電信號(hào)的特征，針對(duì)系統(tǒng)進(jìn)行功能拓展，記錄仿真結(jié)果，并進(jìn)行相應(yīng)的分析。三、設(shè)計(jì)方案過(guò)程及論證（一）matlab部分 1.設(shè)計(jì)流程：開(kāi)始采集原始心電信號(hào)線(xiàn)性插值帶通濾波器低通濾波器帶阻濾波器高通濾波器帶阻濾波器繪圖結(jié)束2.程序（1）m文件%讀取心電信號(hào)并轉(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); %時(shí)間xn=d(:,2); %幅度%線(xiàn)性插值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個(gè)點(diǎn)計(jì)算頻率響應(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); %變直接形式為級(jí)聯(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ū)⒚績(jī)蓚€(gè)極點(diǎn)和兩個(gè)零點(diǎn)組合成二階因子aroots=cplxpair(roots(a); % roots：求多項(xiàng)式的根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心電信號(hào)并轉(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); %時(shí)間xn=e(

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

15、fid);shuru=t1,xn1 figure(1)subplot(2,2,1)plot(t,xn)title('初始信號(hào)時(shí)域波形') %原心電信號(hào)時(shí)域波形subplot(2,2,2)y=fft(xn);plot(abs(y)title('初始信號(hào)頻譜') %原時(shí)域信號(hào)頻譜subplot(2,2,3)plot(t1,xn1)title('插值信號(hào)時(shí)域波形')y1=fft(xn1);subplot(2,2,4)plot(abs(y1)title('插值信號(hào)頻譜')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); %濾波后心電信號(hào)時(shí)域波形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('濾波后心電信號(hào)頻譜')figure(4)subplot(2,1,1)plot(t1,xn1)title('濾波前信號(hào)')subplot(2,1,2)plot(t1,abs(f1)title('濾波后信號(hào)')axis(0,10,0,0.5)圖 1圖 2 圖 3圖 4圖 5%直接通過(guò)帶通濾波器程序后通過(guò)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è)計(jì)h=h1.*h2 y=h1'.*h2'.*abs(

20、fft(xn1); %經(jīng)過(guò)濾波后心電信號(hào)頻譜y=ifft(y); %濾波后心電信號(hào)時(shí)域波形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è)計(jì)')figure(2)plot(f,abs(y)axis(0,100,0,80)title('濾波后心電信號(hào)頻譜')figure(3)subplot(2,1,1)plot(t1,xn1)title(&

21、#39;濾波前信號(hào)')subplot(2,1,2) %6階帶阻濾波器plot(t1,y)title('濾波后信號(hào)') axis(0,10,0,0.5) 圖 6 圖 7圖 8 圖 9%直接通過(guò)帶通濾波器濾波程序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('初始信號(hào)時(shí)域波形')subplot(2,2,2)y=fft(xn);plot(abs(y)title('初始信號(hào)頻譜')subplot(2,2,3)plot(t1,xn1)title('插值信號(hào)時(shí)域波形')y1=fft(xn1);subplot(2,2,4)plot(abs(y1)title('插值信號(hào)頻譜')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)過(guò)濾波后心電信號(hào)頻譜y=ifft(y); %濾波后心電信號(hào)時(shí)域波形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('濾波后心電信號(hào)頻譜')figure(3)subplo

24、t(2,2,1)plot(t1,xn1)title('濾波前信號(hào)')subplot(2,2,2)plot(t1,abs(y)axis(0,10,0,0.5)title('濾波后信號(hào)')subplot(2,2,3)plot(t1,xn1)axis(0,1.5,-1.5,1.5)title('濾波前截取的一部分信號(hào)')subplot(2,2,4)plot(t1,abs(y)axis(0,1.5,-1.5,1.5)title('濾波后截取一部分信號(hào)') 圖 10圖 11%將信號(hào)通過(guò)低通、高通組合成的帶通濾波器程序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); %畫(huà)原始信號(hào)和插值后信號(hào)波形和頻譜 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)過(guò)濾波后心電信號(hào)頻譜 y=ifft(y); %濾波后心電信號(hào)時(shí)域波形 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('濾波后心電信號(hào)頻譜') figure(2) subplot(2,1,1) plot(t1,xn1) title('濾波前信號(hào)') subplot(2,1,2) plot(t1,y) title('濾波后信號(hào)&

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)的零點(diǎn)a1=roots(a) %求出系統(tǒng)的極點(diǎn)zplane(b,a) %zplane函數(shù)畫(huà)出零極點(diǎn)圖title('帶通幅零極點(diǎn)圖') b0,b,a=dir2cas(b,a) p=roots(a1) %求h(z)的極點(diǎn) pm=abs(p); %求h(z)的極點(diǎn)的模 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 運(yùn)行結(jié)果：系統(tǒng)不因果穩(wěn)定 圖 15圖 16（4） 擴(kuò)展題%主函數(shù) %將信號(hào)通過(guò)低通、高通、帶阻濾波器程序t

31、,xn=duqushuju1('117.txt'); %讀取2.5s原心電信號(hào)fid=fopen('xn.txt','wt'); %保存原信號(hào)fprintf(fid,'%gn',xn);fclose(fid);t1,xn1=xianxingchazhi(t,xn); %線(xiàn)性插值fid=fopen('xn1.txt','wt'); %保存插值后信號(hào) 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); %濾波后心電信號(hào)時(shí)域波形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('濾波后心電信號(hào)頻譜')figure(3)subplot(2,1,1)plot(t1,xn1)title('濾波前信號(hào)')subplot(2,1,2)plot(t1,abs(f1)title('濾波后信號(hào)') 圖17 圖18圖 19（4）simulink部分圖 20圖 21圖 22圖 23（5）gui部分圖 24圖 25（2） labview部分1.設(shè)計(jì)流程 設(shè)計(jì)一個(gè)簡(jiǎn)單的心電信號(hào)分析系統(tǒng)。其基本功能包括：輸入原始心電信號(hào)，對(duì)其做一定的數(shù)字信號(hào)處理，進(jìn)行時(shí)域顯示、分析及頻譜分析。其具體流程如下：讀取原始心電信號(hào)對(duì)原始

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