ECG：机器学习之预处理算法与R波定位算法（附部分代码）

ECG：机器学习之预处理算法与R波定位算法（附部分代码）

目前网上有关ECG信号处理的文章越来越多，说明这一块做的人也越来越多，但是网上的资源却很有限，而且很多资源都有一定的缺陷且不全。本文会对网上的资源进行整合，然后结合自己处理ECG信号的经验进行一些补充。这篇文章主要是对ECG信号的预处理算法和R波定位算法进行讨论，如有不足，望指正。

一、ECG信号简介

1. 心电信号是人体心脏的心脏细胞的细胞膜产生的电势差，心房和心室肌在静止的间歇中，由于细胞内外离子(包括K＋，Na＋，Ca2＋，cl-等)浓度差别很大，处于“极化状态”。一旦受到自搏细胞传来的激动，这极化状态便暂时瓦解，在心电图上称为“除极”(有少数学者称为“去极”)，由此产生心电活动。心房肌的除极在心电图上表现为P波，心室肌的除极表现为QRS波群。当然在一次除极后，心肌又会恢复原来的极化状态，此过程称为“复极”。复极过程远较除极缓慢，电活动所产生的振幅也较低。心房的复极在P—R段上，一般很不明显(唯有在右心房扩大时，P—R段轻度压低)。心室肌复极则表现为心电图上的ST段及T波。
2. 完整的一个PQRST波形如下：

（1）P波：
正常心脏的电激动从窦房结开始。由于窦房结位于右心房与上腔静脉的交界处，所以窦房结的激动首先传导到右心房，通过房间束传到左心房，形成心电图上的P波。P波代表了心房的激动，前半部代表右心房激动，后半部代表左心房的激动。P波时限为0.12秒，高度为0.25mv。当心房扩大，两房间传导出现异常时，P波可表现为高尖或双峰的P波。

（2）PR间期：
PR间期代表由窦房结产生的兴奋经由心房、房室交界和房室束到达心室并引起心室肌开始兴奋所需要的时间，故也称为房室传导时间。正常PR间期在0.12～0.20秒。当心房到心室的传导出现阻滞，则表现为PR间期的延长或P波之后心室波消失。

（3）QRS波群：
激动向下经希氏束、左右束枝同步激动左右心室形成QRS波群。QRS波群代表了心室的除极，激动时限小于0.11秒。当出现心脏左右束枝的传导阻滞、心室扩大或肥厚等情况时，QRS波群出现增宽、变形和时限延长。

（4）J点：
QRS波结束，ST段开始的交点。代表心室肌细胞全部除极完毕。

（5）ST段：
心室肌全部除极完成，复极尚未开始的一段时间。此时各部位的心室肌都处于除极状态，细胞之间并没有电位差。因此正常情况下ST段应处于等电位线上。当某部位的心肌出现缺血或坏死的表现，心室在除极完毕后仍存在电位差，此时表现为心电图上ST段发生偏移。

（6）T波：
T波代表了心室的复极。在QRS波主波向上的导联，T波应与QRS主波方向相同。心电图上T波的改变受多种因素的影响。例如心肌缺血时可表现为T波低平倒置。T波的高耸可见于高血钾、急性心肌梗死的超急期等。

（7）U波：
某些导联上T波之后可见U波，目前认为与心室的复极有关。

（8）QT间期：
代表了心室从除极到复极的时间。正常QT间期为0.44秒。由于QT间期受心率的影响，因此引入了矫正的QT间期（QTC）的概念。其中一种计算方法为QTc=QT/√RR。QT间期的延长往往与恶性心律失常的发生相关。

3. ECG信号的预处理
人体的心电信号是一种非平稳、非线性、随机性比较强的微弱生理信号,幅值约为毫伏(mV)级,频率在0.05-100Hz之间。在我们获取到的原始ECG信号中都会有噪声干扰，本节将介绍噪声来源以及消除噪声的方法。

（1）ECG信号噪声种类及来源：

基线漂移：
一般是由呼吸和电极滑动变化所异致的,频率一般低于1Hz，其表现为变化比较缓慢的类正弦曲线，对心电波形中的ST段识别影响较大。基线漂移的频率很低，其范围为0.05Hz至几Hz，主要分量在0.1Hz左右，而心电信号的P波、T波及ST段的频率也很低，其范围为0.5Hz至10Hz，两者的频谱非常接近，在消除噪声的同时，不可避免地对心电信号成分造成一定的损失。（由于ST段的频率也很低，无法使用低通滤波器去除基线漂移。）

肌电干扰：
它是由人体肌肉颤抖产生不规则的高频电分扰所导致的,其频率范围很宽，一般在10-1000Hz之间,严重的肌电干扰信号频率在10～300Hz之间，其频谱特性接近于瞬时发生的高斯零均值带限白噪声。

