algorithm - （心電図）波のパターンを検出する方法は？

Question

心電図から画像を読み取り、その中の主要な波 (P 波、QRS 群、T 波) をそれぞれ検出しようとしています。画像を読み取ってベクトルを取得できます（のように(4.2; 4.4; 4.9; 4.7; ...)）。このベクトルをたどって、これらの各波の開始と終了を検出できるアルゴリズムが必要です。例：

サイズが常に同じだったり、心電図の波数が事前に分かっていれば簡単です。与えられた波:

ベクトルを抽出します。

[0; 0; 20; 20; 20; 19; 18; 17; 17; 17; 17; 17; 16; 16; 16; 16; 16; 16; 16; 17; 17; 18; 19; 20; 21; 22; 23; 23; 23; 25; 25; 23; 22; 20; 19; 17; 16; 16; 14; 13; 14; 13; 13; 12; 12; 12; 12; 12; 11; 11; 10; 12; 16; 22; 31; 38; 45; 51; 47; 41; 33; 26; 21; 17; 17; 16; 16; 15; 16; 17; 17; 18; 18; 17; 18; 18; 18; 18; 18; 18; 18; 17; 17; 18; 19; 18; 18; 19; 19; 19; 19; 20; 20; 19; 20; 22; 24; 24; 25; 26; 27; 28; 29; 30; 31; 31; 31; 32; 32; 32; 31; 29; 28; 26; 24; 22; 20; 20; 19; 18; 18; 17; 17; 16; 16; 15; 15; 16; 15; 15; 15; 15; 15; 15; 15; 15; 15; 14; 15; 16; 16; 16; 16; 16; 16; 16; 16; 16; 15; 16; 15; 15; 15; 16; 16; 16; 16; 16; 16; 16; 16; 15; 16; 16; 16; 16; 16; 15; 15; 15; 15; 15; 16; 16; 17; 18; 18; 19; 19; 19; 20; 21; 22; 22; 22; 22; 21; 20; 18; 17; 17; 15; 15; 14; 14; 13; 13; 14; 13; 13; 13; 12; 12; 12; 12; 13; 18; 23; 30; 38; 47; 51; 44; 39; 31; 24; 18; 16; 15; 15; 15; 15; 15; 15; 16; 16; 16; 17; 16; 16; 17; 17; 16; 17; 17; 17; 17; 18; 18; 18; 18; 19; 19; 20; 20; 20; 20; 21; 22; 22; 24; 25; 26; 27; 28; 29; 30; 31; 32; 33; 32; 33; 33; 33; 32; 30; 28; 26; 24; 23; 23; 22; 20; 19; 19; 18; 17; 17; 18; 17; 18; 18; 17; 18; 17; 18; 18; 17; 17; 17; 17; 16; 17; 17; 17; 18; 18; 17; 17; 18; 18; 18; 19; 18; 18; 17; 18; 18; 17; 17; 17; 17; 17; 18; 17; 17; 18; 17; 17; 17; 17; 17; 17; 17; 18; 17; 17; 18; 18; 18; 20; 20; 21; 21; 22; 23; 24; 23; 23; 21; 21; 20; 18; 18; 17; 16; 14; 13; 13; 13; 13; 13; 13; 13; 13; 13; 12; 12; 12; 16; 19; 28; 36; 47; 51; 46; 40; 32; 24; 20; 18; 16; 16; 16; 16; 15; 16; 16; 16; 17; 17; 17; 18; 17; 17; 18; 18; 18; 18; 19; 18; 18; 19; 20; 20; 20; 20; 20; 21; 21; 22; 22; 23; 25; 26; 27; 29; 29; 30; 31; 32; 33; 33; 33; 34; 35; 35; 35; 0; 0; 0; 0;]

たとえば、次のように検出したいと思います。

• P波で[19 - 37]。
• の QRS コンプレックス[51 - 64]。
• 等

Answer 1

私が最初にすることは、すでにそこにあるものを見ることです。実際、この特定の問題はすでに十分に研究されています。以下は、いくつかの非常に単純なメソッドの簡単な概要です: link。

私は別の答えにも答えなければなりません。信号処理と音楽情報検索の研究をしています。表面的には、この問題はオンセット検出に似ているように見えますが、問題のコンテキストは同じではありません。このタイプの生物学的信号処理、すなわち、P、QRS、および T フェーズの検出は、これらの波形のそれぞれの特定の時間領域特性に関する知識を活用できます。MIR での開始検出は、実際にはそうではありません。（少なくとも確実ではありません。）

