小波时间散射在心电信号分类中的应用

这个示例使用:

打开生活的脚本

该示例显示如何使用小波时间散射和支持向量机（SVM）分类器对人的心电图（ECG）信号进行分类。万博1manbetx在小波散射中，数据通过一系列小波变换，非线性和平均传播，以产生时间序列的低方差表示。小波时间散射产生信号表示不敏感以在输入信号中移位而不牺牲类别辨别性。您必须拥有小波工具箱™和统计和机器学习工具箱™来运行此示例。本例中使用的数据是公开的生理网．在此示例中，您可以找到对此分类问题的深入学习方法利用小波分析和深度学习对时间序列进行分类还有这个例子中的机器学习方法基于小波特征和支持向量机的信号分类万博1manbetx．

数据描述

本例使用了三组或三类人的心电图数据:心律失常者、充血性心力衰竭者和窦性节律正常者。该示例使用了来自三个PhysioNet数据库的162个心电记录:MIT-BIH心律失常数据库[3] [5],MIT-BIH正常窦性节律数据库[3]BIDMC充血性心力衰竭数据库[2][3]。共有96例心律失常患者的记录，30例充血性心力衰竭患者的记录，36例窦性节律正常的患者的记录。目的是训练一个分类器来区分心律失常(ARR)、充血性心力衰竭(CHF)和正常窦性心律(NSR)。

下载数据

第一步是从GitHub库．要下载数据，请单击代码并选择下载zip.．保存文件physionet_ecg_data-main.zip.在有写权限的文件夹中。本例的说明假设您已经将文件下载到临时目录，(tempdir在MATLAB中)。如果您选择在不同的文件夹中下载数据，请修改后续解压缩和加载数据的说明tempdir．

该文件physionet_ecg_data-main.zip.包含

ECGData.zip
README.md

和ECGData.zip包含

ECGData.mat
Modified_physionet_data.txt
License.txt。

ECGData.mat持有s 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中。

load（fullfile（tempdir，“ECGData”那“ECGData.mat”）)

ECGData是一个有两个字段的结构数组:数据和标签．数据是一个162乘65536的矩阵，其中每一行是以128赫兹采样的心电图记录。每个心电时间序列的总持续时间为512秒。标签是一个162×1个单元格阵列诊断标签，每排一个数据．三个诊断类别是:“ARR”(心律失常)，“CHF”(充血性心力衰竭)，和“NSR”(正常窦性心律)。

创建培训和测试数据

将数据随机分成训练数据集和测试数据集。辅助函数helperRandomSplit执行随机分割。helperRandomSplit接受训练数据和所需的分割百分比ECGData．的helperRandomSplit函数输出两个数据集以及每个数据集的一组标签。每一行的训练DATA.和testdata.是一个心电图信号。每个元素trainLabels和testLabels包含数据矩阵的相应行的类标签。在本例中，我们将每个类中70%的数据随机分配给训练集。剩下的30%用于测试(预测)，并分配给测试集。

百分之= 70;[TrainData，TestData，TrainLabels，TestLabels] =．..helperRandomSplit (percent_train ECGData);

共有113项记录训练DATA.设置和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

ctest = countcats（分类（testlabels））./ numel（testlabels）。* 100

ct =3×159.1837 18.3673 22.4490

情节样品

绘制来自四个随机选择的记录的前几千样本ECGData．辅助函数helperPlotRandomRecords这是否。helperPlotRandomRecords接受ECGData和一个随机种子作为输入。初始种子设置为14，这样每个类至少可以绘制一条记录。你可以执行helperPlotRandomRecords借ECGData作为唯一的输入参数，只要你希望得到与每个类相关的各种心电图波形的感觉。您可以在本例末尾的支持函数一节中找到这个函数和所有帮助函数的源代码。万博1manbetx

helperPlotRandomRecords (ECGData 14)

小波时间散射

在小波时间散射网络中需要指定的关键参数是时不变的尺度、小波变换的次数以及每个小波滤波器组中每倍频程的小波数。在许多应用中，两个滤波器组的级联就足以实现良好的性能。在这个例子中，我们构造了一个带有默认滤波器组的小波时间散射网络:第一个滤波器组每倍频程有8个小波，第二个滤波器组每倍频程有1个小波。不变性比例设置为150秒。

N =大小(ECGData.Data 2);sn = waveletScattering (“SignalLength”N“InvarianceScale”, 150,“SamplingFrequency”, 128);

您可以使用以下内容可视化两个过滤器库中的小波滤波器。

