报告.docx

1、报告信号与线性系统课程设计报告（课题二 心电信号分析系统的设计与仿真）成绩：指导教师：刘翠响日期：2016年1月7日心电信号分析系统的设计与仿真摘要：心脏病发病率呈日渐上升趋势，是威胁全人类健康的首要疾病之一。心电信号分析是判断心脏疾病最直接最有效的手段，研究心电信号的处理和分析系统是非常有意义的。本文针对心电信号特征,采用MATLAB和LabVIEW两种软件实现心电信号分析系统的设计与仿真，主要实现对心电信号进行去噪和预处理。MATLAB的设计包括对心电信号进行提取、线性插值、滤波和分析，实现simulink仿真，并且利用GUI做出系统的图形用户界面设计。LabVIEW的设计包括在完成上述要

2、求的基础上使用选项卡等控件做出简洁友好的人机交互界面。关键词：心电信号，预处理，MATLAB，LabVIEW，人机交互界面1课程设计的目的及意义1.1 课程设计的目的1了解MATLAB软件的特点和使用方法，熟悉基于Simulink的动态建模和仿真的步骤和过程；2了解LabVIEW虚拟仪器软件的特点和使用方法，熟悉采用LabVIEW进行信号分析、系统设计及仿真的方法。3了解人体心电信号的时域特征和频谱特征；4通过设计具体的滤波器进一步加深对滤波器性能的理解；5掌握数字心电信号的分析方法，学会系统设计与软件仿真方法；6通过本课题的训练，培养学生运用所学知识分析和解决实际问题的能力。1.2 课程设计

3、的意义心电信号是组成人体心脏的肌细胞电活动在体表的综合表现,在一定程度上客观反映了心脏各部位的生理状况,因此,对心电信号的研究处理具有实用价值和重要意义。但由于心电信号自身特征的复杂性,在对其研究处理中仍存在一些有待于进一步深入和解决的问题。本课题通过使用Matlab和LabVIEW软件，实现简单的心电信号分析系统的设计及仿真，完成对心电信号的预处理及简单分析。2 课程设计任务及要求2.1 基于Matlab的简单心电信号分析系统设计设计一个简单的心电信号分析系统。其基本功能包括：输入原始心电信号，对其做一定的数字信号处理，进行时域显示、分析及频谱分析。采用Matlab软件设计相关程序。对基于M

4、atlab软件的程序设计，要求分别采用两种方式进行仿真，即直接采用Matlab语言编程的静态系统仿真方式、采用Simulink进行动态建模仿真的方式。根据心电信号的具体特性参数设计系统各功能模块的源程序，进行调试。具体要求如下：1.对原始数字心电信号进行读取，由数字信号数据绘制出其时域波形并加以分析。2.对数字信号数据做一次线性插值，使其成为均匀数字信号，以便后面的信号分析。3.根据心电信号的频域特征（自己查阅相关资料），设计相应的滤波器去除噪声。4.绘制进行信号处理前后的频谱，做频谱分析，得出相关结论。5.使用GUI进行系统的图形用户界面设计，（包含以上功能）。2.2 基于LabVIEW虚拟

5、仪器的简单心电信号分析系统设计1.进行心电信号的频谱分析，根据心电信号的频域特征（自己查阅相关资料），设计相应的滤波器去除噪声。要求给出系统的前面板和框图，并记录仿真结果。2.根据心电信号的特征，针对系统进行功能拓展，记录仿真结果，并进行相应的分析。3 设计方案及论证3.1基于Matlab的简单心电信号分析系统设计3.1.1设计流程图图1-1 基于MATLAB的心电信号分析系统设计流程图3.1.2设计原理1心电信号的读取：文件的第一列为采样时间，后两列是不同的导联方式所得到的采样数据，所以只读取前两列数据，并且忽略前两行解释说明文字，分别存放在两个数组中。2插值：把时间分隔成0.001s，添加

