研究报告-基于matlab的操作员心电信号QRS波检测及分析.doc

基于matlab的操作员心电信号QRS波检测及分析1引言众所周知，疲劳是人的一个自然的生理现象，它的产生与一个人的劳动强度、工作环境、本身的健康和情绪等等有关，是人体的一种自我调节和保护功能。当疲劳被人们感觉时就意味着此刻需要好好休息一下。但是随着现在经济的高速发展与生活节奏的加快，许多工作人员因为种种原因，都是顶着疲劳继续工作，这样就使得人体产生的疲劳成为一些事故的根本导火索。本文要探索的是与人疲劳有关的心电信号，通过对人体心电信号的研究，来找出人的疲劳度量，预防人因为过度疲劳工作而导致事故的发生。本项目通过实验方法采集了执行过程控制任务的操作员的心电信号，通过确定阈值求极值确立R波，然后通过R波确立其前后最近的最低点分别为QS波信号，从而计算出操作员的心电QRS波等时间间隔，最后分析心电R波时间间隔和操作员主观疲劳程度、焦虑程度、努力程度和客观任务负荷的相关性。2心电原理 人体内蕴藏着大量丰富的生物电信号，如：脑电信号、心电信号、眼电信号、肌电信号、皮肤电压反应信号等，这些电信号都包含着大量的信息，它们大多都是在各种因素下由复杂生命体发出的不稳定的、幅值微弱的低频生理信号。而我们研究的便是其中的心电信号，心信号较之其他生物电信号更容易检测，便于实时监测。心脏周围的组织和体液都能导电，因此可将人体看成为一个具有长、宽、厚三度空间的容积导体。心脏好比电源，无数心肌细胞动作电位变化的总和可以传导并反映到体表。在体表很多点之间存在着电位差，也有很多点彼此之间无电位差是等电的。心脏在每个心动周期中，由起搏点、心房、心室相继兴奋，伴随着生物电的变化，这些生物电的变化称为心电。现在的科技已经可以通过采取一些有效设备和技术手段捕获人体可用的电信号，并对其进行信号采集与分析，进而有效破译其包含信息，利用这些大量的生命体征信息来对人的生活、工作、医疗等许多方面进行准确而有效的改善和调整。心电信号受到许多因素的影响，如：疲劳、疾病等等。而我们所要探索的便是与人疲劳有关的心电信号，通过对这类心电信号的研究，来找出人的疲劳与心电之间的关系，防止人因为过度疲劳工作而导致事故的发生。3 QRS波的检测方法3.1 以软件为主的方法实现QRS波的检测以软件为主的方法实现QRS波的检测在有大量数据需要检测的时候就显得方便快捷，极度占优。用软件来实现，就需要确定各波的阈值，当心律失常或QRS波幅度变小，阈值设置过高，会导致漏检产生假阴性(FN， falSe negative)结果。由于固定阈值的这些特点，可以以此来判断病人的心脏是否存在疾病以及隐患。3.2 QRS波检测算法Q波和S波通常是低幅高频波，一般Q波位于S波之前，S波位于R波之后 ，由于他们是一般向下的波，所以他们的峰值点和极值是对应的。因次只需要在检测到R波后向左和向右分别搜寻到极值点，对应的就是Q波和S波。 R波的检测利用其为心电信号图形中为最高点的特性，在整个数据中，只需要找出所有的峰值点，将所有峰值点从大到小排序，依据总数据长度以及采样频率来确定取前几位的平均值，确立R波的最大幅度平均值以及最小幅度平均值，确立R波阈值并取大于平均值的最大值为R波，在检测R波的过程中，我们尝试过使用小波变换来平滑数据，但是在寻找到R波并返回到原始数据上的时候，R波会发生偏移，得到的结果并不准确，所以最后在确定R波的过程中，未采用小波变换，而是直接进行逐个判断。Q、S波的检测在R波检测出来之后就容易了许多，根据QRS波的分布特性，只需从R波前后各推出一个最低点即可。4 眼电数据分析过程及结果 为方便观察，所有例举图形只截取部分（1）首先导入心电数据，然后将数据赋值给ECG具体实验步骤如下：function y=cacu(x); %定义函数，方便之后做相关性分析z=textRead(x); %读取原始数据为zECG=z(:,1); %将原始数据保存进变量ECGinput=ECG(1:150000); Rate=ECG(500); %将采样频率Rate设为500Sig=input; %将所采数据保存进变量SiglenSig=length(Sig); %将Sig长度保存为lenSiglenSig=length(Sig);Sig(1:40)=0; %将Sig1前40位置为0Sig(lenSig-39:lenSig)=0; %将Sig1后40位置为0y=Sig; %将Sig赋值给y Sigtemp=y; %将y赋值给SigtempSiglen=length(y); %将y的长度保存为SiglenSigmax=; %建立空矩阵SigmaxfoR i=1:Siglen-2 %令变量i依次从1加图1 部分心电信号图（2）根据心电信号特征求R波 foR i=1:Siglen-2 %令变量i依次从1加 if (y(i+1)=y(i)&y(i+1)y(i+2) Sigmax=Sigmax;abS(Sigtemp(i+1),i+1; end; %将y中的所有极点依次填入Sigmax end; thRtemp=SoRt(Sigmax); %将Sigmax内元素排序赋值给thRtempthRlen=length(Sigmax); %将Sigmax长度保存为thRlenthR=0;foR i=(thRlen-899):thRlen thR=thR+thRtemp(i); end;thRmax=thR/900; %最大幅度平均值（依据输入数据而定）zeRotemp=SoRt(y);zeRovalue=0;foR i=1:100 zeRovalue=zeRovalue+zeRotemp(i); end;zeRovalue=zeRovalue/1000; %最小幅度平均值，对消幅度（依据输入数据而定）thR=(thRmax-zeRovalue)*0.3; %最大、最小幅度的差值的30%为判别R波的阈值Rvalue=; %定位R波 foR i=1:thRlen if Sigmax(i,1)thRmax %寻找大于最大幅度平均值的极点，此处更改为thRmax Rvalue=Rvalue;Sigmax(i,2); end; end;lenRvalue=length(Rvalue); 图2 部分心电信号R波图（3）根据已知R波以及心电信号特征求解Q、S波(此处只例举解Q波的程序，求解S波的程序与其类似)： wtSig2=cwt(Sig,8,mexh); %对数据进行小波变换平滑处理Qvalue=; %设Q波集合foR i=1:lenRvalue foR j=Rvalue(i):-1:(Rvalue(i)-30) if Sig(Rvalue(i)0 if wtSig2(j)wtSig2(j-1)&wtSig2(j)wtSig2(j-1)&wtSig2(j)wtSig2(j+1) tempQvalue=j-9; %确定检测窗的起点 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; lenQvalue=length(Qvalue); %取Q波集合长度end;图3 部分心电信号QRS输出图形（4）得到所有QRS波之后，我们对数据进行了整理，提取出了各类特征值，并以图像形式输出。程序如下:foR i=1:(lenRvalue-1) tR(i)=(Rvalue(i+1)-Rvalue(i)/500;end; %计算R波时间间隔tQ=;foR i=1:(lenQvalue-1) tQ(i)=(Qvalue(i+1)-Qvalue(i)/500;end; %计算Q波时间间隔tS=;foR i=1:(lenSvalue-1) tS(i)=(Svalue(i+1)-Svalue(i)/500;end; %计算S波时间间隔 tQR=;foR i=1:lenRvalue tQR(i)=(Rvalue(i)-Qvalue(i)/500;end; %计算同次QR波时间间隔tRS=;foR i=1:lenSvalue tRS(i)=(Svalue(i)-Rvalue(i)/500;end; %计算同次RS波时间间隔tQS=;foR i=1:lenSvalue tQS(i)=(Svalue(i)-Qvalue(i)/500;end; %计算同次QS波时间间隔%R指标计算avetR=0; %定义RR时间间隔平均值 %定义RR间隔标准差foR i=1:lenRvalue-1 avetR=avetR+tR(i);end;图4 QRS波特征值（5）在求的QRS波的特征值之后，我们对其进行了进一步的计算得出了相邻的R、Q、S、QS、QR、RS波时间间隔的平均值，方差，标准差及标准差与均值的比值，具体程序如下：avetR=avetR/i;SDtR=0;foR i=1:lenRvalue-1 SDtR=SDtR+(tR(i)-avetR)2;end;vaRtR=SQRt(SDtR); %定义RR时间间隔方差SDR=SDtR/avetR; %定义RR时间间隔标准差除以平均值%Q指标计算avetQ=0;foR i=1:lenQvalue-1 avetQ=avetQ+tQ(i);end;avetQ=avetQ/i;SDtQ=0;foR i=1:lenQvalue-1 SDtQ=SDtQ+(tQ(i)-avetQ)2;end;vaRtQ=SQRt(SDtQ);SDQ=SDtQ/avetR;%S指标计算avetS=0;foR i=1:lenSvalue-1 avetS=avetS+tS(i);end;avetS=avetS/i;SDtS=0;foR i=1:lenSvalue-1 SDtS=SDtS+(tS(i)-avetS)2;end;vaRtS=SQRt(SDtS);SDS=SDtS/avetS;%QR指标计算avetQR=0;foR i=1:lenRvalue-1 avetQR=avetQR+tQR(i);end;avetQR=avetQR/i;SDtQR=0;foR i=1:lenRvalue-1 SDtQR=SDtQR+(tQR(i)-avetQR)2;end;vaRtQR=SQRt(SDtQR);SDQR=SDtQR/avetQR; avetRS=0; %RS指标计算foR i=1:lenSvalue avetRS=avetRS+tRS(i);end;avetRS=avetRS/i;SDtRS=0;foR i=1:lenRvalue; SDtRS=SDtRS+(tRS(i)-avetRS)2;end;vaRtRS=SQRt(SDtRS);SDRS=SDtRS/avetRS;%QS指标计算avetQS=0;foR i=1:lenSvalue avetQS=avetQS+tQS(i);end;avetQS=avetQS/i;SDtQS=0;foR i=1:lenSvalue SDtQS=SDtQS+(tQS(i)-avetQS)2;end;vaRtQS=SQRt(SDtQS);SDQS=SDtQS/avetQS;通过8组信号各得到的24组对应的特征值分别与已知任务难度标准0 1 3 4 4 3 1 0，疲劳程度标准0 20 30 40 60 80 60 30，焦虑程度标准0 10 10 10 10 20 10 10，努力程度标准0 20 30 40 50 50 40 10进行相关性分析分别得表1如下：表1 特征与任务难度相关性分析结果特征值ctRctQctSctQRctRSctQS相关性-0.8105-0.8105-0.81040.2913-0.11630.1558特征值cSDtRcSDtQcSDtScSDtQRcSDtRScSDtQS相关性0.44580.44640.44510.48580.64100.5562特征值cvaRRcvaRQcvaRScvaRQRcvaRRScvaRQS相关性0.44860.44960.44780.48860.62840.5569特征值cSDRcSDQcSDScSDQRcSDRScSDQS相关性0.47930.47970.47870.44700.64900.5131表2 特征与疲劳程度相关性分析特征值ctRctQctSctQRctRSctQS相关性-0.7026-0.7025-0.70250.54560.08460.3985特征值cSDtRcSDtQcSDtScSDtQRcSDtRScSDtQS相关性0.25450.25330.25370.25370.35740.0407特征值cvaRRcvaRQcvaRScvaRQRcvaRRScvaRQS相关性0.35590.35550.35460.25090.37570.0456特征值cSDRcSDQcSDScSDQRcSDRScSDQS相关性0.27660.27520.27580.2070.36060.0046表3 特征与焦虑程度相关性分析特征值ctRctQctSctQRctRSctQS相关性-0.6242-0.624-0.62410.46400.05530.3296特征值cSDtRcSDtQcSDtScSDtQRcSDtRScSDtQS相关性0.15010.15000.14940.27410.1437-0.0741特征值cvaRRcvaRQcvaRScvaRQRcvaRRScvaRQS相关性0.22580.22640.22460.25930.1616-0.0718特征值cSDRcSDQcSDScSDQRcSDRScSDQS相关性0.16840.16810.16770.23440.1447-0.098表4 特征与努力程度做相关性分析特征值ctRctQctSctQRctRSctQS相关性-0.8134-0.8134-0.81340.3809-0.05480.2370特征值cSDtRcSDtQcSDtScSDtQRcSDtRScSDtQS相关性0.36150.36090.36080.48740.50690.3302特征值cvaRRcvaRQcvaRScvaRQRcvaRRScvaRQS相关性0.42510.42510.42400.48870.5110.3361特征值cSDRcSDQcSDScSDQRcSDRScSDQS相关性0.39140.39050.39070.44240.51370.2912 （7）分析上表数据可得，所有四组程度标准都与相邻Q、R、S波的时间间隔相关性程度最高，即通过分析操作员当前心电信号的QRS频率可分析出操作员当前的各项数据。5. 总结在本项目中，我们通过编程实现了对QRS波群的检测及定位，并求出了相邻的R、Q、S、QS、QR、RS波时间间隔的平均值，方差，标准差及标准差与均值的比值，将其分别与任务难度系数、操作员的疲劳程度、焦虑程度、努力程度进行了相关性分析。相关性分析的结果表明：任务难度系数、操作员的疲劳程度、焦虑程度、努力程度等都可以从操作员的QRS波的平均时间间隔表现出来。这个结果可以广泛应用到各个领域，能够帮助掌握操作员的当前工作状态，以保障其工作效率最优化。6.项目组对项目的自我评价在落实这个项目的过程中，我们运用了所学的知识，对心电信号进行了分析和处理。在这个分析问题、处理问题的探索过程中，培养了我们的创新意识，增强了我们的动手实践能力，有助于加强团队合作的精神，我们学到了解决实际问题的思路并掌握了MATLAB软件编程的方法。同时本项目激发了我们的学习兴趣，提高了我们的实践能力。我们学到了一个道理，在今后的生活当中，只要用心去做，问题，将不是问题。参考文献：1 董承琅 陶寿淇 陈灏珠. 实用心脏病学 上海科技出版社 19942 周辉 数字信号处理基础及MATLAB实现 中国林业出版社 20053 肖伟 刘忠 MATLAB程序设计与应用 清华大学出版社 20054 姚天任 江太辉 数字信号处理第2版武汉理工大学出版社 20005 黄文梅 熊佳林 杨勇 信号分析与处理MATALB语言及应用 国防科技大学出版社，2000