信号与系统课程设计报告心电信号分析系统的设计与仿真

信号与系统课程设计报告课题二 心电信号分析系统的设计与仿真班级：姓名：学号：组号及同组人：指导教师：日期：2015年1月5日心电信号分析系统的设计与仿真摘要：本文利用MATLAB对MIT-BIH数据库中的心电信号进行分析，利用MATLAB软件、simulink平台、GUI图形用户界面、LABVIEW软件对心电信号进行读取、插值、高通低通滤波等处理。并画出时域、频域波形进行比较分析。同时将滤波器的系统函数进行读取，分析，画出滤波的信号流程图，并画出系统的冲击响应、幅频响应、相位响应和零极点图来判断系统的稳定性。关键词：MATLAB，simulink，心电信号，数字滤波器，GUI，LABVIEWAbstract:Thisarticle makes use of MATLAB to analyze ECG signal of MIT-BIH ECG Database.To ECG signal .we collect it first.then we makelinear interpolation.finally we carry a variable of filter including lowpass andHigh Pass.we will compare differences after painting 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 pole figure.Key words:MATLAB, simulink, ECG signal, digitalfilter, GUI, LABVIEW一、 课程设计目的、意义本设计课题主要研究数字心电信号的初步分析方法及滤波器的设计与应用。通过完成本课题的任务，拟主要达到以下几个目的：1了解MATLAB软件的特点和使用方法，熟悉基于Simulink的动态建模和仿真的步骤和过程；2. 了解LabVIEW虚拟仪器软件的特点和使用方法，熟悉采用LabVIEW进行信号分析、系统设计及仿真的方法。3了解人体心电信号的时域特征和频谱特征；4通过设计具体的滤波器进一步加深对滤波器性能的理解；5掌握数字心电信号的分析方法，学会系统设计与软件仿真方法；6通过本课题的训练，培养学生运用所学知识分析和解决实际问题的能力。二、课程设计任务及要求（一）基于Matlab的简单心电信号分析系统设计1.对原始数字心电信号进行读取，由数字信号数据绘制出其时域波形并加以分析。2.对数字信号数据做一次线性插值，使其成为均匀数字信号，以便后面的信号分析。3.根据心电信号的频域特征（自己查阅相关资料），设计相应的滤波器去除噪声。4.绘制进行信号处理前后的频谱，做频谱分析，得出相关结论。5.使用GUI进行系统的图形用户界面设计，（包含以上功能）。（二）基于LabVIEW虚拟仪器的简单心电信号分析系统设计1.进行心电信号的频谱分析，根据心电信号的频域特征（自己查阅相关资料），设计相应的滤波器去除噪声。要求给出系统的前面板和框图，并记录仿真结果。 2.根据心电信号的特征，针对系统进行功能拓展，记录仿真结果，并进行相应的分析。三、设计方案过程及论证（一）matlab部分 1.设计流程：开始采集原始心电信号线性插值带通滤波器低通滤波器带阻滤波器高通滤波器带阻滤波器绘图结束2.程序（1）M文件%读取心电信号并转化为数组形式function t,Xn=duqushuju(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 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=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=(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=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);%滤波器的幅值响应、相位响应及群延迟响应function db,mag,pha,w=freqz_m(b,a);H,w=freqz(b,a,1000,whole); %在0-2*pi之间选取N个点计算频率响应H=(H(1:501); %频率响应 w=(w(1:501); %频率mag=abs(H); %响应幅度db=20*log10(mag+eps)/max(mag); %增益pha=angle(H); %变直接形式为级联形式 function b0,B,A=dir2cas(b,a)b0=b(1);b=b/b0;a0=a(1);a=a/a0;b0=b0/a0; %以上步骤求出系数b0M=length(b); N=length(a);if NM b=b zeros(1,N-M);elseif MN 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 b=b 0; a=a 0;end broots=cplxpair(roots(b); %以下程序将每两个极点和两个零点组合成二阶因子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心电信号并转化为数组形式function t,Xn=duqushuju1(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(:,2); %幅度（2）（3） 主程序%主函数 %将信号通过低通、高通、带阻滤波器程序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(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;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); %低通和高通和带阻组合的滤波网络Y=H.*abs(fft(Xn1); f1=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,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(fid);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工频陷波器设计H=H1.*H2 Y=H1.*H2.*abs(fft(Xn1); %经过滤波后心电信号频谱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工频陷波器设计)figure(2)plot(f,abs(Y)axis(0,100,0,80)title(滤波后心电信号频谱)figure(3)subplot(2,1,1)plot(t1,Xn1)title(滤波前信号)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);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*pi;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); %经过滤波后心电信号频谱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,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=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(fid); %画原始信号和插值后信号波形和频谱 xy=t1,Xn1; wp=90*2*pi; %低通滤波器滤波 ws=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); H=abs(H1).*abs(H2); %低通和高通组合的带通 Y=H.*abs(fft(Xn1); %经过滤波后心电信号频谱 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(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(滤波后信号) 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(带通阶跃响应)figure(2)db,mag,pha,w=freqz_m(b,a);subplot(1,2,1)plot(mag)axis(0,500,0,2)title(带通幅频特性)subplot(1,2,2)B1=roots(b) %求出系统的零点A1=roots(a) %求出系统的极点zplane(b,a) %zplane函数画出零极点图title(带通幅零极点图) b0,B,A=dir2cas(b,a) p=roots(A1) %求H(z)的极点 pm=abs(p); %求H(z)的极点的模 if max(pm)1 disp(系统因果稳定), else disp(系统不因果稳定), 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 运行结果：系统不因果稳定 图 15图 16（4） 扩展题%主函数 %将信号通过低通、高通、带阻滤波器程序t,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; %低通滤波器滤波ws=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); %低通和高通和带阻组合的滤波网络Y=H.*abs(fft(Xn1); f1=ifft(Y); %滤波后心电信号时域波形figure(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,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（二） LABVIEW部分1.设计流程 设计一个简单的心电信号分析系统。其基本功能包括：输入原始心电信号，对其做一定的数字信号处理，进行时域显示、分析及频谱分析。其具体流程如下：读取原始心电信号对原始信号进行线性插值开始对插值后的信号进行滤波观察及分析各阶段信号时域波形及频谱波形结束 2.整体设计框图及整体前面板展示图 26前面板：图 273.各部分设计原理及程序框图，滤波器参数设计 （一）信号读取 a.设计原理 麻省理工学院提供的MIT-BIH数据库，采用txt格式的数据文件作为我们的原心电信号数据。利用labVIEW提供的文件I/O函数，读取txt数据文件中的信号，并且还原实际波形。此数据库的数据文件的前两行为解释说明文字，不是真正的信号数据，读取信号程序需要自动忽略前两行文字，只读取真正的数字信号数据。另外labvIEW默认的从文本文件中读取的数据都是字符串，因此在使用心电信号数据前需要将其转换为数值才可以。第一列时间数据均为0:00.001这种格式，因此需要将字符串0:00.001先转化为字符串0.001，即去除字符串中冒号以前的部分，然后再将其转为数值。最后利用已经转为数值的分别代表心电信号时间和幅值的两个一维数组，图形化还原原始心电信号波形，图形化时可使用控件XY图。 b.程序框图图 28读取信号子VI展开：图 29（二）线性插值 a.设计原理 由于原始心电信号数据不是通过等间隔采样得到的，也就是说原始的心电数据并不是均匀的，而用Matlab中提供的数字滤波器处理数据时，要求数据是等间隔的。因此设计的系统首先应对原始心电信号做线性插值处理，使其变为等间隔的数字信号，否则直接处理后会出现偏差，根据心电信号的特点,把时间分隔成0.001s。添加的幅值点采用一次线性插值。对二维数据进行插值,相连幅值间数据的插值根据时间进行,运算公式如下:=+=-其中ti是第i个数据时间点，Ai是与之对应的数据，N是两数据之间需要的插值数， AD是需要插值的两点数据差，i,j=1时数组tj,Aj依次排列，即得到了插值后等间隔的新数据。 此步骤的实现主要是基于labvIEW中的数组操作函数（索引数组，数组子集等）来实现，插值方法的思路是：第一步中读取的心电信号数据的时间数据和幅值数据分别存放在一个一维数组中。然后利用For循环结构把所有数据依次读取进来。判断时间数据数组中前后两个相邻的数据间隔是否为0.001s，如果是则判断下一对相邻两个数据；如果间隔大于0.001s则在一个CASE条件结构里面做插值处理。对幅值数据同样做插值处理，时间数据和幅值数据一定是相互对应的。 b.程序框图 图 30线性插值子VI展开：图 31 （三）滤波器原理及滤波器参数设计 1. Butterworth滤波器图 32子VI展开：图 33 2. 低通滤波器图 34子VI展开：图 35 3. InverseChebyshev滤波器图36子VI展开：图 177 4. 60Hz工频陷波器图 38子VI展开：图 39 4.各部分调试结果及分析 （一）读取信号 图 40 图 41 （二）线性插值 图42 （三）巴特沃斯滤波器滤波图 43 （四）InverseChebyshev滤波器滤波 图 44 （五）60Hz工频陷波图 45 四、 实验结论参考文献1高西全，丁美玉编著.数字信号处理（第三版）.西安：西安电子科技大学出版社，2008.8.Page149-230.2张志涌，杨祖樱等编著.MATLAB教程：R2010a.北京：北京航空航天大学出版社，2010.8.Page252-302.3赵海滨等编著.MATLAB应用大全.北京：清华大学出版社， 2012.5.Page446-540.4 北京迪阳正泰科技发展公司.综合通信实验系统信号与系统指导书（第二版）. 2006,65 吴大正. 信号与线性系统分析（第四版）. 高等教育出版社，2005,86 谢嘉奎. 电子线路-线性部分(第四版). 高等教育出版社，2003,27 陈后金. 信号分析与处理实验. 高等教育出版社，2006，8