信号使用小波的特性和支持向量机分类金宝app

数据描述

这个示例使用心电图数据来自三组,或类,人:人与心律失常、人与充血性心力衰竭,而人与正常窦性节律。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单元阵列的诊断标签,一个用于每一行的数据。这三个诊断类别是:“加勒比海盗”(心律失常),瑞士法郎(充血性心力衰竭),和“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),瑞士法郎类占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波形与每个类有关。你能找到这个辅助函数的源代码在支持功能部分结束时这个例子。金宝app

helperPlotRandomRecords (ECGData 14)

特征提取

提取每个信号的信号分类中使用的特性。这个示例使用以下特性提取8块的每个信号大约一分钟时间(8192年样本):

此外,多尺度小波提取的方差估计为每个信号在整个数据长度[6]。使用小波方差的无偏估计。这要求只有水平至少一个小波系数不受边界条件用于方差估计。的信号长度2 ^ 16(65536)和db2的小波这导致14水平。

这些特性选择基于出版研究展示心电图波形分类的有效性。这并不是一个详尽的或优化的特性列表。

每个窗口的AR系数估计使用伯格的方法,arburg。在[8]中,作者使用模型顺序选择的方法来确定,一个基于“增大化现实”技术(4)模型提供最适合心电图波形相似的分类问题。在[5],信息理论测量,香农熵,计算终端节点的小波包树和使用一个随机森林分类器。在这里我们使用nondecimated小波包变换,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作为一个衡量的ECG信号的多重分形性质。我们也使用的第二个累积量比例指数。标度指数是scale-based描述幂律指数行为信号在不同的决议。金宝搏官方网站第二个累积量大致代表了从线性标度指数的离开。

小波变化在整个信号得到使用modwtvar。小波变化信号的方差措施规模,或者说可变性在倍频带信号频率间隔。

