ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

心音図データのウェーブレット時間散乱の分類

この例では、ウェーブレット時間散乱とサポート ベクター マシン (SVM) 分類器を使用して人間の心音図 (PCG) の記録を分類する方法を示します。心音図は、心臓の伸縮と拡張の位相によって作成された音の音響記録です。心臓の聴診は、心臓病患者の健康を評価するうえで重要な診断の役割を継続して担います。ただし、世界の多くの地域では、心臓の聴診について学習した医療関係者の数が足りません。したがって、心音図データを解釈する信頼性の高い自動化された方法を開発する必要があります。

この例では、PCG 分類の特徴抽出器としてウェーブレット散乱を使用します。ウェーブレット散乱では、データは、データの低分散の表現を作成するために一連のウェーブレット変換、非線形性、および平均化全体に伝播されます。これらの低分散の表現は、その後分類器で入力として使用されます。この例は、各 PCG 記録が "正常" または "異常" のいずれかであるバイナリ分類問題です。

データの説明

この例では、心臓機能が正常と異常な人から取得した心音 (PCG) データを使用します。データセットは 3829 個の記録で構成されます。心臓機能が正常な人が 2575 個、心臓機能が異常な人が 1254 個です。各記録は、10,000 サンプルの長さで、2 kHz でサンプリングされます。これは、5 秒の心音データを表します。データセットは、PhysioNet Computing in Cardiology Challenge 2016 で使用された学習データと検証データから構築されます [1][2]。

データのダウンロード

1 番目のステップは、GitHub リポジトリからデータをダウンロードすることです。データをダウンロードするには、[Clone or download] をクリックして [Download ZIP] を選択します。書き込み権限のあるフォルダーに、ファイル physionet_phonocardiogram-master を保存します。この例の手順では、ファイルを一時ディレクトリ (MATLAB™ の tempdir) にダウンロードしているものと仮定します。tempdir とは異なるフォルダーにデータをダウンロードすることを選択した場合は、データの解凍および読み込みに関する後続の手順を変更してください。Git に精通している場合は、最新バージョンのツール (git) をダウンロードし、git clone https://github.com/mathworks/physionet_phonocardiogram を使用してシステム コマンド プロンプトからデータを取得できます。

ファイル physionet_phonocardiogram-master.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 のみです。

データの読み込み

前のセクションのダウンロード手順に従った場合、次のコマンドを入力して 2 つのアーカイブ ファイルを解凍します。

unzip(fullfile(tempdir,'physionet_phonocardiogram-master.zip'),tempdir)
unzip(fullfile(tempdir,'physionet_phonocardiogram-master','PCG_Data.zip'),...
    fullfile(tempdir,'PCG_Data'))

PCG_Data.zip ファイルを解凍したら、データを MATLAB に読み込みます。

load(fullfile(tempdir,'PCG_Data','heartSoundData.mat'))

heartSoundData は 2 つのフィールド (DataClasses) をもつ構造体配列です。Data は 10000 行 3829 列の行列で、各列は PCG 記録です。Classes は 3829 行 1 列の診断ラベルの categorical 配列で、それぞれの列が Data の各列に対応します。これはバイナリ分類問題のため、クラスは "正常" および "異常" です。前述のとおり、2575 個の正常の記録と 1254 個の異常の記録があります。つまり、データ内の例の 67.25% は心臓機能が正常な人のもので、32.75% は心臓機能が異常な人のものです。これは、次を入力して検証できます。

summary(heartSoundData.Classes)
     normal        2575 
     abnormal      1254 
countcats(heartSoundData.Classes)./sum(countcats(heartSoundData.Classes))
ans = 2×1

    0.6725
    0.3275

ウェーブレット散乱フレームワーク

waveletScatteringを使用して、ウェーブレット時間散乱フレームワークを構築します。信号長を一致させるために不変スケールを設定します。既定の散乱フレームワークには、2 つのウェーブレット変換 (フィルター バンク) があります。最初のウェーブレット フィルター バンクにはオクターブあたりに 8 個のウェーブレットがあります。2 つ目のフィルター バンクにはオクターブあたりに 1 個のウェーブレットがあります。

N = 1e4;
sn = waveletScattering('SignalLength',N,'InvarianceScale',N);

学習セットとテスト セットの作成

70% (2680) が学習セット内にあり、1802 個の正常、878 個の異常をもつように、補助関数 partition_heartsounds は 3829 個の観測値を分割します。残りの 1149 個の記録 (773 個の正常と 376 個の異常) は予測のテスト セットに割り当てられます。結果が繰り返し再現されるように、乱数発生器は補助関数内にシードされます。partition_heartsounds とこの例で使用されるその他すべての補助関数のコードは、例の最後にある「サポート関数」の節に示しています。

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