6、的幅值点采用一次线性插值。对二维数据进行插值,相连幅值间数据的插值根据时间进行。对时间数据做插值的同时对幅值数据同样做插值处理，时间数据和幅值数相互对应。3根据心电信号的频域特征，设计相应滤波器：一般正常人的心电信号频率在0.7100HZ范围内，幅度为10(胎儿)5mV(成人)，采集心电数据时，由于人的说话呼吸，常常会混有约为0.1Hz到0.25Hz频段的干扰。方案一：根据指导书中的要求，设置一个巴特沃斯模拟低通滤波器滤除高频信号干扰，设置一个巴特沃斯模拟高通滤波器来滤除低频信号干扰。低频滤波器的设计指标为：通带截止频率为92Hz，阻带起始频率为99Hz，通带最大衰减为1dB，阻带最小衰减为3

7、0dB。高频滤波器的设计指标为：通带起始频率为0.7Hz，阻带起始频率为0.25Hz, 通带最大衰减为1dB，阻带最小衰减为30dB。方案二：根据心电信号频谱范围，考虑到过渡带的衰减，我们设计了一个巴特沃斯数字带通滤波器并通过多次试验确定最合适的设计指标，滤波器的阻带范围为00.25Hz和170Hz及以上，通带范围为1.7Hz90Hz。将99Hz设在大概3dB截止频率处，这样能有效的滤除低频和高频信号干扰，又能避免有用的频带被滤除。4去除50Hz工频干扰：由于电子设备采集到的信号经常会混有电源线干扰。电源线干扰是以50 Hz为中心的窄带噪声，带宽小于1Hz。可以通过合适的带阻滤波器滤除电源线干

8、扰。方案一：从零极点分布对系统频率特性的影响出发,设计一个数字陷波器。系统函数完全由它的全部零点和极点来确定。对于因果稳定系统，极点必须全部在单位圆内。极点位置主要影响频响的峰值位置及尖锐程度，极点越靠近单位圆，峰值越高越尖锐；零点位置主要影响频响的谷点位置及形状，零点越靠近单位圆，谷值越接近零。对进行陷波,则取零点。但零点引起的凹谷会影响滤波器的通带范围，所以配置相应的极点来抵消这种影响，为了保证系统的因果稳定性，极点必须在单位圆内。取极点,应小于1且接近于1。于是,所设计的陷波器的系统函数为对50Hz进行陷波，则)，其中是采样频率，本实验中取，所以为，通过试验当取0.95时，滤波效果最好。

9、方案二：根据工频干扰的频带特征，设计一个巴特沃斯数字带阻滤波器，技术指标为：通带范围为045Hz和55Hz及以上，阻带范围为49.5Hz50.5Hz，通带最大衰减2dB，阻带最小衰减30dB。5截取原始信号的约3个周期：根据实验时选用的具体心电信号数据，截取心电信号中不正常的部分约3个周期。截取相应的时间段，再根据截取部分在时间数组中的位置，截取相应的信号幅值。6二次采样：根据采样频率及信号时域范围计算出采样间隔，根据采样间隔进行采样，就得到二次采样后的数据。所用公式如下： 其中T是时域二次采样间隔，是二次采样频率，是时域范围，N是时域二次采样点数，是第一次采样点数，A是二次采样时数据数组元素

10、间隔点数。3.1.3 仿真步骤、结果及分析1原始信号和插值后信号的时域频域分析图1-2 原始心电信号及插值后信号的波形和频谱分析：原始心电信号的频谱主要分布在0100Hz,插值后相邻两点间的幅值之差减小了，相当于波形变平缓。对应的频谱中可以看出，心电信号的低频分量幅度增加，高频分量幅度相对减少，能量重新分配。2. 相应的巴特沃斯低通、高通滤波器（其中以带通滤波器作为分析系统）图1-3 巴特沃斯模拟低通、高通滤波器的幅频特性图1-4 巴特沃斯带通滤波器的幅频特性级联型系统的常数及分子、分母的系数：b0 = 3.6999e-004B = 1.0000 2.0022 1.0022 1.0000 1.

11、9978 0.9978 1.0000 -1.9977 0.9977 1.0000 -2.0023 1.0023 1.0000 2.0000 1.0000 1.0000 -2.0000 1.0000A = 1.0000 -1.0280 0.2742 1.0000 -1.1331 0.4123 1.0000 -1.3855 0.7361 1.0000 -1.9788 0.9790 1.0000 -1.9991 0.99921.0000 -1.9861 0.9861分析：采用级联型系统，因为级联型调整零极点方便，有利于控制频率响应，运算误差的累积相对直接型小。该系统是一个6阶的因果稳定系统，可以实现

