心音图数据的小波时间散射分类

这个示例使用:

打开生活的脚本

这个例子展示了如何使用小波时间散射和支持向量机分类器对人类心音图(PCG)记录进行分类。金宝app心音图是心脏收缩期和舒张期声音的录音。心脏听诊在评估心脏健康方面继续发挥重要的诊断作用。不幸的是，世界上许多地区缺乏足够数量的接受过心脏听诊培训的医务人员。因此，有必要开发可靠的自动解释心音图数据的方法。

本例使用小波散射作为特征提取器进行PCG分类。在小波散射中，数据通过一系列小波变换、非线性和平均来产生数据的低方差表示。然后将这些低方差表示用作分类器的输入。这个例子是一个二进制分类问题，其中每个PCG记录要么是“正常”，要么是“异常”。

数据描述

本例使用心音图(PCG)数据从心功能正常和异常的人获得。该数据集由3829条记录组成，其中2575条来自心功能正常者，1254条来自心功能异常者。每条记录是10,000个样本长，采样在2千赫。这是5秒的心音图数据。方法中使用的训练和验证数据构造了数据集PhysioNet Computing in Cardiology Challenge 2016[1][2]。

下载数据

第一步是从GitHub库．要下载数据，请单击代码并选择下载ZIP．保存文件physionet_phonocardiogram-main.zip在有写权限的文件夹中。本例的说明假设您已经将文件下载到临时目录，(tempdir在MATLAB™)。如果您选择在不同的文件夹中下载数据，请修改后续解压缩和加载数据的说明tempdir．

该文件physionet_phonocardiogram-main.zip包含

PCG_Data.zip
README.md

和PCG_Data.zip包含

heartSoundData.mat
extrafiles.mat
Modified_physionet_data.txt
License.txt。

heartSoundData.mat保存本例中使用的数据和类标签。txt文件Modified_physionet_data.txt是PhysioNet的复制策略所需要的，它提供了数据的源属性以及每个信号如何进入的描述heartSoundData.mat对应于原始PhysioNet数据中的一个文件。extrafiles.mat还包含源文件属性，并在Modified_physionet_data.txt文件中进行解释。运行示例所需的唯一文件是heartSoundData.mat．

加载数据

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

