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

此示例使用：

打开生活的脚本

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

数据描述

本例使用从三组或三类人群中获得的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中。

加载（完整文件）（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.功能输出两个数据集以及每个数据集。每一排列车数据和测试数据是一个ECG信号。每个元素列车标签和testlabels.包含数据矩阵的相应行的类标签。在此示例中，我们随机将每个类中的70％分配给训练集。剩下的30％被列出用于测试（预测），并分配给测试集。

列车百分比=70；[列车数据、测试数据、列车标签、测试标签]=......helperRandomSplit (percent_train ECGData);

有113条记录列车数据设置和49条记录测试数据。根据设计，培训数据包含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

Ctest=countcats（分类（测试标签））./numel（测试标签）。*100

ct =3×159.1837 18.3673 22.4490

绘制样本

绘制四个随机选择的记录的前几千个样本EcgData..辅助功能HelperplotrandomRecords.做这个。HelperplotrandomRecords.接受EcgData.和一个随机种子作为输入。初始种子设置为14，以便每个类至少绘制一条记录。您可以执行HelperplotrandomRecords.具有EcgData.作为唯一的输入参数，您可以多次了解与每个类相关的各种ECG波形。您可以在本示例末尾的“支持函数”部分中找到此帮助函数和所有帮助函数的源代码。万博1manbetx

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]）网格在头衔('第一个过滤器银行小波过滤器')；子地块（212）图（f，fb{3}.psift）xlim（[0 128]）网格在头衔('第二滤波器银行小波过滤器'); 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.5E-4 1.6E-4]，“k——”); xlabel(“秒”）;ylabel（“振幅”)；wavplt=绘图（t，[real（psiL1）imag（psiL1）]；图例（[shopplt wavplt（1）wavplt（2）]{“缩放函数”,'小波真的部分',“小波虚部”}）;标题({“缩放函数”;“最粗尺度小波第一滤波器组”})举行从

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

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

产量特征矩阵在这种情况下，409×16-113。张量的每个页面，scat_特色列车，是一个信号的散射变换。基于缩放功能的带宽，小波散射变换批判性地下采样。在这种情况下，这导致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-by-16-by-49，因为训练集中有49个ECG波形。对SVM分类器进行重塑后，特征矩阵是784-by-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_功能，......所有标签，......“学习者”模板,......“编码”,“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

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

准确率为99.92%。

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

类别=分类({“啊”,'chf',“签约”}); [ClassVotes，ClassCounts]=helperMajorityVote（预标签，[trainLabels；testLabels]，类）；

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

CVAccurance=总和（等式（类投票，分类（[trainLabels；testLabels]））/162*100；fprintf('真正的交叉验证精度为%2.2f%。\n'，准确度）；

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

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

mvconfmatcv =4×496 0 0 0 0 30 0 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;fprintf('测试精度为%2.2f%\不，testaguracy）;

测试准确率为97.96%。

confusionchart(分类(testLabels)、TestVotes)

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

总结

此示例使用小波时间散射和SVM分类器将ECG波形分类为三个诊断类中的一个。小波散射被证明是一个强大的特征提取器，它只需要最小的用户指定参数集，以产生一组用于分类的强大功能。将此与示例进行比较基于小波特征和支持向量机的信号分类万博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.用口服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，S.，2012。群不变散射。纯数学和应用数学中的通信，65,10，第1331-1398页。
Moody GB，Mark RG.《麻省理工学院波黑心律失常数据库的影响》。医学和生物学IEEE Eng 20（3）：45-50（2001年5-6月）。（PMID:11446209）

万博1manbetx支持功能

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

函数helperPlotRandomRecords (ECGData randomSeed)%此函数仅用于支持XpwWaveletMLExample。它可能万博1manbetx%更改或在未来的版本中删除。如果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)})终止终止

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

函数[classvots，ClassCounts]=helperMajorityVote（预标签，或标签，类）%此函数支持ECGWaveletTimeSc万博1manbetxatteringExample。它可能%更改或在未来的版本中删除。如果标签尚未分类，％制作分类阵列predLabels =分类(predLabels);origLabels =分类(origLabels);%预期PredLabel和OrigLabel都是分类向量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）=分类({“错误”})；类别投票=类别投票（：）；％--------------------------------------------------------------------函数modecnt=modecount（ClassCounts，mxcount）modecnt=Inf（size（ClassCounts，2），1）；对于nc=1:size（ClassCounts，2）modecnt（nc）=histc（ClassCounts（：，nc），mxcount（nc））；终止终止%EOF终止