基于耳机接口的动态心电监护——心电信号自动分析

在手机端采集到心电信号后需要进行心电信号的分析来帮助医生更好的做出判断

心电自动分析主要是指：计算机自动分析和诊断，主要依靠心电信号特征点标定和特征参数测量，主要包括特征点的标定、间期测量、幅度测量以及波形形状描述等。典型的心电信号心动周期包含的特征点主要有：P波起始点和终止点、QRS波群起始点个终止点，T波起始点和终止点等。

心电信号特征点标定后就可以进行各种指标的测量了，首先是间期或者时限测量，下表是临床常用的基本时限的参考指标。

因为是实时的心电监护，所以要求算法的实时性要强，鉴于我们比赛deadline所以这里用的都是比较简单的算法

心率计算与R波定位

R-R间隔

根据以下公式

按照这样的思路，如果知道每秒钟R-R间隔就可以得到每分钟的心率，这样就可以做到每秒钟刷新一次。

在该项目中的应用，8s获取到8000个数据，也就是1s钟有1000个数据，可以看成是1ms一个数据。从数据首开始检测R波，计算相邻两个R波的间隔，就可以得到这一秒的心率值。 在具体应用中，用相邻两个（或两个以上的）R波下标求出奇数下标的平均心率，再用相邻两个的奇数下标的心率求出偶数下标的心率，相当于进行了两次滑动平均求值。

目前R波检测的方法有很多,其中最简单地就是阈值判断法[12]，主要的算法如下：

1. 计算数据的导数，计算公式：DIFF(i)=f(i+1)-f(i-1)+2*f(i+2)-2*f(i-2)（差分算子中的一种）

2. 设定阈值为导数最大值的0.375倍（经验值）

判断大于阈值的点且前后导数相乘小于零的点即为R波的顶点

3. 计算相邻两个R波的间隔，换算为心率

用matlab编程，实现以上算法，定位出Ｒ波，结果如下（代码见附录1）：

这是8s的数据应该产生8个心率，但是从图上可以看出有9个R-R间隔，所以产生的心率是9个，结果如下：

表格4  实验结果2心率值

时间s	1	2	3	4	5	6	7	8	9
心率	74	74	76	76	76	73	92	60	71

可以与图像对比一下，发现刚开始的6个R-R间隔较均匀的，所以得到的心率值差距不大，但是第七个R-R间隔很小，所以产生很大的心率值，而第八个R-R间隔较大，所以产生心率值很小，上述说明这个算法的结果是符合人为判断的。

    在实际应用中应该将得到的结果四舍五入得到整数，下面验证一下，肉眼可观察到8s内有10个QRS波，那么可以换算得到1min内的QRS波数=10/8*60=75,然后计算我们的结果的平均值为这1min内的QRS波数=74.125，所以可以看出，两者的差距很小，从而验证了该算法的准确度和可行性。

 QRS起止点定位算法

在上一小节中通过改进的微分阈值法准确识别出R波后，接着要检测的是QRS波的起始点与终止点。对于起止点的检测，这里采用的是局域变换法，用心电信号减去一段斜率不变的直线，以将心电信号中斜率变化最大的部分凸现出来^[13]，然后选取绝对值最大的点为所检测的点。如图所示，其中y(n)如下式计算得出。图 为局域变换法在QRS起止点检测中的应用。

QRS波群起始位置检测的具体实现步骤：

（1）  定位R波后，在其波前80ms（M点）到R波之间（M，Q1）求绝对值最大的点Q1；

（2）  然后在区间（M,Q1）上求绝对值最大的点Q2；

（3）  计算Q2与R的间期，若间期小于QR时限最大值QR_timelimit，并且 x(Q2)>y(Q2)，则R波前有负向波，以Q2为QRS波的起点，否则，R波前无负向波，以Q1为QRS波起点。

（4）  在QRS波起点Q1或Q2附近，小范围内搜索绝对值最小的点作为QRS波群修正后的起点。

采用同样的方法对QRS波终点检测，只是M点有所不同，其定义在R波波峰之后的110ms处，RS波时限RS_timelimit为80ms。

QRS波起止点检测结果如下图所示：

Matlab代码实现心率计算
clc
clear
data=load('NewECG8K.dat');
%%
fs=1000;
samptime=8;
samples=fs*samptime;
threshFactor=0.375;
%%
differ=zeros(1,samples);
for i=3:1:samples-2 %求导数
    differ(i)=data(i+1)-data(i-1)+2*data(i+2)-2*data(i-2);
end   
thresh=max(differ(:))*threshFactor
k=zeros(1,samptime*2);
j=1;
for i=1:samples %找出data中大于阈值且前后导数相乘为0的点，即极大值点也就是R波
    if  (data(i)>thresh)&&(differ(i)*differ(i+1)<0)
        k(j)=i;
        j=j+1;
    end
end
DR_R=zeros(1,samptime*2);
for i=1:samptime+1  
    DR_R(i)=60*fs/(k(i+1)-k(i));%计算心率
end