解压缩(fullfile (tempdir,“physionet_phonocardiogram-main.zip”), tempdir)解压缩(fullfile (tempdir,“physionet_phonocardiogram-main”，“PCG_Data.zip”），．..fullfile (tempdir“PCG_Data”)）

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

负载(fullfile (tempdir“PCG_Data”，“heartSoundData.mat”)）

heartSoundData是一个有两个字段的结构数组:数据和类．数据是一个10000 × 3829矩阵，其中每一列是一个PCG记录。类是否有一个3829 × 1的诊断标签分类数组，每列对应一个数据．因为这是一个二元分类问题，分类分为“正常”和“异常”。如前所述，正常记录2575条，异常记录1254条。同样，67.25%的样本来自心功能正常的人，32.75%来自心功能异常的人。你可以输入:

总结(heartSoundData.Classes)

正常2575不正常1254

countcats (heartSoundData.Classes)。/笔(countcats (heartSoundData.Classes))

ans =2×10.6725 - 0.3275

小波散射网络

使用waveletScattering构造小波时间散射网络。设置不变比例以匹配信号长度。默认的散射网络有两个小波变换(滤波器组)。第一个小波滤波器组每八度有八个小波。第二组滤波器每八度有一个小波。设置“OptimizePath”财产真正的．

N = 1 e4;sn = waveletScattering (“SignalLength”N“InvarianceScale”N“OptimizePath”,真正的);

创建培训和测试集

辅助功能,partition_heartsounds，对3829个观测数据进行分区，使70%(2680)位于训练集中，其中1802个为正态，878个为异常。剩余的1149条记录(773条正常，376条异常)保留在测试集中进行预测。随机数生成器被放置在helper函数中，因此结果是可重复的。的代码partition_heartsounds本例中使用的所有其他辅助函数都在示例最后的支持函数一节中给出。金宝app

[trainData, testData, trainLabels, testLabels] =．..partition_heartsounds(70年,heartSoundData.Data heartSoundData.Classes);

您可以检查训练集和测试集中每个类的数量。

总结(trainLabels)

正常1802不正常878

总结(testLabels)

正常773不正常376

注意，训练集和测试集已经被划分，因此训练集和测试集中的“正常”和“异常”记录的比例与它们在整体数据中的比例相同。你可以用下面的内容来确认。

countcats (trainLabels)。/笔(countcats (trainLabels))

ans =2×10.6724 - 0.3276

countcats (testLabels)。/笔(countcats (testLabels))

ans =2×10.6728 - 0.3272

散射特性

获取训练集中所有2680条记录的散射变换。对于多变量时间序列，散射变换假设每一列都是一个单独的信号。使用“日志”选择获得散射系数的自然对数。

trainData scat_features_train = featureMatrix (sn,“转换”，“日志”）;

对于给定的散射参数，scat_features_train是一个279 × 5 × 2680的矩阵。在2680个信号中，每个信号有279个散射路径和5个散射窗口。为了将其传递给SVM分类器，将张量重塑为13400 × 279矩阵，其中每一行代表279散射路径上的单个散射窗口。总行数等于5和2680(训练数据中的记录数)的乘积。

Nseq =大小(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), []);

对测试数据重复上述过程。

testData scat_features_test = featureMatrix (sn,“转换”，“日志”）;Scat_features_test = permute(Scat_features_test，[2 3 1]);scat_features_test =重塑(scat_features_test,．..大小(scat_features_test 1) *大小(scat_features_test, 2), []);

这里我们复制标签，这样我们就有了每个散射时间窗口的标签。

[sequence_labels_train, sequence_labels_test] =．..createSequenceLabels_heartsounds (Nseq trainLabels testLabels);

将支持向量机与训练数据进行拟合。在这个例子中，我们使用三次多项式核。在将支持向量机与训练数据拟合后，我们执行5倍交叉验证来估计训练数据上的泛化误差。在这里，每个散射窗口是分开分类的。

rng默认的；classificationSVM = fitcsvm (．..scat_features_train,．..sequence_labels_train,．..“KernelFunction”，多项式的，．..“PolynomialOrder”3,．..“KernelScale”，“汽车”，．..“BoxConstraint”, 1．..“标准化”,真的,．..“类名”分类({“正常”，“异常”}));kfoldmodel = crossval (classificationSVM,“KFold”5);

以百分比计算损失并显示混淆矩阵。

predLabels = kfoldPredict (kfoldmodel);损失= kfoldLoss (kfoldmodel) * 100;流(损失是%2.2f % n'、损失);

损失是0.96%

精度= 100 -损失;流('准确度为%2.2f % \n'、准确性);

准确度为99.04%

confmatCV = confusionchart (sequence_labels_train predLabels);

请注意，当每个时间窗口单独分类时，散射网络的准确率约为99%。然而，性能实际上比这个值更好，因为我们每个记录有5个散射窗口，而99%的准确率是基于对所有窗口分别分类。在这种情况下，使用多数投票来获得每个记录的单个类分配。类投票对应于五个窗口的投票模式。如果没有找到唯一模式，helperMajorityVote将散射窗的集合分类为“NoUniqueMode”表示分类错误。这将导致混淆矩阵中多出一列。

类=分类({“异常”，“正常”}）;ClassVotes = helperMajorityVote (predLabels、trainLabels、类);CVaccuracy =总和(eq (ClassVotes trainLabels))。/元素个数(trainLabels) * 100;流(真正的交叉验证准确率是%2.2f %。, CVaccuracy);

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

显示多数投票分类的混淆矩阵。

cmCV = confusionchart (trainLabels ClassVotes);

训练数据的交叉验证准确率实际上是99.89%。有两个正常记录，被误分类为异常。一个异常记录被归类为正常记录。

使用适合训练数据的SVM模型对保留的测试数据进行类预测。

predTestLabels =预测(classificationSVM scat_features_test);

使用多数投票来确定测试集上预测的准确性。

ClassVotes = helperMajorityVote (predTestLabels、testLabels、类);testaccuracy =总和(eq (ClassVotes testLabels))。/元素个数(testLabels) * 100;流(“测试准确率为%2.2f %。”, testaccuracy);

测试准确率为91.82%。

cmt = confusionchart (testLabels ClassVotes);

在1149条测试记录中，大约92%被正确归类为“正常”或“异常”。在测试集中773条正常PCG记录中，732条被正确分类。测试集376条异常记录中，有323条分类正确。

精确度，回忆和F1分数

在一个分类任务中，一个类的精度是正确的正结果数除以正结果数。换句话说，在分类器分配给给定标签的所有记录中，实际属于类的比例是多少。Recall定义为正确的标签数除以给定类的标签数。具体来说，在属于一个类的所有记录中，分类器将其标记为那个类的比例是多少?在判断分类器的准确性时，理想情况下，您希望在精度和回忆两方面都做得很好。例如，假设我们有一个分类器，将每个PCG记录标记为异常。那么我们对异常类的召回率将是100%。所有属于异常类的记录将被标记为异常。然而，精确度会很低。因为我们的分类器将所有记录标记为异常，在这种情况下，精度为1254/3829，即32.75%，将有2575个假阳性。 The F1 score is the harmonic mean of precision and recall and provides a single metric that summarizes the classifier performance in terms of both recall and precision. The helper function,helperF1heartSounds，计算测试集中分类结果的精度、召回率和F1分数，并将这些结果返回到表中。

PRTable = helperF1heartSounds (cmTest.NormalizedValues);disp (PRTable)

精确召回F1_Score _________ ______ ________异常88.736 85.904 87.297正常93.248 94.696 93.967

在这种情况下，对于异常组和正常组的F1分数证实了我们的模型具有良好的精度和召回率。在二元分类中，直接从混淆矩阵中确定精度和召回率是很简单的。为了理解这个，为了方便，再次绘制混淆矩阵。

cmt = confusionchart (testLabels ClassVotes);

异常类的召回率是被识别为异常的异常记录的数量，即混淆矩阵第二行和第二列的条目除以第二行条目的和。异常类的精度是真异常记录占分类器识别出的异常总数的比例。它对应于混淆矩阵第二行第二列的元素除以第二列元素的和。F1分数是两者的调和平均值。

RecallAbnormal = cmTest.NormalizedValues(2, 2) /笔(cmTest.NormalizedValues (2:));PrecisionAbnormal = cmTest.NormalizedValues(2, 2) /笔(cmTest.NormalizedValues (:, 2));F1Abnormal = harm ([RecallAbnormal PrecisionAbnormal]);流('RecallAbnormal = %2.3f\nPrecisionAbnormal = %2.3f\nF1Abnormal = %2.3f\n'，．..100 * 100 * RecallAbnormal PrecisionAbnormal, 100 * F1Abnormal);

RecallAbnormal = 85.904 PrecisionAbnormal = 88.736 F1Abnormal = 87.297

在正常课堂上重复上述步骤。

RecallNormal = cmTest.NormalizedValues(1,1) /笔(cmTest.NormalizedValues (1:));PrecisionNormal = cmTest.NormalizedValues(1,1) /笔(cmTest.NormalizedValues (: 1));F1Normal = harm ([RecallNormal PrecisionNormal]);流('RecallNormal = %2.3f\nPrecisionNormal = %2.3f\nF1Normal = %2.3f\n'，．..100 * 100 * RecallNormal PrecisionNormal, 100 * F1Normal);

RecallNormal = 94.696 PrecisionNormal = 93.248

总结

本例在二值分类问题中使用小波时间散射稳健地识别人类心音图记录的正常或异常。小波散射只需要指定一个参数，即尺度不变式的长度，就可以产生PCG数据的低方差表示，从而使支持向量机分类器能够准确地建模两组数据之间的差异。金宝app尽管训练金宝app集和测试集的正常和异常PCG记录数量显著不平衡，但小波散射支持向量机分类器在精确度和召回率方面都取得了优异的成绩。

参考文献

戈德伯格，a.l.， l.a. N. Amaral, L. Glass, J. M. Hausdorff, P. Ch. Ivanov, R. G. Mark, J. E. miietus, G. B. Moody, c - k。彭，还有h·e·斯坦利。“PhysioBank, PhysioToolkit和PhysioNet:复杂生理信号新研究资源的组成部分”。循环．第101卷，第23卷，2000年6月13日，页e215-e220。http://circ.ahajournals.org/content/101/23/e215.full
刘等人。“一个用于评估心音算法的开放获取数据库”。生理测量．2016年11月21日，第37卷第12期，第2181-2213页。https://www.ncbi.nlm.nih.gov/pubmed/27869105

金宝app支持功能

partition_heartsounds创建由指定比例的数据组成的训练和测试集。该功能还保留了异常和正常PCG记录在每集的比例。

函数[trainData, testData, trainLabels, testLabels] = partition_heartsounds(percent_train_split,Data,Labels)%此函数仅支持小波时间散射金宝app心音图数据分类实例。它可能会改变%将在未来的版本中移除。心音数据中的标签不是连续的。percent_train_split = percent_train_split / 100;%每一列是一个观察值NormalData = Data(:，标签==“正常”）;AbnormalData = Data(:，Labels == .“异常”）;LabelsNormal =标签(标签== . .“正常”）;LabelsAbnormal = Labels(Labels == . .“异常”）;Nnormal =大小(NormalData, 2);Nabnormal =大小(AbnormalData, 2);num_train_normal =圆(percent_train_split * Nnormal);num_train_abnormal =圆(percent_train_split * Nabnormal);rng默认的；Pnormal = randperm (Nnormal num_train_normal);Pabnormal = randperm (Nabnormal num_train_abnormal);notPnormal = setdiff (1: Nnormal, Pnormal);notPabnormal = setdiff (1: Nabnormal, Pabnormal);trainNormalData = NormalData (:, Pnormal);trainNormalLabels = LabelsNormal (Pnormal);trainAbnormalData = AbnormalData (:, Pabnormal);trainAbnormalLabels = LabelsAbnormal (Pabnormal);testNormalData = NormalData (:, notPnormal); testNormalLabels = LabelsNormal(notPnormal); testAbnormalData = AbnormalData(:,notPabnormal); testAbnormalLabels = LabelsAbnormal(notPabnormal); trainData = [trainNormalData trainAbnormalData]; trainData = (trainData-mean(trainData))./std(trainData,1); trainLabels = [trainNormalLabels; trainAbnormalLabels]; testData = [testNormalData testAbnormalData]; testData = (testData-mean(testData))./std(testData,1); testLabels = [testNormalLabels; testAbnormalLabels];结束

createSequenceLabels_heartsounds为小波时间散射序列创建类标签。

函数[sequence_labels_train, sequence_labels_test] = createSequenceLabels_heartsounds (Nseq、trainLabels testLabels)%此函数仅支持小波时间散射金宝app心音图数据分类实例。它可能会改变%将在未来的版本中移除。Ntrain =元素个数(trainLabels);Nseq trainLabels = repmat (trainLabels ', 1);sequence_labels_train =重塑(trainLabels Nseq * Ntrain 1);元=元素个数(testLabels);Nseq testLabels = repmat (testLabels ', 1);sequence_labels_test =重塑(testLabels nt * Nseq 1);结束

helperMajorityVote基于模式实现对分类的多数投票。如果没有唯一的模式，则投票NoUniqueMode返回，以确保记录了分类错误。

函数[ClassVotes, ClassCounts] = helperMajorityVote (predLabels、origLabels类)这个函数支持小波工具箱的例子。金宝app它可能%更改或在未来的版本中删除。%如果标签不是分类的，则创建分类数组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) =分类({“NoUniqueMode”}）;ClassVotes = ClassVotes (:);%-------------------------------------------------------------------------函数modecnt = Inf(size(ClassCounts,2)，1); / / mxcount = mdecnt为c = 1:size(ClassCounts,2) modecnt(nc) = histc(ClassCounts(:，nc)，mxcount(nc));结束结束% EOF结束

helperF1heartSounds计算分类器结果的精度、召回率和F1分数。

函数PRTable = helperF1heartSounds (confmat)%此函数仅支持小波时间散射金宝app心音图数据分类实例。它可能会改变%将在未来的版本中移除。precisionAB = confmat(2, 2) /笔(confmat (:, 2)) * 100;precisionNR = ref (vol,1)/sum(vol,1) *100;recallAB = confmat(2, 2) /笔(confmat (2:)) * 100;recallNR = confmat(1,1) /笔(confmat (1:)) * 100;F1AB = 2 * (precisionAB * recallAB) / (precisionAB + recallAB);F1NR = 2 * (precisionNR * recallNR) / (precisionNR + recallNR);%构造一个MATLAB表来显示结果。PRTable = array2table([precisionAB recallAB F1AB;．..precisionNR recallNR F1NR),．..“VariableNames”, {“精度”，“回忆”，“F1_Score”}，“RowNames”，．..｛“异常”，“正常”}）;结束