1、信号与系统课程设计报告课题二 心电信号分析系统的设计与仿真班级：姓名：学号：组号及同组人：指导教师：xx日期：2015 年 1 月 5 日心电信号分析系统的设计与仿真摘要：本文利用 MATLAB 对 MIT-BIH 数据库中的心电信号进行分析，利用MATLAB 软件、simulink 平台、GUI 图形用户界面、LABVIEW 软件对心电信号进行读取、插值、高通低通滤波等处理。并画出时域、频域波形进行比较分析。同时将滤波器的系统函数进行读取，分析，画出滤波的信号流程图，并画出系统的冲击响应、幅频响应、相位响应和零极点图来判断系统的稳定性。关键词：MATLAB，simulink，心电信号，数字滤
2、波器，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
3、 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:
4、MATLAB, simulink, ECG signal, digital filter, GUI, LABVIEW1、课程设计目的、意义本设计课题主要研究数字心电信号的初步分析方法及滤波器的设计与应用。通过完成本课题的任务，拟主要达到以下几个目的：1了解 MATLAB 软件的特点和使用方法，熟悉基于 Simulink 的动态建模和仿真的步骤和过程；2. 了解 LabVIEW 虚拟仪器软件的特点和使用方法，熟悉采用 LabVIEW 进行信号分析、系统设计及仿真的方法。3了解人体心电信号的时域特征和频谱特征；4通过设计具体的滤波器进一步加深对滤波器性能的理解；5掌握数字心电信号的分析方法，学会系
5、统设计与软件仿真方法；6通过本课题的训练，培养学生运用所学知识分析和解决实际问题的能力。二、课程设计任务及要求（一）基于 Matlab 的简单心电信号分析系统设计1.对原始数字心电信号进行读取，由数字信号数据绘制出其时域波形并加以分析。2.对数字信号数据做一次线性插值，使其成为均匀数字信号，以便后面的信号分析。3.根据心电信号的频域特征（自己查阅相关资料），设计相应的滤波器去除噪声。4.绘制进行信号处理前后的频谱，做频谱分析，得出相关结论。5.使用 GUI 进行系统的图形用户界面设计，（包含以上功能）。（二）基于 LabVIEW 虚拟仪器的简单心电信号分析系统设计1.进行心电信号的频谱分析，根
6、据心电信号的频域特征（自己查阅相关资料），设计相应的滤波器去除噪声。要求给出系统的前面板和框图，并记录仿真结果。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;
7、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)
8、;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,W
9、c=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
10、,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 NMb=b
11、 zeros(1,N-M);elseif MNa=a zeros(1,M-N);elseNM=0;endK=floor(N/2); B=zeros(K,3); A=zeros(K,3);if K*2=Nb=b 0; a=a 0;end broots=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:
12、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: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,:)
13、;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(
14、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;
15、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(H
16、2)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
17、,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;w
18、s=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)titl
19、e(滤波后心电信号频谱 )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);
20、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);
21、 %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(滤波后信号 )sub
22、plot(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
23、,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,
24、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
25、(滤波后信号 )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,
26、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
27、.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(
28、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; %
29、带阻滤波器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)
30、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（2）LABVI
31、EW 部分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 展开：