QRS 検出に適した方法の 1 つ (ノート オンセット検出には必ずしも有効ではありません) は、ダイナミック タイム ワーピングです。時間領域の特性が変わらない場合、DTW は非常にうまく機能します。この問題に DTW を使用する IEEE の短い論文を次に示します: link。

これは、多くの方法を比較する素敵な IEEE 雑誌の記事です:リンク。多くの一般的な信号処理モデルが試されていることがわかります。紙をざっと読んで、基本的なレベルで理解しているものを試してください。

編集: これらの記事を閲覧した後、ウェーブレット ベースのアプローチが最も直感的に思えます。DTW もうまく機能し、DTW モジュールも存在しますが、私にはウェーブレット アプローチが最適と思われます。他の誰かが、信号の派生物を利用して答えました。私の最初のリンクは、それを行う 1990 年以前の方法を調べていますが、それらは最新の方法ほど堅牢ではないのではないかと思います。

編集:機会があれば簡単な解決策を提示しようと思いますが、ここでウェーブレットが適していると思う理由は、時間や振幅のスケーリングに関係なく、さまざまな形状をパラメーター化するのに役立つからです。言い換えれば、同じ時間的形状が繰り返されているが、時間スケールと振幅が異なる信号がある場合でも、ウェーブレット解析ではこれらの形状が類似していると認識できます (大まかに言えば)。また、フィルタ バンクをこのカテゴリにまとめていることに注意してください。似たようなもの。

Answer 2

このパズルのピースは「開始検出」であり、この問題を解決するために多数の複雑なアルゴリズムが作成されています。オンセットの詳細については、こちらをご覧ください。

次のピースはハミング距離です。このアルゴリズムを使用すると、ファジー比較を行うことができます。入力は 2 つの配列で、出力は整数の「距離」または 2 つのデータ セット間の差です。数字が小さいほど、2 つが似ています。これは必要なものに非常に近いですが、正確ではありません。私は先に進み、新しい距離を計算するためにハミング距離アルゴリズムにいくつかの変更を加えました。おそらく名前がありますが、それが何であるかはわかりません。基本的に、配列内の各要素間の絶対距離を合計し、合計を返します。Pythonでのコードは次のとおりです。

import math

def absolute_distance(a1, a2, length):
       total_distance=0
       for x in range(0,length):
               total_distance+=math.fabs(a1[x]-a2[x])
       return total_distance

print(absolute_distance([1,3,9,10],[1,3,8,11],4))

このスクリプトは、これら 2 つの配列間の距離である 2 を出力します。

さて、これらのピースを組み立てます。オンセット検出を使用して、データ セット内のすべての波の始まりを見つけることができます。次に、各波をサンプル P 波と比較して、これらの場所をループできます。QRS コンプレックスに当たると、距離が最大になります。別の P 波に当たった場合、数値はゼロにはなりませんが、はるかに小さくなります。P-Wave と T-Wave の間の距離はかなり小さくなりますが、次の仮定を行う場合、これは問題ではありません。

The distance between any p-wave and any other p-wave will be smaller than the distance between any p-wave and any t-wave.

シリーズは次のようになります: pQtpQtpQt... p 波と t 波は隣り合っていますが、このシーケンスは予測可能であるため、読みやすくなります。

そうでない面では、おそらくこの問題に対する微積分ベースの解決策があります。しかし、私の考えでは、カーブフィッティングと積分はこの問題をより混乱させます。私が書いた距離関数は、両方の曲線の積分を差し引いて非常によく似 た面積の差を見つけます。

一度に 1 ポイントずつ反復して O(n) 距離計算を実行することを優先して、開始計算を犠牲にすることが可能かもしれません。ここで、n はグラフ内のポイントの数です。これらの距離計算のすべてのリストがあり、50 の pQt シーケンスがどこにあるかを知っている場合、p 波のすべての場所で 重複しない50 の最短距離がわかります。ビンゴ！ 簡単にするためにどうですか？ただし、トレードオフは、距離計算の数が増えることによる効率の低下です。

Answer 3