学習セットとテスト セットに含まれる各クラスの数を調べることができます。

summary(trainLabels)
     normal        1802 
     abnormal       878 
summary(testLabels)
     normal        773 
     abnormal      376 

学習セットとテスト セットの "正常" と "異常" の記録の比率が全体のデータの比率と同じになるように、学習セットとテスト セットは分割されています。以下を使用してこれを確認することができます。

countcats(trainLabels)./sum(countcats(trainLabels))
ans = 2×1

    0.6724
    0.3276

countcats(testLabels)./sum(countcats(testLabels))
ans = 2×1

    0.6728
    0.3272

散乱特徴量

学習セットの 2680 個のすべての記録の散乱変換を取得します。多変数の時系列の場合、散乱変換は各列が個別の信号であると仮定しています。'log' オプションを使用して、散乱係数の自然対数を取得します。

scat_features_train = featureMatrix(sn,trainData,'Transform','log');

指定された散乱パラメーターの場合、scat_features_train は 340 行 5 列の 2680 行列です。2680 個の信号のそれぞれに 340 個の散乱パスと 5 つの散乱時間枠があります。SVM 分類器にこれを渡すには、各行が 340 個の散乱パス間で単一の散乱時間枠を表す 13400 行 340 列の行列にテンソルの形状を変更します。合計の行数は、5 と 2680 の積 (学習データの記録の数) に等しくなります。

Nseq = size(scat_features_train,2);
scat_features_train = permute(scat_features_train,[2 3 1]);
scat_features_train = reshape(scat_features_train,...
    size(scat_features_train,1)*size(scat_features_train,2),[]);

テスト データに対して手順を繰り返します。

scat_features_test = featureMatrix(sn,testData,'Transform','log');
scat_features_test = permute(scat_features_test,[2 3 1]);
scat_features_test = reshape(scat_features_test,...
    size(scat_features_test,1)*size(scat_features_test,2),[]);

ここで、各散乱時間枠にラベルをもつようにラベルを複製します。

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

SVM を学習データに適合させます。この例では、3 次多項式カーネルを使用します。学習データに SVM を適合させた後、5 分割の交差検証を実行して、学習データで汎化誤差を推定します。

rng default;
classificationSVM = fitcsvm(...
    scat_features_train, ...
    sequence_labels_train , ...
    'KernelFunction', 'polynomial', ...
    'PolynomialOrder', 3, ...
    'KernelScale', 'auto', ...
    'BoxConstraint', 1, ...
    'Standardize', true, ...
    'ClassNames', categorical({'normal','abnormal'}));
kfoldmodel = crossval(classificationSVM, 'KFold', 5);

割合として損失を計算して、混同行列を表示します。

predLabels = kfoldPredict(kfoldmodel);
loss = kfoldLoss(kfoldmodel)*100;
fprintf('Loss is %2.2f percent\n',loss);
Loss is 0.79 percent
accuracy = 100-loss;
fprintf('Accuracy is %2.2f percent\n',accuracy);
Accuracy is 99.21 percent
confmatCV = confusionchart(sequence_labels_train,predLabels);

各時間枠が個別に分類されるときに散乱フレームワークの精度は 99.21% になります。ただし、記録ごとに 5 つの散乱時間枠があり、99.21% の精度はすべての枠を個別に分類することに基づいているため、性能はこの値より実際には良くなります。この場合、多数決を使用して、記録ごとに単一のクラス割り当てを取得します。クラスの票は、5 つの枠の投票のモードに対応します。固有のモードが見つからない場合、helperMajorityVote は、分類誤差を示すために散乱時間枠のセットを 'NoUniqueMode' として分類します。これは、混同行列の余分の列になります。

classes = categorical({'abnormal','normal'});
ClassVotes = helperMajorityVote(predLabels,trainLabels,classes);
CVaccuracy = sum(eq(ClassVotes,trainLabels))./numel(trainLabels)*100;
fprintf('The true cross-validation accuracy is %2.2f percent.\n',CVaccuracy);
The true cross-validation accuracy is 99.93 percent.

多数決の分類に対する混同行列を表示します。

cmCV = confusionchart(ClassVotes,trainLabels);

学習データの交差検証の精度は実際には 99.93% です。2 つの正常の記録がありますが、異常として誤分類されます。

学習データに適合させる SVM モデルを使用して、ホールドアウトされたテスト データでクラス予測を行います。

predTestLabels = predict(classificationSVM,scat_features_test);

多数決を使用して、テスト セットで予測の精度を判断します。

ClassVotes = helperMajorityVote(predTestLabels,testLabels,classes);
testaccuracy = sum(eq(ClassVotes,testLabels))./numel(testLabels)*100;
fprintf('The test accuracy is %2.2f percent.\n',testaccuracy);
The test accuracy is 92.25 percent.
cmTest = confusionchart(ClassVotes,testLabels);

