小波时间散射用于ECG信号分类

此示例使用：

打开直播脚本

该示例显示如何使用小波时间散射和支持向量机（SVM）分类器对人的心电图（ECG）信号进行分类。金宝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中。

load（fullfile（tempdir，“埃格达塔”那'ecgdata.mat'）））

EcgData.是一个具有两个字段的结构阵列：资料和标签。资料是一个162×65536矩阵，其中每行是在128赫兹采样的心电图录制。每个ECG时间序列都有512秒的总持续时间。标签是一个162×1个单元格阵列诊断标签，每排一个资料. 这个three diagnostic categories are: 'ARR' (arrhythmia), 'CHF' (congestive heart failure), and 'NSR' (normal sinus rhythm).

创建培训和测试数据

随机将数据分成两组 - 培训和测试数据集。辅助功能Helperrandomsplit.执行随机拆分。Helperrandomsplit.接受培训数据的所需分割百分比EcgData.. 这个Helperrandomsplit.功能输出两个数据集以及每个数据集。每一排训练DATA.和测试数据是一个心电图信号。每个元素列车标签和testlabels.包含数据矩阵的相应行的类标签。在此示例中，我们随机将每个类中的70％分配给训练集。剩下的30％被列出用于测试（预测），并分配给测试集。

百分之= 70;[TrainData，TestData，TrainLabels，TestLabels] =......Helperrandomsplit（百分号，EcgData）;

有113条记录训练DATA.设置和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（分类（testlabels））./ numel（testlabels）。* 100

ctest =3×159.1837 18.3673 22.4490

绘制样本

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

HelperplotrandomRecords（Ecgdata，14）

小波时间散射

在小波时间散射网络中指定的关键参数是时间不变的比例，小波变换的数量，以及每个小波滤波器组中的每个八度音程的数量。在许多应用中，两个过滤器组的级联足以实现良好的性能。在该示例中，我们构建一个小波时间散射网络，其中默认滤波器组：在第二滤波器组中的第一个滤波器组中的每个倍频度的8个小波，在第二滤波器组中每个倍频值。不变性刻度设置为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秒的边界。另外，从第一滤波器组绘制粗略级（最低频率）小波的实部和虚部。请注意，粗略级小波不会超过由缩放功能的时间支持确定的不变尺度。金宝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.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 = plot（t，[真实（psil1）imag（psil1）]）;图例（[scalplt wavplt（1）波普（2）]，{'缩放功能'那'小波真的部分'那“小波虚部”}）;标题（{'缩放功能';“最粗尺度小波第一滤波器组”}） 抓住离开

构造散射网络后，将训练数据的散射系数作为矩阵进行计算。当你跑的时候Featurematrix.对于多个信号，每列被视为一个信号。

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

输出Featurematrix.在这种情况下，409×16-113。张量的每个页面，scat_features_train.，是一个信号的散射变换。基于缩放功能的带宽，小波散射变换批判性地下采样。在这种情况下，这导致409个散射路径中的每一个的16个时间窗口。

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

nwin = size（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×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）;

交叉验证

对于分类，进行了两次分析。第一种方法使用所有散射数据，并拟合具有二次核的多类支持向量机。整个数据集中总共有2592个散射序列，162个信号中的每一个都有16个。使用5倍交叉验证估计错误率或损失。

scat_features = [scat_features_train;scat_features_test];alllabels_scat = [sequence_labels_train;sequence_labels_test];RNG（1）;template = templatesvm（......'骨箱'那'多项式'那......'polynomialOrder'2，......“内核尺度”那'汽车'那......'boxconstraint'，1，......“标准化”， 真的）;Classificationsvm = fitcecoc（......scat_功能，......AllLabels_Scat，......'学习者'， 模板，......'编码'那'Onevsone'那......“类名”，{'arr';'chf';'nsr'}); kfoldmodel=crossval（分类VM，'kfold'5）;

计算丢失和混淆矩阵。显示精度。

predLabels=kfoldPredict（kfoldmodel）；损失=kfoldLoss（kfoldmodel）*100；confmatCV=confusionmat（所有标签、预标签）

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

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

准确率为99.92%。

准确性为99.88％，这是相当不错的，但实际结果更好地考虑到这里每次窗口分别分类。每个信号有16个单独的分类。使用简单的多数投票来获取每个散射表示的单级预测。

classes = patporical（{'arr'那'chf'那'nsr'}）;[classvotes，classcounts] = HelpermajorityVote（Predlabels，[TrainLabels; Testlabels]，课程）;

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

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

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

MVconfmatCV=confusionmat（分类（[trainLabels；测试标签），类别投票）；MVconfmatCV

mvconfmatcv =4×4.96 0 0 0 0 0 0 0 0 0 0 36 0 0 0 0 0

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

支持向量机分类

对于下一次分析，我们仅适用于培训数据的多级二次SVM（70％），然后使用该模型来对所删除的30％的数据进行预测。测试集中有49个数据记录。在各个散射窗口上使用大多数投票。

Model = FitCecoc（......scat_features_train，......sequence_labels_train，......'学习者'， 模板，......'编码'那'Onevsone'那......“类名”，{'arr'那'chf'那'nsr'}）;predlabels =预测（模型，scat_features_test）;[testvotes，testcounts] = HelpermajorityVote（Predlabels，Testlabels，课程）;testAccuracy = sum（eq（testvotes，分类（testlabels））/ numel（testlabels）* 100;fprintf（'测试精度为％2.2f％。\ n'，testaguracy）;

测试精度为97.96％。

ConfusionChart（分类（testlabels），testvotes）

测试数据集的分类准确性约为98％。困惑矩阵表明，一个CHF记录被错误分类为ARR。所有48个其他信号都是正确分类的。

总结

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

工具书类

Anden，J.，Mallat，S. 2014.深度散射谱，信号处理的IEEE交易，62,16，PP。4114-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,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）

金宝app支持功能

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

功能HelperplotrandomRecords（Ecgdata，随机）%此函数仅用于支持XpwWaveletMLExample。可能金宝app％更改或在将来的释放中删除。如果nargin==2 rng（随机种子）结尾M=尺寸（ECGData.Data，1）；idxsel=randperm（M，4）；为了numplot = 1：4子图（2,2，numplot）绘图（EcgData.data（Idxsel（Numplot），1：3000））Ylabel（“伏特”）如果numplot> 2 xlabel（'样品'）结尾标题（EcgData.Labels {idxsel（numplot）}）结尾结尾

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

功能[classvots，ClassCounts]=helperMajorityVote（预标签，或标签，类）%此函数支持ECGWaveletTimeSc金宝appatteringExample。可能％更改或在将来的释放中删除。如果标签尚未分类，％制作分类阵列predlabels =分类（predlabels）;origlabels =分类（origlabels）;％期望Predlabels和Origlabels成为分类向量Npred=numel（预标签）；Norig=努梅尔（origLabels）；Nwin=Npred/Norig；predLabels=重塑（predLabels、Nwin、Norig）；ClassCounts=COUNTCTS（预标签）[mxcount，idx]=最大（类计数）；类别投票=类别（idx）；％检查最大值中的任何关系，并确保它们标记为%如果模式出现多次，则出错modecnt=modecount（ClassCounts，mxcount）；类别投票（modecnt>1）=分类({“错误”}）;classvotes = classvotes（:);％--------------------------------------------------------------------功能modecnt = modecount（classcounts，mxcount）modecnt = Inf（大小（classcounts，2），1）;为了nc=1:size（ClassCounts，2）modecnt（nc）=histc（ClassCounts（：，nc），mxcount（nc））；结尾结尾％eof.结尾