相互相関を使用できます。各パターンのモデル サンプルを取得し、信号と関連付けます。相関が高いピークが得られます。qrs波とt波を抽出するこの手法で良い結果が得られると期待しています。その後、qrs の前にある相関信号のピークを探すことで、p 波を抽出できます。

相互相関は、非常に簡単に実装できるアルゴリズムです。基本的：

x is array with your signal of length Lx
y is an array containing a sample of the signal you want to recognize of length Ly
r is the resulting correlation

for (i=0; i<Lx - Ly; i++){
  r[i] = 0;
  for (j=0; j<Ly ; j++){
    r[i] += x[i+j]*y[j];
  }
}

そして、r のピークを探します (たとえば、しきい値を超える値)

Answer 4

最初に行うことは、データを単純化することです。

絶対データを分析する代わりに、あるデータ ポイントから次のデータ ポイントへの変化量を分析します。

;これは、分離されたデータを入力として受け取り、そのデータのデルタを出力する簡単なワンライナーです。

perl -0x3b -ple'( $last, $_ ) = ( $_, $_-$last )' < test.in > test.out

提供したデータで実行すると、次の出力が得られます。

0;0;20;0;0;-1;-1;-1;0;0;0;0;-1;0;0;0;0;0;0;1;0;1;1; 1;1;1;1;0;0;2;0;-2;-1;-2;-1;-2;-1;0;-2;-1;1;-1;0;- 1;0;0;0; 0;-1;0;-1;2;4;6;9;7;7;6;-4;-6;-8;-7;-5;-4;0;-1;0;- 1;1;1;0;1;0;-1;1;0;0;0;0;0;0;-1;0;1;1;-1;0;1;0;0;0 ;1;0;-1;1; 2;2;0;1;1;1;1;1;1;0;0;1;0;0;-1;-2;-1;-2;-2;-2;-2 ;0;-1;-1;0;-1;0;-1;0;-1;0;1;-1;0;0;0;0;0;0;0;0;-1; 1;1;0;0;0; 0;0;0;0;0;-1;1;-1;0;0;1;0;0;0;0;0;0;0;-1;1;0;0;0;0 ;-1;0;0;0;0;1;0;1;1;0;1;0;0;1;1;1;0;0;0;-1;-1;-2;- 1;0;-2;0; -1;0;-1;0;1;-1;0;0;-1;0;0;0;1;5;5;7;8;9;4;-7;-5;-8 ;-7;-6;-2;-1;0;0;0;0;0;1;0;0;1;-1;0;1;0;-1;1;0;0;0 ;1;0;0;0; 1;0;1;0;0;0;1;1;0;2;1;1;1;1;1;1;1;1;-1;1;0;0;-1; -2;-2;-2;-2;-1;0;-1;-2;-1;0;-1;-1;0;1;-1;1;0;-1;1; -1;1;0;-1; 0;0;0;-1;1;0;0;1;0;-1;0;1;0;0;1;-1;0;-1;1;0;-1;0;0 ;0;0;1;-1;0;1;-1;0;0;0;0;0;0;1;-1;0;1;0;0;2;0;1;0; 1;1;1;-1; 0;-2;0;-1;-2;0;-1;-1;-2;-1;0;0;0;0;0;0;0;0;-1;0;0; 4;3;9;8;11;4;-5;-6;-8; -8;-4;-2;-2;0;0;0;-1;1;0;0;1;0;0;1;-1; 0;1;0;0;0;1;-1;0;1;1;0;0;0;0;1;0;1;0;1;2;1;1;2;0;1 ;1;1;1;0;0;1;1;0;0;-35;0;0;0;

_{上記のテキストには、元々出力に存在しない改行が挿入されています。}

---

それが終わったら、qrs コンプレックスを見つけるのは簡単です。

perl -F';' -ane'@F = map { abs($_) > 2 and $_ } @F; print join ";", @F'< test.out

;;20;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;4;6;9;7;7;6;-4;-6;-8;-7;-5;-4;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;5;5;7;8;9;4;-7;-5;-8;-7;-6
;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;4;3;9;8;11;4;-5;-6;-8;-8;-4;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;-35 ;;;

20およびデータ ポイントは-35、 で始まり で終わる元のデータから生成され0ます。

他のデータ ポイントを見つけるには、パターン マッチングに頼る必要があります。

---

最初の p 波を見ると、はっきりとパターンが見えます。

