心电R波检测之差分奇点法和局部极值排除法

心电R波检测之差分奇点法和局部极值排除法

心电图(Electrocardiogram，ECG)中R波检测是心律失常分析、心率变异性分析、心率趋势分析的必要基础。心电R波检测的算法众多，有很多文献可以参考。实测的心电图有各种各异的表现，在基本算法原理的基础上需要针对具体信号情况和应用场景对算法进行优化。这里介绍两种简单直接的方法：差分奇点法和局部极值排除法。把心电信号中R波的出现看作是信号中的奇异点，是小概率事件，据此设计一种简单的心电R波检测方法，称之为差分奇点法，该方法计算简单，速度快，具有自适应跟踪能力，适合于实时心率计算场合。把心电R波顶点看成是局部极大值，并依据R波的特征参数来确认该极大值是否为R波，这就是局部极值排除法的基本出发点。依据局部极值排除法，设计了实时心电R波检测的算法流程。

一、差分奇点法：一种心电R波快速自适应检测算法

心率检测实际上是心电信号R波间期的检测，R波是心电信号中最尖的波峰，利用这一点就可以实现R-R间期的检测。设心电信号为x(t)，为了消除基线漂移，先对信号进行差分：d(t) = x(t+1)-x(t)，d(t)在R波波峰附近的值之绝对值较大，假设d(t)服从正态分布，则d(t)在R波附近的大数值在整个序列中出现的概率较小，属于奇异数值，由于d(t)是差分后的值，可以认为它的均值为零，其方差var(t)可用递推公式计算：

    var(t) = var(t-1) + [d(t)*d(t) - var(t-1)]/t

按上式计算的方差随着时间的推移趋于常数，为了让其具有一定的自适应跟踪能力，可改用下式计算：

    var(t) = var(t-1) + A*(d(t)*d(t) - var(t-1) )

其中A是学习率，且0<A<1，A越大对过去数据遗忘越快，学习能力越强，var(t)随时间的波动也越大；A越小对过去数据遗忘越慢，学习能力越差，var(t)随时间的波动越小。如果信号的采样速率为Fs，则我们可取A=1/(Fs*5)。

在R波波峰附近d(t)的绝对值比较大，而出现的概率相对于整个信号长度来说却很小，假设d(t)序列的分布是正态的，则这些大值点属于奇异点。假设R波是波峰朝上的，令s(t)=[var(t)]1/2，近似为t时刻信号的标准差，则如果d(t)> B*s(t) (B>1是判决门限, 例如B=2.5)，则可认为R波波峰即将出现，记下此时刻为T1，继续向前搜索，如果再次遇到d(t)>B*s(t)并不改变T1之值，直到检测到- B*s(t)<d(t)<0记下T2为止，在T1和T2之间有一个最大的值对应的点即为R波波峰的位置。

下图是原始心电及心电差分波形的例子。

在检测心率时，实际上并不要求检测到精确的R波波峰位置，而只需要R-R间隔，这样我们只需要得到T1即可，因而可大大减少算法的复杂度。根据以上原理，可将实时心率检测算法描述如下：

初始化全局或静态变量
oldx=0; //保存前一点样值
var=0; //保存方差
Hrmax=200; //心率最大可能值
Hrmin=35; //心率最小可能值
LearnR=1.0/(5*Fs); //学习率，Fs是心电信号的采样速率
Th=2.5; //判决门限
T=0; //间期计数器

心率检测函数实现。每输入一个采样点之值x，进行一点运算，并返回计算结果。返回0表明正常迭代运算，未检测到R波；返回值大于零，则表明检测到了R波，并返回心率；返回小于零，则表明检测到了非正常的R波，因为此时的心率小于Hrmin或大于Hrmax。

float DetectHR(float x)
{
    float d;
    float T1,T2,HR;

    T1=60.0/Hrmax*Fs;
    T2=60.0/Hrmin*Fs;
    d=x-oldx;
    oldx=x;
    var+=LearnR*(d*d-var);
    if(d>Th*sqrt(var) && T>T1){
        HR=60.0/(T/Fs);
        If(HR>Hrmax || HR<Hrmin) return –1.0;
        return HR;
        T=0;
    }
    T++;
    return 0.0;
}

由上述算法可知，只有在经过了几秒钟的学习后，方可得到正确的检测结果。为了让该算法在面对实际信号的复杂情况具有适用性，需要对信号中存在的高频干扰（如工频干扰和肌电干扰）和基线漂移做预先处理，即平滑和去基线漂移，如平滑窗口长度为Fs/50，去基线窗口长度为Fs。需要依据信号的具体情况，设置合适的学习率和门限。该算法简单、速度快、逐点计算、具有一定的幅度变化跟踪能力，适合于一边采集一边计算的实时心率计算场合。

二、心电R波检测的局部极值排除法

心电图R波检测局部极值排除法，是在信号中建立一个滑动窗口，窗口长度大致为QRS波宽度，窗口每滑动一点，判断窗口的中点位置的信号值是否为窗口内信号的最大值，如果是最大值，则计算斜率、幅度、间期、峰宽等参数来判断该最大值是否R波。详细描述可参考文献：成奇明 等. 生理信号波峰检测的局部极值排除方法[J]. 医疗卫生装备，2008, 29(9): 256-257.

三、心电R波检测流程

1、信号预处理

（1）消除基线

自适应基线检测，并祛除基线漂移。

（2）消除50Hz工频干扰

对于采样率为200Hz的信号，可采用4点滑动平均来抑制50Hz工频干扰。

2、干扰分类检测

判断一小段信号是否是干扰，段的持续时间可以设置成1秒左右。

