数据描述

这个例子使用了从三组人(或三类人)获得的心电图数据:心律失常的人、充血性心力衰竭的人和窦性心律正常的人。该示例使用了来自三个PhysioNet数据库的162个心电图记录:MIT-BIH心律失常数据库[3] [7],MIT-BIH正常窦性节律数据库[3]充血性心力衰竭数据库[1][3]。总共有96个心律失常患者的录音，30个充血性心力衰竭患者的录音，36个窦性心律正常患者的录音。目的是训练分类器区分心律失常(ARR)、充血性心力衰竭(CHF)和正常窦性心律(NSR)。

下载数据

第一步是从。下载数据GitHub库．如需下载数据，请单击代码并选择下载ZIP．保存文件physionet_ECG_data-main.zip在您有写权限的文件夹中。本例的说明假设您已经将文件下载到临时目录(tempdir在MATLAB中)。如果选择从不同的文件夹下载数据，请修改后续关于解压缩和加载数据的说明tempdir．

该文件physionet_ECG_data-main.zip包含

和ECGData.zip包含

ECGData.matholds the data used in this example. The .txt file, Modified_physionet_data.txt, is required by PhysioNet's copying policy and provides the source attributions for the data as well as a description of the pre-processing steps applied to each ECG recording.

加载文件

如果您按照上一节中的下载说明操作，请输入以下命令解压这两个存档文件。