0;0;0;0;0;0;1;0;1;1;1;1;1;1;0;0;2;0;-2;-1;-2;-1;-2;-1;0;-2;-1;1;-1;0;-1;0;0;0;0;
#           \________ up _______/   \________ down _________/

ただし、2番目のp波のパターンを見るのは簡単ではありません. これは、2番目のものがさらに広がっているためです

0;0;0;1;0;1;1;0;1;0;0;1;1;1;0;0;0;-1;-1;-2;-1;0;-2;0;-1;0;-1;0;1;-1;0;0;-1;0;0;0;
#     \________ up _______/       \________________ down ________________/

3 番目の p 波は、他の 2 つよりも少し不安定です。

0;0;0;0;0;1;-1;0;1;0;0;2;0;1;0;1;1;1;-1;0;-2;0;-1;-2;0;-1;-1;-2;-1;0;0;0;0;0;
#                \_______ up ______/  \__________ down __________/

t 波は、p 波と同様の方法で見つけることができます。主な違いは、それらがいつ発生するかです。

---

これは、開始するのに十分な情報です。

_{2 つのワンライナーは単なる例であり、日常的な使用はお勧めしません。}

Answer 5

他の 2 つの鋭い山と谷も qrs 複合体ですか?

頭のてっぺんから、あなたがする必要があるのは、このグラフの各ポイントでの傾きを計算することだと思います。次に、勾配がどのくらい速く変化するかも確認する必要があります (2 次導関数???)。急激な変化がある場合は、ある種の鋭いピークに達したことがわかります。もちろん、変化の検出を制限したいので、グラフの小さな隆起を拾わないように、「時間間隔 T にわたって勾配が X だけ変化した場合」のようなことをしたいかもしれません。

数学をやったのは久しぶりです...これは数学の質問のようです;)ああ、信号分析も行っていません:)。

もう一点追加するだけです。信号平均化を試すこともできると思います。たとえば、最後の 3 つまたは 4 つのデータ ポイントを平均化します。急激な変化もそのように検出できると思います。

Answer 6

私はこの特定の問題の専門家ではありませんが、より一般的な知識から頭のてっぺんから外れているだけです: QRS コンプレックス (または他の機能の 1 つですが、この例では QRS コンプレックスを使用します) を知っているとしましょう。は、おおよそ長さ L の固定された期間に発生します。これを分類問題として次のように扱うことができるでしょうか。

1. 信号を長さ L の重なり合うウィンドウに分割します。各ウィンドウには、完全な QRS 群が含まれているか含まれていないかのいずれかです。
2. 各ウィンドウをフーリエ変換します。あなたの特徴は、各周波数での信号強度です。
3. 手作業で注釈を付けたデータで決定木、サポート ベクター マシンなどをトレーニングします。

Answer 7

これは素晴らしい質問です！私はいくつかの考えを持っています：

ここでは、動的タイムワーピングが興味深いツールになる可能性があります。3つのクラスの「テンプレート」を確立し、DTWを使用すると、テンプレートと信号の「チャンク」との相関関係を確認できます（信号をたとえば.5秒ビット、つまり0〜.5に分割します）。 1-.6 .2-.7 ...）。私は加速度計データを使用した歩行分析に似たものを使用しましたが、それはかなりうまく機能しました。

もう1つのオプションは、信号処理と機械学習の組み合わせアルゴリズムです。信号を再び「チャンク」に分割します。もう一度「テンプレート」を作成し（クラスごとに1ダース程度必要）、各チャンク/テンプレートのFFTを取得し、ナイーブベイズ分類器（または別のML分類器、ただしNBはそれをカットする必要があります）を使用して、それぞれを分類します。あなたの3つのクラス。歩行データでもこれを試しましたが、98％以上の精度を取得し、比較的複雑な信号で再現率を上げることができました。これがどのように機能するか教えてください、それは非常にエキサイティングな問題です。

Answer 8

良い結果が得られる可能性が非常に高い 1 つのアプローチは、カーブ フィッティングです。

• 連続波を間隔に分割します (おそらく、間隔の境界を qrs 群の鋭いピークの間の約半分にするのが最善です)。一度に 1 つの間隔のみを考慮してください。

