数据描述

本例使用从三组或三类人群获得的心电图数据:心律失常患者、充血性心力衰竭患者和窦性心律正常患者。该示例使用来自三个PhysioNet数据库的162条ECG记录:MIT-BIH心律失常数据库[3] [7],麻省理工学院- 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波形的多样性。您可以在本示例末尾的“支持函数”一节中找到这个辅助函数的源代码。金宝app

helperPlotRandomRecords (ECGData 14)

特征提取

对每个信号提取用于信号分类的特征。本例使用在每个信号的8个块上提取以下特征，持续时间约为一分钟(8192个样本):

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

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

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

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

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 = 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应用程序快速评估大量分类器。在这个例子中，使用了一个二次核的多类支持向量机。进行了两项分析。首先，我们使用整个数据集(训练集和测试集)，并使用5倍交叉验证估计误分类率和混淆矩阵。

features = [trainFeatures;testFeatures];rng(1) template = templateSVM()…“KernelFunction”，多项式的，…“PolynomialOrder”2,…“KernelScale”，“汽车”，…“BoxConstraint”, 1…“标准化”,真正的);模型= fitcecoc(…的特性,…(trainLabels; testLabels),…“学习者”模板,…“编码”，“onevsone”，…“类名”, {“加勒比海盗”，瑞士法郎的，“签约”}）;Kfoldmodel = crossval(模型，“KFold”5);classLabels = kfoldPredict(kfoldmodel);loss = 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通过检查support Functions部分中的代码，根据混淆矩阵计算精度、召回率和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类别的查全率明显较低。

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

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

使用下面的公式来确定正确预测的数量并获得混淆矩阵。

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

testAccuracy = 97.9592

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

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

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

testTable = helperPrecisionRecall(conmattest);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);loss = kfoldLoss(kfoldmodel)*100

损失= 33.3333

[confmatCVraw,grouporder] = confusionmat([trainLabels;testLabels]，classLabels);rawTable = helperPrecisionRecall(conmatcvraw);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);loss = kfoldLoss(kfoldmodel)*100

损失= 19.1358

[confmatCVDFT,grouporder] = confusionmat([trainLabels;testLabels]，classLabels);dftTable = helperPrecisionRecall(conmatcvdft);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默认的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人。

总结

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

参考文献

金宝app支持功能

helperPlotRandomRecords绘制随机选择的四个心电信号ECGData．

函数helperPlotRandomRecords (ECGData randomSeed)此函数仅用于支持xpwwavetmlexample。金宝app它可能%更改或在将来的版本中删除。如果输入参数个数= = 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 > 2 xlabel(“样本”）结束标题(ECGData.Labels {idxsel (numplot)})结束结束

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

函数[trainFeatures, testFeatures,featureindices] = helperExtractFeatures(trainData,testData,T,AR_order,level)此函数仅支持xpwwavetmlexample。金宝app它可能会改变%将在以后的版本中删除。trainFeatures = [];testFeatures = [];为idx =1:size(trainData,1) x = trainData(idx，:);X = dettrend (X,0);arcoefs = blockAR(x,AR_order,T);se = shannonEntropy(x,T,level);[cp,rh] = leader (x,T);Wvar = modwtvar(modwt(x);“db2”)，“db2”）;trainFeatures = [trainFeatures;[参考译文];% #好< AGROW >结束为idx =1:size(testData,1);X1 = dettrend (X1,0);arcoefs = blockAR(x1,AR_order,T);se = shannonEntropy(x1,T,level);[cp,rh] = leaders(x1,T);Wvar = modwtvar(modwt(x1，“db2”)，“db2”）;testFeatures = [testFeatures;arcoefs se . rwvar '];% #好< 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) numwindow = nummel (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*numwindow,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*numwindow,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) = range(h);结束结束

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

函数PRTable = helperPrecisionRecall(确认)此函数仅支持xpwwavetmlexample。金宝app它可能会改变%将在以后的版本中删除。precisionARR = conmat (1,1)/sum(conmat (:，1))*100;precisionCHF = conmat (2,2)/sum(conmat (:，2))*100;precisionNSR = conmat (:，3) /sum(conmat (:，3))*100;recallARR = conmat (1,1)/sum(conmat (1，:))*100;recallCHF = conmat (2,2)/sum(conmat (2，:))*100;recallNSR = conmat (3,3)/sum(conmat (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;…precision;召回量;precisionNSR recallNSR…F1NSR),“VariableNames”, {“精度”，“回忆”，“F1_Score”}，“RowNames”，…{“加勒比海盗”，瑞士法郎的，“签约”}）;结束