数据描述

本例使用了三组或三类人的心电图数据:心律失常者、充血性心力衰竭者和窦性节律正常者。该示例使用了来自三个PhysioNet数据库的162个心电记录:MIT-BIH心律失常数据库[3] [7],MIT-BIH正常窦性节律数据库[3]BIDMC充血性心力衰竭数据库[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作为唯一的输入参数，只要你希望得到与每个类相关的各种心电图波形的感觉。您可以在本例末尾的支持函数一节中找到这个帮助函数的源代码。万博1manbetx

helperPlotRandomRecords (ECGData 14)

特征提取

提取每个信号用于信号分类的特征。这个例子使用了从每个信号的8个块中提取的以下特征，大约持续一分钟(8192个样本):

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

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

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

给出[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函数计算这些特征，并将它们连接成每个信号的特征向量。您可以在本例末尾的支持函数一节中找到这个帮助函数的源代码。万博1manbetx

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．您可以使用这些索引来按组研究特性。作为一个例子，检查霍尔德指数的范围在奇异谱的第一次窗口。绘制整个数据集的数据。

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.0687 0.0218 0.1764 2.0000 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, NSR组中有1例误诊为ARR。

精确度，回忆和F1分数

在一个分类任务中，一个类的精度是正确的正结果数除以正结果数。换句话说，在分类器分配给给定标签的所有记录中，实际属于类的比例是多少。Recall定义为正确的标签数除以给定类的标签数。具体来说，在属于一个类的所有记录中，分类器将其标记为那个类的比例是多少?在判断你的机器学习系统的准确性时，理想情况下，你希望在精度和回忆两方面都做得很好。例如，假设我们有一个分类器，它将每一条记录标记为ARR。然后我们对ARR类的回忆将是1(100%)。所有属于ARR类的记录都被标记为ARR。然而，精确度会很低。因为我们的分类器将所有记录标记为ARR，在这种情况下，对于96/162或0.5926的精度，将有66个假阳性。 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分数。万博1manbetx

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 CHF 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 CHF 100 26.667 42.105 NSR 93.548 80.556 86.567

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

为了回答关于分类器的作用的问题，尝试仅使用特征向量对数据进行聚类。使用k-means聚类和间隙统计量来确定最佳的聚类数量和聚类分配。允许数据有1到6个集群的可能性。

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

eva = GapEvaluation with properties: NumObservations: 162 InspectedK: [1 2 3 4 5 6] criteria values: [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类之间的差异，这在交叉验证结果和SVM分类器在测试集上的性能上得到了证明。该例子进一步证明，将支持向量机分类器应用于原始数据会导致较差的性能，就像不使用分类器聚类特征向量一样。无论是分类器还是特性本身都不足以分离类。然而，在使用分类器之前，将特征提取作为数据约简步骤，这三个类被很好地分离了。

参考文献

万博1manbetx支持功能

helperPlotRandomRecords图中随机选取的四个心电信号ECGData．

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

helperExtractFeatures提取特定大小的数据块的小波特征和AR系数。特征被连接成特征向量。

函数[trainFeatures, testFeatures,featureindices] = helperExtractFeatures(trainData,testData,T,AR_order,level)%此函数仅支持XpwWaveletMLExample。万博1manbetx它可能会改变%将在未来的版本中被删除。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 >结束为idx =1:size(testData,1) 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;结束函数numwindows = nummel (x)/numbuffer;缓冲(x, y = numbuffer);se = 0(2 ^的级别,大小(y, 2));为Kk = 1: WPT = modwpt(y(:， Kk)，水平);跨越时间总和E = (wpt。^ 2,2)总和;j = wpt。^ 2. / E;%以下是eps(1)se (: kk) =总和(j。*日志(j + eps), 2);结束se =重塑(se、2 ^ * numwindows水平,1);se = se”;结束函数numwindows = nummel (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。万博1manbetx它可能会改变%将在未来的版本中被删除。precisionARR = confmat(1,1) /笔(confmat (: 1)) * 100;precisionCHF = confmat(2,2)/sum(confmat(: 2))*100;precisionNSR = confmat(3,3)/sum(confmat(: 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”，．..｛“加勒比海盗”，瑞士法郎的，“签约”}）;结束