• 心電図曲線のすべての可能なバリエーションを近似するために使用できるモデル関数を定義します。これは、最初に思われるほど難しくありません。モデル関数は、各波の原点 (t_)、振幅 (a_)、および幅 (w_) のパラメーターを持つ 3 つの関数の合計として構築できます。

     f_model(t) = a_p   *  f_p  ((t-t_p  )/w_p) + 
                  a_qrs *  f_qrs((t-t_qrs)/w_qrs) +
                  a_t   *  f_t  ((t-t_t  )/w_t)
  

  関数f_p(t)、f_qrs(t)、f_t(t)は、3 つの波のそれぞれをモデル化するために使用できる単純な関数です。

• フィッティング アルゴリズム (Levenberg-Marquardt-Algorithm http://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithmなど) を使用して、フィッティング パラメータ a_p、t_p、w_p、a_qrs、t_qrs、w_qrs、a_t を決定します。 、t_t、各間隔のデータセットの w_t。

  パラメーター t_p、t_qrs、および t_p は、関心のあるパラメーターです。

Answer 9

まず、標準的な心電図波のさまざまな成分が、任意のプロットから欠落している可能性があります。このようなプロットは一般的に異常であり、通常は何らかの問題を示していますが、それらが存在することを約束することはできません。

第二に、それらを認識することは、特に何かがうまくいかない場合には、科学と同じくらい芸術です。

私のアプローチは、コンポーネントを識別するためにニューラルネットワークをトレーニングしようとすることかもしれません。過去30秒間のデータを正規化して、最低点が0、最高点が1.0になるように指定すると、11個の出力が得られます。異常評価ではなかった出力は、最後の10秒間の重み付けになります。0.0は現在から-10秒であり、1.0は現在を意味します。出力は次のようになります。

1. 最新のP波が始まった場所
2. 最新のP波が終わった場所
3. 最新のP波の異常評価。1つの極値が「不在」です。
4. 最新のQRS群が始まった場所
5. 最新のQRS群のQ部分がR部分に変わったところ。
6. 最新のQRS群のR部分がS部分に変わったところ。
7. 最新のQRS群が終了した場所。
8. 最新のQRS群の異常評価で、極端なものが1つ「存在しない」。
9. 最新のT波が始まった場所。
10. 最新のT波が終わった場所。
11. 最新のT波の異常評価で、1つの極値が「不在」です。

人々が提案した他の種類の分析でこれを再確認するか、ニューラルネットワークの出力と一緒にそれらの他の種類の分析を使用してあなたの答えを与えるかもしれません。

もちろん、ニューラルネットワークのこの詳細な説明を規範的なものと見なすべきではありません。たとえば、必ずしも最適な出力を選択したわけではないと確信しています。たとえば、それらが何であるかについていくつかのアイデアを投げかけただけです。

Answer 10

「ウェーブレット変換」は関連キーワードかもしれません。私はかつて、この技術を使用してノイズの多い心電図のさまざまな心拍フェーズを検出した誰かのプレゼンテーションに参加したことがあります。

私の限られた理解に関する限り、それはフーリエ変換に似ていますが、あなたの場合は心拍の形をしたパルスの（スケーリングされた）コピーを使用しています。

Answer 11

ウェーブレットは、ピークが「異なるサイズ」であるこのタイプのデータのピークを特定するための最適なツールであることが示されています。ウェーブレットのスケーリング プロパティにより、このタイプのマルチスケール ピーク検出に理想的なツールになります。これは非定常信号のように見えるため、一部の人が示唆しているように DFT を使用することは適切なツールではありませんが、これが探索的プロジェクトである場合は、信号のスペクトルを使用して調べることができます (本質的に自己相関の FFT を使用して推定されます)。シグナル。）

これは、いくつかのピーク検出方法をレビューしている優れた論文です。これは、開始するのに適した場所です。

-ポール

Answer 12

私はお互いの回答を完全に読んだことはありませんが、それらをスキャンしたところ、フーリエ変換を見てこれらの波をセグメント化することを誰も勧めていないことに気付きました.

私には、数学における高調波解析の明確なアプリケーションのように思えます。私が見逃しているかもしれないいくつかの微妙な点があるかもしれません。

離散フーリエ変換係数は、離散時間信号を構成するさまざまな正弦波成分の振幅と位相を示します。これは、基本的に、問題が見つけたい状態です。

ここで何かが足りないかもしれませんが...