主要内容

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

这个例子使用了:

打开实时脚本

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

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

术语说明:在小波散射的背景下，术语“时间窗口”指的是对平滑操作的输出进行下采样后获得的样本数量。有关更多信息，请参见时间窗口．

数据描述

本例使用心音图(PCG)数据，这些数据来自心功能正常和异常的人。该数据集包括3829份记录，其中2575份来自心功能正常的人，1254份来自心功能异常的人。每个录音有10,000个样本长，并以2 kHz采样。这表示5秒的心音图数据。数据集是由训练和验证数据构建的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是具有两个字段的结构数组:数据而且类．数据是一个10,000 × 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 = 1e4;小波散射(“SignalLength”N“InvarianceScale”N“OptimizePath”,真正的);

创建训练和测试集

辅助函数，partition_heartsounds，将3829个观测值进行划分，使70%(2680)处于训练集中，其中1802个正常，878个异常。剩余的1149条记录(773条正常记录和376条异常记录)被保留在测试集中用于预测。随机数生成器被植入helper函数内部，因此结果是可重复的。的代码partition_heartsounds和本例中使用的所有其他辅助函数在示例末尾的支持函数部分中给出。万博1manbetx

[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个录音的散射变换。对于多元时间序列，散射变换假设每一列是一个单独的信号。使用“日志”选项，以获得散射系数的自然对数。

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

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

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

对测试数据重复此过程。

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), []);

在这里，我们复制标签，以便每个散射时间窗口都有一个标签。

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

将SVM与训练数据进行拟合。在这个例子中，我们使用了一个三次多项式核。在将SVM拟合到训练数据后，我们执行5倍交叉验证来估计训练数据的泛化误差。这里对每个散射窗口分别进行分类。

rng默认的；分类svm = 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-loss;流('精度为%2.2f % \n'、准确性);
准确率为99.04%
confmatCV = confusichart (sequence_labels_train,predLabels);

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

Classes = categorical({“异常”，“正常”});ClassVotes = helperMajorityVote(predLabels,trainLabels,classes);CVaccuracy = sum(eq(ClassVotes,trainLabels))./ nummel (trainLabels)*100;流(“真正的交叉验证准确率是%2.2f %。”, CVaccuracy);
真正的交叉验证准确率为99.89%。

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

cmCV = confusionchart(trainLabels,ClassVotes);

训练数据的交叉验证准确率实际上为99.89%。有两条正常的记录被错误的归类为异常。一条异常记录被判定为正常记录。

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

predTestLabels = predict(classificationSVM,scat_features_test);

使用多数票确定测试集中预测的准确性。

ClassVotes = helperMajorityVote(predTestLabels,testLabels,classes);testaccuracy = sum(eq(ClassVotes,testLabels))./numel(testLabels)*100;流('测试精度为%2.2f %。', testaccuracy);
测试准确率为91.82%。
cmTest = confusionchart(testLabels,ClassVotes);

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

精确度，召回率和F1分数

在分类任务中，一个类的精度是正确的阳性结果数除以阳性结果数。换句话说，在分类器分配给给定标签的所有记录中，实际属于该类的比例是多少。召回被定义为正确标签的数量除以给定类的标签数量。具体来说，在属于一个类的所有记录中，分类器标记为该类的比例是多少。在判断分类器的准确性时，理想情况下，您希望在精度和召回率方面都做得很好。例如，假设我们有一个分类器，将每个PCG记录标记为异常。那么我们对异常类的召回率将是100%。属于异常类的所有记录将被标记为异常。然而，精度会很低。由于我们的分类器将所有记录标记为异常，因此在这种情况下会有2575个假阳性，精度为1254/3829，即32.75%。 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分数证实了我们的模型具有良好的精度和召回率。在二元分类中，直接从混淆矩阵中确定精度和查全率是很简单的。为了方便起见，再次绘制混淆矩阵。

cmTest = confusionchart(testLabels,ClassVotes);

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

