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

此示例使用：

打开生活的脚本

这个例子展示了如何使用小波时间散射和支持向量机分类器对人体心电图信号进行分类。金宝app在小波散射中，数据通过一系列小波变换、非线性和平均来产生时间序列的低方差表示。小波时间散射产生的信号表示对输入信号的移位不敏感，而不牺牲类的可辨别性。要运行此示例，必须有小波工具箱™和统计和机器学习工具箱™。本例中使用的数据可以从以下网站公开获得生理网．在这个例子中，你可以找到一种深度学习的方法来解决这个分类问题使用小波分析和深度学习分类时间序列和这个例子中的机器学习方法基于小波特征和支持向量机的信号分类金宝app．

数据描述

本例使用从三组或三类人群中获得的ECG数据：心律失常患者、充血性心力衰竭患者和窦性心律正常患者。该示例使用三个PhysioNet数据库中的162个ECG记录：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.xt.
License.txt。

EcgData.mat保存此示例中使用的数据。Physionet的复制策略需要.txt文件modified_physionet_data.txt，并为数据提供了数据的源属性以及应用于每个ECG录制的预处理步骤的描述。

加载文件

如果您按照前一节中的下载说明操作，请输入以下命令来解压缩这两个归档文件。

解压（完整文件）（tempdir，'physionet_ecg_data-main.zip')，tempdir）解压缩（完整文件（tempdir，“physionet\u ECG\u data-main”，“ECGData.zip”），......fullfile (tempdir“埃格达塔”)）

解压ECGData.zip文件后，将数据加载到MATLAB中。

负载(fullfile (tempdir“埃格达塔”，'ecgdata.mat')）

EcgData.是一个具有两个字段的结构阵列：数据和标签．数据是一个162乘65536的矩阵，其中每一行是以128赫兹采样的心电图记录。每个心电时间序列的总持续时间为512秒。标签是一个162 × 1的细胞阵列诊断标签，每一行一个数据这个three diagnostic categories are: 'ARR' (arrhythmia), 'CHF' (congestive heart failure), and 'NSR' (normal sinus rhythm).

创建培训和测试数据

随机将数据分成两组 - 培训和测试数据集。辅助功能Helperrandomsplit.执行随机分割。Helperrandomsplit.接受培训数据的所需分割百分比EcgData.这个Helperrandomsplit.功能输出两个数据集以及每个数据集。每一排trainData和testData是心电信号。的每个元素列车标签和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（分类（列车标签））/numel（列车标签）。*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个小波，第二个滤波器组每倍频程有1个小波。不变性比例设置为150秒。

n =大小（EcgData.data，2）;Sn =小波示踪剂（“SignalLength”，n，'invariarcescale',150,'采样频率', 128);

你可以用下面的方法在两个滤波器组中可视化小波滤波器。

