ECG ×AI: 机器/深度学习的ECG应用入门（4）

Aiwiscal

传统机器学习:特征工程+分类器
1.引言
经过前面的工作，我们已经解决了数据来源和识别对象问题。那么接下来，我们就要进行机器/深度学习算法的应用了。由于本人写这些博文的目的不在于讲解机器/深度学习理论，所以，涉及到这部分的内容基本不会做太深的原理探讨，不懂的请自行百度，网上的资料有的是。
这一部分，我们先从传统机器学习框架开始，几乎所有利用传统机器学习算法进行分类的问题都遵循“特征工程+分类器”的思路，即先人工构造若干能体现样本特性的特征，然后后端选择分类器，例如支持向量机（SVM），k最近邻分类器（k-NN），决策树，贝叶斯分类器等等单个或集成的算法。然后，再不断进行优化，包括特征的选优和分类器参数的调整，达到理想的效果。一般说来，后端的分类器算法已经成型，想实质改进它们需要很深的理论基础，在应用领域很少对它们进行大的改动。所以，我们已经清楚了，若要采用传统机器学习算法框架，设计好特征是关键。那么对于我们的ECG信号来说，又该怎样设计特征呢？

2.特征设计与分类
1）对于传统机器学习，不设计特征行不行？
我们之前已经提取了心拍，每个长度为250个采样点，那么能不能直接把这些采样点当作特征呢？
乍一想，这么做好像也没什么问题，还能省下特征提取的步骤。但其实这么做确实不是个好的选择，如果我们不做任何处理，那么特征维度就是特征点个数，为250维。对于这么多的特征，机器学习算法往往都会遭遇“维数灾难”问题（机器学习领域经典理论，可自行查阅相关资料），从而达不到有效学习的目的。况且，这些特征点的冗余度很大，除去信号的快变区域，相邻的几个点表征的信息几乎没有什么区别，这样做无疑也会加重分类器的负担。因此，对于ECG信号，设计有效的特征对于后端的分类来说，是很有必要的。

2）常用ECG特征
（1）ECG形态特征
通过查阅一些ECG领域的医用手册，我们可以知道，医生在通过ECG诊断心血管疾病时，其实关注的是各个波形的变化情况，例如，当QRS波变大变宽时，可能发生了室性早搏；ST段抬高时，可能发生了心肌梗死。这样，通过最直观的波形变化，结合医生的经验，可以进行疾病的诊断。所以，这就引出来了第一类常用的特征：形态特征。常用的形态特征有：
① P 波振幅 ②QRS波振幅 ③T 波振幅 ④PR间期 ⑤QRS间期 ⑥QT间期 ⑦ST间期 ⑧PR段水平 ⑨ST段水平 ⑩RR间期，如下图所示：


以上这些都是可以直接从ECG信号中提取到的“一级”特征，由这些特征，还可以组合为各种丰富的“二级”特征，例如例如QRS波面积（QRS间期×QRS振幅）等。对于形态特征，它的优点在于直观，可解释性强，但是缺点也很明显。那就是需要对心拍的各个波段进行精细的定位，而做到这一点是相当困难的，前面已经说过，目前除了QRS波定位检测算法已经足够成熟可靠外，其他波段的定位算法可靠性都不高。如果不能对这些波段进行精准定位，那对于后端分类器的影响是很大的。因此，在近几年的文献中，单纯使用形态特征的论文已经不多见了。
（2）ECG变换系数特征
上面直接提取形态特征的方法虽然直接，但存在严重问题。那能不能通过一些间接的方式得到有效的，较少的，能够有效表示心拍的特征呢？目前大多数采用传统机器学习框架的论文都采用了这样一种方法：使用一些数学变换处理ECG，得到较少的系数，用这些系数来表征心拍。最常见的就是小波变换。利用小波变换能提取多尺度特征的特性，得到有效的小波系数，来表征心拍：



小波变换的基础理论相当复杂，涉及到了大量难懂的数学概念和运算。不过从应用角度来说，我们可以把小波变化看作一个个级联的低通和高通滤波器，然后还有下采样操作。上图中，G和H分别表示不同的高通和低通滤波器，后接2倍下采样，a和d表示了不同尺度下的近似系数和细节系数，这些系数是跟ECG信号的内在特点相关的，通过直接使用或进一步处理这些系数，我们可以得到丰富的特征用于分类。
另外，还有部分文献采用了多项式拟合的思路。即对一个心拍进行多项式拟合，得到各阶的系数作为特征。总体来说，这些使用ECG变换系数的特征，避免了对一些特征点的直接定位，从提取难度上来说小了很多。但是，付出的代价是特征的含义不再直观，这使得对于这些系数特征进行选择时变得困难，往往比较盲目。
以上两种常用特征各有利弊，不过从目前的趋势来看，提取难度更小的第二种特征更受青睐。而且这些特征基于数学变换，也具有更多的理论价值和研究价值。这可能是目前此类论文较多的原因。
3）常用分类器
这一点其实不用多说了。基本上所有经典的分类器都可以用，只要能够得出好的效果。而目前各种分类器算法的工具包支持也很完善了，如果从应用角度来看，我们可以直接套用高效率的工具包，而不用自己重造轮子。废话不多说，下面直接开始示例。

