小波时间散射在心电信号分类中的应用
这个例子展示了如何使用小波时间散射和支持向量机分类器对人体心电图信号进行分类。在小波散射中，数据通过一系列小波变换、非线性和平均来传播，以产生时间序列的低方差表示。小波时间散射产生对输入信号位移不敏感的信号表示，而不牺牲类的区别性。您必须拥有小波工具箱™和统计和机器学习工具箱™才能运行此示例。本例中使用的数据可从生理网．您可以在本例中找到解决这个分类问题的深度学习方法用小波分析和深度学习对时间序列进行分类在这个例子中是机器学习的方法基于小波特征和支持向量机的信号分类．
术语说明:在小波散射的背景下，“时间窗”一词是指对平滑操作的输出进行下采样后得到的样本数量。有关更多信息，请参见时间窗口．
数据描述
这个例子使用了从三组人(或三类人)获得的心电图数据:心律失常的人、充血性心力衰竭的人和窦性心律正常的人。该示例使用了来自三个PhysioNet数据库的162个心电图记录:MIT-BIH心律失常数据库[3] [5],MIT-BIH正常窦性节律数据库[3]充血性心力衰竭数据库[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中。
负载(fullfile (tempdir“ECGData”，“ECGData.mat”）)
ECGData是一个包含两个字段的结构数组:数据而且标签．数据是一个162 × 65536的矩阵，其中每一行都是128赫兹采样的心电图记录。每个心电图时间序列的总持续时间为512秒。标签是一个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(分类(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波形的多样性的感觉。您可以在本示例末尾的支持函数一节中找到这个和所有帮助函数的源代码。
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])网格在标题(“第一滤波器组小波滤波器”）;Subplot (212) plot(f,fb{3}.psift) xlim([0 128])网格在标题(“第二滤波器组小波滤波器”）;包含(“赫兹”）;
为了证明不变性尺度，得到尺度函数的傅里叶反变换，并以0秒为中心。两条垂直的黑线标记了-75秒和75秒的边界。此外，绘制出来自第一个滤波器组的最粗尺度(最低频率)小波的实部和虚部。注意，尺度最粗的小波不超过由尺度函数的时间支持决定的不变尺度。这是小波时间散射的一个重要性质。
图;φ= ifftshift(传输线(神奇动物{1}.phift));psiL1 = ifftshift(传输线(神奇动物{2}.psift (:,)));t = (2 ^ 15:2 ^ 15 - 1) * 1/128;scalplt =情节(t,φ);持有在网格在ylim(1.6[-1.5军医的军医]);情节(-75[-75],[-1.5 1.6的军医]的军医,“k——”）;情节(75[75],[-1.5 1.6的军医]的军医,“k——”）;包含(“秒”）;ylabel (“振幅”）;wavplt = plot(t，[real(psiL1) imag(psiL1)]);传奇([scalplt wavplt (1) wavplt (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_test是409 × 16 × 49因为在训练集中有49个心电图波形。对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”，“汽车”，.．.“BoxConstraint”, 1.．.“标准化”,真正的);classificationSVM = fitcecoc (.．.scat_features,.．.allLabels_scat,.．.“学习者”模板,.．.“编码”，“onevsone”，.．.“类名”, {“加勒比海盗”；瑞士法郎的；“签约”});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 0 576
流(“精度为%2.2f %。\n”100 -损失);
准确率为99.92%。
准确率为99.88%，非常好，但考虑到这里每个时间窗口都是单独分类的，实际结果更好。每个信号有16个单独的分类。使用简单多数投票获得每个散射表示的单一类预测。
类=分类({“加勒比海盗”，瑞士法郎的，“签约”});[ClassVotes, ClassCounts] = helperMajorityVote (predLabels [trainLabels;testLabels)、类);
根据每组散射时间窗的类预测模式确定实际交叉验证精度。如果给定集合的模式不是唯一的，helperMajorityVote返回由指示的分类错误“NoUniqueMode”．这导致在混淆矩阵中增加了一列，在这种情况下是全零，因为每一组散射预测都存在唯一的模式。
CVaccuracy =总和(eq (ClassVotes分类([trainLabels;testLabels]))) / 162 * 100;流('真正的交叉验证精度为%2.2f %。\n', CVaccuracy);
真正的交叉验证精度为100.00%。
MVconfmatCV = confusionmat(分类([trainLabels;testLabels]), ClassVotes);MVconfmatCV
MVconfmatCV =4×496 0 0 0 30 0 0 0 0 36 0 0 0 0 0 0 0 0
散射对交叉验证模型中的所有信号进行了正确的分类。如果你检查ClassCounts，你会看到2个错误分类的时间窗口confmatCV归因于两个信号，其中16个散射窗口中的15个被正确分类。
支持向量机分类
对于下一个分析，我们只对训练数据(70%)拟合一个多类二次支持向量机，然后使用该模型对30%的数据进行预测，以供测试。在测试集中有49个数据记录。使用多数投票的个别散射窗口。
模型= fitcecoc (.．.scat_features_train,.．.sequence_labels_train,.．.“学习者”模板,.．.“编码”，“onevsone”，.．.“类名”, {“加勒比海盗”，瑞士法郎的，“签约”});predLabels =预测(模型、scat_features_test);[TestVotes, TestCounts] = helperMajorityVote (predLabels、testLabels、类);testaccuracy =总和(eq (TestVotes分类(testLabels))) /元素个数(testLabels) * 100;流(“测试准确度为%2.2f %。\ n”, testaccuracy);
测试正确率为97.96%。
confusionchart(分类(testLabels)、TestVotes)
测试数据集上的分类准确率约为98%。混淆矩阵显示有一条CHF记录被错误分类为ARR。其他48个信号都被正确分类了。
总结
本例使用小波时间散射和支持向量机分类器将心电图波形分类为三个诊断类之一。小波散射被证明是一种强大的特征提取器，它只需要用户指定的最小参数集，就可以生成一组用于分类的健壮特征。将其与示例进行比较基于小波特征和支持向量机的信号分类这需要大量的专业知识来手工制作分类中使用的特征。使用小波时间散射，您只需要指定时不变的尺度，滤波器组(或小波变换)的数量，以及每个八度的小波数量。小波散射变换和支持向量机分类器的组合在交叉验证模型上产生了100%的分类和98%的正确分类，当将支持向量机应用到保留测试集的散射变换时。
参考文献
Anden, J.， Mallat, S. 2014。深度散射光谱，IEEE信号处理学报，62,16，第4114-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, Mietus JE, Moody GB, Peng C-K, Stanley HE。PhysioBank, PhysioToolkit和PhysioNet:复杂生理信号新研究资源的组成部分。循环．第101卷第23期，2000年6月13日，页e215-e220。
http://circ.ahajournals.org/content/101/23/e215.fullMallat。2012。集团不变的散射。通信在纯粹数学与应用数学中的应用，65,10,pp. 1331-1398。
穆迪·GB，马克·RG。MIT-BIH心律失常数据库的影响。医学与生物工程20(3):45-50(2001年5- 6月)。(PMID: 11446209)
支持功能
helperPlotRandomRecords绘制4个随机选取的心电信号ECGData．
函数helperPlotRandomRecords (ECGData randomSeed)此函数仅用于支持XpwWaveletMLExample。它可能%更改或在未来版本中删除。如果输入参数个数= = 2 rng (randomSeed)结束M =大小(ECGData.Data, 1);idxsel = randperm (M, 4);为subplot(2,2,numplot) plot(ECGData.Data(idxsel(numplot)，1:300))“伏”）如果2 xlabel(“样本”）结束标题(ECGData.Labels {idxsel (numplot)})结束结束
helperMajorityVote找到每一组散射时间窗的预测类标签中的模式。该函数返回每组散射时间窗口的类标签模式和类预测的数量。如果没有唯一模态，helperMajorityVote返回一个“error”的类标签，以指示散射窗口集是一个分类错误。
函数[ClassVotes, ClassCounts] = helperMajorityVote (predLabels、origLabels类)此函数支持ecgwaveettimescatteringexample。它可能%更改或在未来版本中删除。如果标签还不是分类的，则创建分类数组。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(size(ClassCounts,2)，1);为nc = 1:size(ClassCounts,2) modecnt(nc) = histc(ClassCounts(:，nc)，mxcount(nc));结束结束% EOF结束
另请参阅
相关的话题
您也可以从以下列表中选择网站: