这个例子展示了如何使用小波时间散射和支持向量机分类器对人类心音图(PCG)记录进行分类。万博1manbetx心音图是心脏收缩期和舒张期声音的录音。心脏听诊在评估心脏健康方面继续发挥重要的诊断作用。不幸的是,世界上许多地区缺乏足够数量的接受过心脏听诊培训的医务人员。因此,有必要开发可靠的自动解释心音图数据的方法。
本例使用小波散射作为特征提取器进行PCG分类。在小波散射中,数据通过一系列小波变换、非线性和平均来产生数据的低方差表示。然后将这些低方差表示用作分类器的输入。这个例子是一个二进制分类问题,其中每个PCG记录要么是“正常”,要么是“异常”。
本例使用心音图(PCG)数据从心功能正常和异常的人获得。该数据集由3829条记录组成,其中2575条来自心功能正常者,1254条来自心功能异常者。每条记录是10,000个样本长,采样在2千赫。这是5秒的心音图数据。方法中使用的训练和验证数据构造了数据集PhysioNet Computing in Cardiology Challenge 2016[1][2]。
第一步是从github存储库.要下载数据,请单击代码
并选择下载ZIP
.保存文件physionet_phonocardinogram-main.zip.
在您具有写入权限的文件夹中。此示例的说明假定您已将文件下载到临时目录(tempdir
在matlab™中)。修改后续说明解压缩并加载数据如果您选择在不同的文件夹中下载数据tempdir
.
该文件physionet_phonocardinogram-main.zip.
包含
pcg_data.zip.zip.zip.
readme.md.
和PCG_Data.zip包含
heartSoundData.mat
Uttrafiles.mat
Modified_physionet_data.txt
License.txt。
heartSoundData.mat
保存本例中使用的数据和类标签。txt文件Modified_physionet_data.txt是PhysioNet的复制策略所需要的,它提供了数据的源属性以及每个信号如何进入的描述heartSoundData.mat
对应于原始PhysioNet数据中的一个文件。Uttrafiles.mat
还包含源文件属性,并在Modified_physionet_data.txt文件中进行解释。运行示例所需的唯一文件是heartSoundData.mat
.
如果您遵循上一节中的下载说明,请输入以下命令以解压缩两个存档文件。
解压缩(fullfile (tempdir,“physionet_phonocardiogram-main.zip”), tempdir)解压缩(fullfile (tempdir,'physoonet_phonocardinogram-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个异常)在测试集中举出预测。随机数发生器接种在辅助功能的内部,因此结果是可重复的。代码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)./ sum(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), []);
对测试数据重复上述过程。
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,...'骨箱','多项式',...“PolynomialOrder”3,...“KernelScale”,“汽车”,...“BoxConstraint”, 1...“标准化”,真的,...“类名”分类({“正常”,“异常”}));kfoldmodel = crossval(分类vm,“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%。\ n', CVaccuracy);
真正的交叉验证准确率为99.89%。
显示多数投票分类的混淆矩阵。
cmCV = confusionchart (trainLabels ClassVotes);
训练数据的交叉验证准确率实际上是99.89%。有两个正常记录,被误分类为异常。一个异常记录被归类为正常记录。
使用SVM模型适合培训数据,使CALL-OUT测试数据进行类预测。
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条分类正确。
在一个分类任务中,一个类的精度是正确的正结果数除以正结果数。换句话说,在分类器分配给给定标签的所有记录中,实际属于类的比例是多少。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.NORMALIZEVALUES);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]);流('recallabnnormal =%2.3f \ n precisionabnormal =%2.3f \ nf1abnnormal =%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数据的低方差表示,从而使支持向量机分类器能够准确地建模两组数据之间的差异。万博1manbetx尽管训练万博1manbetx集和测试集的正常和异常PCG记录数量显著不平衡,但小波散射支持向量机分类器在精确度和召回率方面都取得了优异的成绩。
Goldberger,A.L.,L. A. N.Maral,L.玻璃,J.M.Hausdorff,P. Ch。Ivanov,R.G.Mark,J.E.Mietus,G. B.Coody,C. -K。彭和H. E. Stanley。“Physiobank,PhysioLoolkit和PhysioIoneet:复杂生理信号的新研究资源的组成部分”。循环.第101卷,第23卷,2000年6月13日,页e215-e220。http://circ.ahajournals.org/content/101/23/e215.full
刘等人。“一个用于评估心音算法的开放获取数据库”。生理测量.卷。37,2016年11月21日,第12届,第2181-2213号。https://www.ncbi.nlm.nih.gov/pubmed/27869105
partition_heartsounds
创建由指定的数据的比例组成的培训和测试集。该功能还保留了每个集合中的异常和正常PCG录制的比例。
功能[trainData, testData, trainLabels, testLabels] = partition_heartsounds(percent_train_split,Data,Labels)%此功能仅支持小波时间散射万博1manbetx心音图数据分类实例。它可能会改变%将在未来的版本中移除。心音数据中的标签不是连续的。percient_train_split = percent_train_split / 100;%每一列是一个观察值NormalData = Data(:,标签==“正常”);AbnormalData = Data(:,Labels == .“异常”);LabelSnormal =标签(标签==“正常”);LabelsAbnormal = Labels(Labels == . .“异常”);nnormal = size(normaldata,2);Nabnormal =尺寸(异常异常,2);num_train_normal = round(percient_train_split * nnormal);num_train_abnnormal = round(percient_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)%此功能仅支持小波时间散射万博1manbetx心音图数据分类实例。它可能会改变%将在未来的版本中移除。ntrain = numel(trainlabels);Trainlabels = Repmat(Trainlabels',NSeq,1);sequence_labels_train = Rehape(TrainLabels,NSeq * Ntrain,1);ntest = numel(testlabels);testlabels = repmat(testlabels',nseq,1);sequence_labels_test =重塑(testlabels,ntest * nseq,1);结束
helperMajorityVote
基于模式实现对分类的多数投票。如果没有唯一的模式,则投票NoUniqueMode
返回,以确保记录了分类错误。
功能[ClassVotes, ClassCounts] = helperMajorityVote (predLabels、origLabels类)%此功能支持小波工具箱示例。万博1manbetx它可能%更改或在将来的释放中删除。%如果标签不是分类的,则创建分类数组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)%此功能仅支持小波时间散射万博1manbetx心音图数据分类实例。它可能会改变%将在未来的版本中移除。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”,...{“异常”,“正常”});结束