RecallAbnormal = cmTest.NormalizedValues(2,2)/sum(cmTest.NormalizedValues(2，:));PrecisionAbnormal = cmTest.NormalizedValues(2,2)/sum(cmTest.NormalizedValues(:，2));F1Abnormal = harmmean([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)/sum(cmTest.NormalizedValues(1，:));PrecisionNormal = cmTest.NormalizedValues(1,1)/sum(cmTest.NormalizedValues(:，1));F1Normal = harmmean([RecallNormal PrecisionNormal]);流('RecallNormal = %2.3f\nPrecisionNormal = %2.3f\nF1Normal = %2.3f\n'，.．.100 * 100 * RecallNormal PrecisionNormal, 100 * F1Normal);
RecallNormal = 94.696 PrecisionNormal = 93.248 F1Normal = 93.967

总结

本例使用小波时间散射在二元分类问题中稳健地识别人类心音图记录为正常或异常。小波散射只需要指定一个参数，即尺度不变的长度，以产生PCG数据的低方差表示，从而使支持向量机分类器能够准确地模拟两组之间的差异。万博1manbetx小波散射万博1manbetx支持向量机分类器在训练集和测试集正常和异常PCG记录数量显著不平衡的情况下，在精度和查全率方面都能获得更好的表现。

参考文献

  1. 戈德伯格，a.l.， L. A. N.阿马拉尔，L.格拉斯，J. M.豪斯多夫，P. Ch.伊万诺夫，R. G.马克，J. E.米耶图斯，G. B.穆迪，c . k。彭先生和斯坦利先生。PhysioBank, PhysioToolkit，和PhysioNet:复杂生理信号新研究资源的组成部分循环．第101卷，第23期，2000年6月13日，第e215-e220页。http://circ.ahajournals.org/content/101/23/e215.full

  2. Liu等。“用于评估心音算法的开放访问数据库”。生理测量．第37卷，第12期，2016年11月21日，第2181-2213页。https://www.ncbi.nlm.nih.gov/pubmed/27869105

万博1manbetx支持功能

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

函数[trainData, testData, trainLabels, testLabels] = partition_heartsounds(percent_train_split,Data,Labels)此函数仅支持小波时间散射万博1manbetx%心音图数据分类示例。它可能改变，也可能是%将在将来的版本中删除。心音数据中的标签不是连续的。Percent_train_split = Percent_train_split /100;每一列都是一个观察值NormalData = Data(:，Labels ==“正常”）;AbnormalData = Data(:，Labels ==“异常”）;LabelsNormal =标签(标签==“正常”）;LabelsAbnormal = Labels(标签==“异常”）;Nnormal = size(NormalData,2);Nabnormal = size(AbnormalData,2);num_train_normal = round(percent_train_split*Nnormal);num_train_abnormal = round(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 =标签正常(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)此函数仅支持小波时间散射万博1manbetx%心音图数据分类示例。它可能改变，也可能是%将在将来的版本中删除。Ntrain = number (trainLabels);trainLabels = repmat(trainLabels'，Nseq,1);sequence_labels_train =重塑(trainLabels,Nseq*Ntrain,1);Ntest =数字(testLabels);testLabels = repmat(testLabels'，Nseq,1);sequence_labels_test =重塑(testLabels,Ntest*Nseq,1);结束

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

函数[ClassVotes,ClassCounts] = helperMajorityVote(predLabels,origLabels,classes)此函数支持小波工具箱示例。万博1manbetx它可能%更改或在将来的版本中删除。如果标签不是分类的，则创建分类数组predLabels = categorical(predLabels);origLabels = categorical(origLabels);%期望predLabels和origLabels都是分类向量Npred =数字(predLabels);Norig =数字(origLabels);Nwin = Npred/Norig;predLabels =重塑(predLabels,Nwin,Norig);ClassCounts = countcats(predLabels);[mxcount,idx] = max(ClassCounts);ClassVotes = classes(idx);检查最大值中是否有任何联系，并确保它们被标记为如果该模式出现多次，则错误%modecnt = modecount(ClassCounts,mxcount);ClassVotes(modecnt>1) = categorical({“NoUniqueMode”});ClassVotes = ClassVotes(:);%-------------------------------------------------------------------------函数modecnt = modecount(ClassCounts,mxcount) modecnt = Inf(size(ClassCounts,2)，1);为nc = 1:size(ClassCounts,2) modecnt(nc) = histc(ClassCounts(:，nc)，mxcount(nc));结束结束% EOF结束

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

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

另请参阅

相关的例子

更多关于