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

此示例使用：

打开实时脚本

此示例演示如何使用小波时间散射和支持向量机（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
自述文件

而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.

加载文件

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

解压（完整文件）（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赫兹采样的ECG记录。每个ECG时间序列的总持续时间为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；[列车数据、测试数据、列车标签、测试标签]=．..HelperAndSplit（列车百分比，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

测试=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 = waveletScattering (“信号长度”N“InvarianceScale”,150,“SamplingFrequency”,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

图;φ= 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]，-1.5e-4 1.6e-4]，“k——”); xlabel(“秒”）;ylabel (“振幅”)；wavplt=绘图（t，[real（psiL1）imag（psiL1）]；图例（[shopplt wavplt（1）wavplt（2）]{“缩放函数”,“Wavelet-Real部分”,“小波虚部”}); 头衔({“缩放函数”;“最粗尺度小波第一滤波器组”})持有关

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

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

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

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

Nwin=尺寸（scat\u特征\u序列，2）；scat_特征_序列=排列（scat_特征_序列[2 3 1]）；scat_特征_序列=重塑（scat_特征_序列，．..尺寸（scat_特征_序列，1）*尺寸（scat_特征_序列，2），[]）；

对测试数据重复该过程。首先，scat_特性_测试是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] = createSequenceLabels (Nwin、trainLabels testLabels);

交叉验证

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

scat_功能=[scat_功能_训练；scat_功能_测试]；所有标签_scat=[sequence_标签_训练；sequence_标签_测试]；rng（1）；模板=模板SVM(．..“内核函数”,“多项式”,．..“PolynomialOrder”2,．..“内核尺度”,“自动”,．..“BoxConstraint”,1．..“标准化”，true）；classificationSVM=FitCecc(．..scat_功能，．..所有标签，．..“学习者”样板．..“编码”,“onevsone”,．..“类名”, {“啊”;瑞士法郎的;“NSR”}); kfoldmodel=crossval（分类VM，“KFold”, 5);

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

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

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

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

准确率为99.92%。

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

类别=分类({“啊”,瑞士法郎的,“NSR”}); [ClassVotes，ClassCounts]=helperMajorityVote（预标签，[trainLabels；testLabels]，类）；

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

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个信号，其中16个散射窗口中有15个分类正确。

支持向量机分类

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

模型= fitcecoc (．..scat_features_train,．..序列标签列车，．..“学习者”样板．..“编码”,“onevsone”,．..“类名”, {“啊”,瑞士法郎的,“NSR”}); predLabels=预测（模型、scat\u特征\u测试）；[TestVotes，TestCounts]=HelperMajorityNote（预标签、测试标签、类）；testaccuracy=sum（eq（测试投票，分类（测试标签）））/numel（测试标签）*100；fprintf('测试精度为%2.2f%\不, testaccuracy);

测试准确率为97.96%。

混淆图（分类（测试标签），测试投票）

测试数据集的分类精度约为98%。混淆矩阵显示，一个CHF记录被错误分类为ARR。所有其他48个信号都被正确分类。

总结

这个例子使用小波时间散射和支持向量机分类器将心电波形分类为三个诊断类之一。小波散射被证明是一种强大的特征提取方法，它只需要用户指定的最小参数集就可以产生一组鲁棒的特征用于分类。将这个与示例进行比较基于小波特征和支持向量机的信号分类金宝app，这需要大量的专业知识来手工制作用于分类的特征。使用小波时间散射，您只需要指定时间不变性的比例、滤波器组（或小波变换）的数量，以及每倍频程的小波数。将小波散射变换与SVM分类器相结合，可在交叉验证模型上实现100%的分类，并在将SVM应用于保持测试集的散射变换时实现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, miietus 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.full
Mallat，S.，2012。群不变散射。纯数学和应用数学中的通信，65,10，第1331-1398页。
Moody GB，Mark RG.《麻省理工学院波黑心律失常数据库的影响》。医学和生物学IEEE Eng 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子地块（2,2，numplot）图（ECGData.Data（idxsel（numplot），1:3000））ylabel(“伏特”)如果Numplot > 2 xlabel(“样本”)终止标题（ECGData.Labels{idxsel（numplot）}）终止终止

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

作用[classvots，ClassCounts]=helperMajorityVote（预标签，或标签，类）%此函数支持ECGWaveletTimeSc金宝appatteringExample。它可能%在将来的版本中更改或删除。%如果标签不是分类的，则创建分类数组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终止