matlab心电信号的QRS波检测与分析_-精品资料.doc

燕山大学课 程 设 计 说 明 书 题目：基于matlab的心电信号QRS波检测与分析 学院（系）： 电气工程学院 年级专业： 09医疗仪器 学 号： 学生姓名： 指导教师： 孟辉 赵勇 教师职称： 讲师 讲师 燕山大学课程设计（论文）任务书院（系）：电气工程学院 基层教学单位：生物医学工程系 学 号学生姓名专业（班级）09医疗仪器设计题目基于matlab的心电信号QRS波检测与分析设计技术参数阈值=0.6相对幅值的差采样点数 N=256采样频率 f=100 Hz设计要求 能够准确的提取与分析QRS波绘制图表来说明检测结果工作量学会使用matlab软件了解人体生理信号的特征对心电信号有个初步的认识与掌握搜集一组正常人体心电信号数据，编写matlab程序，并对它进行分析与检测根据设计结果写出报告工作计划1.查资料，确认题目 2.编写任务书及审定 3.编写程序及调试4.编写说明书 5.完成任务书及设计要求 参考资料1 周辉 数字信号处理基础及Matlab实现 中国林业出版社 20052 肖伟 刘忠 Matlab程序设计与应用 清华大学出版社 20053 钱同惠编著 数字信号处理.北京：机械工业出版社，2004指导教师签字 孟辉 赵勇基层教学单位主任签字徐永红说明：此表一式四份，学生、指导教师、基层教学单位、系部各一份2012年12 月10 日 燕山大学课程设计说明书目录一、 MATLAB软件介绍2二、 概述3三、 ECG特征参数及分析4 3.1 心电信号的特点4 3.2 心电信号的特征参数4四、 QRS波得检测与分析5 4.1 以软件为主的方法实现QRS波的检测5 4.2 QRS波检测方法与程序 5五、 心得体会 14六、 参考文献 14一、 MATLAB软件介绍MATLAB是矩阵实验室的简称，是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。MATLAB是由美国MathWorks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。MATLAB的基本数据单元是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解决问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点。在新的版本中也加入了对C，FORTRAN，C+，JAVA的支持。可以直接调用，用户也可以将自己编写的使用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。二、 概述2.1 ECG处理的意义生物医学信号属于强噪声背景下的低频微弱信号，它是由复杂的生命体发出的不稳定的自然信号。作为一种对判断人体生命状况极其重要的生理信号，处理心电信号就显得很有必要，尤其在临床诊断上的应用。应用计算机分析心电信号，已经越来越广泛的用于心脏功能检查(Holter系统)、心电监护等方面，而心电分析中的首要的关键问题是QRS波的检测可靠的检测不仅是诊断心律失常的重要依据，而且只有在QRS波确定之后，有可能计算心率并进行心率变异分析，才能检测ST段的参数和分析心电的其它细节信息。进行全面综合分析，才能对心脏的功能结构做出正确的判断。在临床上，分析心电信号，可以确诊心肌梗塞及急性冠状动脉供血不足，协助诊断慢性冠脉供血不足、心肌炎、心肌病及心包炎，判定有无心房、心室肥大，从而协助某些心脏病的病因学诊断，例如风湿性、肺源性、高血压性和先天性心脏病等，观察某些药物对心肌的影响，包括治疗心血管疾病的药物（如洋地黄、抗心律失常药物）及可能对心肌有损害的药物。此外，对某些电解质紊乱（如血钾、血钙的过高或过低），心电信号不仅有助于诊断，还可以对指导治疗有重要参考价值。本文主要介绍心电信号的预处理和QRS复波检测的方法及演示结果。 三、 ECG特征参数及分析3.1心电信号的特点心电信号属于强噪声背景下的低频微弱信号，它是由复杂的生命体发出的不稳定的自然信号，由于受到人体诸多因素的影响，因而有着一般信号所没有的特点。信号弱。例如从母体腹部取到的胎儿心电信号仅为10v,成人的心电信号范围也仅为5mv.噪声强。由于人体自身信号弱，加之人体又是一个复杂的整体，因此信号易受噪声的干扰。如胎儿心电混有很强噪声，它一方面来自肌电、工频等干扰，另一方面，在胎儿心电中不可避免地含有母亲心电，母亲心电相对我们要提取的胎儿心电则变成了噪声。随机性强。心电信号信号不但是随机的，而且是非平稳的。正是因为生物医学信号的这些特点，使得心电信号处理成为当代信号处理技术最可发挥其威力的一个重要领域。3.2 心电信号的特征参数图2.1如图2.1为完整的心电信号波形图，分别由P波、PR段、PR间期、QRS复合波、ST段、T波和U波组成。本文将重点讲诉QRS波。并且对 它进行相关分析。QRS复波。代表两个心室兴奋传播过程的电位变化。由窦房结发生的兴奋波经传导系统首先到达室间隔的左侧面，以后按一定路线和方向，并由内层向外层依次传播。随着心室各部位先后去极化形成多个瞬间综合心电向量，在额面的导联轴上的投影，便是心电图肢体导联的QRS复合波。典型的QRS复合波包括三个相连的波动。第一个向下的波为Q波，继Q波后一个狭高向上的波为R波，与R波相连接的又一个向下的波为S波。由于这三个波紧密相连且总时间不超过0.10秒，故合称QRS复合波。QRS复合波所占时间代表心室肌兴奋传播所需时间，正常人在0.060.10秒之间 四、 QRS波得检测与分析4.1 以软件为主的方法实现QRS波的检测以软件为主的方法实现QRS波的检测滤波之后的信号一般经过一些变换以提高QRS波的份量，进而采用一系列阈值进行判别，这些阈值有固定阈值法，也有可变阈值法。前者由于可能的干扰或高P、高T波的存在，若其滤波后超过其阈值便会产生假阳性(FP，falsepositive)结果；另外，当心律失常或QRS波幅度变小，阈值设置过高，会导致漏检产生假阴性(FN，falsenegative)结果。由于固定阈值的这些缺点，有研究者提出了用可变阈值检测，以提高检测的精确率，所采用的可变阈值包括幅度阈值、斜率阈值和时间间隔阈值等。4.2 QRS波检测方法与程序Q波和S波通常是低幅高频波，一般Q波位于S波之前，S波位于R波之后 ，由于他们是一般向下的波，所以他们的峰值点和极值是对应的。因次在检测到R波向左和向右分别搜寻到极值点，对应的就是Q波和S波。具体程序如下：clear all;clc;z=textread('ECG.txt');ECG=z(:,1);input=ECG(1:256);rate=ECG(100);sig=input;lensig=length(sig);wtsig1=cwt(sig,6,'mexh');lenwtsig1=length(wtsig1);wtsig1(1:20)=0;wtsig1(lenwtsig1-20:lenwtsig1)=0;y=wtsig1;yabs=abs(y); %?sigtemp=y;siglen=length(y);sigmax=;for i=1:siglen-2 if (y(i+1)>y(i)&y(i+1)>y(i+2)|(y(i+1)<y(i)&y(i+1)<y(i+2) sigmax=sigmax;abs(sigtemp(i+1),i+1; end;end;%打印原信号及变换信号figure(1);subplot(2,1,1),plot(sig);subplot(2,1,2),plot(wtsig1);%取阈值,阈值为相对幅值的差的60%thrtemp=sort(sigmax);thrlen=length(sigmax);thr=0;for i=(thrlen-7):thrlen thr=thr+thrtemp(i);end;thrmax=thr/8; %最大幅度平均值，8个最大幅值点的平均值zerotemp=sort(y);zerovalue=0;for i=1:100 zerovalue=zerovalue+zerotemp(i);end;zerovalue=zerovalue/100; %最小幅度平均值，对消幅度，100个最小幅值点的平均值thr=(thrmax-zerovalue)*0.3; %最大、最小幅度的差值的30%为判别R波的阈值 %定位R波rvalue=;for i=1:thrlen if sigmax(i,1)>thr rvalue=rvalue;sigmax(i,2); end;end;rvalue_1=rvalue;%排除误检，如果相邻两个极大值间距小于0.4，则去掉幅度较小的一个lenvalue=length(rvalue);i=2;while i<=lenvalue if (rvalue(i)-rvalue(i-1)*rate<0.4 if yabs(rvalue(i)>yabs(rvalue(i-1) rvalue(i-1)=; else rvalue(i)=; end; lenvalue=length(rvalue); i=i-1; end; i=i+1;end; lenvalue=length(rvalue);%在原信号上精确校准for i=1:lenvalue if (wtsig1(rvalue(i)>0) k=(rvalue(i)-5):(rvalue(i)+5); a,b=max(sig(k); rvalue(i)=rvalue(i)-6+b; else k=(rvalue(i)-5):(rvalue(i)+5); a,b=min(sig(k); rvalue(i)=rvalue(i)-6+b; end;end;%打印纠正及校准前后的R波信号figure(2);subplot(2,1,1),plot(1:lensig,wtsig1,rvalue_1,wtsig1(rvalue_1),'r.');subplot(2,1,2),plot(1:lensig,sig,rvalue,sig(rvalue),'r.');%检测Q波wtsig2=cwt(sig,8,'mexh');lenrvalue=length(rvalue);qvalue=;for i=1:lenrvalue for j=rvalue(i):-1:(rvalue(i)-30) if wtsig1(rvalue(i)>0 if wtsig2(j)<wtsig2(j-1)&wtsig2(j)<wtsig2(j+1) tempqvalue=j-10; %确定检测窗的起点 break; %正向波，取第一个负极大值 end; else if wtsig2(j)>wtsig2(j-1)&wtsig2(j)>wtsig2(j+1) tempqvalue=j-10; %确定检测窗的起点 break; %倒置R波，取第一个正极大值 end; end; end; x1=tempqvalue; y1=sig(tempqvalue); x2=rvalue(i); y2=sig(rvalue(i); a0=(y2-y1)/(x2-x1); b0=-1; c0=-a0*x1+y1; %求直线公式参数ax+by+c=0 dist=; for k=tempqvalue:rvalue(i) tempdist=(abs(a0*k+b0*sig(k)+c0)/sqrt(a02+b02); dist=dist;tempdist; end; %求点到直线距离 a,b=max(dist); %找到距离最大值，Q波就在附近 tempqvalue=tempqvalue+b-1;% l=(tempqvalue-5):rvalue(i);% c,d=min(sig(l);% tempqvalue=tempqvalue-6+d; %在最大值附近修正Q波，得到结果 qvalue=qvalue;tempqvalue; end;%检测S波svalue=;for i=1:lenrvalue-1 for j=rvalue(i):1:(rvalue(i)+100) if wtsig1(rvalue(i)>0 if (wtsig2(j)<wtsig2(j-1)&(wtsig2(j)<wtsig2(j+1) tempsvalue=j+10; %在小波变换域从R波开始向后寻找第一个极小值 break; end; else if (wtsig2(j)>wtsig2(j-1)&(wtsig2(j)>wtsig2(j+1) tempsvalue=j+10; %在小波变换域从R波开始向后寻找第一个极大值 break; end;end; end; x1=tempsvalue; y1=sig(tempsvalue); x2=rvalue(i); y2=sig(rvalue(i); a0=(y2-y1)/(x2-x1); b0=-1; c0=-a0*x1+y1; %求直线公式参数ax+by+c=0 dist=; for k=rvalue(i):tempsvalue tempdist=(abs(a0*k+b0*sig(k)+c0)/sqrt(a02+b02); dist=dist;tempdist; end; %求点到直线距离 a,b=max(dist); %找到距离最大值，S波就在附近 tempsvalue=rvalue(i)+b-1;% l=rvalue(i):(tempsvalue+10);% c,d=min(sig(l);% tempsvalue=rvalue(i)+d-1; %在最大值附近修正S波，得到结果 svalue=svalue;tempsvalue; end;%检测QRS起点start=;for i=1:lenrvalue for j=qvalue(i):-1:(qvalue(i)-100) if wtsig1(j)>0 start=start;j; break; end; end;end;%打印Q,S波信号qrvalue=qvalue;rvalue;qrvalue=sort(qrvalue);qrsvalue=qvalue;rvalue;svalue;start;qrsvalue=sort(qrsvalue);figure(3);subplot(2,1,1),plot(1:lensig,sig,qrvalue,sig(qrvalue),'r.');subplot(2,1,2),plot(1:lensig,sig,qrsvalue,sig(qrsvalue),'r.');运行结果图如下：原信号及变换信号纠正及校准前后的R波信号Q,S波信号五、心得体会回顾起此课程设计，至今我仍感慨颇多，从理论到实践，在这段日子里，可以说得是苦多于甜，但是可以学到很多很多的东西，同时不仅可以巩固了以前所学过的知识，而且学到了很多在书本上所没有学到过的知识。通过这次课程设计使我懂得了理论与实际相结合是很重要的，只有理论知识是远远不够的，只有把所学的理论知识与实践相结合起来，从理论中得出结论，才能真正为社会服务，从而提高自己的实际动手能力和独立思考的能力。在设计的过程中遇到问题，可以说得是困难重重，但可喜的是最终都得到了解决。 最后我要感谢老师对我的耐心指导和同学对我的热心帮助。 六、参考文献1 周辉 数字信号处理基础及Matlab实现 中国林业出版社 20052 肖伟 刘忠 Matlab程序设计与应用 清华大学出版社 20053 钱同惠编著.数字信号处理.北京：机械工业出版社，20044 姚天任，江太辉编著.数字信号处理.第2版.武汉：武汉理工大学出版社，20005 黄文梅，熊佳林，杨勇编著.信号分析与处理MATALB言 及应用.长沙：国防科技大学出版社，2000- 15 - 燕山大学课程设计评审意见表指导教师评语：该生学习态度 （认真 较认真 不认真） 该生迟到、早退现象 （有 无）该生依赖他人进行设计情况 （有 无）平时成绩： 指导教师签字： 2012 年 12月15 日图面及其它成绩：答辩小组评语：设计巧妙，实现设计要求，并有所创新。 设计合理，实现设计要求。 实现了大部分设计要求。 没有完成设计要求，或者只实现了一小部分的设计要求。 答辩成绩： 组长签字： 2012 年12月15 日课程设计综合成绩：答辩小组成员签字： 2012年12月 15 日附录资料：matlab画二次曲面一、螺旋线1.静态螺旋线a=0:0.1:20*pi;h=plot3(a.*cos(a),a.*sin(a),2.*a,'b','linewidth',2);axis(-50,50,-50,50,0,150);grid onset(h,'erasemode','none','markersize',22);xlabel('x轴');ylabel('y轴');zlabel('z轴');title('静态螺旋线'); 2.动态螺旋线t=0:0.1:10*pi;i=1;h=plot3(sin(t(i),cos(t(i),t(i),'*','erasemode','none');grid onaxis(-2 2 -2 2 0 35)for i=2:length(t) set(h,'xdata',sin(t(i),'ydata',cos(t(i),'zdata',t(i); drawnow pause(0.01)endtitle('动态螺旋线');(图略) 3.圆柱螺旋线t=0:0.1:10*pi;x=r.*cos(t);y=r.*sin(t);z=t;plot3(x,y,z,'h','linewidth',2);grid onaxis('square')xlabel('x轴');ylabel('y轴');zlabel('z轴');title('圆柱螺旋线') 二、旋转抛物面b=0:0.2:2*pi;X,Y=meshgrid(-6:0.1:6);Z=(X.2+Y.2)./4;meshc(X,Y,Z);axis('square')xlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('旋转抛物面')或直接用：ezsurfc('(X.2+Y.2)./4') 三、椭圆柱面load clownezsurf('(2*cos(u)','4*sin(u)','v',0,2*pi,0,2*pi)view(-105,40) %视角处理shading interp %灯光处理colormap(map) %颜色处理grid on %添加网格线axis equal %使x,y轴比例一致xlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('椭圆柱面') %添加标题四、椭圆抛物面b=0:0.2:2*pi;X,Y=meshgrid(-6:0.1:6);Z=X.2./9+Y.2./4;meshc(X,Y,Z);axis('square')xlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('椭圆抛物面')或直接用：ezsurfc('X.2./9+Y.2./4')五、'双叶双曲面ezsurf('8*tan(u)*cos(v)','8.*tan(u)*sin(v)','2.*sec(u)',-pi./2,3*pi./2,0,2*pi)axis equalgrid onaxis squarexlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('双叶双曲面')六、双曲柱面load clownezsurf('2*sec(u)','2*tan(u)','v',-pi/2,pi/2,-3*pi,3*pi)hold on %在原来的图上继续作图ezsurf('2*sec(u)','2*tan(u)','v',pi/2,3*pi/2,-3*pi,3*pi)colormap(map)shading interpview(-15,30)axis equalgrid onaxis equalxlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('双曲柱面')七、双曲抛物面（马鞍面）X,Y=meshgrid(-7:0.1:7);Z=X.2./8-Y.2./6;meshc(X,Y,Z);view(85,20)axis('square')xlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('双曲抛物面')或直接用：ezsurfc('X.2./8-Y.2./6') 八、抛物柱面X,Y=meshgrid(-7:0.1:7);Z=Y.2./8;h=mesh(Z);rotate(h,1 0 1,180) %旋转处理%axis(-8,8,-8,8,-2,6);axis('square')xlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('抛物柱面')或直接用：ezsurfc('Y.2./8') 九、环面ezmesh('(5+2*cos(u)*cos(v)','(5+2*cos(u)*sin(v)','2*sin(u)',0,2*pi,0,2*pi)axis equalgrid onxlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('环面')十、椭球ezsurfc('(5*cos(u)*sin(v)','(3*sin(u)*sin(v)','4*cos(v)',0,2*pi,0,2*pi)axis equalgrid onxlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('椭球')十一、单叶双曲面ezsurf('4*sec(u)*cos(v)','2.*sec(u)*sin(v)','3.*tan(u)',-pi./2,pi./2,0,2*pi)axis equalgrid onxlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('单叶双曲面')十二、旋转单叶双曲面load clownezsurf('8*sec(u)*cos(v)','8.*sec(u)*sin(v)','2.*tan(u)',-pi./2,pi./2,0,2*pi)colormap(map)view(-175,30)%alpha(.2) %透明处理axis equalgrid onaxis squarexlabel('x轴');ylabel('y轴');zlabel('z轴');shading flat;title('旋转单叶双曲面')十三、圆柱面subplot(1,2,1)ezsurf('(2*cos(u)','2*sin(u)','v',0,2*pi,0,2*pi)grid onshading interpaxis equalxlabel('x轴');ylabel('y轴');zlabel('z轴');title('圆柱面')subplot(1,2,2)cylinder(30)shading interpaxis squaretitle('调用cylinder函数所得圆柱面')十四、二次锥面clc,clear;P=1,0,0; 0,cos(45*pi/180),sin(45*pi/180); 0,-sin(45*pi/180),cos(45*pi/180);for k2 = 1:31 for k1 = 1:31 x(k1,k2) = (k2-1)*cos ( (k1-1)*12*pi/180); y(k1,k2) = (k2-1)*sin ( (k1-1)*12*pi/180); z(k1,k2) = sqrt(x(k1,k2)2+y(k1,k2)2); Vxyz = P*x(k1,k2),y(k1,k2),z(k1,k2)' x1(k1,k2)=Vxyz(1); y1(k1,k2)=Vxyz(2); z1(k1,k2)=Vxyz(3); endendsurf(x,y,z)hold on;surf(x1,y1,z1);shading flat;