[fb，f，filterparams]=filterbank（sn）；图子地块（211）图（f，fb{2}.psift）xlim（[0 128]）网格在头衔('第一个过滤器银行小波过滤器'）;Subplot (212) plot(f,fb{3}.psift) xlim([0 128]) grid在头衔('第二滤波器银行小波过滤器'）;包含(“赫兹”）;

为了演示不变性尺度，获取尺度函数的傅里叶反变换，并将其居中于0秒时间。两条垂直的黑线分别是-75秒和75秒的分界线。另外，从第一个滤波器组中绘制出最粗尺度(最低频率)小波的实部和虚部。注意，最粗尺度小波不会超过由尺度函数的时间支持所决定的不变尺度。金宝app这是小波时间散射的一个重要性质。

数字;PHI = IFFTSHIFT（IFFT（FB {1} .FLIFT））;psil1 = ifftshift（ifft（fb {2} .psift（：，结束）））;t =（-2 ^ 15：2 ^ 15-1）。* 1/128;scalplt = plot（t，phi）;抓住在网格在ylim(1.6[-1.5军医的军医]);情节(-75[-75],[-1.5 1.6的军医]的军医,“k——”）;图（[75 75]，[ -  1.5E-4 1.6E-4]，“k——”）;包含(“秒”）;ylabel（“振幅”）;wavplt = plot(t，[real(psiL1) imag(psiL1)]);传奇([scalplt wavplt (1) wavplt (2)), {“扩展功能”，'小波真的部分'，“小波虚部”}）;标题({“扩展功能”；“最粗尺度小波第一滤波器组”})举行从

构建散射网络后，获取训练数据的散射系数作为矩阵featureMatrix对于多个信号，每一列都被视为单个信号。

scat_特征_列车=特征矩阵（序列号，列车数据’）；

的输出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_特征_序列，1）*尺寸（scat_特征_序列，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] = createSquenceLabels（nwin，trainlabels，testlabels）;

交叉验证

对于分类，执行了两种分析。第一种分析使用所有散射数据并拟合具有二次核的多类SVM。整个数据集中总共有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，......“内核尺度”，“汽车”，......'boxconstraint'，1，......“标准化”,真正的);classificationSVM = fitcecoc (......scat_功能，......allLabels_scat,......“学习者”模板,......“编码”，“onevsone”，......“类名”, {“加勒比海盗”；'chf'；“签约”}); kfoldmodel=crossval（分类VM，'kfold'5);

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

predLabels=kfoldPredict（kfoldmodel）；loss=kfoldLoss（kfoldmodel）*100；confmatCV=confusionmat（所有标签，predLabels）

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

流('准确率为%2.2f%。\n'，损失100元）；

准确率为99.92%。

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

类=分类({“加勒比海盗”，'chf'，“签约”}）;[ClassVotes, ClassCounts] = helperMajorityVote (predLabels [trainLabels;testLabels)、类);

根据每组散射时间窗口的类别预测模式确定实际交叉验证精度。如果该模式对于给定集不是唯一的，则helperMajorityVote返回由指示的分类错误'nounquemode'．这将导致混淆矩阵中额外的一列，在这种情况下，这一列都是零，因为对于每一组散射预测，都存在一个独特的模式。

CVaccuracy =总和(eq (ClassVotes分类([trainLabels;testLabels]))) / 162 * 100;流(“真正的交叉验证准确率为%2.2f %。”, CVaccuracy);

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

MVconfmatCV=confusionmat（分类（[trainLabels；testLabels]），类投票；MVconfmatCV

mvconfmatcv =4×496 0 0 0 30 0 0 0 36 0 0 0 0 0

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

支持向量机分类

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

Model = FitCecoc（......scat_features_train,......sequence_labels_train，......“学习者”模板,......“编码”，“onevsone”，......“类名”, {“加勒比海盗”，'chf'，“签约”}）;predLabels =预测(模型、scat_features_test);[TestVotes, TestCounts] = helperMajorityVote (predLabels、testLabels、类);testaccuracy =总和(eq (TestVotes分类(testLabels))) /元素个数(testLabels) * 100;流(测试准确率为%2.2f %。\ n”，testaguracy）;

试验准确率为97.96%。

confusionchart(分类(testLabels)、TestVotes)

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

总结

此示例使用小波时间散射和SVM分类器将ECG波形分类为三个诊断类中的一个。小波散射被证明是一个强大的特征提取器，它只需要最小的用户指定参数集，以产生一组用于分类的强大功能。将此与示例进行比较基于小波特征和支持向量机的信号分类金宝app，这需要大量的专业知识来手工制作特征，以用于分类。对于小波时间散射，您只需要指定时不变性的规模，滤波器组(或小波变换)的数量，以及每个八度的小波数量。小波散射变换和支持向量机分类器的结合在交叉验证模型上产生了100%的分类结果，当将支持向量机应用到持度测试集的散射变换上时，分类的正确率达到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.用口服Milrinone治疗严重充血性心力衰竭患者的生存。J美国心脏病学院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，PhysioLoolkit和PhysioIoneet：复杂生理信号的新研究资源组件。循环．第101卷，第23卷，2000年6月13日，页e215-e220。http://circ.ahajournals.org/content/101/23/e215.full.
Mallat。2012。集团不变的散射。《纯粹与应用数学通讯》，65,10,1331-1398页。
穆迪GB，马克RG。MIT-BIH心律失常数据库的影响。IEEE医学与生物学工程20(3):45-50(2001年5- 6月)。(PMID: 11446209)

金宝app支持功能

HelperplotrandomRecords.绘制从中随机选择的四个ECG信号EcgData.．

函数helperPlotRandomRecords (ECGData randomSeed)%此函数仅用于支持XpwWaveletMLExample。它可能金宝app%更改或在未来的版本中删除。如果nargin==2 rng（随机种子）结束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返回一个类标签“error”，表示散射窗口集是一个分类错误。

函数[classvots，ClassCounts]=helperMajorityVote（预标签，或标签，类）%此函数支持ECGWaveletTimeSc金宝appatteringExample。它可能%更改或在未来的版本中删除。如果标签尚未分类，％制作分类阵列predLabels =分类(predLabels);origLabels =分类(origLabels);%要求predLabels和origLabels都是分类向量Npred=numel（predLabels）；Norig=numel（origLabels）；Nwin=Npred/Norig；predLabels=Reformate（predLabels，Nwin，Norig）；ClassCounts=COUNTCTS（predLabels）；[mxcount，idx]=max（ClassCounts）；CLASSVOUTS=classes（idx）；％检查最大值中的任何关系，并确保它们标记为%如果模式出现多次，则出错modecnt=modecount（类计数，mxcount）；类投票（modecnt>1）=分类({“错误”}）;ClassVotes = ClassVotes (:);％--------------------------------------------------------------------函数modecnt = Inf(size(ClassCounts,2)，1); / / mxcount = mdecnt为nc=1:size（ClassCounts，2）modecnt（nc）=histc（ClassCounts（：，nc），mxcount（nc））；结束结束% EOF结束