12、带通滤波功能。系统的信号流程图如下：图1-5 系统信号流图图1-6 系统的冲激响应和阶跃响应分析：在系统的时域特性中，冲激响应最终趋近于0，阶跃响应最终趋近于常数-0.38，说明系统是稳定的。图1-7 系统的零极点图分析：根据零极点图能定性分析系统函数的幅频响应。在z=1处，部分零点与极点抵消，剩下的零点使系统函数的频响在=0处幅度为零，其它极点决定系统函数的频响的峰值位置，极点越靠近单位圆峰值越高越尖锐，在z=-1处，只有零点，使系统函数的频响在=处幅度为零。3心电信号滤噪图1-8通过组合滤波器和带通滤波器滤噪后的信号波形和频谱图分析：根据指导书要求设计的滤波器，最终有效的滤除了0.25Hz

13、以下及99Hz以上的干扰，但靠近99Hz的部分有用信号也被滤除了。自行设计的带通滤波器将99Hz大概设在3dB截止频率处，既滤除了高频干扰，又保留了有用的信号。所以巴特沃斯数字带通滤波器的滤波效果较好。4滤除50Hz工频干扰图1-9 通过设置零极点设计的50Hz陷波器的幅频特性图1-10 50Hz陷波器的系统零极点图图1-11 巴特沃斯带阻滤波器的幅频特性图1-12 巴特沃斯带阻滤波器的零极点图图1-13 滤除50Hz工频干扰后的信号时域波形及频谱分析：采用设置零极点方法设计的陷波器只在60Hz处滤噪，阶数低，只有2阶，易实现，但是过渡带较宽，会影响其他频率分量的幅度。巴特沃斯带阻滤波器设计为

14、3阶，过渡带窄，避免了对其它频率分量的影响，但阻带有一定范围，会滤掉60Hz附近的部分频率分量。滤噪后，心电信号的时域波形更平滑。5截取信号的3个周期图1-14 截取3个周期后的原始及插值之后的心电信号波形及频谱图1-15 3个周期的心电信号滤噪后的波形及频谱分析：截取了心电信号中不正常的约3个周期的部分后，信号的时域范围缩小了，由可知，谱分辨率减小。滤噪后，信号波形更平滑。6二次采样图1-16 二次采样后的信号波形及频谱分析：取第二次采样频率为200Hz，二次采样后谱分辨率不变，谱分析范围缩小。3.1.4 GUI设计及实现图1-17 GUI初始界面点击“进入”，进入心电信号处理界面，如下：图

15、1-18 心电信号分析系统的使用界面导入数据后，点击“线性插值”、“截取三个周期”，选择巴特沃斯带通滤波器和50Hz陷波器分别进行滤波后，得到如下图形：图1-19 心电信号分析系统的GUI演示点击“保存图片”，保存成功后，点击“退出使用”，就可以退出GUI界面。保存的图片如下：图1-20 GUI界面中保存的心电信号图片3.2基于simulink的简单心电信号分析系统设计3.2.1仿真步骤及参数设置1数字带通滤波器和50Hz陷波器图2-1 建模方框图图2-2 数字带通滤波器参数设定及幅频特性图2-3 50Hz陷波器参数设定及幅频特性2数字低通滤波器+数字高通滤波器+50Hz陷波器图2-4 数字低

16、通+数字高通+50Hz陷波器建模方框图图2-5数字低通滤波器参数设定及幅频特性图2-6数字高通滤波器参数设定及幅频特性3. Simulink线性插值模块图2-7线性插值模块的建模方框图3.2.2仿真结果及分析 图2-7例231.txt插值后的波形及频谱1.数字带通滤波器和50Hz陷波器图2-8滤波后的波形图频谱图2数字低通滤波器+数字高通滤波器+60Hz陷波器图2-9滤波之后的波形图和频谱图3Simulink线性插值之后的波形图频谱图无明显变化。分析：经过多次修改参数，设计参数为：ws1=0.25Hz wp1=1.7Hz wp2=90Hz ws2=170Hz Rp=1 Rs=30 的数字带通滤