解压缩(fullfile (tempdir,“physionet_ECG_data-main.zip”), tempdir)解压缩(fullfile (tempdir,“physionet_ECG_data-main”，“ECGData.zip”)，.．.fullfile (tempdir“ECGData”）)

解压ECGData.zip文件后，将数据加载到MATLAB中。

负载(fullfile (tempdir“ECGData”，“ECGData.mat”）)

ECGData是一个包含两个字段的结构数组:数据而且标签．数据是一个162 × 65536的矩阵，其中每一行都是128赫兹采样的心电图记录。标签是一个162 × 1的诊断标签单元阵列，每一行一个数据．三个诊断类别是:“ARR”(心律失常)、“CHF”(充血性心力衰竭)和“NSR”(正常窦性心律)。

创建培训和测试数据

将数据随机分为两组——训练数据集和测试数据集。辅助函数helperRandomSplit执行随机拆分。helperRandomSplit接受训练数据所需的分割百分比ECGData．的helperRandomSplit函数输出两个数据集以及每个数据集的一组标签。每一行的trainData而且testData是心电信号。的每个元素trainLabels而且testLabels包含数据矩阵对应行的类标签。在本例中，我们将每个类的70%的数据随机分配到训练集。剩下的30%用于测试(预测)并分配给测试集。

percent_train = 70;[trainData, testData trainLabels testLabels] =.．.helperRandomSplit (percent_train ECGData);

有113条记录trainData集和49条记录于一体testData．根据设计，训练数据包含69.75%(113/162)的数据。回想一下，ARR类占数据的59.26% (96/162)，CHF类占18.52% (30/162)，NSR类占22.22%(36/162)。检查每个类在训练和测试集中的百分比。每个类中的百分比与数据集中的总体类百分比一致。

Ctrain = countcats(分类(trainLabels)。/元素个数(trainLabels)。* 100

Ctrain =3×159.2920 18.5841 22.1239

ct = countcats(分类(testLabels)。/元素个数(testLabels)。* 100

ct =3×159.1837 18.3673 22.4490

情节样品

绘制四个随机选择的记录的前几千个样本ECGData．辅助函数helperPlotRandomRecords这是否。helperPlotRandomRecords接受ECGData和一个随机种子作为输入。初始种子设置为14，以便每个类中至少有一个记录被绘制出来。你可以执行helperPlotRandomRecords与ECGData作为唯一的输入参数，只要你想要得到与每个类相关的ECG波形的多样性的感觉。您可以在本示例末尾的支持函数一节中找到该助手函数的源代码。

helperPlotRandomRecords (ECGData 14)

特征提取

提取每个信号分类中使用的特征。本例使用从每个信号的8个块上提取的特征，时长约为一分钟(8192个样本):

此外，在整个数据长度[6]上提取每个信号的多尺度小波方差估计。采用了小波方差的无偏估计。这要求在方差估计中只使用至少一个小波系数不受边界条件影响的水平。对于2^16(65,536)的信号长度和“db2”小波，结果是14级。

这些特征的选择基于已发表的研究，证明了它们在心电波形分类中的有效性。这并不是一个详尽的或优化的功能列表。

每个窗口的AR系数使用Burg方法估计，arburg．在[8]中，作者使用模型顺序选择方法来确定AR(4)模型为类似分类问题中的心电波形提供了最佳拟合。在[5]中，在小波包树的终端节点上计算信息论度量香农熵，并将其与随机森林分类器一起使用。这里我们使用非抽取小波包变换，modwpt，降至第4层。

给出了下面[5]的未模拟小波包变换香农熵的定义: $$\mathrm{年代}{E}_j＝-\underset{k＝1}{\overset{N}{\sum }}{p}_{j，k}＊\mathrm{日志}{p}_{j，k}$$在哪里 $$N$$第j个节点对应系数的个数和 $$p_{j，k}$$为第j个终端节点的小波包系数归一化平方。

用小波估计的两个分形测度作为特征。下面[4]，我们使用由得到的奇异谱的宽度dwtleader用来测量心电信号的多重分形特性。我们也使用缩放指数的第二个累积。标度指数是描述信号在不同分辨率下幂律行为的基于标度的指数。第二个累积量广泛地表示比例指数偏离线性。

得到整个信号的小波方差modwtvar．小波方差测量信号的变异性的尺度，或等效的变异性的信号在八度频带的频率间隔。

总共有190个特征:32个AR特征(每个块4个系数)，128个Shannon熵值(每个块16个值)，16个分形估计(每个块2个)，14个小波方差估计。

的helperExtractFeatures函数计算这些特征，并将它们连接到每个信号的特征向量。您可以在本示例末尾的支持函数一节中找到该助手函数的源代码。

timeWindow = 8192;ARorder = 4;MODWPTlevel = 4;[trainFeatures, testFeatures featureindices] =.．.helperExtractFeatures (trainData、testData timeWindow、ARorder MODWPTlevel);

trainFeatures而且testFeatures分别是113 × 190和49 × 190矩阵。这些矩阵的每一行都是对应心电数据的特征向量trainData而且testData,分别。在创建特征向量时，数据从65536个样本减少到190个元素向量。这大大减少了数据量，但目标不仅仅是减少数据量。目标是将数据减少到更小的特征集，从而捕获类之间的差异，以便分类器能够准确地分离信号。构成这两者的特征的索引trainFeatures而且testFeatures包含在结构数组中，featureindices．您可以使用这些索引按组探索特性。作为一个例子，考察了第一个时间窗的奇异谱中的Holder指数的范围。绘制整个数据集的数据。

allFeatures = [trainFeatures; testFeatures];allLabels = [trainLabels; testLabels];图箱线图(allFeatures (:, featureindices.HRfeatures (1)), allLabels,“缺口”，“上”) ylabel (“霍尔德指数范围”)标题(奇异谱群范围(第一时间窗)网格)在

您可以对该特征进行单向方差分析，并确认箱线图中出现的情况，即ARR和NSR组的范围明显大于CHF组。

[p anovatab st] = anova1 (allFeatures (:, featureindices.HRfeatures (1)),.．.allLabels);

c = multcompare(圣“显示”，“关闭”）

c =3×61.0000 2.0000 0.0176 0.1144 0.2112 0.0155 1.0000 3.0000 -0.1591 -0.0687 0.0218 0.1764 2.000 3.0000 -0.2975 -0.1831 -0.0687 0.0005

作为一个额外的例子，考虑三组的第二低频率(第二大尺度)小波子带的方差差异。

箱线图(allFeatures (:, featureindices.WVARfeatures (end-1)), allLabels,“缺口”，“上”) ylabel (小波变化的)标题(“分组小波方差”网格)在

如果你对这个特征进行方差分析，你会发现NSR组在这个小波子带的方差明显低于ARR组和CHF组。这些示例只是为了说明单个特性如何用于分离类。虽然单独一个特性是不够的，但我们的目标是获得足够丰富的特性集，以使分类器能够分离所有三个类。

信号分类

既然数据已经被简化为每个信号的特征向量，下一步就是使用这些特征向量对心电信号进行分类。你可以使用分类学习者应用程序快速评估大量的分类器。在本例中，使用了一个具有二次核的多类支持向量机。进行了两种分析。首先，我们使用整个数据集(训练集和测试集)，并使用5倍交叉验证估计误分类率和混淆矩阵。

特点= [trainFeatures;testFeatures];rng(1) template = templateSVM.．.“KernelFunction”，多项式的，.．.“PolynomialOrder”2,.．.“KernelScale”，“汽车”，.．.“BoxConstraint”, 1.．.“标准化”,真正的);模型= fitcecoc (.．.的特性,.．.(trainLabels; testLabels),.．.“学习者”模板,.．.“编码”，“onevsone”，.．.“类名”, {“加勒比海盗”，瑞士法郎的，“签约”});kfoldmodel = crossval(模型,“KFold”5);classLabels = kfoldPredict (kfoldmodel);损失= kfoldLoss (kfoldmodel) * 100

损失= 8.0247

[confmatCV, grouporder] = confusionmat ([trainLabels; testLabels], classLabels);

5倍分类误差为8.02%(正确率为91.98%)。混淆矩阵,confmatCV，显示了哪些记录分类错误。grouporder给出了组的顺序。ARR组2例误诊为CHF, CHF组8例误诊为ARR, 1例误诊为NSR, NSR组2例误诊为ARR。

精确度，召回率和F1分数

在分类任务中，一个类的精度是正确的积极结果的数量除以积极结果的数量。换句话说，在分类器分配给给定标签的所有记录中，实际属于该类的比例是多少。召回率定义为正确标签的数量除以给定类的标签数量。具体地说，在属于一个类的所有记录中，分类器将该类标记为该类的比例是多少。在判断机器学习系统的准确性时，理想情况下，您希望在精确度和召回率两方面都做得好。例如，假设我们有一个分类器，它将每个记录标记为ARR。那么我们对ARR类的召回率将是1(100%)。属于ARR类的所有记录将被标记为ARR。然而，精确度会很低。因为我们的分类器将所有记录标记为ARR，所以在本例中会有66个假阳性，精度为96/162，或0.5926。 The F1 score is the harmonic mean of precision and recall and therefore provides a single metric that summarizes the classifier performance in terms of both recall and precision. The following helper function computes the precision, recall, and F1 scores for the three classes. You can see howhelperPrecisionRecall通过检查支持函数部分中的代码，计算精度、召回率和基于混淆矩阵的F1得分。

CVTable = helperPrecisionRecall (confmatCV);

可以显示返回的表helperPrecisionRecall使用以下命令。

disp (CVTable)

精度召回F1_Score _________ ______ ________ ARR 90.385 97.917 94 CHF 91.304 70 79.245 NSR 97.143 94.444 95.775

ARR和NSR类的精确度和召回率都很好，而CHF类的召回率明显较低。

对于下一个分析，我们只对训练数据(70%)拟合一个多类二次支持向量机，然后使用该模型对30%的数据进行预测，以供测试。在测试集中有49个数据记录。

模型= fitcecoc (.．.trainFeatures,.．.trainLabels,.．.“学习者”模板,.．.“编码”，“onevsone”，.．.“类名”, {“加勒比海盗”，瑞士法郎的，“签约”});predLabels =预测(模型、testFeatures);

使用以下方法确定正确预测的个数，得到混淆矩阵。

correctPredictions = strcmp (predLabels testLabels);testAccuracy = (correctPredictions) /长度总和(testLabels) * 100

testAccuracy = 97.9592

[confmatTest, grouporder] = confusionmat (testLabels predLabels);

测试数据集上的分类准确率约为98%，混淆矩阵显示有一个CHF记录被错误分类为NSR。

与交叉验证分析中所做的类似，获得测试集的精度、召回率和F1分数。

testTable = helperPrecisionRecall (confmatTest);disp (testTable)

精度召回F1_Score _________ ______ ________ ARR 100 100 100瑞士法郎100 88.889 94.118 NSR 91.667 100 95.652

原始数据分类与聚类

从前面的分析中自然产生了两个问题。为了获得良好的分类结果，特征提取是必要的吗?是否需要分类器，或者这些特性可以在没有分类器的情况下将组分隔开?要解决第一个问题，请对原始时间序列数据重复交叉验证结果。注意，下面的步骤计算量很大，因为我们将SVM应用于162 × 65536的矩阵。如果您不希望自己运行此步骤，结果将在下一段中描述。

rawData = [trainData; testData];标签= [trainLabels; testLabels];rng(1) template = templateSVM.．.“KernelFunction”，多项式的，.．.“PolynomialOrder”2,.．.“KernelScale”，“汽车”，.．.“BoxConstraint”, 1.．.“标准化”,真正的);模型= fitcecoc (.．.rawData,.．.(trainLabels; testLabels),.．.“学习者”模板,.．.“编码”，“onevsone”，.．.“类名”, {“加勒比海盗”，瑞士法郎的，“签约”});kfoldmodel = crossval(模型,“KFold”5);classLabels = kfoldPredict (kfoldmodel);损失= kfoldLoss (kfoldmodel) * 100

损失= 33.3333

[confmatCVraw, grouporder] = confusionmat ([trainLabels; testLabels], classLabels);rawTable = helperPrecisionRecall (confmatCVraw);disp (rawTable)

精度召回F1_Score _________ ______ ________ ARR 64 100 78.049 CHF 100 13.333 23.529 NSR 100 22.222 36.364

原始时间序列数据的误分类率为33.3%。重复精度、回忆率和F1得分分析显示，CHF组(23.52)和NSR组(36.36)的F1得分都非常低。获取每个信号的幅值离散傅里叶变换(DFT)系数，进行频域分析。由于数据是实值的，我们可以利用傅立叶幅度是偶函数这一事实，利用DFT实现一些数据约简。

rawDataDFT = abs (fft (rawData [], 2));rawDataDFT = rawDataDFT (:, 1:2 ^ 16/2 + 1);rng(1) template = templateSVM.．.“KernelFunction”，多项式的，.．.“PolynomialOrder”2,.．.“KernelScale”，“汽车”，.．.“BoxConstraint”, 1.．.“标准化”,真正的);模型= fitcecoc (.．.rawDataDFT,.．.(trainLabels; testLabels),.．.“学习者”模板,.．.“编码”，“onevsone”，.．.“类名”, {“加勒比海盗”，瑞士法郎的，“签约”});kfoldmodel = crossval(模型,“KFold”5);classLabels = kfoldPredict (kfoldmodel);损失= kfoldLoss (kfoldmodel) * 100

损失= 19.1358

[confmatCVDFT, grouporder] = confusionmat ([trainLabels; testLabels], classLabels);dftTable = helperPrecisionRecall (confmatCVDFT);disp (dftTable)

精度召回F1_Score _________ ______ ________ ARR 76.423 97.917 85.845瑞士法郎100 26.667 42.105 NSR 93.548 80.556 86.567

使用DFT量级将误分类率降低到19.13%，但这仍然是使用我们的190个特征获得的错误率的两倍多。这些分析表明，分类器得益于对特征的仔细选择。

要回答关于分类器角色的问题，尝试只使用特征向量对数据进行聚类。使用k-均值聚类和差距统计量来确定最佳的聚类数量和聚类分配。允许数据有1到6个集群的可能性。

rng默认的伊娃= evalclusters(特性,“kmeans”，“差距”，“中”[1:6]);伊娃

eva = GapEvaluation with properties: NumObservations: 162 InspectedK: [1 2 3 4 5 6] CriterionValues: [1.2777 1.3539 1.3644 1.3570 1.3591 1.3752] OptimalK: 3

间隙统计量表明最佳集群数量为3个。然而，如果您查看三个聚类中的每一个记录的数量，您会发现基于特征向量的k-means聚类在分离三个诊断类别方面做得很差。

countcats(分类(eva.OptimalY))

ans =3×161 74 27

回想一下，有96人属于“抗辐射”级别，30人属于“抗辐射”级别，36人属于“噪音敏感”级别。

总结

本例利用信号处理技术从心电信号中提取小波特征，并将心电信号分为三类。特征提取不仅减少了大量的数据，它还捕获了ARR、CHF和NSR类之间的差异，交叉验证结果和支持向量机分类器在测试集上的性能证明了这一点。该示例进一步说明，将支持向量机分类器应用于原始数据会导致较差的性能，不使用分类器聚类特征向量也是如此。单独的分类器和特征都不足以区分类。然而，当在使用分类器之前使用特征提取作为数据缩减步骤时，这三个类被很好地分离了。

参考文献

支持功能

helperPlotRandomRecords绘制4个随机选取的心电信号ECGData．

函数helperPlotRandomRecords (ECGData randomSeed)此函数仅用于支持XpwWaveletMLExample。它可能%更改或在未来版本中删除。如果输入参数个数= = 2 rng (randomSeed)结束M =大小(ECGData.Data, 1);idxsel = randperm (M, 4);为subplot(2,2,numplot) plot(ECGData.Data(idxsel(numplot)，1:300))“伏”）如果2 xlabel(“样本”）结束标题(ECGData.Labels {idxsel (numplot)})结束结束

helperExtractFeatures提取指定大小的数据块的小波特征和AR系数。这些特征被串联成特征向量。

函数[trainData,testData, featureindexes] = helperExtractFeatures(trainData,testData,T,AR_order,level)此函数只支持XpwWaveletMLExample。它可能会改变%将在未来的版本中删除。trainFeatures = [];testFeatures = [];为idx =1:size(trainData,1) x = trainData(idx，:);x =去趋势(x, 0);AR_order arcoefs = blockAR (x, T);se = shannonEntropy (x T级);(cp, rh) =领导人(x, T);wvar = modwtvar (modwt (x,“db2”)，“db2”）;trainFeatures = [trainFeatures;Arcoefs se cp rh wvar'];% #好< AGROW >结束为x1 = testData(idx，:);x1 =去趋势(x1, 0);arcoefs = blockAR (x1, AR_order, T);se = shannonEntropy (x1, T,水平);(cp, rh) =领导人(x1, T);wvar = modwtvar (modwt (x1,“db2”)，“db2”）;testFeatures = [testFeatures;arcoefs se cp rh wvar'];% #好< AGROW >结束featureindices =结构();% 4 * 8featureindices。ARfeatures = 1:32;startidx = 33;endidx = 33 + (16 * 8) 1;featureindices。年代Efeatures = startidx:endidx; startidx = endidx+1; endidx = startidx+7; featureindices.CP2features = startidx:endidx; startidx = endidx+1; endidx = startidx+7; featureindices.HRfeatures = startidx:endidx; startidx = endidx+1; endidx = startidx+13; featureindices.WVARfeatures = startidx:endidx;结束函数se = shanonentropy (x,numbuffer,level) numwindows = numel(x)/numbuffer;缓冲(x, y = numbuffer);se = 0(2 ^的级别,大小(y, 2));为Kk = 1:size(y,2) WPT = modwpt(y(:， Kk)，level);%跨时间总和E = (wpt。^ 2,2)总和;j = wpt。^ 2. / E;以下是eps(1)se (: kk) =总和(j。*日志(j + eps), 2);结束se =重塑(se、2 ^ * numwindows水平,1);se = se”;结束函数arcfs = blockAR(x,order,numbuffer) numwindows = numl (x)/numbuffer;缓冲(x, y = numbuffer);arcfs = 0(顺序,大小(y, 2));为Kk = 1:size(y,2) artmp = arburg(y(:， Kk)，order);arcfs (:, kk) = artmp(2:结束);结束arcfs =重塑(arcfs订单* numwindows 1);arcfs = arcfs ';结束函数[cp,rh] = leaders(x,numbuffer) y = buffer(x,numbuffer);cp = 0(1、大小(y, 2));rh = 0(1、大小(y, 2));为kk = 1:尺寸(y, 2) [~ h, cptmp] = dwtleader (y (:, kk));cp (kk) = cptmp (2);rh (kk) = (h);结束结束

helperPrecisionRecall返回基于混淆矩阵的精度、召回率和F1分数。将结果输出为MATLAB表。

函数PRTable = helperPrecisionRecall (confmat)此函数只支持XpwWaveletMLExample。它可能会改变%将在未来的版本中删除。precisionARR = confmat(1,1) /笔(confmat (: 1)) * 100;precisionCHF = confat (2,2)/sum(confat (:，2))*100;precisionNSR = confat (3,3)/sum(confat (:，3))*100;recallARR = confmat(1,1) /笔(confmat (1:)) * 100;recallCHF = confmat(2, 2) /笔(confmat (2:)) * 100;recallNSR = confmat(3,3) /笔(confmat (3:)) * 100;F1ARR = 2 * precisionARR * recallARR / (precisionARR + recallARR);F1CHF = 2 * precisionCHF * recallCHF / (precisionCHF + recallCHF);F1NSR = 2 * precisionNSR * recallNSR / (precisionNSR + recallNSR);构建一个MATLAB表来显示结果。PRTable = array2table([precisionARR recallARR F1ARR;.．.precisionCHF recallCHF F1CHF;precisionNSR recallNSR.．.F1NSR),“VariableNames”, {“精度”，“回忆”，“F1_Score”}，“RowNames”，.．.｛“加勒比海盗”，瑞士法郎的，“签约”});结束