通信工程课程设计——信号与线性系统课程设计

信号与系统课程设计报告课题二心电信号分析系统的设计与仿真班级：姓名：学号：组号及同组人：指导教师：王宝珠日期：2015年1月5日心电信号分析系统的设计与仿真摘要：本文利用MATLAB对MIT-BIH数据库中的心电信号进行分析，利用MATLAB软件、simulink平台、GUI图形用户界面、LABVIEW软件对心电信号进行读取、插值、高通低通滤波等处理。

并画出时域、频域波形进行比较分析。

同时将滤波器的系统函数进行读取，分析，画出滤波的信号流程图，并画出系统的冲击响应、幅频响应、相位响应和零极点图来判断系统的稳定性。

关键词：MATLAB，simulink，心电信号，数字滤波器，GUI，LABVIEWAbstract: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 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, digital filter, 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=C{2};b=C{1};k=length(b);for i=1:kc(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;[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 N>Mb=[b zeros(1,N-M)];elseif M>Na=[a zeros(1,M-N)];elseendK=floor(N/2); B=zeros(K,3); A=zeros(K,3);if K*2==Nb=[b 0]; a=[a 0];endbroots=cplxpair(roots(b)); %以下程序将每两个极点和两个零点组合成二阶因子aroots=cplxpair(roots(a)); % roots：求多项式的根for i=1:2:2*KBrow=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=C{2};b=C{1};k=length(b);for i=1:kc(i)=strread(b(i,:),'%*s %f','delimiter',':');endc=c';d=[c,a];for i=1:kif 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,'%g\n',Xn); fclose(fid);[t1,Xn1]=xianxingchazhi(t,X n); %线性插值fid=fopen('Xn1.txt','wt');%保存插值后信号fprintf(fid,'%g\n',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,'%g\n',Xn);fclose(fid);[t1,Xn1]=xianxingchazhi(t,Xn);fid=fopen('Xn1.txt','wt');fprintf(fid,'%g\n',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.*H2Y=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,'%g\n',Xn); fclose(fid);[t1,Xn1]=xianxingchazhi(t,X n);fid=fopen('Xn1.txt','wt'); fprintf(fid,'%g\n',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,'%g\n',Xn);fclose(fid);[t1,Xn1]=xianxingchazhi(t,X n);fid=fopen('Xn1.txt','wt'); fprintf(fid,'%g\n',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./100 0,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('系统不因果稳定'), endb0 =0.1672B = 1.0000 1.0560 1.00001.0000 0.6125 1.00001.0000 -1.9997 1.00001.0000 -1.9998 1.00001.0000 0.0000 -1.0000A = 1.0000 0.3973 0.93701.0000 0.0433 0.61821.0000 -1.9883 0.98891.0000 -1.9983 0.99871.0000 -1.2649 0.2960运行结果：系统不因果稳定图15图16（4）扩展题%主函数%将信号通过低通、高通、带阻滤波器程序[t,Xn]=duqushuju1('117.txt' ); %读取2.5S原心电信号fid=fopen('Xn.txt','wt');%保存原信号fprintf(fid,'%g\n',Xn); fclose(fid);[t1,Xn1]=xianxingchazhi(t,X n); %线性插值fid=fopen('Xn1.txt','wt');%保存插值后信号fprintf(fid,'%g\n',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.设计流程设计一个简单的心电信号分析系统。