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

基于matlab的操作员心电信号QRS波检测及分析1引言众所周知，疲劳是人的一个自然的生理现象，它的产生与一个人的劳动强度、工作环境、本身的健康和情绪等等有关，是人体的一种自我调节和保护功能。当疲劳被人们感觉时就意味着此刻需要好好休息一下。但是随着现在经济的高速发展与生活节奏的加快，许多工作人员因为种种原因，都是顶着疲劳继续工作，这样就使得人体产生的疲劳成为一些事故的根本导火索。本文要探索的是与人疲劳有关的心电信号，通过对人体心电信号的研究，来找出人的疲劳度量，预防人因为过度疲劳工作而导致事故的发生。本项目通过实验方法采集了执行过程控制任务的操作员的心电信号，通过确定阈值求极值确立R波，然后通过R波确立其前后最近的最低点分别为QS波信号，从而计算出操作员的心电QRS波等时间间隔，最后分析心电R波时间间隔和操作员主观疲劳程度、焦虑程度、努力程度和客观任务负荷的相关性。2心电原理人体内蕴藏着大量丰富的生物电信号，如：脑电信号、心电信号、眼电信号、肌电信号、皮肤电压反应信号等，这些电信号都包含着大量的信息，它们大多都是在各种因素下由复杂生命体发出的不稳定的、幅值微弱的低频生理信号。而我们研究的便是其中的心电信号，心信号较之其他生物电信号更容易检测，便于实时监测。心脏周围的组织和体液都能导电，因此可将人体看成为一个具有长、宽、厚三度空间的容积导体。心脏好比电源，无数心肌细胞动作电位变化的总和可以传导并反映到体表。在体表很多点之间存在着电位差，也有很多点彼此之间无电位差是等电的。心脏在每个心动周期中，由起搏点、心房、心室相继兴奋，伴随着生物电的变化，这些生物电的变化称为心电。现在的科技已经可以通过采取一些有效设备和技术手段捕获人体可用的电信号，并对其进行信号采集与分析，进而有效破译其包含信息，利用这些大量的生命体征信息来对人的生活、工作、医疗等许多方面进行准确而有效的改善和调整。心电信号受到许多因素的影响，如：疲劳、疾病等等。而我们所要探索的便是与人疲劳有关的心电信号，通过对这类心电信号的研究，来找出人的疲劳与心电之间的关系，防止人因为过度疲劳工作而导致事故的发生。3QRS波的检测方法3.1以软件为主的方法实现QRS波的检测以软件为主的方法实现QRS波的检测在有大量数据需要检测的时候就显得方便快捷，极度占优。用软件来实现，就需要确定各波的阈值，当心律失常或QRS波幅度变小，阈值设置过高，会导致漏检产生假阴性(FN，falSenegative)结果。由于固定阈值的这些特点，可以以此来判断病人的心脏是否存在疾病以及隐患。3.2QRS波检测算法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具体实验步骤如下：functiony=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赋值给ySigtemp=y;%将y赋值给SigtempSiglen=length(y);%将y的长度保存为SiglenSigmax=[];%建立空矩阵SigmaxfoRi=1:Siglen-2%令变量i依次从1加40004500500055006000650070007500-500050010001500图1部分心电信号图（2）根据心电信号特征求R波foRi=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中的所有极点依次填入Sigmaxend;thRtemp=SoRt(Sigmax);%将Sigmax内元素排序赋值给thRtempthRlen=length(Sigmax);%将Sigmax长度保存为thRlenthR=0;foRi=(thRlen-899):thRlenthR=thR+thRtemp(i);end;thRmax=thR/900;%最大幅度平均值（依据输入数据而定）zeRotemp=SoRt(y);zeRovalue=0;foRi=1:100zeRovalue=zeRovalue+zeRotemp(i);end;zeRovalue=zeRovalue/1000;%最小幅度平均值，对消幅度（依据输入数据而定）thR=(thRmax-zeRovalue)*0.3;%最大、最小幅度的差值的30%为判别R波的阈值Rvalue=[];%定位R波foRi=1:thRlenifSigmax(i,1)thRmax%寻找大于最大幅度平均值的极点，此处更改为thRmaxRvalue=[Rvalue;Sigmax(i,2)];end;end;lenRvalue=length(Rvalue);40060080010001200140016001800200022002400-1500-1000-500050010001500图2部分心电信号R波图（3）根据已知R波以及心电信号特征求解Q、S波(此处只例举解Q波的程序，求解S波的程序与其类似)：wtSig2=cwt(Sig,8,'mexh');%对数据进行小波变换平滑处理Qvalue=[];%设Q波集合foRi=1:lenRvaluefoRj=Rvalue(i):-1:(Rvalue(i)-30)ifSig(Rvalue(i))0ifwtSig2(j)wtSig2(j-1)&wtSig2(j)wtSig2(j+1)tempQvalue=j-9;%确定检测窗的起点bReak;%正向波，取第一个负极大值end;elSeifwtSig2(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=0diSt=[];foRk=tempQvalue:Rvalue(i)tempdiSt=(abS(a0*k+b0*Sig(k)+c0))/SQRt(a0^2+b0^2);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;30003100320033003400350036003700380039004000-1000-500050010001500图3部分心电信号QRS输出图形（4）得到所有QRS波之后，我们对数据进行了整理，提取出了各类特征值，并以图像形式输出。程序如下:foRi=1:(lenRvalue-1)tR(i)=(Rvalue(i+1)-Rvalue(i))/500;end;%计算R波时间间隔tQ=[];foRi=1:(lenQvalue-1)tQ(i)=(Qvalue(i+1)-Qvalue(i))/500;end;%计算Q波时间间隔tS=[];foRi=1:(lenSvalue-1)tS(i)=(Svalue(i+1)-Svalue(i))/500;end;%计算S波时间间隔tQR=[];foRi=1:lenRvaluetQR(i)=(Rvalue(i)-Qvalue(i))/500;end;%计算同次QR波时间间隔tRS=[];foRi=1:lenSvaluetRS(i)=(Svalue(i)-Rvalue(i))/500;end;%计算同次RS波时间间隔tQS=[];foRi=1:lenSvaluetQS(i)=(Svalue(i)-Qvalue(i))/500;end;%计算同次QS波时间间隔%R指标计算avetR=0;%定义RR时间间隔平均值%定义RR间隔标准差foRi=1:lenRvalue-1avetR=avetR+tR(i);end;01002003004000.70.80.91相邻r波时间间隔曲线01002003004000.70.80.9相邻q波时间间隔曲线01002003004000.70.80.9相邻s波时间间隔曲线01002003004000.0250.030.0350.04每次心跳qr波时间间隔曲线01002003004000.0220.0240.026每次心跳rs波时间间隔曲线01002003004000.040.050.06每次心跳qs波时间间隔曲线图4QRS波特征值（5）在求的QRS波的特征值之后，我们对其进行了进一步的计算得出了相邻的R、Q、S、QS、QR、RS波时间间隔的平均值，方差，标准差及标准差与均值的比值，具体程序如下：avetR=avetR/i;SDtR=0;foRi=1:lenRvalue-1SDtR=SDtR+(tR(i)-avetR)^2;end;vaRtR=SQRt(SDtR);%定义RR时间间隔方差SDR=SDtR/avetR;%定义RR时间间隔标准差除以平均值%Q指标计算avetQ=0;foRi=1:lenQvalue-1avetQ=avetQ+tQ(i);end;avetQ=avetQ/i;SDtQ=0;foRi=1:lenQvalue-1SDtQ=SDtQ+(tQ(i)-avetQ)^2;end;vaRtQ=SQRt(SDtQ);SDQ=SDtQ/avetR;%S指标计算avetS=0;foRi=1:lenSvalue-1avetS=avetS+tS(i);end;avetS=avetS/i;SDtS=0;foRi=1:lenSvalue-1SDtS=SDtS+(tS(i)-avetS)^2;end;vaRtS=SQRt(SDtS);SDS=SDtS/avetS;%QR指标计算avetQR=0;foRi=1:lenRvalue-1avetQR=avetQR+tQR(i);end;avetQR=avetQR/i;SDtQR=0;foRi=1:lenRvalue-1SDtQR=SDtQR+(tQR(i)-avetQR)^2;end;vaRtQR=SQRt(S