Main Content

このページは前リリースの情報です。該当の英語のページはこのリリースで削除されています。

深層学習を使用した波形セグメンテーション

この例は、再帰型深層学習ネットワークと時間-周波数解析を使用して、人間の心電図 (ECG) 信号をセグメント化する方法を説明します。

はじめに

人間の心臓の電気的活動は、ベースライン信号から離れた一連の振幅として測定できます。単一の正常な心拍サイクルの場合、ECG 信号は次の心拍形態に分けられます [1]。

  • P 波 — QRS 群が心房性脱分極を表す前の小さいたわみ

  • QRS 群 — 心拍の最大振幅部分

  • T 波 — QRS 群が心室再分極を表した後の小さいたわみ

ECG 波形のこれらの領域のセグメンテーションでは、人間の心臓全体の健康状態と異常の存在を評価するのに役立つ測定値の基準を提供できます [2]。手作業で ECG 信号の各領域に注釈を付けると、手間と時間がかかるタスクになる可能性があります。信号処理と深層学習の手法を使用すると、関心領域の注釈を効率的かつ自動的に付けられるようになる可能性があります。

この例では、一般に公開されている QT データベースからの ECG 信号を使用します [3] [4]。データは、サンプル レート 250 Hz で合計 105 人の患者から測定された、約 15 分間の ECG の記録で構成されています。各記録を取得するために、検査員は患者の胸部の異なる場所に 2 つの電極を配置して、2 チャネル信号にします。データベースは、自動化されたエキスパート システムによって生成される信号領域ラベルを提供します [2]。この例は、深層学習の解を使用して、サンプルがある領域に応じて各 ECG 信号のサンプルにラベル付けすることを目的としています。信号に対して関心領域をラベル付けするこのプロセスは、多くの場合、波形セグメンテーションと呼ばれます。

信号領域を分類するよう深層ニューラル ネットワークに学習させるために、長短期記憶 (LSTM) ネットワークを使用できます。この例は、 LSTM セグメンテーション性能を向上させるために信号前処理技術と時間-周波数解析をどのように使用できるかを示します。特に、例は、フーリエ シンクロスクイーズド変換を使用して、ECG 信号の非定常動作を表しています。

データのダウンロードと準備

105 個の 2 チャネルの ECG 信号から成る各チャネルは、自動化されたエキスパート システムによってそれぞれラベル付けされ、210 個の MAT ファイルに領域ラベルと共に保存された合計 210 個の ECG 信号に対して個別に扱われます。このファイルは、https://www.mathworks.com/supportfiles/SPT/data/QTDatabaseECGData.zip で入手できます。

MATLAB® の tempdir コマンドによって指定された場所にある一時ディレクトリに、データ ファイルをダウンロードします。tempdir と異なるフォルダーにデータ ファイルを置きたい場合、これ以降の手順でディレクトリ名を変更してください。

% Download the data
dataURL = 'https://www.mathworks.com/supportfiles/SPT/data/QTDatabaseECGData.zip';
datasetFolder = fullfile(tempdir,'QTDataset');
zipFile = fullfile(tempdir,'QTDatabaseECGData.zip');
if ~exist(datasetFolder,'dir')
     websave(zipFile,dataURL);
end

% Unzip the data
unzip(zipFile,tempdir)

unzip 操作を実行すると、一時ディレクトリにフォルダー QTDatabaseECGData が作成され、このフォルダーに 210 個の MAT ファイルが格納されます。各ファイルでは、変数 ecgSignal に ECG 信号が格納され、変数 signalRegionLabels に領域ラベルのテーブルが格納されます。また、各ファイルでは、変数 Fs に信号のサンプル レートが格納されます。この例では、信号のサンプル レートはすべて 250 Hz です。

ファイル内のデータにアクセスするための信号データストアを作成します。この例では、一時ディレクトリの QTDatabaseECGData フォルダーにデータセットが保存されていると仮定します。そうでない場合は、以下のコードでデータのパスを変更してください。SignalVariableNames パラメーターを使用して、各ファイルから読み取りたい信号変数名を指定します。

sds = signalDatastore(datasetFolder,'SignalVariableNames',["ecgSignal","signalRegionLabels"])
sds = 
  signalDatastore with properties:

                       Files:{
                             '/tmp/QTDataset/ecg1.mat';
                             '/tmp/QTDataset/ecg10.mat';
                             '/tmp/QTDataset/ecg100.mat'
                              ... and 207 more
                             }
    AlternateFileSystemRoots: [0×0 string]
                    ReadSize: 1
         SignalVariableNames: ["ecgSignal"    "signalRegionLabels"]