工频干扰：
主要来源于工频电源以及器件周围环境中的传输线辐射出的电磁场，频率为50Hz或60Hz,在ECG上出现为周期性的细小波纹，其频率成分主要为工频频率及其谐波 。（可以通过硬件滤波和软件滤波实现，使用50Hz和60Hz的陷波器。）

（2）小波变换去除基线漂移

对于一般肌电干扰和工频干扰一般通过滤波器可以去除掉，麻烦一点的就是基线漂移，在这里将介绍小波变换去除基线漂移。

1）小波变换的原理：
小波变换的公式为（觉得复杂可以不看）：

小波变换有两个变量，一个是伸缩 a、一个是平移τ，伸缩反映的是频率情况，平移反映的是时间，假定用小波不同尺度（伸缩）下都与源信号进行运算，那最终能得到的是任意位置（时刻）的频率情况，即时频图，这样我们就可以知道每个位置（时刻）的频率成分了，这就是非常著名的小波变换了，小波变换就像信号处理里面的一个显微镜，能很好的分析出信号在时间、频率上的信息。（至于算法的实现部分这里不做详细说明，直接说怎么用）

2）小波变换应用到ECG信号中

在网上找到一个ECG信号采样率为250HZ时使用小波变换对其进行分解时的图，图中包含了不同分解层中含有的不同频率的信号分量：

在MATLAB中，可以利用函数wavedec来进行小波分解，利用函数waverec进行小波重构，具体用法可以在MATLAB帮助文档中查看。

3）实例：
这里我们使用一个三分钟左右采样频率为512HZ的ECG信号进行演示。

使用小波分解与重构之前的信号如下图所示：

使用小波分解与重构之后的ECG信号如下图所示：

使用的MATLAB代码如下：

[c,l]=wavedec(x,5,'dmey');%分解
 c(1:l(1))=0;
 x=waverec(c,l,'dmey'); %重构

根据实际的信号，发现5层分解的效果就很好了。对于预处理好的ECG信号，下一步就可以进行R波的定位了。

4. 常用的R波定位算法

常用的R波定位算法有两种，一种是Pan-Tompkins算法还有一种是自适应滑动窗算法。

（1）Pan-Tompkins算法简介：

Pan-Tompkins法检测R波峰值的具体步骤如下：
1）将信号通过给定的滤波器；
2）对滤波后的信号求一阶导数；
3）对求导之后的信号进行平方运算；
4）将信号通过滑动窗口进行积分；
5）使用阈值法检测经过处理之后的R波峰值。

算法流程图如下：
在实际的处理中没有用到过，所以理解的不够透彻，代码也有，代码链接：R波定位算法

（2）自适应滑动窗口定位R波

1）首先，运用离散小波变换将获得的心电数据去除基线漂移干扰；其次，把信号根据时间窗长进行分割，最开始的时间窗长定义为经验窗长，经验 RR 间期长度T1设置为最慢的正常 RR 间期值，长度设置为 N =t×f （一般情况 t 为 1.2秒，运动情况下 t 为 0.8 秒；f 为采样频率，设置为 512Hz）；最后，通过 max()函数计算出时间窗内的最大值点，该点就是 R 波的波峰点。
2）当检测到两个 R 波波峰时，计算出两个相邻 R 波波峰间距离，计算心率值，为下一个时间窗长的设定做准备。顾名思义，自适应时间窗是指时间窗长度的设置是根据心率值的变化而变化。当心率值低于 0.45 时，时间窗长值设置为 1.2倍的 RR 间期序列；当心率值超过 0.5 时，时间窗长值设置为 1.5 倍的 RR 间期序列。
3）当检测到三个 R 波波峰时，计算出前两个 RR 间期序列的相邻差，设置RR 间期序列相邻差的经验阈值 T2 为 0.35。当 RR 间期的相邻差超出规定的阈值范围时，说明前面的 R 波峰出现错误定位，需要重新设置时间窗长，随后跳转到步骤（1），对 R 波峰重新进行检测。

R 波峰的相邻差决定了时间窗长值是不断改变的过程。因此当 R 波峰出现错误定位时，并不会影响后续的 R 波峰定位。按照上述描述的 R 波峰定位的流程，依次往后定位 R 波峰，直至心电数据段中的所有 R波波峰被全部定位。自适应滑动时间窗 R 波波峰定位的算法流程图如下图所示。
这个算法参考的是一篇硕士学位论文，题目为《运动干扰焦虑实时检测的自主神经机制研究》，实际应用上这个算法定位还是很准确的，代码就不放了，需要的在评论区留下联系方式。

定位R波之后就是提取RR间期了，后续操作会在后面的文章进行介绍。