3.示例：心律失常的SVM识别
这一部分中，我们采用大名鼎鼎的支持向量机SVM作为后端的分类器。仍然是在Matlab环境下进行实验，相应的工具包使用了著名的libsvm。由于libsvm不是Matlab自带的工具包，其安装过程略微复杂，需要外部C++编译器的辅助。但是这里我直接将编译好的工具包共享出来，放在本人的github上:libsvm-3.21.zip。解压后，这个工具包是一个文件夹（libsvm-3.21），想在Matlab中使用，需要在Matlab中设置路径（Set path），将这个文件夹中名为“matlab”子文件夹添加到Matlab可检索的目录中，就可以调用了。建议将libsvm文件夹放在Matlab安装目录下的toolbox文件夹中。另外，由于我是在64位平台上编译的，因此这个工具包只能在64位版本上运行。如果是32位，可以参考书籍《MATLAB神经网络43个案例分析》中有关SVM的章节，里面有比较详细的安装过程。
在进行以下步骤之前，需要使用load函数把我们之前存储的心拍数据（.mat）导入到工作空间中.这时我们可以发现,正常类的心拍远远多于其他三类非正常心拍, 即存在数据分布不平衡问题.这里我们为了减少数据对之后分类器训练和测试的影响,每一类只选取5000个进行实验,也就是数据集规模为20000.以下所有步骤的Matlab代码可在本人github上下载到：Classification-SVM.m,可对照代码理解并上手.
1)特征提取
。这里我们采用小波变换的系数作为我们的特征. 原因除了在于小波变换的多分辨率属性外,还在于小波变换中还包含了下采样操作,可以理解为内嵌了降维的过程.从实现角度来讲,Matlab平台中提供了相应的内建函数供我们使用,大大降低了我们的实现难度.
具体的,我们对每个心拍进行5阶的小波分解,小波函数利用db6小波,使用Matlab中的wavedec函数（函数具体功能建议查看Matlab帮助文档）可以轻松完成:

        [C,L]=wavedec(sig,5,'db6');
        
其中，C包含了各阶小波变换后的系数，sig在这里表示我们所提取的250点心拍。经过5阶小波分解和2倍下采样后，我们取小波变换系数中原信号的“近似”系数，即5阶分解后的a系数，根据L中的信息可以知道是前25个点。提取“近似”系数的动机是这些系数去除了一些细节，对更一般性的规律体现更好。这一点在疾病的诊断中非常关键，因为不同的病人或相同的病人在不同时间，由于身体状况的改变可能造成ECG的细节变化。这些细节变化是带有特殊性的，往往不可以作为疾病的一般性特征。注重一般性规律，是有利于防止机器学习中的“过拟合”问题的。因此，我们这里选取这25个系数为代表每个心拍的特征。
从这里看起来，提取特征貌似也没什么难度。但是，我们这里是最简单的示例，不一定保证这些特征就是最优的。真正好的特征需要非常多的经验和实验，同样以小波变换系数为例，分解阶数，特征个数，小波函数等等的选择，都会影响特征的质量。并且这只是直接提取的系数，在此基础上还可以做各种处理，形成二级特征。而这些因素在真正的研究中都需要认真考虑和探究。
2）数据集选取与划分
随机划分训练集和测试集,各占一半,即训练集和测试集各有10000个样本. 随机选取可通过Matlab内建randperm函数随机打乱样本索引,然后截取前10000个索引对应的样本作为训练集,剩余的作为测试集来实现.
3）特征归一化
在机器学习任务中，特征归一化不是必须的步骤。这里为了加快SVM收敛，进行了特征归一化。这里的归一化指的是同一特征在不同样本上取值的归一化,借助Matlab内建归一化函数mapminmax实现. 不过这里要注意的是mapminmax函数对于归一化对象,即样本特征矩阵的行列含义的要求,必要时要对矩阵进行转置. 实际操作时,先使用mapminmax函数的正常模式对训练集特征归一化到0,1之间,得到归一化后的训练集和归一化信息; 然后使用mapminmax的apply模式,将训练集归一化得到的归一化信息应用至测试集,完成测试集的归一化:

