数据描述

本例使用了三组或三类人的ECG数据:心律失常者、充血性心力衰竭者和窦性心律正常者。该示例使用了来自三个PhysioNet数据库的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单元阵列的诊断标签，每行一个数据．三种诊断类型是:“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(categorical(trainLabels))./numel(trainLabels).*100

Ctrain =3×159.2920 18.5841 22.1239

Ctest = countcats(categorical(testLabels))./numel(testLabels).*100

ct =3×159.1837 18.3673 22.4490

情节样品

绘制4个随机选择的记录的前几千个样本ECGData．辅助函数helperPlotRandomRecords这是否。helperPlotRandomRecords接受ECGData和一个随机种子作为输入。初始种子设置为14，以便绘制每个类的至少一条记录。你可以执行helperPlotRandomRecords与ECGData作为唯一的输入参数，你可以多次获得与每个类别相关的ECG波形的多样性。您可以在本示例末尾的支持函数一节中找到这个helper函数的源代码。万博1manbetx

helperPlotRandomRecords (ECGData 14)

特征提取

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

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

这些特征是根据已发表的研究选择的，这些研究证明了它们在分类ECG波形方面的有效性。这并不是一个详尽或优化的特性列表。

每个窗口的AR系数用Burg方法估计，arburg．在[8]中，作者使用模型顺序选择方法来确定AR(4)模型在类似的分类问题中为ECG波形提供了最佳拟合。在[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函数计算这些特征，并将它们连接到每个信号的特征向量中。您可以在本示例末尾的支持函数一节中找到这个helper函数的源代码。万博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．您可以使用这些索引按组探索特性。作为一个例子，在奇异谱的第一时间窗口检查Holder指数的范围。为整个数据集绘制数据图。

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

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

[p,anovatab,st] = anova1(allFeatures(:， featureindexes . hrfeatures (1)))，.．.allLabels);

C = multicompare (st，“显示”，“关闭”）

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.0000 3.0000 -0.2975 -0.1831 -0.0687 0.0005

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

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

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

信号分类

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

features =[火车特征;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分数。万博1manbetx

CVTable = helpprecisionrecall (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%)拟合了一个多类二次SVM，然后使用该模型对用于测试的30%的数据进行预测。在测试集中有49个数据记录。

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

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

correctforecasts = strcmp(predLabels,testLabels);testAccuracy = sum(correctforecasts)/length(testLabels)*100

testAccuracy = 97.9592

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

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

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

testTable = helpprecisionrecall (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-均值聚类和差距统计量来确定最优的聚类数量和聚类分配。允许数据有1到6个集群。

rng默认的Eva = 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

回想一下，ARR类有96人，CHF类有30人，NSR类有36人。

总结

本例采用信号处理的方法从心电信号中提取小波特征，并利用小波特征将心电信号分为三类。特征提取不仅导致了大量的数据缩减，它还捕获了ARR、CHF和NSR类之间的差异，这由交叉验证结果和SVM分类器在测试集上的性能所证明。该示例进一步表明，将SVM分类器应用于原始数据会导致性能较差，就像不使用分类器对特征向量进行聚类一样。单独的分类器和特征都不足以将类分开。然而，当在使用分类器之前使用特征提取作为数据约简步骤时，这三个类被很好地分开了。

参考文献

万博1manbetx支持功能

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

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

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

函数[trainFeatures, testFeatures,featureindices] = helperExtractFeatures(trainData,testData,T,AR_order,level)此函数仅支持xpwwavetmlexample。万博1manbetx它可能会改变%将在未来的版本中删除。trainFeatures = [];testFeatures = [];为idx =1:size(trainData,1) x = trainData(idx，:);X =趋势(X,0);arcoefs = blockAR(x,AR_order,T);se = shannonEntropy(x,T,level);[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,level);[cp,rh] =导联(x1,T);Wvar = modwtvar(modwt(x1，“db2”)，“db2”）;testFeatures = [testFeatures;arcoefs se cp rh wvar'];% #好< AGROW >结束Featureindices = struct();% 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 = shannonEntropy(x,numbuffer,level) numwindows = numel(x)/numbuffer;Y = buffer(x,numbuffer);Se = 0 (2^level,size(y,2));为Kk = 1:size(y,2) WPT = modwpt(y(:， Kk)，level);%时间总和E = sum(wpt.^2,2);Pij = wpt.^2./E;以下为eps(1)se(:，kk) = -sum(Pij.*log(Pij+eps)，2);结束Se =重塑(Se,2^level*numwindows,1);Se = Se ';结束函数arcfs = blockAR(x,order,numbuffer) numwindows = nummel (x)/numbuffer;Y = buffer(x,numbuffer);Arcfs = 0 (order,size(y,2));为Kk = 1:size(y,2) artmp = arburg(y(:， Kk)，order);Arcfs (:，kk) = artmp(2:end);结束Arcfs =重塑(Arcfs,order*numwindows,1);Arcfs = Arcfs ';结束函数[cp,rh] = leaders(x,numbuffer) y = buffer(x,numbuffer);Cp = 0 (1,size(y,2));Rh = 0 (1,size(y,2));为kk = 1:尺寸(y, 2) [~ h, cptmp] = dwtleader (y (:, kk));Cp (kk) = cptmp(2);Rh (kk) =范围(h);结束结束

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

函数PRTable = helperPrecisionRecall(confmat)此函数仅支持xpwwavetmlexample。万博1manbetx它可能会改变%将在未来的版本中删除。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)/sum(confmat(1，:))*100;recallCHF = confmat(2,2)/sum(confmat(2，:))*100;recallNSR = confmat(3,3)/sum(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”，.．.｛“加勒比海盗”，瑞士法郎的，“签约”});结束