関数 read を呼び出すたびに、ECG 信号および領域ラベルのテーブルが格納された 2 要素 cell 配列がデータストアによって返されます。データストアの関数 preview を使用して、最初のファイルの内容が、長さ 225,000 サンプルの ECG 信号、および 3385 個の領域ラベルを含むテーブルであることを確認します。

data = preview(sds)
data=2×1 cell array
    {225000×1 double}
    {  3385×2 table }

領域ラベル テーブルの最初の数行を確認し、領域範囲のインデックスと領域クラスの値 (P、T、または QRS) が各行に含まれていることを確認します。

head(data{2})
ans=8×2 table
    ROILimits     Value
    __________    _____

     83    117     P   
    130    153     QRS 
    201    246     T   
    285    319     P   
    332    357     QRS 
    412    457     T   
    477    507     P   
    524    547     QRS 

displayWaveformLabels 補助関数を使用して、最初の 1000 個のサンプルのラベルを可視化します。

displayWaveformLabels(data,true,1000)

通常の機械学習の分類手順は次のとおりです。

  1. データベースを学習データセットとテスト データセットに分割します。

  2. 学習データセットを使用してネットワークに学習させます。

  3. 学習済みのネットワークを使用して、テスト データセットについて予測を行います。

ネットワークをデータの 70% で学習し、残りの 30% でテストします。

再現可能な結果が必要な場合は、乱数発生器をリセットします。関数 dividerand を使用して、ファイルをシャッフルするためのランダムなインデックスを取得し、signalDatastore の関数 subset を使用して、データを学習とテストのデータストアに分割します。

rng default
[trainIdx,~,testIdx] = dividerand(numel(sds.Files),0.7,0,0.3);
trainDs = subset(sds,trainIdx);
testDs = subset(sds,testIdx);

このセグメンテーション問題では、ECG 信号が LSTM ネットワークに入力され、入力信号と同じ長さをもつラベルのシーケンスまたはマスクが LSTM ネットワークから出力されます。このネットワークのタスクは、各信号サンプルに対し、その信号が属する領域の名前をラベルとして付けることです。そのため、データセット上の領域ラベルを、信号サンプルごとにラベルが 1 つずつ付けられたシーケンスに変換する必要があります。変換されたデータストアと補助関数 roi2mask を使用して、ラベルを変換します。変換されたデータストアをプレビューし、このデータストアによって同じ長さの信号ベクトルとラベル ベクトルが返されることを確認します。カテゴリカル マスク ベクトルの先頭から 1000 個の要素をプロットします。関数 roi2mask は、いずれの関心領域にも属さないラベル サンプルにラベル カテゴリ n/a を追加します。

trainDs = transform(trainDs, @roi2mask);
testDs = transform(testDs, @roi2mask);

transformedData = preview(trainDs)
transformedData=1×2 cell array
    {224993×1 double}    {224993×1 categorical}

plot(transformedData{2}(1:1000))

非常に長い入力信号を LSTM ネットワークに渡すと、推定精度の低下と過剰なメモリ消費を招く可能性があります。これらの影響を回避するには、変換されたデータストアおよび補助関数 resizeData を使用して、ECG 信号およびそれらに対応するラベル マスクを分割します。この補助関数は、5000 サンプルのセグメントをできるだけ多く作成し、残りのサンプルを破棄します。変換されたデータストアの出力のプレビューを見ると、最初の ECG 信号とそのラベル マスクが 5000 サンプルのセグメントに分割されているのが分かります。変換されたデータストアのプレビューには、この cell 配列の最初の 8 要素のみが示されていることに注意してください。この cell 配列の要素数は、本来、floor(224993/5000) = 44 であり、データストアの関数 read を呼び出すと、この数の要素をもつ cell 配列が返されます。

trainDs = transform(trainDs,@resizeData);
testDs = transform(testDs,@resizeData);
preview(trainDs)
ans=8×2 cell array
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}
    {1×5000 double}    {1×5000 categorical}

生の ECG 信号の LSTM ネットワークへの直接入力

最初に、学習データセットから取得した生の ECG 信号を使用して、LSTM ネットワークに学習させます。

学習前のネットワーク アーキテクチャを定義します。サイズ 1 の sequenceInputLayer を指定して、1 次元の時系列を受け入れます。'sequence' 出力モードで LSTM 層を指定して、信号の各サンプルに対して分類できるようにします。最適なパフォーマンスを得るために 200 個の非表示ノードを使用します。出力サイズが 4 の fullyConnectedLayer を波形クラスごとに 1 つ指定します。softmaxLayerclassificationLayer を追加して、推定ラベルを出力します。