[fb, f, filterparams] = filterbank (sn);Subplot (211) plot(f,fb{2}.psift) xlim([0 128]) grid在标题(“第一个滤波器组小波滤波器”）;

Subplot (212) plot(f,fb{3}.psift) xlim([0 128]) grid在标题(第二滤波器组小波滤波器）;Xlabel（“赫兹”）;

为了演示不变性的缩放，获得缩放函数的逆傅里叶变换，并在0秒钟内将其中心。两个垂直黑线标记为-75和75秒的边界。另外，从第一滤波器组绘制粗略级（最低频率）小波的实部和虚部。请注意，粗略级小波不会超过由缩放功能的时间支持确定的不变尺度。万博1manbetx这是小波时间散射的重要特性。

图;φ= ifftshift(传输线(神奇动物{1}.phift));psiL1 = ifftshift(传输线(神奇动物{2}.psift (:,)));t = (2 ^ 15:2 ^ 15 - 1) * 1/128;scalplt =情节(t,φ);持有在网格在ylim（[ -  1.5e-4 1.6e-4]）;绘图（[ -  75 -75]，[ -  1.5E-4 1.6E-4]，'k-'）;情节(75[75],[-1.5 1.6的军医]的军医,'k-'）;Xlabel（'秒'）;ylabel (“振幅”）;Wavplt = plot（t，[真实（psil1）imag（psil1）]）;图例（[scalplt wavplt（1）波普（2）]，{'缩放功能'那“Wavelet-Real部分”那“Wavelet-Imaginary部分”}）;标题({'缩放功能';“最粗尺度小波第一滤波器组”})举行从

构建散射网络后，将训练数据的散射系数以矩阵形式表示。当您运行Featurematrix.对于多个信号，每一列都被处理为单个信号。

scat_features_train = featureMatrix (sn, trainData ');

输出Featurematrix.这里是409乘16乘113。张量的每一页，scat_features_train.，为一个信号的散射变换。基于尺度函数的带宽，对小波散射变换进行了严格的时间下采样。在这种情况下，409条散射路径中的每条都有16个时间窗口。

为了得到一个与SVM分类器兼容的矩阵，将多信号散射变换重构为一个矩阵，其中每一列对应一个散射路径，每一行为散射时间窗。在本例中，您将获得1808行，因为训练数据中的113个信号每个都有16个时间窗口。

Nwin =大小(scat_features_train, 2);Scat_features_train = permute(Scat_features_train，[2 3 1]);scat_features_train =重塑(scat_features_train,．..大小(scat_features_train 1) *大小(scat_features_train, 2), []);

对测试数据重复上述过程。最初,scat_features_test409-by-16-by-49，因为训练集中有49个ECG波形。在为SVM分类器重塑后，特征矩阵为784×416。

scat_features_test = featureMatrix (sn, testData ');Scat_features_test = permute(Scat_features_test，[2 3 1]);scat_features_test =重塑(scat_features_test,．..大小(scat_features_test 1) *大小(scat_features_test, 2), []);

因为对于每个信号，我们获得了16个散射窗口，我们需要创建标签来匹配窗口的数量。辅助函数createSequenceLabels这是基于窗口的数量。

[sequence_labels_train, sequence_labels_test] = createSequenceLabels (Nwin、trainLabels testLabels);

交叉验证

为了分类，进行了两项分析。首先利用所有散射数据，用二次核拟合多类支持向量机;整个数据集共有2592个散射序列，其中162个信号分别对应16个散射序列。错误率，或损失，估计使用5倍交叉验证。

scat_features = [scat_features_train;scat_features_test];allLabels_scat = [sequence_labels_train;sequence_labels_test];rng (1);模板= templateSVM (．..“KernelFunction”那多项式的那．..“PolynomialOrder”2,．..“KernelScale”那'auto'那．..“BoxConstraint”, 1．..“标准化”,真正的);classificationSVM = fitcecoc (．..scat_features,．..Alllabels_Scat，．..'学习者'模板,．..“编码”那“onevsone”那．..“类名”, {'arr';瑞士法郎的;“签约”}）;kfoldmodel = crossval (classificationSVM,“KFold”5）;

计算损失和混淆矩阵。显示的准确性。

predLabels = kfoldPredict (kfoldmodel);损失= kfoldLoss (kfoldmodel) * 100;confmatCV = confusionmat (allLabels_scat predLabels)

confmatCV =3×31535 0 1 1 479 0 0 576

fprintf（'准确度为%2.2f % .\n'100 -损失);

准确率为99.92%。

准确率为99.88%，相当不错，但考虑到这里每个时间窗口都是单独分类的，实际结果更好。每个信号有16个单独的分类。使用简单多数投票获得每个散射表示的单个类预测。

classes = patporical（{'arr'那瑞士法郎的那“签约”}）;[classvotes，classcounts] = HelpermajorityVote（Predlabels，[TrainLabels; Testlabels]，课程）;

根据每组散射时间窗口的类预测模式，确定实际的交叉验证精度。如果模式对于给定的集合不是唯一的，HelpermajorityVote.返回由指示的分类错误“NoUniqueMode”．这导致混淆矩阵中的额外列，在这种情况下是所有零，因为每组散射预测存在唯一模式。

cvaccuracy = sum（eq（classvotes，分类（[trainlabels; testlabels]））/ 162 * 100;fprintf（'真正的交叉验证精度为％2.2f％。\ n'，cvaccuracy）;

真实交叉验证准确率为100.00%。

MVconfmatCV = confusionmat(分类([trainLabels;testLabels]), ClassVotes);MVconfmatCV

MVconfmatCV =4×496 0 0 0 0 0 0 0 0 0 0 36 0 0 0 0 0

散射已正确分类交叉验证模型中的所有信号。如果你检查ClassCounts，你看到2个错误分类的时间窗口confmatCV归因于2个信号，其中15个散射窗被正确分类。

支持向量机分类

在接下来的分析中，我们只对训练数据拟合了一个多类二次支持向量机(70%)，然后使用该模型对等待测试的30%数据进行预测。测试集中有49条数据记录。使用多数投票个别散射窗口。

模型= fitcecoc (．..scat_features_train,．..sequence_labels_train,．..'学习者'模板,．..“编码”那“onevsone”那．..“类名”, {'arr'那瑞士法郎的那“签约”}）;predLabels =预测(模型、scat_features_test);[TestVotes, TestCounts] = helperMajorityVote (predLabels、testLabels、类);testaccuracy =总和(eq (TestVotes分类(testLabels))) /元素个数(testLabels) * 100;fprintf（'测试精度为％2.2f％。\ n', testaccuracy);

测试精度为97.96％。

confusionchart(分类(testLabels)、TestVotes)

在测试数据集上的分类准确率约为98%。混淆矩阵显示一个CHF记录被误分类为ARR。其他48个信号都被正确分类了。

总结

这个例子使用小波时间散射和支持向量机分类器将心电波形分类为三个诊断类之一。小波散射被证明是一种强大的特征提取方法，它只需要用户指定的最小参数集就可以产生一组鲁棒的特征用于分类。将这个与示例进行比较基于小波特征和支持向量机的信号分类万博1manbetx，这需要大量的专业知识来进行分类的手工特征。通过小波时间散射，您只需要指定时间不变性的比例，滤波器库（或小波变换）的数量以及每个八度音程的小波次数。小波散射变换和SVM分类器的组合在交叉验证模型上产生了100％的分类，并且在将SVM应用于保持测试集的散射变换时，98％正确的分类。

参考文献

Anden J.， Mallat, S. 2014。深散射谱，IEEE信号处理汇刊，62,16,414 -4128页。
Baim DS, Colucci WS, Monrad ES, Smith HS, Wright RF, Lanoue A, Gauthier DF, Ransil BJ, Grossman W, Braunwald E.口服米力农治疗严重充血性心力衰竭患者的生存率。美国心脏病学会1986年3月;7(3): 661 - 670。
Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh, Mark RG, miietus JE, Moody GB, Peng C-K, Stanley HE。PhysioBank, PhysioToolkit和PhysioNet:复杂生理信号新研究资源的组成部分。循环．卷。101,23,33，2000年6月13日，PP。E215-E220。http://circ.ahajournals.org/content/101/23/e215.full
Mallat，S.，2012年。集团不变散散射。纯净和应用数学的通信，65,10，pp。1331-1398。
穆迪GB，Mark RG。MIT-BIH心律失常数据库的影响。IEEE Eng在Med和Biol 20（3）：45-50（2001年5月）。（PMID：11446209）

万博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)})结束结束

HelpermajorityVote.为每一组散射时间窗口查找预测类标签中的模式。该函数返回每个散射时间窗口集的类标签模式和类预测数量。如果没有唯一的模式，HelpermajorityVote.返回“错误”的类标签以指示散射窗口集是一个分类错误。

函数[ClassVotes, ClassCounts] = helperMajorityVote (predLabels、origLabels类)%该函数支持ecgwavettimescat万博1manbetxteringexample。它可能%更改或在未来的版本中删除。%如果标签不是分类的，则创建分类数组predLabels =分类(predLabels);origLabels =分类(origLabels);％期望Predlabels和Origlabels成为分类向量npr =元素个数(predLabels);Norig =元素个数(origLabels);Nwin = npr / Norig;predLabels =重塑(predLabels Nwin Norig);ClassCounts = countcats (predLabels);[mxcount, idx] = max (ClassCounts);ClassVotes =类(idx);%检查最大值中的任何领带，并确保它们被标记为%错误，如果模式出现多次modecnt = modecount (ClassCounts mxcount);ClassVotes (modecnt > 1) =分类({“错误”}）;classvotes = classvotes（:);%-------------------------------------------------------------------------函数modecnt = modecount（classcounts，mxcount）modecnt = Inf（大小（classcounts，2），1）;对于c = 1:size(ClassCounts,2) modecnt(nc) = histc(ClassCounts(:，nc)，mxcount(nc));结束结束％eof.结束