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

信号与系统课程设计报告课题二心电信号分析系统的设计与仿真班级：姓名：学号：组号及同组人：指导教师：王宝珠日期：2015年1月5日心电信号分析系统的设计与仿真摘要：本文利用MATLAB对MIT-BIH数据库中的心电信号进行分析，利用MATLAB软件、simulink平台、GUI图形用户界面、LABVIEW软件对心电信号进行读取、插值、高通低通滤波等处理。并画出时域、频域波形进行比较分析。同时将滤波器的系统函数进行读取，分析，画出滤波的信号流程图，并画出系统的冲击响应、幅频响应、相位响应和零极点图来判断系统的稳定性。关键词：MATLAB，simulink，心电信号，数字滤波器，GUI，LABVIEWAbstract:ThisarticlemakesuseofMATLABtoanalyzeECGsignalofMIT-BIHECGDatabase.ToECGsignal.wecollectitfirst.thenwemakelinearinterpolation.finallywecarryavariableoffilterincludinglowpassandHighPass.wewillcomparedifferencesafterpaintingthetimedomainandfrequencydomainwaveform.atthesametimewereadandanalyzethesystemfunctionoffilterwithpaintingitstheflowchartofthesignal.fanallywepaintsystemshockresponsealongwithamplitude-frequencyresponseandphaseresponse.wejudgesystemstabilitybyZeropolefigure.Keywords:MATLAB,simulink,ECGsignal,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=C{2};b=C{1};k=length(b);fori=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;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);ifNMb=[bzeros(1,N-M)];elseifMNa=[azeros(1,M-N)];elseNM=0;endK=floor(N/2);B=zeros(K,3);A=zeros(K,3);ifK*2==Nb=[b0];a=[a0];endbroots=cplxpair(roots(b));%以下程序将每两个极点和两个零点组合成二阶因子aroots=cplxpair(roots(a));%roots：求多项式的根fori=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);fori=1:kc(i)=strread(b(i,:),'%*s%f','delimiter',':');endc=c';d=[c,a];fori=1:kifc(i)=2.5;%读取2.5秒e(i,:)=d(i,:);elsebreak;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,Xn);%线性插值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