layers = [ ...
    sequenceInputLayer(1)
    lstmLayer(200,'OutputMode','sequence')
    fullyConnectedLayer(4)
    softmaxLayer
    classificationLayer];

適切なネットワーク性能を確認する学習プロセスに対するオプションを選択します。各パラメーターの説明については、trainingOptions (Deep Learning Toolbox)のドキュメンテーションを参照してください。

options = trainingOptions('adam', ...
    'MaxEpochs',10, ...
    'MiniBatchSize',50, ...
    'InitialLearnRate',0.01, ...
    'LearnRateDropPeriod',3, ...
    'LearnRateSchedule','piecewise', ...
    'GradientThreshold',1, ...
    'Plots','training-progress',...
    'shuffle','every-epoch',...
    'Verbose',0,...
    'DispatchInBackground',true);

この学習データセットはメモリにすべて収まるため、Parallel Computing Toolbox™ を使用できる場合、データストアの関数 tall を使用してデータを並列で変換してから、変換したデータをワークスペースに集めることができます。ニューラル ネットワークの学習は反復的です。データストアでは、反復のたびにファイルからデータを読み取り、そのデータを変換してからネットワークの係数を更新します。データがコンピューターのメモリに収まる場合、データをワークスペースにインポートしておけば、データの読み取りと変換を 1 回のみ実行すれば済むため、学習が高速になります。データがメモリに収まらない場合、学習エポックごとに学習関数にデータストアを渡して変換を実行しなければならないことに注意してください。

学習セットとテスト セットの両方に tall 配列を作成します。システムに応じて、MATLAB が作成する並列プール内のワーカー数は異なる場合があります。

tallTrainSet = tall(trainDs);
Starting parallel pool (parpool) using the 'local' profile ...
Connected to the parallel pool (number of workers: 8).
tallTestSet = tall(testDs);

次に、tall 配列の関数 gather を呼び出してデータセット全体の変換を計算し、学習信号、テスト信号、およびラベルが格納された cell 配列を取得します。

 trainData = gather(tallTrainSet);
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: 0% complete
Evaluation 0% complete
- Pass 1 of 1: Completed in 28 sec
Evaluation completed in 28 sec
 trainData(1,:)
ans=1×2 cell array
    {1×5000 double}    {1×5000 categorical}

 testData = gather(tallTestSet);
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 9.8 sec
Evaluation completed in 10 sec

ネットワークの学習

trainNetwork コマンドを使用して LSTM ネットワークを学習させます。データセットのサイズが大きいため、このプロセスには数分かかる場合があります。マシンに GPU と Parallel Computing Toolbox™ がある場合は、MATLAB は学習に GPU を自動で使用します。そうでない場合は CPU が使用されます。

図の学習の精度と損失のサブプロットは、すべての反復回数間の学習の進行状況を追跡します。生の信号データを使用して、ネットワークは、P 波、QRS 群、T 波、または N/A に属するものとして、サンプルの約 77% を正しく分類します。

net = trainNetwork(trainData(:,1),trainData(:,2),layers,options);

テスト データの分類

学習済み LSTM ネットワークと classify コマンドを使用してテスト データを分類します。学習オプションを一致させるためにミニバッチのサイズ 50 を指定します。

predTest = classify(net,testData(:,1),'MiniBatchSize',50);

混同行列では、分類性能を可視化するための直感的で有益な方法を提供しています。confusionchart コマンドを使用して、テスト データ予測に対する全体の分類精度を計算します。各入力に対して、categorical ラベルの cell 配列を行ベクトルに変換します。列を正視化した表示を指定して、各クラスのサンプルの割合として結果を表示します。

confusionchart([predTest{:}],[testData{:,2}],'Normalization','column-normalized');

ネットワークへの入力として生の ECG 信号を使用すると、T 波のサンプルの約 60%、P 波のサンプルの 40%、QRS 群のサンプルの 60% のみが正しくなりました。パフォーマンスを向上させるには、深層学習のネットワークに入力する前に ECG 信号の特性の情報を適用します。たとえば、患者の呼吸の動きによって引き起こされる基線変動です。

フィルター処理方法の適用による基線変動と高周波ノイズの除去

3 つの心拍形態は、異なる周波数バンドを占めます。QRS 群のスペクトルには一般的に、10 ~ 25 Hz あたりの中心周波数があり、そのコンポーネントは 40 Hz 未満にあります。P 波と T 波は低周波でも起こります。P 波のコンポーネントは 20 Hz 未満で、T 波のコンポーネントは 10 Hz 未満です [5]。