1149 個のテスト記録の 92.25% は "正常" または "異常" として正確に分類されます。テスト セットの 773 個の正常の PCG 記録のうち、728 個は正確に分類されます。テスト セットの 376 個の異常の記録のうち、332 個は正確に分類されます。

適合率、再現率および F1 スコア

分類タスクにおいて、クラスの適合率とは正しい正の結果の数を正の結果の数で除算したものです。言い換えれば、分類器が特定のラベルを割り当てるすべてのレコードのうち、実際にクラスに属している割合のことです。再現率は、正しいラベルの数を特定のクラスのラベル数で除算したものと定義されます。具体的には、クラスに属しているすべてのレコードのうち、分類器によってそのクラスのラベルが付けられた割合です。分類器の精度を判断する場合に、理想的には適合率と再現率の両方で良い結果を得る必要があります。たとえば、1 つ 1 つの PCG の記録に異常のラベルを付けた分類器があると仮定します。異常クラスの再現率は 100% になります。異常クラスに属しているすべての記録が異常とラベル付けされます。しかし、適合率は低くなります。分類器によってすべての記録が異常とラベル付けされたため、この場合は 2575 個の偽陽性があり、適合率は 1254/3829 (つまり 32.75%) になります。F1 スコアは適合率および再現率の調和平均であり、再現率と適合率の両方の観点から分類器の性能を要約する 1 つのメトリクスを提供します。補助関数 helperF1heartSounds は、テスト セットに分類結果に関する適合率、再現率、F1 スコアを計算し、テーブルにそれらの結果を返します。

PRTable = helperF1heartSounds(cmTest.NormalizedValues);
disp(PRTable)
                Precision    Recall    F1_Score
                _________    ______    ________

    Abnormal     88.298      88.064     88.181 
    Normal       94.179      94.301     94.239 

この場合、異常グループと正常グループに関する F1 スコアは、モデルに良好な適合率と再現率が両方あることを確認します。バイナリ分類では、適合率と再現率の両方を混同行列から直接判断することは簡単です。これを確認するには、便宜上、混同行列を再度プロットします。

cmTest = confusionchart(ClassVotes,testLabels);

異常クラスの再現率は異常として識別された異常の記録の数です。これは、混同行列の最初の行と最初の列のエントリを最初の行のエントリの合計で除算したものです。異常クラスの適合率は、分類器によって異常として識別された合計数での、真の異常の記録の比率です。これは、混同行列最初の行と最初の列のエントリを最初の列のエントリの合計で除算したものに対応します。F1 スコアは、2 つの調和平均です。

RecallAbnormal = cmTest.NormalizedValues(1,1)/sum(cmTest.NormalizedValues(1,:));
PrecisionAbnormal = cmTest.NormalizedValues(1,1)/sum(cmTest.NormalizedValues(:,1));
F1Abnormal = harmmean([RecallAbnormal PrecisionAbnormal]);
fprintf('RecallAbnormal = %2.3f\nPrecisionAbnormal = %2.3f\nF1Abnormal = %2.3f\n',...
    100*RecallAbnormal,100*PrecisionAbnormal,100*F1Abnormal);
RecallAbnormal = 88.064
PrecisionAbnormal = 88.298
F1Abnormal = 88.181

正常クラスについて上記の内容を繰り返します。

RecallNormal = cmTest.NormalizedValues(2,2)/sum(cmTest.NormalizedValues(2,:));
PrecisionNormal = cmTest.NormalizedValues(2,2)/sum(cmTest.NormalizedValues(:,2));
F1Normal = harmmean([RecallNormal PrecisionNormal]);
fprintf('RecallNormal = %2.3f\nPrecisionNormal = %2.3f\nF1Normal = %2.3f\n',...
    100*RecallNormal,100*PrecisionNormal,100*F1Normal);
RecallNormal = 94.301
PrecisionNormal = 94.179
F1Normal = 94.239

まとめ

この例では、バイナリ分類問題の正常または異常として人間の心音図の記録を確実に識別するためにウェーブレット時間散乱を使用しました。サポート ベクター マシン分類器が 2 つのグループ間での差異を正確にモデル化できるよう PCG データの低分散表現を作成するために、ウェーブレット散乱には単一のパラメーター、スケール不変の長さの仕様のみが必要でした。ウェーブレット散乱を使用するサポート ベクター マシン分類器は、学習セットとテスト セットの両方の正常と異常の PCG 記録の数がかなり不均衡であるにもかかわらず、両方のグループの適合率と再現率の双方において優れたパフォーマンスの実現が可能でした。