[train_x,ps]=mapminmax(train_x',0,1); %利用mapminmax内建函数特征归一化到0，1之间；
test_x=mapminmax('apply',test_x',ps);
4)模型训练与测试
调用libsvmtrain函数和libsvmpredict函数训练和测试模型,注意两个函数对于数据顺序和格式的要求:
model=libsvmtrain(train_y,train_x,'-c 2 -g 1'); %模型训练；
[ptest,~,~]=libsvmpredict(test_y,test_x,model); %模型预测；
libsvm训练的默认核函数为RBF核函数,需要人为设定两个超参数c和g. 这里我们设定为2和1. 注意,不同取值的c和g可能会导致差异比较大的结果.如果想要取得更好的效果,必须进行"调参",减轻欠拟合与过拟合问题. 而这一部分也是所有机器/深度学习问题中最依赖经验和最耗费时间的. 常用的调参方法除了手动调整,还有网格搜索,启发式寻优等方法.

matlab的旋律

5)结果统计与展示
这一部分主要是根据分类器的实际输出结果(上一步中libsvmpredict输出的ptest向量),与期望结果(标签)对比, 统计各类别的准确率等.详细的代码可参见本人github. 例如,根据以上的特征提取和参数选择, 我们可以得到如下关于准确率的信息(由于这里每次会选取不同的训练和测试样本,因此每次运行结果会有所不同):

在本次实验中,我们的SVM模型总体预测准确率为96.69%, 其中四类目标类型的准确率分别为正常(N):99.68%, 室性早搏(V): 90.90%, 右束支阻滞（R）:97.58%, 左束支阻滞（L）:98.49%. 代码中其实还给出了一种常见的衡量多分类的方法：混淆矩阵。具体定义参见：https://baike.baidu.com/item/混淆矩阵/10087822?fr=aladdin 工作空间中打开名为“Conf_Mat”的矩阵，可见：

类别1，2，3，4分别表示N，V，R，L。从混淆矩阵中，我们可以看出预测结果的分布情况。对角线上是各类别正确预测的个数，其余均为错误预测个数。例如，第2行第4列，有128个早搏心拍被错误预测为左束支阻滞，是所有错误中最多的。可以猜测，早搏心拍与左束支阻滞心拍可能在表现上更相似一些；而第4行第3列错误预测数为0，即没有左束支阻滞心拍被错误预测为右束支阻滞，说明可能两者差异更大一些。通过分析混淆矩阵，可以对分类器的性能和数据分布做更精细的分析，根据经验采取相应措施优化。

综上，我们以经典的SVM为例，完成了一个传统机器学习框架的最基本ECG应用。从特征提取到模型训练与评估，最终得到了一个看起来还不错的效果（这里还存在一些问题，之后的内容中会有所提及）。以上强烈建议结合代码理解和上手。当然，这里面还大有文章可做，包括怎么样提取质量更高的特征，怎么样选取更优的超参数，以及选取更适合的分类器等等。但万变不离其宗，总体上，在传统机器学习框架中，我们主要的精力应该还是集中于“好”特征和参数的选择。
4.小结          这一部分将传统机器学习框架的ECG分类应用大致走了一遍流程。看似简单有效，但其实特征的选取和选择是一件很困难的事情。想要设计出好的特征，必须具有相当的经验和大量的实验。下面推荐了几篇文献，可以更进一步体会。

** 相关代码文件下载：https://github.com/Aiwiscal/ECG-ML-DL-Algorithm-Matlab
** 推荐文献：[1] Yeh Y C, Wang W J, Chiou C W. Feature selection algorithm for ECG signals using range-overlaps method[J]. Expert Systems with Applications, 2010, 37(4): 3499-3512.                 
              [2]   Raj S, Ray K C. ECG signal analysis using DCT-based DOST and PSO optimized SVM[J]. IEEE Transactions on Instrumentation and Measurement, 2017, 66(3): 470-478.                       
              [3]  Liu T, Si Y, Wen D, et al. Dictionary learning for VQ feature extraction in ECG beats classification[J]. Expert Systems with Applications, 2016, 53: 129-137.
                        [4]  Kamath C. ECG beat classification using features extracted from Teager energy functions in time and frequency domains[J]. IET signal processing, 2011, 5(6): 575-581.
                        [5] Rad A B, Eftestøl T, Engan K, et al. ECG-based classification of resuscitation cardiac rhythms for retrospective data analysis[J]. IEEE Transactions on Biomedical Engineering, 2017, 64(10): 2411-2418.