总共有190个特点:32基于“增大化现实”技术的特点(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由- 190和113 - 49 -到- 190矩阵,分别。这些矩阵的每一行是一个特征向量对应的心电图数据trainData和testData,分别。在创建特征向量,样本数据从65536减少到190元素向量。这是一个显著减少数据,但是我们的目标不仅仅是减少数据。目的是为了减少数据小得多的特性捕获一个分类器类之间的差异,以便能够准确地分离信号。的指数特性,两者组成trainFeatures和testFeatures结构中包含数组,featureindices。您可以使用这些指标来探索功能组。作为一个例子,考察的范围霍尔德指数在首次奇点光谱窗口中。绘制整个数据集的数据。

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

可以执行单向方差分析在此功能并确认箱线图中出现,即加勒比海盗和NSR组比瑞士法郎组范围明显增大。

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

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

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

如果你在这个特性进行方差分析,你会发现NSR集团已在这个小波子带方差显著低于加勒比海盗和瑞士法郎组。这些例子只是想说明个人特性提供单独的类。虽然独自一个特性是不够的,我们的目标是获得一个足够丰富的特性集,使分类器分离三个类。

信号分类

现在数据已经减少为每个信号特征向量,下一步是使用这些特征向量分类ECG信号。您可以使用分类学习者应用快速评估大量的分类器。在这个例子中,多支持向量机的二次使用内核。两个执行分析。首先我们使用整个数据集(训练集和测试集)和估计误分类率和混淆矩阵使用5倍交叉验证。

特点= [trainFeatures;testFeatures];rng(1)模板= 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给组的顺序。两个加勒比海盗集团并被错误地归类为瑞士法郎,瑞士法郎的八组是不是加勒比海盗,一个NSR,和两个从NSR组并被错误地归类为加勒比海盗。

精度、召回和F1的分数

在分类任务中,一个类的精度是正确的积极成果的数量除以数量的积极成果。换句话说,所有的记录,分类器分配给定的标签,比例实际上属于类。召回的定义是正确的标签的数量除以标签对于一个给定的类的数量。具体来说,属于一个类的所有记录,什么比例我们作为这个类的分类标签。机器学习系统在判断的准确性,你最好要做的精度和召回。例如,假设我们有一个分类器,每个记录ARR的标签。我们的回忆的ARR类将是1 (100%)。所有记录属于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得分基于混淆矩阵通过检查支持函数中的代码部分。金宝app

CVTable = helperPrecisionRecall (confmatCV);

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

disp (CVTable)

精密召回F1_Score _____ ______月______ ARR 90.385 - 97.917 94 91.304 70 79.245 NSR 97.143 94.444 95.775瑞士法郎

精度和召回都有利于ARR和NSR类,而回忆是瑞士法郎类显著降低。

接下来的分析中,我们只适合多层次二次SVM训练数据(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%,混淆矩阵显示一个记录是瑞士法郎并被错误地归类为“新丝路”。

相似是什么做的交叉验证分析,得到精度,召回,F1测试集的分数。

testTable = helperPrecisionRecall (confmatTest);disp (testTable)

精密召回F1_Score _____ ______月______ ARR 100 100 100 95.652 100 88.889 94.118 91.667 NSR 100瑞士法郎

在原始数据分类和聚类

从前面分析两个自然的问题出现。特征提取是必要的为了达到好的分类结果吗?分类器是必要的,这些特性也可以单独的组没有一个标识符?为了解决第一个问题,对原始时间序列数据重复交叉验证结果。注意,下面是一个计算代价高昂,因为我们将支持向量机应用于162 - 65536矩阵。如果你不希望这一步自己运行,结果被描述在接下来的段落。

rawData = [trainData; testData];标签= [trainLabels; testLabels];rng(1)模板= 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 100 13.333 23.529 NSR 100 22.222 - 36.364瑞士法郎

原始时间序列数据的误分类率是33.3%。重复精度、召回和F1得分分析显示非常贫穷的F1得分(23.52)瑞士法郎和NSR组(36.36)。获得级离散傅里叶变换(DFT)系数为每个信号进行频域分析。因为数据是实值,我们可以使用DFT实现减少一些数据利用傅里叶级数是偶函数。

rawDataDFT = abs (fft (rawData [], 2));rawDataDFT = rawDataDFT (:, 1:2 ^ 16/2 + 1);rng(1)模板= 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 - means聚类以及统计的差距确定最优数量的集群和集群任务。允许1到6集群的数据的可能性。

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

伊娃= GapEvaluation属性:NumObservations: 162 InspectedK: (1 2 3 4 5 6) CriterionValues: [1.2777 1.3539 1.3644 1.3570 1.3591 1.3752] OptimalK: 3

统计表明,集群的最优数量的差距是三。然而,如果你看看记录的数量在每个三个集群,可以看到,基于特征向量的k - means聚类做了很好地分离三个诊断类别。

countcats(分类(eva.OptimalY))

ans =3×161 74 27

回想一下,有96人ARR类,30在瑞士法郎类,36 NSR类。

总结

这个例子使用信号处理从心电信号中提取小波的特性和使用这些特性对ECG信号进行分类为三个类。特征提取的结果不仅在减少大量的数据,它还捕获的加勒比海盗之间的差异,瑞士法郎,NSR类证明了交叉验证结果和SVM分类器的性能测试集。进一步证明了应用支持向量机分类器的原始数据聚类特征向量一样导致表现不佳不使用一个标识符。分类和功能就有充分的理由去单独的类。然而,当特征提取被用作数据简化步骤使用分类器之前,三个类的分离。

引用

金宝app支持功能

helperPlotRandomRecords情节四ECG信号随机选择ECGData。

函数helperPlotRandomRecords (ECGData randomSeed)%这个函数仅仅是为了支持XpwWaveletMLExample。金宝app它可能变化百分比或在将来的版本中被删除。如果输入参数个数= = 2 rng (randomSeed)结束M =大小(ECGData.Data, 1);idxsel = randperm (M, 4);为numplot = 1:4次要情节(2,2,numplot)情节(ECGData.Data (idxsel (numplot) 1:3000)) ylabel (“伏”)如果numplot > 2包含(“样本”)结束标题(ECGData.Labels {idxsel (numplot)})结束结束

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

函数[trainFeatures, testFeatures featureindices] = helperExtractFeatures (trainData testData T、AR_order级别)%这个函数只支持XpwWaveletMLExample金宝app。它可能改变或%在将来的版本中被删除。trainFeatures = [];testFeatures = [];为idx = 1:尺寸(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:尺寸(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;结束函数se = shannonEntropy (x, numbuffer级别)numwindows =元素个数(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 = blockAR (x,秩序,numbuffer) numwindows =元素个数(x) / numbuffer;缓冲(x, y = numbuffer);arcfs = 0(顺序,大小(y, 2));为kk = 1:尺寸(y, 2) artmp = arburg (y (:, kk),顺序);arcfs (:, kk) = artmp(2:结束);结束arcfs =重塑(arcfs订单* numwindows 1);arcfs = arcfs ';结束函数(cp, rh) =领导人(x, numbuffer) y =缓冲区(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) /笔(confmat (: 1)) * 100;precisionCHF = confmat(2, 2) /笔(confmat (:, 2)) * 100;precisionNSR = confmat(3,3) /笔(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”,…{“加勒比海盗”,瑞士法郎的,“签约”});结束