数据描述

该示例使用从三组或课程中获得的心电图数据：心脏心律失常的人，具有充血性心力衰竭的人，以及具有正常窦性心律的人。该示例使用来自三个PhysioIoneT数据库的162个ECG录制：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的细胞阵列诊断标签，每一行一个资料. 这个three diagnostic categories are: 'ARR' (arrhythmia), 'CHF' (congestive heart failure), and 'NSR' (normal sinus rhythm).

创建培训和测试数据

随机将数据分成两组 - 培训和测试数据集。辅助功能Helperrandomsplit.执行随机分割。Helperrandomsplit.接受培训数据的所需分割百分比EcgData.. 这个Helperrandomsplit.功能输出两个数据集以及每个数据集。每一排trainData和testData是心电信号。的每个元素列车标签和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））./ NOWEL（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波形的感觉。您可以在本例末尾的支持函数一节中找到这个帮助函数的源代码。金宝app

helperPlotRandomRecords (ECGData 14)

特征提取

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

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

根据已发布的研究选择这些功能，证明了它们在分类ECG波形方面的有效性。这并非旨在是穷举或优化的功能列表。

每个窗口的AR系数用Burg方法估计，阿尔堡. 在[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函数计算这些特征，并将它们连接成每个信号的特征向量。您可以在本例末尾的支持函数一节中找到这个帮助函数的源代码。金宝app

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

训练费和testfeatures.分别是113 × 190和49 × 190矩阵。这些矩阵的每一行都是相应心电数据的特征向量trainData和testData,分别。在创建特征向量时，将数据从65536个样本减少到190个元素向量。这是数据的显著减少，但目标不仅仅是数据的减少。目标是将数据减少到更小的特征集，以捕获类之间的差异，以便分类器能够准确地分离信号。这些特征的指数，它们构成了两者训练费和testfeatures.包含在结构数组中，featureindices。您可以使用这些索引按组探索特征。例如，第一次检查奇点光谱中的Holder指数范围。绘制整个数据集的数据。

AllFeatures = [训练性; Testfeatures];AllLabels = [TrainLabels; Testlabels];Figure Boxplot（AllFeatures（：，feathindices.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）模板=模板SVM(…“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分数

在分类任务中，类的精度是正确的阳性结果数除以阳性结果数。换句话说，在分类器分配给给定标签的所有记录中，实际属于该类的比例是多少。召回定义为正确标签的数量除以给定类别的标签数量。具体来说，在属于一个类的所有记录中，我们的分类器标记为该类的比例是多少。在判断机器学习系统的准确性时，理想情况下，您希望在精确度和召回率方面都做得很好。例如，假设我们有一个分类器，将每个记录标记为ARR，那么我们对ARR类的召回率将为1（100%）。属于ARR类的所有记录都将标记为ARR。但是，精度较低。因为我们的分类器将所有记录标记为ARR，所以在这种情况下，如果精度为96/162或0.5926，将有66个误报。F1分数是精确性和召回率的调和平均值，因此提供了一个单一的指标，总结了分类器在召回率和召回率方面的性能。以下帮助函数计算三个类的精度、召回率和F1分数。你可以看到HelperPrecisionRecall.通过检查支持功能部分中的代码来计算基于混淆矩阵的精度，调用和F1分数。金宝app

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类的查全率明显较低。

对于下一次分析，我们仅适用于培训数据的多级二次SVM（70％），然后使用该模型来对所删除的30％的数据进行预测。测试集中有49个数据记录。

模型= fitcecoc (…训练疗法，…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 x 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”，'差距'，“KList”[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

差距统计表明群集的最佳数量是三个。但是，如果您查看三个集群中的每一个中的记录数，则会看到基于特征向量的K-Means群集已经做出了分离三个诊断类别的糟糕工作。

countcats（分类（eva.optimality））

ans =3×161 74 27.

回想一下，ARR类中有96人，CHF等级中的30人，NSR课程中有36个。

总结

该示例使用信号处理来从ECG信号中提取小波特征，并使用这些功能将ECG信号分为三类。该特征提取不仅导致大量数据减少，它还捕获了ARR，CHF和NSR类之间的差异，如跨验证结果所示以及SVM分类器在测试集上的性能。该示例进一步说明了将SVM分类器应用于原始数据导致性能不佳，因为在不使用分类器的情况下群集特征向量的性能不佳。分类器和单独的特征都没有足以分离类。然而，当在使用分类器之前使用特征提取作为数据还原步骤时，三类分开很好。

工具书类

金宝app支持功能

HelperplotrandomRecords.绘制四个心电图信号随机选择EcgData.．

函数helperPlotRandomRecords (ECGData randomSeed)％此函数仅旨在支持XPWWaveletMlexample。金宝app它可能%更改或在未来的版本中删除。如果nargin == 2 rng（随机）结束m = size（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。金宝app它可能会改变％在将来的释放中删除。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：大小（y，2）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”;结束函数ARCFS = SLOBLAR（x，订单，numbuffer）numwindows = numel（x）/ numbuffer;缓冲(x, y = numbuffer);ARCFS =零（顺序，大小（y，2））;为Kk = 1:size(y,2) artmp = arburg(y(:， Kk)，order);arcfs (:, kk) = artmp(2:结束);结束arcfs=重塑（arcfs，order*numindows，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。金宝app它可能会改变％在将来的释放中删除。PrecisionArr = Confmat（1,1）/ Sum（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;Revallchf = Confmat（2,2）/总和（Confmat（2，：））* 100;Recallnsr = Confmat（3,3）/总和（Confmat（3，：））* 100;F1ARR = 2 * precisionARR * RECALLARR /（PrecisionArr + Recallarr）;f1chf = 2 * precisionchf * recallchf /（Precisionchf + Revallchf）;f1nsr = 2 * precisionnsr * Recallnsr /（Precisionnsr + Recallnsr）;％构造MATLAB表以显示结果。PRTable = array2table([precisionARR recallARR F1ARR;…precisionCHF recallCHF F1CHF;precisionNSR recallNSR…F1nsr]，'variablenames'，{“精度”，“回忆”，“F1_Score”}，'rownames'，…｛“加勒比海盗”，瑞士法郎的，“签约”}）;结束