基線変動は、患者の呼吸の動きによって引き起こされる低周波数 (< 0.5 Hz) の振動です。この発振は、心拍形態から独立していて、有効な情報を提供しません [6]。

[0.5, 40] Hz の通過帯域周波数幅を持つバンドパス フィルターを設計して、変動と高周波ノイズを削除します。これらのコンポーネントを削除すると、ネットワークが関係のない特徴を学習しないため LSTM 学習が向上します。tall データの cell 配列に対して cellfun を使用し、データセットを並列でフィルター処理します。

% Bandpass filter design
hFilt = designfilt('bandpassiir', 'StopbandFrequency1',0.4215,'PassbandFrequency1', 0.5, ...
    'PassbandFrequency2',40,'StopbandFrequency2',53.345,...
    'StopbandAttenuation1',60,'PassbandRipple',0.1,'StopbandAttenuation2',60,...
    'SampleRate',250,'DesignMethod','ellip');

% Create tall arrays from the transformed datastores and filter the signals
tallTrainSet = tall(trainDs);
tallTestSet = tall(testDs);

filteredTrainSignals = gather(cellfun(@(x)filter(hFilt,x),tallTrainSet(:,1),'UniformOutput',false));
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: 0% complete
Evaluation 0% complete
- Pass 1 of 1: Completed in 19 sec
Evaluation completed in 19 sec
trainLabels = gather(tallTrainSet(:,2));
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 19 sec
Evaluation completed in 19 sec
filteredTestSignals = gather(cellfun(@(x)filter(hFilt,x),tallTestSet(:,1),'UniformOutput',false));
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 8.5 sec
Evaluation completed in 8.6 sec
testLabels = gather(tallTestSet(:,2));
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 8.3 sec
Evaluation completed in 8.5 sec

一般的なケースの生の信号とフィルター処理した信号をプロットします。

trainData = gather(tallTrainSet);
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 19 sec
Evaluation completed in 19 sec
figure
subplot(2,1,1)
plot(trainData{95,1}(2001:3000))
title('Raw')
grid
subplot(2,1,2)
plot(filteredTrainSignals{95}(2001:3000))
title('Filtered')
grid

フィルター済み ECG 信号を使用したネットワークの学習

前と同様に、同じネットワークのアーキテクチャを使用してフィルター済み ECG 信号で LSTM ネットワークに学習させます。

filteredNet = trainNetwork(filteredTrainSignals,trainLabels,layers,options);

信号の前処理は 80% 以上まで学習精度を向上させます。

フィルター処理された ECG 信号の分類

前処理済みテスト データを更新した LSTM ネットワークで分類します。

predFilteredTest = classify(filteredNet,filteredTestSignals,'MiniBatchSize',50);

混同行列として分類性能を可視化します。

figure
confusionchart([predFilteredTest{:}],[testLabels{:}],'Normalization','column-normalized');

シンプルな前処理で、T 波の分類が約 15%、QRS 群の分類が 10% 向上します。P 波の分類精度は変わりませんでした。

ECG 信号の時間-周波数の表現

時系列データの正常な分類に対する共通アプローチは、時間-周波数の特徴を抽出し、それらを元のデータの代わりにネットワークに供給することです。その後、ネットワークは、時間と周波数間で同時にパターンを学習します [7]。

フーリエ シンクロスクイーズド変換 (FSST) は、信号サンプルごとに周波数スペクトルを計算します。そのため、これは元の信号と同じ時間分解能を維持する必要があるこのセグメンテーション問題に最適です。fsst 関数を使用して、学習信号の 1 つの変換を検査します。適切な周波数分解能を与えるために、長さが 128 のカイザー ウィンドウを指定します。

data =  preview(trainDs);
figure
fsst(data{1,1},250,kaiser(128),'yaxis')

対象となる周波数範囲 [0.5, 40] Hz に対し、学習データセットの各信号の FSST を計算します。個別の機能として FSST の実数部と虚数部を取り扱い、ネットワークに両方のコンポーネントを供給します。さらに、平均を減算し、標準偏差で除算することで、学習用の特徴を標準化します。変換されたデータストア、補助関数 extractFSSTFeatures、および関数 tall を使用して、データを並列で処理します。

fsstTrainDs = transform(trainDs,@(x)extractFSSTFeatures(x,250));
fsstTallTrainSet = tall(fsstTrainDs);
fsstTrainData = gather(fsstTallTrainSet);
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: 0% complete
Evaluation 0% complete
- Pass 1 of 1: Completed in 2 min 35 sec
Evaluation completed in 2 min 35 sec

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