17、波器滤波效果最为明显。利用simulink进行动态仿真时我们设计并加上了50Hz工频陷波器，该陷波器是一个数字带阻滤波器设计参数为:wp1=48Hz ws1=49Hz ws2=51Hz wp2=53Hz Rp=1 Rs=30可以有效地滤除电源线干扰。要求显示原始波形和频谱以及滤完波之后的波形和频谱，可以加上Averaging Power Spectral Density显示模块可同时显示波形和频谱。3.3基于LabVIEW虚拟仪器的简单心电信号分析系统3.3.1设计原理 1心电信号的读取本部分内容在以子vi的形式在主程序中调用。图3-1读取心电信号.vi的程序流程图图3-2读取心电信号.vi的

18、程序框图2心电信号的线性插值处理得到双精度型的时间和幅值序列之后，使用索引函数得到时间序列的第一个和最后一个数据，再使用斜坡信号函数得到间隔为0.001的时间序列，最后使用一维插值函数对幅度序列做线性插值。图3-3心电信号的线性插值程序框图3设计相应的数字滤波器利用枚举输入控件，while循环结构和条件结构实现对带通滤波器和带阻滤波器种类和参数的实时改变。图3-4设计不同数字滤波器的程序框图图3-5选择不同的数字滤波器的前面板4.频谱分析使用XY图显示控件对原始信号，插值信号显示。使用波形图显示控件显示滤波前和滤波后的时域波形图和频谱图。图3-6使用XY图显示控件显示原始信号和插值信号的程序框

19、图图3-7 XY图显示原始信号和插值后信号前面板图3-8使用波形图显示控件显示时域波形图和频谱图的前面板5.使用选项卡控件设计简洁友好的人机交互界面图3-9 开始运行后的主界面图3-10 导入心电信号后的主界面图3-11 导入心电信号后的滤波前界面图3-12 导入心电信号的滤波后界面3.3.2 仿真调试步骤1导入心电信号文件选择想要分析的心电信号文件，然后导入，并把插值后的幅值序列保存。在导入文件后，程序运行，得到初始的时间幅值序列、插值后的时间幅值序列，插值前后的心电图以及滤波前后的波形图频谱图。2改变数字滤波器参数对于带通滤波器和带阻滤波器，分别设计了Butterworth、Chebysh

20、ev、Inverse Chebyshev三种滤波器。对于Butterworth滤波器可以更改滤波器阶数和3dB截止频率，对于Chebyshev和Inverse Chebyshev滤波器可以设置滤波器阶数、3dB截止频率和波纹衰减。 3.3.3分析相同阶数下，Chebyshev和Inverse Chebyshev滤波器要比Butterworth滤波器的过渡带更窄，滤波效果更好；对于Chebyshev滤波器通带波纹衰减越小，滤波器性能越好；对于Inverse Chebyshev滤波器阻带波纹衰减越大，滤波器性能越好。参考文献:1王璨，章佳荣编著.LabVIEW2011程序设计与案例分析M.北京：北

21、京航空航天大学出版社，20132罗华飞.MATLAB GUI设计学习手记.北京航空航天大学出版社，2009.83丁玉美.数字信号处理（第二版）.西安电子科技大学出版社，20014吴大正.信号与线性系统分析（第四版）. 高等教育出版社，2005.85丁亦农.Simulink与信号处理.北京航空航天大学出版社，2010.86王立会，潘东明.一种消除心电信号中工频干扰的陷波器设计. 2007.7附录：1心电信号的读取和插值程序（main.m文件）：clear,clct,xn=duqu(231.txt); %读取心电信号数据subplot(2,2,1),plot(t,xn)xlabel(t/s);yl

