信号与线性系统课程设计 心电信号分析系统的设计与仿真.doc

信号与线性系统课程设计报告课题二 心电信号分析系统的设计与仿真班级:通信C114姓名：胡伟学号：115665成绩：指导教师：日期：2013.1.5摘要： 本课题设计了一个简单的心电信号分析系统。直接采用Matlab语言编程的静态仿真方式、采用Simulink进行动态建模和仿真的方式，对输入的原始心电信号，进行线性插值处理，并通过matlab语言编程设计对其进行时域和频域的波形频谱分析，根据具体设计要求完成系统的程序编写、调试及功能测试。得出一定的结论。 关键字：matlab、心电信号提取、线性插值、滤波、simulink仿真。本课题的目的本设计课题主要研究数字心电信号的初步分析方法及滤波器的设计与应用。通过完成本课题的任务，拟主要达到以下几个目的：1了解MATLAB软件的特点和使用方法，熟悉基于Simulink的动态建模和仿真的步骤和过程；2. 了解LabVIEW虚拟仪器软件的特点和使用方法，熟悉采用LabVIEW进行信号分析、系统设计及仿真的方法。3了解人体心电信号的时域特征和频谱特征；4通过设计具体的滤波器进一步加深对滤波器性能的理解；5掌握数字心电信号的分析方法，学会系统设计与软件仿真方法；6通过本课题的训练，培养学生运用所学知识分析和解决实际问题的能力。2 设计任务及技术指标设计一个简单的心电信号分析系统。其基本功能包括：输入原始心电信号，对其做一定的数字信号处理，进行时域显示、分析及频谱分析。采用Matlab软件（或LabVIEW虚拟仪器软件）设计相关程序。对基于Matlab软件的程序设计，要求分别采用两种方式进行仿真，即直接采用Matlab语言编程的静态系统仿真方式、采用Simulink进行动态建模仿真的方式。根据心电信号的具体特性参数设计系统各功能模块的源程序，进行调试。1.对原始数字心电信号进行读取，由数字信号数据绘制出其时域波形并加以分析。2.对数字信号数据做一次线性插值，使其成为均匀数字信号，以便后面的信号分析。3.根据心电信号的频域特征（自己查阅相关资料），设计相应的滤波器去除噪声。4.绘制进行信号处理前后的频谱，做频谱分析，得出相关结论。5.对系统功能进行综合测试，整理数据，撰写设计报告。3主要设备和软件1PC机一台。2. MATLAB6.5以上版本软件，一套。4 设计内容以及实验结果与分析读取文件的function函数function t,Xn=duquexinhao1(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);保存文件1function baocun1(t,Xn)fid = fopen(t.txt,wt);fprintf(fid,%gn,t); fclose(fid);fid = fopen(Xn.txt,wt);fprintf(fid,%gn,Xn); fclose(fid);保存文件2function baocun2(t1,Xn1)fid = fopen(t1.txt,wt);fprintf(fid,%gn,t1); fclose(fid);fid = fopen(Xn1.txt,wt);fprintf(fid,%gn,Xn1); fclose(fid);插值文件function t2,Xn2=chazhi1(t,Xn)n=0;y=0;t=t.*1000; m=length(t);for i=1:me(i)=round (t(i);endfor i=1:(length(t)-1) if(e(i+1)-e(i)=1) N=(e(i+1)-e(i)/1; A=(Xn(i+1)-Xn(i)/N;for j=1:N z(y+j,1)=e(i)+(j-1)*1; z(y+j,2)=Xn(i)+(j-1)*A n=n+1; j=j+1; end y=n; endi=i+1;endz(y+1,2)=Xn(i);z(y+1,1)=t(i);t2=z(:,1);t2=t2./1000;Xn2=z(:,2);afd_butt文件function b,a = afd_butt(Wp,Ws,Rp,As); if Wp = 0 error(Passband edge must be larger than 0) end if Ws = Wp error(Stopband edge must be larger than Passband edge) end if (Rp = 0) | (As 0) error(PB ripple and/or SB attenuation ust be larger than 0) end N = ceil(log10(10(Rp/10)-1)/(10(As/10)-1)/(2*log10(Wp/Ws); fprintf(n* Butterworth Filter Order = %2.0f n,N) OmegaC = Wp/(10(Rp/10)-1)(1/(2*N); b,a=u_buttap(N,OmegaC);freqz_m文件functiondb,mag,pha,grd,w=freqz_m(b,a)H,w=freqz(b,a,1000,whole); H=(H(1:1:501);w=(w(1:501); mag=abs(H); db=20*log10(mag+eps)/max(mag); pha=angle(H); grd=grpdelay(b,a,w); u_buttap文件function b,a = u_buttap(N,Omegac);% Unnormalized Butterworth Analog Lowpass Filter Prototype% -% b,a = u_buttap(N,Omegac);% b = numerator polynomial coefficients of Ha(s)% a = denominator polynomial coefficients of Ha(s)% N = Order of the Butterworth Filter% Omegac = Cutoff frequency in radians/sec%z,p,k = buttap(N); p = p*Omegac; k = k*OmegacN; B = real(poly(z); b0 = k; b = k*B; a = real(poly(p);一.提取txt格式心电信号：fid=fopen(1.txt);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;plot(c,a)二.对原数据进行插值：tict,Xn=duquexinhao1(1.txt)baocun1(t,Xn);figure(1)plot(t,Xn)t2,Xn2=chazhi1(t,Xn);baocun2(t2,Xn2);figure(2)plot(t2,Xn2)toc 三.插值前后波形对比：subplot(2,2,1) plot(t,Xn) title(初始信号时域波形) axis(0 2.5 -1 2) subplot(2,2,2) fs=1000; N=length(Xn) n=1:N; f1=n*fs/N; Y1=fft(Xn); plot(f1,abs(Y1) title(初始信号频谱) axis(0 1000 0 200) subplot(2,2,3) plot(t2,Xn2) title(插值后信号时域波形) axis(0 2.5 -1 2) M=length(Xn2); m=1:M; f2=m*fs/M; Y2=fft(Xn2); subplot(2,2,4) plot(f2,abs(Y2) title(插值后信号频谱) axis(0 1000 0 200) tt=t2,Xn2四 .把 数 据 读 到txt中：fid = fopen(1.txt,wt); fprintf(fid,%gn,t); fclose(fid); fid = fopen(F.txt,wt); fprintf(fid,%gn,Xn2); fclose(fid); 五.模拟低通滤波器：wp=50*2*pi;ws=99*2*pi;Rp=1;As=40;%设置滤波器参数N,wc=buttord(wp,ws,Rp,As,s) %计算滤波器阶数N和3dB截止频率wc B,A=butter(N,wc,s) %计算滤波器系数函数分子分母多项式k=0:511;fk=0:1000/512:1000;wk=2*pi*fk; Hk=freqs(B,A,wk); plot(fk,20*log10(abs(Hk); axis(0,200,-40,5);结果N =8wc =349.7984B =1.0e+020 *Columns 1 through 6 0 0 0 0 0 0Columns 7 through 9 0 0 2.2415A =1.0e+020 *Columns 1 through 6 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000Columns 7 through 9 0.0002 0.0328 2.2415六.模拟高通滤波器： wp=0.6*2*pi;ws=0.25*2*pi;Rp=0.1;As=40; N,wc=buttord(wp,ws,Rp,As,s) B0,A0=butter(N,wc,s); wph=2*pi*0.25; BH,AH=lp2hp(B0,A0,wph); h,w=freqs(BH,AH); plot(w,20*log10(abs(h); axis(0,1,-80,5); 结果N = 8wc =2.7933七.滤波前后图形对比：fid=fopen(1.txt);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;figure(1)subplot(3,1,1); plot(c,a) wp=0.6*2*pi;ws=0.25*2*pi;Rp=0.1;As=40;T=1; N,wc=buttord(wp,ws,Rp,As,s) B,A=butter(N,wc,s); b,a=impinvar(B,A,T) db,mag,pha,w=freqz_m(b,a); y1=filter(b,a,Xn2); subplot(3,1,2);plot(y1); title(高通滤波后) wp1=2*pi*50;ws1=2*pi*99;Rp1=0.1;As1=40;T1=1000; OmegaP1=wp1/T1;OmegaS1=ws1/T1; cs1,ds1=afd_butt(OmegaP1,OmegaS1,Rp1,As1); b1,a1=impinvar(cs1,ds1,T) db1,mag1,pha1,w1=freqz_m(b1,a1); y2=filter(b1,a1,y1); subplot(3,1,3);plot(y2); title(低通滤波后) M=length(Xn2); m=1:M; fs=1000; f2=m*fs/M; F1=fft(Xn2); Y1=fft(y1); Y2=fft(y2) figure(2) subplot(3,1,1) plot(f2,abs(F1) axis(0,1000,0,200) title(原始信号频谱_9.997) subplot(3,1,2) plot(f2,abs(Y1) axis(0,1000,0,200) title(高通滤波后信号频谱_9.997) subplot(3,1,3) plot(f2,abs(Y2) axis(0,1000,0,200) title(低通滤波后信号频谱_9.997)八.系统函数的级联、冲击响应、幅度响应、相位响应： figure(3) H1=impz(b,a); H2=impz(b1,a1); Hn=conv(H1,H2); %将系统一的系统函数与系统二的函数相级联 subplot(3,1,1); plot(H1);title(高通系统函数冲击响应曲线); %系统函数的冲击响应曲线subplot(3,1,2); plot(H2);title(低通系统函数冲击响应曲线); subplot(3,1,3); plot(Hn);title(级联后系统函数冲击响应曲线); figure(4) fs=1000; H,f=freqz(B,A,256,fs); %系统一的幅频特性曲线mag=abs(H); %幅度响应ph=angle(H); %相位响应ph=ph*180/pi; subplot(2,2,1),plot(f,mag);grid %h1的幅度响应title(h1幅度响应) %xlabel(frequency(Hz);ylabel(magnitude); subplot(2,2,2);plot(f,ph);grid %h1的相位响应 title(h1相位响应 ) fs=1000; H,f=freqz(BH,AH,256,fs); mag=abs(H); ph=angle(H); ph=ph*180/pi; subplot(2,2,3),plot(f,mag);grid title(h2幅度响应) subplot(2,2,4);plot(f,ph);grid title(h2相位响应) figure(5) zr=roots(B) %系统一的零点图 pk=roots(A) %系统一的极点图zplane(B,A); %zplane函数画出系统一的零极点图 figure(6) zr1=roots(cs1) %系统二的零点图pk1=roots(ds1) %系统二的极点图zplane(cs1,ds1); %zplane函数画出系统一的零极点图 9. Simulink仿真以及仿真结果图 高通滤波器：低通滤波器：带阻滤波器：原信号波形：滤波后信号波形：滤波后 滤波前简要分析：5设计方案论证1.心电信号读取设计方案论证与选择美国麻省理工学院提供的MIT-BIH数据库是一个权威性的国际心电图检测标准库，近年来应用广泛，为我国的医学工程界所重视。MIT-BIH数据库共有48个病例，每个病例数据长30min，总计约有116000多个心拍，包含有正常心拍和各种异常心拍，内容丰富完整。为了读取简单方便，采用其txt格式的数据文件作为我们的原心电信号数据。利用Matlab提供的文件textread或textscan函数，读取txt数据文件中的信号，并且还原实际波形。（或者利用labVIEW提供的文件I/O函数，读取txt数据文件中的信号，并且还原实际波形。）由于不熟悉Labview提供的文件I/O函数，所以采用老师提供的Matlab提供的文件textscan。2.线性插值的设计方案与选择由于原始心电信号数据不是通过等间隔采样得到的，也就是说原始的心电数据并不是均匀的，而用Matlab中提供的数字滤波器处理数据时，要求数据是等间隔的。因此设计的系统首先应对原始心电信号做线性插值处理，使其变为等间隔的数字信号，否则直接处理后会出现偏差，根据心电信号的特点, 把时间分隔成0.001s。添加的幅值点采用一次线性插值。对二维数据进行插值,相连幅值间数据的插值根据时间进行,运算公式如下: ，其中是第i个数据时间点，Ai是与之对应的数据，N是两数据之间需要的插值数，是需要插值的两点数据差，时数组依次排列,即得到了插值后等间隔的新数据。虽然Labview提供的插值法有很多，但是由于自身对Labview不熟悉，所以采用教师提供的线性插值法（即上述所列公式）。3.滤波器的设计方案论证与选择一般正常人的心电信号频率在0.7100HZ范围内，幅度为(胎儿)5mV(成人)。人体心电信号微弱，信噪比小，因此，在采集心电信号时，易受到仪器、人体活动等因素的影响，而且所采集的心电信号常伴有干扰。采集心电数据时，由于人的说话呼吸，常常会混有约为0.1Hz到0.25Hz频段的干扰，对于这些低频干扰，可以让信号通过一个高频滤波器，低截止频率设置为0.25，来滤除低频信号，对于高频信号干扰，可以让信号再通过一个低频滤波器，其中截止频率设置为99Hz。所以我选择butterworth滤波器。 6总结7参考文献1 北京迪阳正泰科技发展公司.综合通信实验系统信号与系统指导书（第二版）. 2006,62 丁玉美.数字信号处理（第二版）.西安电子科技大学出版社，20013 吴大正. 信号与线性系统分析（第四版）. 高等教育出版社，2005,84 谢嘉奎. 电子线路-线性部分(第四版). 高等教育出版社，2003,25 陈后金. 信号分析与处理实验. 高等教育出版社，2006，86陈锡辉，张银鸿编著.LabVIEW 8.20 程序设计从入门到精通M.北京：清华大学出版社，2007. 8 附录1.心电信号的读取txt格式的数据文件内容及格式如图1-1所示（以100.txt为例）。 图1-1 txt格式的心电数据文件其中文件的第一列为采样时间，第二列是在以MLII这种导联方式所得到的采样数据，第三列是以V5这种导联方式所得到的采样数据，全文件记录了约为10s的心电数据，3600个采样数据，每一行数据之间用Tab符分隔。由于数据文件中后两列数据是对同一种心电信号进行不同的导联方式所得到的采样数据，所以可以只采用其中的一种采样数据，摒弃另外一种，即可完成对此心电信号的分析。全部的心电文件记录时间约为10s，共计12个左右周期的心电信号。实际设计心电信号数据文件时应注意：2.心电信号的线性插值处理根据上文中提到的插值公式，以此为原理，设计Matlab程序，对心电信号数据做线性插值处理。插值完以后的数据应该是时间均匀的、以0.001秒为间隔的。此步骤的实现主要是基于Matlab中的数组操作函数来实现，建议一定首先熟悉并掌握Matlab中的所有数组操作函数的作用和操作方法。其中一种插值方法的思路是：将第一步中读取的心电信号数据的时间数据和幅值数据分别存放在一个一维数组中。然后利用for循环结构把所有数据依次读取进来。判断时间数据数组中前后两个相邻的数据间隔是否为0.001s，如果是则判断下一对相邻两个数据；如果间隔大于0.001s则进行一维插值处理。注意对时间数据做插值的同时一定不要忘记对幅值数据同样做插值处理，时间数据和幅值数据一定是相互对应的。3. 滤波器设计3.1 模拟滤波器设计原理（1）模拟巴特沃思滤波器原理 巴特沃斯滤波器具有单调下降的幅频特性：在小于截止频率的范围内，具有最平幅度的响应，而在后，幅频响应迅速下降。巴特沃思低通滤波器幅度平方函数为： (3-1)式中N为滤波器阶数，为3dB截止角频率。将幅度平方函数写成s的函数形式： (3-2) 该幅度平方函数有2N个等间隔分布在半径为的圆上的极点， 为了形成稳定的滤波器，取左半平面的N个极点构成，即： (3-3) 为使设计统一，将频率归一化，得到归一化极点，相应的归一化系统函数为： (3-4)多项式形式为： (3-5) （2）模拟切比雪夫滤