（1）直线干扰，信号饱和；

x(i)-x(i-1)=0 持续一定时间长度(直线持续限制lineTH)，如0.1s。

例1

例2

例3

例4

（2）突变干扰，电极脱落；

|x(i)-x(i-1)|>th(突变门限deltaTH)。正常信号变化是连续的，不可能超出该门限。

例1

例2

（3）低幅度宽带噪声干扰；

max-min<a(低幅度门限lowTH)，过零率大于门限(过零率门限zrTH)。

（4）趋势项干扰；

信号趋势性检验。

（5）抖动干扰，运动引起干扰；

标准差大于门限（标准差上限maxStdTH），并且幅度差也大于门限(幅度上限maxTH)。

（6）低幅度信号。

幅度差小于门限（幅度下限minTH），或者标准差小于门限（标准差下限minStdTH）。

建立如下门限：

幅度下限minTH，幅度上限maxTH，低幅度门限lowTH，突变门限deltaTH，标准差上限maxStdTH，标准差下限minStdTH，直线持续限制lineTH，过零率门限zrTH。

（7）R波向下

3、心电R波检测

（1）设置数据缓存，长度大于等于3秒。m≥3*Fs。B[m]，m是奇数。

（2）关注点：t0=(m-1)/2，是缓存的中间点。

（3）假设R波波峰是波信号波形在关注点t0局部范围的顶点，也是最大值，局部范围长度为L，设置L为0.1秒，L=0.1*Fs。

（4）判断t0点是否局部最大值。

（5）判断t0点是否局部波峰顶点。

       x(t0-L+1)≤x0(t0-L+2)≤…≤x(t0-1)≤x(t0)≤x(t0+1)≤…≤x(t0+L-2)≤x(t0+L-1)。

（6）如果是顶点，则顶点的幅度需要大于一定的门限，避免小幅度毛刺影响。

       Amin=min{x(t0-L+1)…x(t0)…x(t0+L-1)}

   顶点幅度为：A=x(t0)-Amin, A>ath(门限)。

（7）RR间期计算

t0当前位置顶点，右边两个顶点位置为p1,p2，左边两个顶点位置为n1,n2，那么

右间期：r1=p1-t0, r2=p2-p1, r3=(p2-t0)/2，

左间期：l1=t0-n1, l2=n1-n2, l3=(t0-n2)/2

在实时算法中，右边顶点尚未求出来，不能使用。

（8）间期异常判别

根据心率变异性范围，确定当前顶点是否可疑。

|L1-L2|<thms, |L3-L1|<thms, |L3-L2|<thms

minRR<L1<maxRR, minRR<L2<maxRR

thms设置为：150ms

如果L1<minRR，或者|L1-L2|>thms，|L1-L3|>thms则认为间期异常，考虑是否抛弃该顶点的选择。

（9）如果幅度也变异很大，可考虑抛弃该顶点。

|x(t0)-x(n1)|>thv

|x(t0)-x(n2)|>thv

thv是幅度变异门限

（10）信号特征参数自适应学习

  局部幅度，判断为峰值点附近的幅度：

   A = (1-a)*A + a*(Xmax-Xmin)

   a是学习率，例如a=0.01

  平均心动周期：RR = (1-a)*RR + a*L1

  平均左、右斜率：SlopeL, SlopeR

  平均间期变异：RRV = (1-a)*RRV + a*|L1-L2|

  平均幅度（基线）：base = (1-a)*base+a*x

 4、实时心率计算流程

（1）初始化，包括设置缓存大小，分配内存空间。（2）输入信号，在信号输入之前，对信号进行滤波预处理。（3）逐点更新缓存。（4）判断缓存中固定点附近是否干扰信号。（5）判断缓存中固定点是否是波峰。（6）判断心动间期。

四、心电R波检测实例

1、差分奇点法

多通道信号分析软件菜单操作：《分析》→《波峰检测》→《心跳检测》，出现如下对话框，选择方法4：差分奇点法，输入心率范围，对指定的通道信号进行R波检测。

下图是差分奇点法心电图R波检测结果例子。

2、局部极值排除法

多通道信号分析软件菜单操作：《分析》→《波峰检测》→《心跳检测》，出现对话框，选择方法0：直接极值排除法，对指定的通道信号进行R波检测。

下图是局部极值排除法检测结果例子。

参考文献：
[1] 赵文哲，方滨，沈毅，王普. 心电信号中R波检测方法的比较研究[J]. 生物医学工程学杂志，2009，26(1): 55-58.
[2] 成奇明 等. 生理信号波峰检测的局部极值排除方法[J]. 医疗卫生装备，2008, 29(9): 256-257.

联系作者：chengbowork@163.com

相关文章：

多通道信号分析软件系统
多通信号分析软件(MUSIA)功能介绍
生物医学信号处理与分析软件系统设计
曲线拟合软件
离散小波变换用于信号滤波
多维特征参数机器学习算法
多维特征参数机器学习软件
Kohonen自组织特征映射神经网络(环形和球面形网络)
矩阵的三维图形显示软件
图片浏览软件工具
主成分分析(K-L变换)与信号的分解与合成(滤波)
信号的样本熵序列计算
信号的双谱分析
信号的经验模态分解(EMD)
希尔伯特(Hilbert)变换信号瞬时频率计算
给定概率分布的随机变量仿真
信号的特征参数计算
数据分布点纹图
多通道信号数据压缩存储
图像处理基本算法及软件简介
信号功率谱估计
特殊函数计算器
滤波器设计、自适应滤波、匹配滤波方法
离散正交变换及其应用
短时傅里叶变换时变功率谱分析
数值矩阵的图形表示