22、abel(xn);title(原始心电信号波形)Fs=1000;Xk=fft(xn,length(xn);n=0:length(xn)-1;f=Fs*n/length(xn);subplot(2,2,2),plot(f,abs(Xk)xlabel(f/Hz);ylabel(X(k);title(原始心电信号频谱)axis(-10,150,0,1500)t1,xn1=chazhi(t,xn); %对心电信号数据进行插值baocun(t1,xn1);subplot(2,2,3),plot(t1,xn1)xlabel(t/s);ylabel(xn1);title(插值后心电信号波形)Xk1=fft(

23、xn1,length(xn1);n1=0:length(xn1)-1;f1=Fs*n1/length(xn1);subplot(2,2,4),plot(f1,abs(Xk1)xlabel(f/Hz);ylabel(X(k)1);title(插值后心电信号频谱)axis(-10,150,0,1500)%simulink数据输入端simchazhi=t1,xn1;心电信号读取函数程序如下（duqu.m文件）：function t,Xn=duqu(w)fid=fopen(w); %打开文档wC=textscan(fid,%8c %f %*f,headerlines,2); %读出前两列数据且除去前两

24、行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);信号线性插值函数程序如下（chazhi.m文件）：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+

25、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);数据保存函数程序如下（baocun.m文件）：function baocun(t,Xn)fid = fopen(t.txt,wt);fprintf(fid,%gn,t); %将t数组写入文档t.txt 中fclose(fid);fid

26、 = fopen(Xn.txt,wt); fprintf(fid,%gn,Xn); %将t数组写入文档t.txt中 fclose(fid);2滤波器的设计及使用2.1 巴特沃斯低通滤波器和高通滤波器（lowhigh.m文件）%根据要求中给的设计指标，设计模拟巴特沃斯低通滤波器和高通滤波器clear,clct,xn=duqu(231.txt); %读取心电信号数据t1,xn1=chazhi(t,xn); %对心电信号数据进行插值Fs=1000;wp=92*2*pi;ws=99*2*pi;Rp=1;Rs=30;N,wc=buttord(wp,ws,Rp,Rs,s);b,a=butter(N,wc,

27、s);n=0:length(xn1)-1;f=n*Fs/length(xn1);w=f*2*pi;h=freqs(b,a,w);subplot(2,1,1),plot(f,abs(h)xlabel(f/Hz);ylabel(H(z);title(巴特沃斯模拟低通滤波器)axis(0,110,0,1.2)wph=0.7*2*pi;wsh=0.25*2*pi;Rp=1;Rs=30;Nh,wch=buttord(wph,wsh,Rp,Rs,s);bh,ah=butter(Nh,wch,high,s);hh=freqs(bh,ah,w);subplot(2,1,2),plot(f,abs(hh)xla

28、bel(f/Hz);ylabel(H(z);title(巴特沃斯模拟高通滤波器)axis(0,10,0,1.2)2.2 巴特沃斯带通滤波器（bandpass.m文件）%巴特沃斯数字带通滤波器的分析clear,clcwp=1.7*2*pi/1000,90*2*pi/1000;ws=0.25*2*pi/1000,70*2*pi/1000;Rp=1;Rs=30;N,wc=buttord(wp/pi,ws/pi,Rp,Rs);Nb,a=butter(N,wc);h,w=freqz(b,a);plot(w/max(w),abs(h); %幅频特性xlabel(w/pi);ylabel(H(k);titl

29、e(系统幅频特性);axis(-0.05,1,0,1.2);%直接型转换成级联型b0,B,A=dir2cas(b,a)%分析系统的时域特性y1=dimpulse(b,a,100); %冲击响应figuresubplot(2,1,1),plot(y1);xlabel(时间);ylabel(幅值);title(系统冲激响应);gridaxis(0,100,-0.2,0.4);y2=dstep(b,a,100); %阶跃响应subplot(2,1,2),plot(y2);xlabel(时间);ylabel(幅值);title(系统阶跃响应);grid%分析系统零极点figurelingjdt(a,b

30、);其中绘制零极点图函数程序如下（lingjdt.m文件）：function lingjdt(A,B)p=roots(A) ; % polesq=roots(B); % zerosp=p;q=q;x=max(abs(p,q);x=x+0.1;y=x;clfhold onaxis(-x x -y y);w=0:pi/300:2*pi;t=exp(i*w);plot(t); % 画圆axis(square);plot(-x x,0 0);plot(0 0,-y y);text(0.1,x,jIMz);text(y,1/10,Rez);plot(real(p),imag(p),x); % polesplot(real(q),imag(q),o); % zerosti