参考文献

  1. Goldberger, A. L., L. A. N. Amaral, L. Glass, J. M. Hausdorff, P. Ch. Ivanov, R. G. Mark, J. E. Mietus, G. B. Moody, C.-K.Peng, and H. E. Stanley."PhysioBank, PhysioToolkit, and PhysioNet: Components of a New Research Resource for Complex Physiologic Signals".Circulation.Vol. 101, No. 23, 13 June 2000, pp. e215-e220. http://circ.ahajournals.org/content/101/23/e215.full

  2. Liu et al. "An open access database for the evaluation of heart sound algorithms".Physiological Measurement.Vol. 37, No. 12, 21 November 2016, pp. 2181-2213. https://www.ncbi.nlm.nih.gov/pubmed/27869105

サポート関数

partition_heartsounds は、データの指定したプロパティを構成する学習セットとテスト セットを作成します。さらに、関数は、各セットの異常と正常の PCG 記録の比率を保持します。

function [trainData, testData, trainLabels, testLabels] = partition_heartsounds(percent_train_split,Data,Labels)
% This function is only in support of the Wavelet Time Scattering 
% Classification of Phonocardiogram Data example. It may change or be
% removed in a future release.

% Labels in heart sound data are not sequential.
percent_train_split = percent_train_split/100;
% Each column is an observation
NormalData = Data(:,Labels == 'normal');
AbnormalData = Data(:,Labels == 'abnormal');
LabelsNormal = Labels(Labels == 'normal');
LabelsAbnormal = Labels(Labels == 'abnormal');
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 default;
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];
   
end

createSequenceLabels_heartsounds は、ウェーブレット時間散乱シーケンスに対するクラス ラベルを作成します。

function [sequence_labels_train,sequence_labels_test] = createSequenceLabels_heartsounds(Nseq,trainLabels,testLabels)
% This function is only in support of the Wavelet Time Scattering 
% Classification of Phonocardiogram Data example. It may change or be
% removed in a future release.
Ntrain = numel(trainLabels);
trainLabels = repmat(trainLabels',Nseq,1);
sequence_labels_train = reshape(trainLabels,Nseq*Ntrain,1);
Ntest = numel(testLabels);
testLabels = repmat(testLabels',Nseq,1);
sequence_labels_test = reshape(testLabels,Ntest*Nseq,1);
end

helperMajorityVote は、モードを基にした分類に多数決を実装します。固有モードが存在しない場合、NoUniqueMode の投票は分類誤差が記録されていることを確認するために返されます。

function [ClassVotes,ClassCounts] = helperMajorityVote(predLabels,origLabels,classes)
% This function is in support of ECGWaveletTimeScatteringExample. It may
% change or be removed in a future release.

% Make categorical arrays if the labels are not already categorical
predLabels = categorical(predLabels);
origLabels = categorical(origLabels);
% Expects both predLabels and origLabels to be categorical vectors
Npred = numel(predLabels);
Norig = numel(origLabels);
Nwin = Npred/Norig;
predLabels = reshape(predLabels,Nwin,Norig);
ClassCounts = countcats(predLabels);
[mxcount,idx] = max(ClassCounts);
ClassVotes = classes(idx);
% Check for any ties in the maximum values and ensure they are marked as
% error if the mode occurs more than once
modecnt = modecount(ClassCounts,mxcount);
ClassVotes(modecnt>1) = categorical({'NoUniqueMode'});
ClassVotes = ClassVotes(:);

%-------------------------------------------------------------------------
function modecnt = modecount(ClassCounts,mxcount)
modecnt = Inf(size(ClassCounts,2),1);
for nc = 1:size(ClassCounts,2)
    modecnt(nc) = histc(ClassCounts(:,nc),mxcount(nc));
end

end

% EOF
end

helperF1heartSounds は、分類器の結果に対する適合率、再現率、および F1 スコアを計算します。

function PRTable = helperF1heartSounds(confmat)
% This function is only in support of the Wavelet Time Scattering 
% Classification of Phonocardiogram Data example. It may change or be
% removed in a future release.

precisionAB = confmat(1,1)/sum(confmat(:,1))*100;
precisionNR = confmat(2,2)/sum(confmat(:,2))*100 ;
recallAB = confmat(1,1)/sum(confmat(1,:))*100;
recallNR = confmat(2,2)/sum(confmat(2,:))*100;
F1AB = 2*(precisionAB*recallAB)/(precisionAB+recallAB);
F1NR = 2*(precisionNR*recallNR)/(precisionNR+recallNR);
% Construct a MATLAB Table to display the results.
PRTable = array2table([precisionAB recallAB F1AB;...
    precisionNR recallNR F1NR],...
    'VariableNames',{'Precision','Recall','F1_Score'},'RowNames',...
    {'Abnormal','Normal'});

end