fsstTTestDs = transform(testDs,@(x)extractFSSTFeatures(x,250));
fsstTallTestSet = tall(fsstTTestDs);
fsstTestData = gather(fsstTallTestSet);
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 1 min 23 sec
Evaluation completed in 1 min 23 sec

ネットワーク アーキテクチャの調整

ネットワークが単一の値の代わりに各サンプルに周波数スペクトルを受け入れるように LSTM アーキテクチャを変更します。FSST のサイズを検査し、周波数の数値を確認します。

size(fsstTrainData{1,1})
ans = 1×2

          40        5000

入力の特徴を 40 として sequenceInputLayer に指定します。残りのネットワーク パラメーターは変更しないで保持します。

layers = [ ...
    sequenceInputLayer(40)
    lstmLayer(200,'OutputMode','sequence')
    fullyConnectedLayer(4)
    softmaxLayer
    classificationLayer];

ECG 信号の FSST を使用したネットワークの学習

変換されたデータセットを使用して更新された LSTM ネットワークに学習をさせます。

fsstNet = trainNetwork(fsstTrainData(:,1),fsstTrainData(:,2),layers,options);

時間-周波数の特徴を使用すると学習精度が改善します。現在、90% を超えています。

FSST を使用したテスト データの分類

更新された LSTM ネットワークと抽出された FSST の特徴を使用して、テスト データを分類します。

predFsstTest = classify(fsstNet,fsstTestData(:,1),'MiniBatchSize',50);

混同行列として分類性能を可視化します。

confusionchart([predFsstTest{:}],[fsstTestData{:,2}],'Normalization','column-normalized');

時間-周波数表現を使用すると、生データの結果と比べて、T 波の分類が約 20%、P 波の分類が約 40%、および QRS 群の分類が 30% 向上します。

補助関数 displayWaveformLabels を使用して、ネットワーク予測を単一の ECG 信号のグラウンド トゥルース ラベルと比較します。

testData = gather(tall(testDs));
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 8.3 sec
Evaluation completed in 8.5 sec
figure
subplot(2,1,1)
displayWaveformLabels({testData{1,1}(3000:4000),testData{1,2}(3000:4000)},false)
title('Ground Truth')

subplot(2,1,2)
displayWaveformLabels({testData{1,1}(3000:4000),predFsstTest{1}(3000:4000)},false)
title('Predicted')

まとめ

この例は、信号の前処理と時間-周波数解析が LSTM 波形セグメンテーション性能を向上させることができる方法を説明しています。バンドパス フィルターとフーリエ ベースのシンクロスクイージングによって、すべての出力クラス間の平均が 55% から約 85% まで改善されます。

参考文献

[1] McSharry, Patrick E., et al. "A dynamical model for generating synthetic electrocardiogram signals." IEEE® Transactions on Biomedical Engineering.Vol. 50, No. 3, 2003, pp. 289–294.

[2] Laguna, Pablo, Raimon Jané, and Pere Caminal."Automatic detection of wave boundaries in multilead ECG signals: Validation with the CSE database." Computers and Biomedical Research.Vol. 27, No. 1, 1994, pp. 45–60.

[3] Goldberger, Ary L., Luis A. N. Amaral, Leon Glass, Jeffery M. Hausdorff, Plamen Ch. Ivanov, Roger G. Mark, Joseph E. Mietus, George B. Moody, Chung-Kang Peng, and H. Eugene Stanley. "PhysioBank, PhysioToolkit, and PhysioNet: Components of a New Research Resource for Complex Physiologic Signals." Circulation. Vol. 101, No. 23, 2000, pp. e215–e220. [Circulation Electronic Pages; http://circ.ahajournals.org/content/101/23/e215.full].

[4] Laguna, Pablo, Roger G. Mark, Ary L. Goldberger, and George B. Moody."A Database for Evaluation of Algorithms for Measurement of QT and Other Waveform Intervals in the ECG."Computers in Cardiology.Vol.24, 1997, pp. 673–676.

[5] Sörnmo, Leif, and Pablo Laguna."Electrocardiogram (ECG) signal processing."Wiley Encyclopedia of Biomedical Engineering, 2006.

[6] Kohler, B-U., Carsten Hennig, and Reinhold Orglmeister."The principles of software QRS detection."IEEE Engineering in Medicine and Biology Magazine.Vol. 21, No. 1, 2002, pp. 42–57.

[7] Salamon, Justin, and Juan Pablo Bello."Deep convolutional neural networks and data augmentation for environmental sound classification."IEEE Signal Processing Letters.Vol. 24, No. 3, 2017, pp. 279–283.

参考

関数

関連するトピック