JP3735774B2 - Biological signal processing method - Google Patents

Description

【０００１】
【発明の属する技術分野】
本発明は、生体信号処理方法、詳しくは、ICG（Impedance Cardio Graphy）や脈波などの生体信号を処理する生体信号処理方法に関する。
【０００２】
【従来の技術】
心拍出量を測定する方法として、血流中に既知の熱量を加えて血流温度変化から血流量を測定する熱希釈法があるが、この方法は侵襲的で患者に大きな負担をかける方法といわれている。患者の負担が軽く非侵襲的に測定する方法としては、超音波法とインピーダンス法が知られている。しかし、超音波法は、測定に熟練した検査技師の技術を要し、安静状態での短時間測定しかできないという問題がある。一方、インピーダンス法は、胸部に電極を装着するだけで簡便に測定できる方法である。
【０００３】
インピーダンス法は、以下の原理に基づいた測定方法である。血液は他の組織に比べて電気的に良導体であり、血液が心臓から送り出されると胸部のインピーダンスに変化が生じる。このため、胸部インピーダンス波形は心臓から送り出される血流量を反映したものとなる。ＩＣＧは、この胸部インピーダンス波形の一次微分波のことをいう。インピーダンス波形は、胸部に交流定電流を加えることにより電圧応答波形として測定することができ、また、ＩＣＧから左室駆出時間などの心周期パラメータや心拍出量を求めることができる。このように、インピーダンス法は、非侵襲的で簡便な測定方法である。
【０００４】
【発明が解決しようとする課題】
しかしながら、実際の胸部インピーダンス波形には心拍成分以外に、呼吸、体動、筋電、ハム等の様々な雑音成分が含まれており、必要とする心拍成分がこれらのノイズ成分に埋もれてしまうことが多い。
【０００５】
実際のインピーダンス測定においては、安静仰臥状態で呼吸を停止させるなど、測定条件を規定している。また一方では、フィルタを用いて雑音除去をする試みもなされているが、必要な成分と不要な成分のスペクトラムが重なっていて、うまく分離、除去できないのが現状である。
【０００６】
また、心周期分画においても、雑音の混入やフィルタによる波形歪のために正しい区分点認識ができず、オフライン処理によって平均波形を求めてから手動計測によって各パラメータを測定しようとしているのが現状である。
【０００７】
本発明は、以上のような現状の打開策としてなされたもので、胸部のインピーダンス測定において必要な心拍成分が呼吸や体動による雑音成分に埋もれてしまい本来のＩＣＧが得られずまた心周期に関する正確な測定ができなくなるという問題を解決することを目的とする。
【０００８】
【課題を解決するための手段】
本発明の生体信号処理方法は、生体信号の原波形に近似させるよう心拍間隔及び心臓の収縮時間を基に線形結合モデルによる推定波形を作成し、前記原波形から前記推定波形を減算し、この減算結果をバンドパスフィルタに通して過渡応答成分を抽出し、前記推定波形に前記過渡応答成分を加算して出力波形を作成することを特徴とする。
【０００９】
前記出力波形から心周期の基点を設定し、該心周期の基点から心周期パラメータを算出する。
【００１４】
【発明の実施の形態】
以下、本発明の一実施形態を図面に基づいて説明する。
【００１５】
図１は、一実施形態に係るICG処理方法（以下、 ICG 測定方法ともいう。）の手順を示すフローチャートである。
【００１６】
本実施形態に係るＩＣＧ測定方法は、基本的には、ＩＣＧの原波形ＩＣＧ_measured(i) （以下、ＩＣＧ_m(i)と記述する。なお、ｉは、１心拍内での時間を示すサンプル番号である。）を算出するとともに、各サンプル毎に線形結合モデルによるＩＣＧの推定波形ＩＣＧ_estimated(i)（以下、ＩＣＧ_e(i)と記述する。）を作成し、さらに、過渡応答成分ＩＣＧ_transient(i)（以下、ＩＣＧ_t(i)と記述する。）を抽出し、ＩＣＧの出力波形ＩＣＧ(i) を作成するものである。さらに、出力波形ＩＣＧ(i) における区分点認識を行い、心周期パラメータを算出するものである。以下、原波形ＩＣＧ_m(i)の算出方法、推定波形ＩＣＧ_e(i)の作成方法、過渡応答成分ＩＣＧ_t(i)の抽出方法、出力波形ＩＣＧ(i) の作成方法、区分点認識方法、心周期パラメータの算出方法について順に説明する。
【００１７】
▲１▼原波形ＩＣＧ_m(i)の算出方法
胸部インピーダンス（胸部のベース抵抗Ｚ₀ 、ｄＺ）を測定し（ステップＳ１）、この測定結果に対して微分処理を行って原波形ＩＣＧ_m(i)（＝ｄＺ／ｄｔ(i) ）を算出する（ステップＳ２）。
【００１８】
▲２▼推定波形ＩＣＧ_e(i)の作成方法
推定波形ＩＣＧ_e(i)は、心拍間隔を心電図測定（ステップＳ３）による心電図（図２参照）のＲ−Ｒ間隔（ＲＲ）から求めた場合（ステップＳ４）、線形結合モデルによって下記式(1) で表され、この線形結合モデルによってＩＣＧが推定される（ステップＳ５）。
【００１９】

ここで、ｍは、心拍間隔ＲＲを基本波とする高調波の次数を表し、Ｘ_cos ( ｍ，ｉ）及びＸ_sin ( ｍ，ｉ）は、下記式(2) 、(3) で示すように、それぞれ高調波の余弦波成分及び正弦波成分である。なお、Ｘ_cos ( ｍ，ｉ）及びＸ_sin ( ｍ，ｉ）は、ｃｏｓ，ｓｉｎテーブルに書き込まれ、心拍が検出されるたびにＲＲにしたがって更新される（ステップＳ６）。
【００２０】
Ｘ_cos ( ｍ，ｉ）＝ｃｏｓ｛２πｍｉ／ＲＲ｝ (2)
Ｘ_sin ( ｍ，ｉ）＝ｓｉｎ｛２πｍｉ／ＲＲ｝ (3)
また、Ｗ_cos ( ｍ，ｉ）及びＷ_sin ( ｍ，ｉ）は、各高調波の余弦波成分及び正弦波成分の大きさを表す重み係数である。
【００２１】
重み係数Ｗ_cos ( ｍ，ｉ）、Ｗ_sin ( ｍ，ｉ）は、原波形ＩＣＧ_m(i)から推定波形ＩＣＧ_e(i)を減算する処理（ステップＳ７）によって得た差分ｅ(i) （＝ＩＣＧ_m(i)−ＩＣＧ_e(i)）の大きさが小さくなるように、ＬＭＳ法（最急降下法）を用いて下記式(4) 、(5) にしたがって毎サンプル更新される（ステップＳ８）。
【００２２】

ここで、μ( ｍ）は、重み係数Ｗ_cos ( ｍ，ｉ）、Ｗ_sin ( ｍ，ｉ）の修正量を決める修正係数であり、オンラインもしくはオフライン処理により得られる１心拍分の原波形ＩＣＧ_m の振幅スペクトルに応じて下記式(6) で表される。
【００２３】
μ( ｍ）＝μ0 ・Amplitude ｛ＩＣＧ_m(ｍ) ｝ (6)
（μ0 ＞０）
ここで、μは、正の固定値であるμ0 そのものを用いてもよいが、上記のように１心拍分の原波形ＩＣＧ_m の振幅スペクトルに応じて各高調波ごとに設定することにより、収束速度、収束後の安定性が向上する。
【００２４】
▲２▼’心拍間隔と収縮時間を用いた推定波形ＩＣＧ_e(i)の作成方法
▲２▼の推定波形の作成では、心拍間隔に基づく線形結合モデルを用いているため、日内変動のように長いスパンで起こる心拍間隔（心拍数）の変動に対しても追従できる能力を有している。しかしながら、心拍間隔の短縮、延長に比例して波形も短縮、延長することを仮定しているこのモデルでは、運動や呼吸によって短いスパンで心拍間隔に大きな変動が起きた場合は追従できない。なぜならば、実際の生体信号においては心臓の収縮時間や拡張時間が心拍間隔に比例しないからである。心電図を例にとって見れば分かり易いが、心拍間隔の短縮、延長に対して、収縮時間（心電図の場合、Ｑ波からＴ波までの時間）の短縮、延長の度合は小さく、拡張時間の短縮、延長の度合は逆に大きいのが一般的で、ＩＣＧや脈波についても全く同様である。
【００２５】
つまり、心拍間隔の変動に対して実際の波形が時間軸方向に単純に伸び縮みさせた形にならないため、１心拍ごとの心拍間隔を用いたとしても、変動が大きい場合にはうまく適応できず、波形の急峻な部分がうまく推定できなくなるのである。
【００２６】
これに対し、心拍間隔だけでなく、心臓の収縮時間を用い、線形結合モデルの入力信号である正弦波ならびに余弦波の位相子部分を前半部と後半部の二つに分けることによって心拍変動に伴う波形の短縮、延長をより実際の生体信号に即した形になるようにし、波形の急峻な部分も再生可能にする。
【００２７】
以下、具体的に説明する。
【００２８】
このＩＣＧ測定方法では、線形結合モデルの入力信号として用いていた基本波とその高調波からなる
Ｘ_cos ( ｍ，ｉ）＝ｃｏｓ｛２πｍｉ／ＲＲ｝
Ｘ_sin ( ｍ，ｉ）＝ｓｉｎ｛２πｍｉ／ＲＲ｝
の位相子部分について、次のように１周期を前半部（収縮期部分）と後半部（拡張期部分）の二つに分ける（ステップＳ６’参照）。
【００２９】

ここで、心電図を用いて前半部と後半部を分ける場合、ＲＲｆ＝ＱＴ×α （１＜α＜２）、ＲＲａ＝ＲＲ−ＲＲｆとする。なお、ＱＴは、図２に示すように、心電図のＱ波形からＴ波形の終了部分までの時間である。
【００３０】
その他の構成については、上述した実施形態と同様であるため説明は省略する。
【００３１】
▲３▼過渡応答成分ＩＣＧ_t(i)の抽出方法
過渡応答成分ＩＣＧ_t(i)は、下記式(7) で示すように、原波形ＩＣＧ_m(i)から推定波形ＩＣＧ_e(i)を減算した差分ｅ(i) の信号をバンドパスフィルタに通すことによって得る（ステップＳ９）。
【００３２】
ＩＣＧ_t(i)＝ＢＰＦ｛ＩＣＧ_m(i)−ＩＣＧ_e(i)｝ (7)
心拍間隔を用いた線形結合モデルは、心拍が検出されるたびに基本波、高調波を更新するため、時間方向に対する応答性に優れているが、心拍同期成分が繰返波形であることを仮定したモデルであるため、不整脈のような突発的な波形の変化に対する応答性が良くない。そのため、原波形ＩＣＧ_m(i)の信号を直接バンドパスフィルタに通すことが考えられるが、このようにすると波形歪が生じることになる。そこで、原波形ＩＣＧ_m(i)から推定波形ＩＣＧ_e(i)を減算することによって原波形ＩＣＧ_m(i)から心拍同期成分を分離しておき、その残留成分である過渡応答成分と雑音成分とをバンドパスフィルタに通すことによって、過渡応答成分を抽出するようにしている。
【００３３】
▲４▼出力波形ＩＣＧ(i) の作成方法
出力波形ＩＣＧ(i) は、下記式(8) で示すように、心拍同期成分である推定波形ＩＣＧ_e(i)に過渡応答成分ＩＣＧ_t(i)を加算することによって得る（ステップＳ１０）。
【００３４】
ＩＣＧ(i) ＝ＩＣＧ_e(i)＋ＩＣＧ_t(i) (8)
▲５▼区分点認識方法
出力波形ＩＣＧ(i) から、図３に示すように、心周期パラメータとしての前駆出期ＰＥＰ、左室駆出時間ＬＶＥＴ、等容性拡張期ＩＶＲＰ及び充満時間ＦＴを求めるために、僧帽弁の閉鎖を示すＡポイント、大動脈弁の開放を示すＢポイント、ＩＣＧ(i) の負極性の最大値を示す（ｄＺ／ｄｔ）max ポイント、大動脈弁の閉鎖を示すＸポイント及び僧帽弁の開放を示すＯポイントをそれぞれ求める（ステップＳ１１）。
【００３５】
各ポイントの認識方法を以下に示す。
【００３６】
(i) （ｄＺ／ｄｔ）max ポイント
まず、前心拍の負極性最大値（ｄＺ／ｄｔ）max を参照してその７０〜８０％をトリガレベルに設定し、このトリガレベルを上から下に横切るポイントを検索して心拍検出を行う。次に、この心拍検出ポイントから時間軸の正方向に負極性最大値を検索し、この負極性最大値をもつポイントを（ｄＺ／ｄｔ）max ポイントとする。
【００３７】
(ii)Ｂポイント
検出した負極性最大値（ｄＺ／ｄｔ）max の１０％をトリガレベルに設定し、（ｄＺ／ｄｔ）max ポイントから時間軸の負方向にＢポイントを検索し、トリガレベルを下から上に横切るポイントをＢポイントとする。
【００３８】
(iii) Ｘポイント
（ｄＺ／ｄｔ）max ポイントから時間軸の正方向に正極性最大値を検索し、最大値の更新が一定時間されなくなったら、それまでに保持している最大値のポイントをＸポイントとする。この一定時間は、雑音による変曲点（極大値）を最大値と誤って認識しないように設けるものであり、４０〜５０ｍｓに設定する。なお、ここでの最大値とは、１心拍中の最大値ではなく、（ｄＺ／ｄｔ）max ポイントから時間軸の正方向に対して進んでいった場合の最大値をいう。
【００３９】
(iv)Ａポイント
（ｄＺ／ｄｔ）max ポイントから時間軸の負方向に最大傾斜を検索し、最大傾斜の更新が一定時間されなくなったら、それまでに保持している最大傾斜ポイントの値を最大傾斜値とする。この一定時間は、上記Ｘポイント認識と同様、雑音による誤認識防止のためのものであり、４０〜５０ｍｓに設定する。また、ここでの最大傾斜値とは、１心拍中の最大傾斜値ではなく、（ｄＺ／ｄｔ）max ポイントから時間軸の負方向に対して進んでいった場合の最大傾斜ポイントの傾斜値をいう。次に、上記Ｂポイントから時間軸の負方向にＡポイントを検索し、傾斜値が最大傾斜値の１０％を下回るポイントをＡポイントとする。
【００４０】
(v) Ｏポイント
Ｘポイントから時間軸の正方向にマスク期間をおき、そこから時間軸の正方向に最小値を検索して最小値の更新が一定時間されなくなったら、それまで保持している最小値のポイントをＯポイントとする。この一定時間は、上記Ｘポイント認識、Ａポイント認識と同様、雑音による誤認識防止のためのものであり、４０〜５０ｍｓに設定する。また、ここでの最小値とは、１心拍中の最小値ではなく、Ｘポイントからマスク期間をおき、そこから時間軸の負方向に対して進んでいった場合の極小値をいう。マスク期間は、一定時間に設定されているものではなく、心拍数が高ければ短く、心拍数が低ければ長くなるよう、心拍数に応じて設定する。このマスク期間Ｔmaskを決定する一例を下記式(9) で示す。
【００４１】
Ｔmask＝１００−ＨＲ／β (9)
（β＝３〜４）
ここで、ＨＲは心拍数（bpm)を表す。
【００４２】
▲６▼心周期パラメータの算出方法
区分点認識によって得られた各ポイントから下記式(10)、(11)、(12)、(13) により心周期パラメータを算出する（ステップＳ１２）。
【００４３】
前駆出期：ＰＥＰ＝Ｔ［Ｂ］−Ｔ［Ａ］ (10)
左室駆出時間：ＬＶＥＴ＝Ｔ［Ｘ］−Ｔ［Ｂ］ (11)
等容性拡張期：ＩＶＲＰ＝Ｔ［Ｏ］−Ｔ［Ｘ］ (12)
充満時間：ＦＴ＝Ｔ［Ａ’］−Ｔ［Ｏ］ (13)
ここで、Ｔ［］は［］内のポイントの時刻を表し、式(13)におけるＡ’は次の心拍のＡポイントを表す。
【００４４】
さらに、ＬＶＥＴ、（ｄＺ／ｄｔ）max からKubicek の心拍出量の導出式を用いることにより、下記式(14)、(15) に示すように１心拍の心拍出量ＳＶ(Stroke Volume)(ｍｌ）、１分間の心拍出量ＣＯ(Cardiac Output)（ｍｌ）が導出できる。
【００４５】
ＳＶ＝ρ（Ｌ／Ｚ₀ ）² ・ＬＶＥＴ・（ｄＺ／ｄｔ）max (14)
ＣＯ＝ＳＶ・ＨＲ (15)
ここで、ρは血液の粘性抵抗（Ω・ｃｍ）、Ｌは測定電極間距離（ｃｍ）、Ｚ₀ は上記した胸部のベース抵抗（Ω）を表す。なお、ＬＶＥＴの単位は（ｓ）で、（ｄＺ／ｄｔ）max の単位は（Ω／ｓ）である。
【００４６】
次に、本実施形態のＩＣＧ測定方法による効果を以下に説明する。
【００４７】
(1) 運動負荷試験への適用例
図４に、本発明の比較例としてエルゴメータを用いて運動負荷試験を行いながら測定したＩＣＧを示す。波形はサンプリング周波数２００Ｈｚにて収録したものである。測定された波形には、呼吸や体動による雑音のほか高周波雑音が混入している。
【００４８】
このＩＣＧの信号をバンドパスフィルタに通した結果を図５〜図７に示す。図５〜図７は、中心周波数ｆ₀ が３Ｈｚ、Ｑ値が０．２〜０．５のバンドパスフィルタを用いた結果である。Ｑ値を大きくして帯域を狭めていくにしたがい雑音成分が抑制されていくが、逆に波形の減衰や歪み（Ａ波レベルの低下、Ｘ波レベルの上昇）が大きくなっていく様子が確認できる。このように、従来のフィルタ処理では雑音抑制と信号成分保持とは相反する関係となるため、両者を両立させることは非常に困難であった。
【００４９】
これに対し、本実施形態による線形結合モデルを用いてＩＣＧの推定を行った結果を図８に示す。線形結合モデルの次数は１６次（Ｍ＝１６）とし、１６倍高調波までの推定を行った。このような心拍同期式の推定を行うことにより、心拍とは非同期である雑音成分が排除でき、かつ波形の減衰や歪みもなく安定した心拍同期成分が推定、抽出できていることが分かる。
【００５０】
(2) 心拍間隔（心拍数）変動が大きい場合への適用例
短時間で大きく心拍間隔変動が起きているＩＣＧを図９に、心拍間隔だけを用いた場合の推定波形を図１０に、心拍間隔と収縮時間を用いた場合の推定波形を図１１に示す。
【００５１】
図９に示すＩＣＧは短時間で心拍間隔の短縮（▲１▼〜▲３▼）と延長（▲３▼〜▲５▼）が起きている。このような場合、心拍間隔だけを用いて推定を行うと、雑音成分の抑制はできるものの、図１０に示すように一部の信号成分が減衰し、波形全体も鈍りや横幅の広がりが発生して再現性が悪くなってしまう。一方、心拍間隔と収縮時間を用いて推定を行うと、図１１に示すように図１０で減衰していた部分も保持でき、波形全体の再現性が非常に高いことが分かる。
【００５２】
(3) 不整脈を含む波形への適用例
図１２に、不整脈を含む心電図を示す。左から４拍目で早期期外収縮が起き、リズムとは異なった心電図が現れている。この心電図と同時に測定したＩＣＧを図１３に示す。ＩＣＧも心電図と同様に不整脈のときの波形が全く異なったものとなっているのが確認できる。
【００５３】
図１３の測定波形から線形結合モデルにより推定したＩＣＧを図１４に示す。線形結合モデルは波形の繰返を仮定したモデルであるため、不整脈のような波形の急激な変化には追従できず、４拍目の推定結果は測定波形と全く異なったものになっている。
【００５４】
これに対し、図１３の測定波形から線形結合モデルの出力を減算し、中心周波数が３Ｈｚ、Ｑ値が０．３のバンドパスフィルタを通して過渡応答成分を抽出し、再び線形結合モデルの出力を加算した結果を図１５に示す。過渡応答成分の抽出により、不整脈にも対応しながら安定したＩＣＧが得られることが分かる。
【００５５】
(4) ＩＣＧの自動認識、測定例
運動負荷試験におけるＩＣＧの測定波形からＩＣＧの推定を行った後、区分点認識によって心周期の各基点（Ａポイント等）を自動認識させた結果を図１６に示す。図１６に示すように、運動負荷中においても心周期分画ができ、心臓の収縮期、拡張期の様子が観測できることが分かる。
【００５６】
区分点認識によって求めた各基点から左室駆出時間を自動測定し、運動負荷中の左室駆出時間のトレンドグラフが図１７のように得られた。また、左室駆出時間、（ｄＺ／ｄｔ）max 、心拍数からKubicek の導出式を用いて心拍出量を算出した結果、図１８のようになった。運動負荷試験は、負荷上昇率２０Ｗ／分で直線的に負荷を上昇させていくランプ負荷をかけて行ったものであるが、被験者の心機能がそれに追従していく様子が図１７、１８から確認でき、運動負荷中においても心機能を測定できるようになっている。
【００５７】
なお、上述した実施形態において、心拍間隔については、心拍に起因する生体信号波形であれば心電図以外から求めても構わない。また、ＩＣＧ以外の生体信号が得られない場合は、原波形ＩＣＧ_m(i)そのものをバンドパスフィルタに通すことによって心拍検出を行い、心拍間隔を求めることもできる。この場合に用いるバンドパスフィルタは、心拍を検出できれば十分であるため波形歪を生じるようなフィルタであっても一向に構わない。
【００５８】
以上説明したように、上述した実施形態によると、生体信号（ＩＣＧ）の原波形ＩＣＧ_m(i)に近似させるよう心拍間隔ＲＲを基に線形結合モデルによる推定波形ＩＣＧ_e(i)を作成するようにしたため、例えば胸部のインピーダンス測定において必要な心拍成分が呼吸や体動による雑音成分に埋もれてしまう不具合を解消可能になる。
【００５９】
また、推定波形ＩＣＧ_e(i)は、心拍間隔ＲＲ以外に収縮時間ＱＴに基づいて作成されることから、心拍変動に伴う推定波形ＩＣＧ_e(i)の延長、短縮がより実際の生体信号波形ＩＣＧ_m(i)の延長、短縮に即した形になり、波形の急峻な部分も再生可能になる。
【００６０】
また、原波形ＩＣＧ_m(i)から推定波形ＩＣＧ_e(i)を減算し、この減算結果をバンドパスフィルタに通して過渡応答成分ＩＣＧ_t(i)を抽出することから、呼吸や体動などによる雑音成分を除去可能になる。
【００６１】
また、推定波形ＩＣＧ_e(i)に過渡応答成分ＩＣＧ_t(i)を加算して出力波形ＩＣＧ(i) を作成することから、不整脈にも対応しながら安定した出力波形ＩＣＧ(i) を得ることができる。
【００６２】
また、出力波形ＩＣＧ(i) が安定したものとなるため、出力波形ＩＣＧ(i) から心周期の基点（Ａポイント等）を認識することができ、さらに、心周期の基点から心周期パラメータ（前駆出期等）を算出することができる。
【００６３】
【発明の効果】
本発明によると、生体信号測定例えば胸部のインピーダンス測定において必要な心拍成分が呼吸や体動による雑音成分に埋もれる不具合を解消しまた心周期に関する正確な測定が可能になる。
【図面の簡単な説明】
【図１】本発明の一実施形態に係るＩＣＧ測定方法の手順を示すフローチャートである。
【図２】心電図におけるＲ−Ｒ間隔とＱ−Ｔ間隔の説明図である。
【図３】ＩＣＧにおける心周期の各基点と心周期パラメータの説明図である。
【図４】運動負荷試験中に測定したＩＣＧを表す。
【図５】図４図示のＩＣＧを中心周波数３Ｈｚ、Ｑ＝０．２のバンドパスフィルタで処理した結果の信号波形を表す。
【図６】図４図示のＩＣＧを中心周波数３Ｈｚ、Ｑ＝０．３のバンドパスフィルタで処理した結果の信号波形を表す。
【図７】図４図示のＩＣＧを中心周波数３Ｈｚ、Ｑ＝０．５のバンドパスフィルタで処理した結果の信号波形を表す。
【図８】図４図示のＩＣＧから線形結合モデルによって推定されたＩＣＧを表す。
【図９】短時間で大きく心拍間隔変動が起きているＩＣＧを表す。
【図１０】心拍間隔だけを用いた場合の推定波形を示す。
【図１１】心拍間隔と収縮時間を用いた場合の推定波形を示す。
【図１２】不整脈（早期期外収縮）発生時の心電図を表す。
【図１３】不整脈発生時のＩＣＧの測定波形を表す。
【図１４】不整脈発生時の線形結合モデルによるＩＣＧの推定結果の信号波形を表す。
【図１５】過渡応答成分抽出法の適用結果である出力波形を表す。
【図１６】運動負荷試験におけるＩＣＧの各基点の自動認識結果を表す。
【図１７】左室駆出時間のトレンドグラフである。
【図１８】心拍出量のトレンドグラフである。
【符号の説明】
ＩＣＧ 生体信号
ＩＣＧ_m(i) 原波形
ＲＲ 心拍間隔
ＩＣＧ_e(i) 推定波形
ＱＴ 収縮時間
ＩＣＧ_t(i) 過渡応答成分
ＩＣＧ(i) 出力波形
Ａ、Ｂ（ｄＺ／ｄｔ）max 、Ｘ、Ｏポイント 心周期の基点
ＰＥＰ、ＬＶＥＴ、ＩＶＲＰ、ＦＴ 心周期パラメータ[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a biological signal processing method, and more particularly to a biological signal processing method for processing biological signals such as ICG (Impedance Cardio Graphy) and pulse waves.
[0002]
[Prior art]
As a method of measuring cardiac output, there is a thermodilution method that measures blood flow from changes in blood flow temperature by adding a known amount of heat into the bloodstream, but this method is invasive and places a heavy burden on the patient It is said that. As a method for lightly and non-invasively measuring a patient, an ultrasonic method and an impedance method are known. However, the ultrasonic method requires the skill of an inspection engineer skilled in measurement, and there is a problem that only a short time measurement can be performed in a resting state. On the other hand, the impedance method is a method that can be measured simply by attaching electrodes to the chest.
[0003]
The impedance method is a measurement method based on the following principle. Blood is a better electrical conductor than other tissues, and when blood is pumped from the heart, the impedance of the chest changes. For this reason, the chest impedance waveform reflects the blood flow volume delivered from the heart. ICG refers to the first derivative wave of this chest impedance waveform. The impedance waveform can be measured as a voltage response waveform by applying an alternating constant current to the chest, and cardiac cycle parameters such as left ventricular ejection time and cardiac output can be obtained from ICG. Thus, the impedance method is a noninvasive and simple measurement method.
[0004]
[Problems to be solved by the invention]
However, the actual chest impedance waveform contains various noise components such as breathing, body movement, myoelectricity, and ham in addition to the heart rate component, and the required heart rate component is buried in these noise components. There are many.
[0005]
In actual impedance measurement, the measurement conditions are defined such as stopping breathing in a resting and supine state. On the other hand, attempts have been made to remove noise by using a filter, but the current situation is that the spectrum of necessary components and unnecessary components overlap and cannot be separated and removed well.
[0006]
Moreover, even in cardiac cycle fractionation, correct classification point recognition cannot be performed due to noise mixing or waveform distortion due to filters, and it is currently attempted to measure each parameter by manual measurement after obtaining an average waveform by offline processing. It is.
[0007]
The present invention has been made as a solution to the present situation as described above, and the heartbeat component necessary for measuring the impedance of the chest is buried in a noise component due to breathing and body movement, so that the original ICG cannot be obtained, and the cardiac cycle is concerned. The purpose is to solve the problem of inaccurate measurement.
[0008]
[Means for Solving the Problems]
The biological signal processing method of the present invention creates an estimated waveform based on a linear combination model based on the heartbeat interval and the heart contraction time so as to approximate the original waveform of the biological signal, and subtracts the estimated waveform from the original waveform. The subtraction result is passed through a band pass filter to extract a transient response component, and the transient response component is added to the estimated waveform to create an output waveform .
[0009]
A base point of the cardiac cycle is set from the output waveform, and a cardiac cycle parameter is calculated from the base point of the cardiac cycle.
[0014]
DETAILED DESCRIPTION OF THE INVENTION
Hereinafter, an embodiment of the present invention will be described with reference to the drawings.
[0015]
FIG. 1 is a flowchart showing a procedure of an ICG processing method (hereinafter also referred to as an ICG measurement method) according to an embodiment.
[0016]
The ICG measurement method according to the present embodiment is basically described as an ICG original waveform ICG _measured (i) (hereinafter referred to as ICG _m (i). Here, i is a sample indicating time within one heartbeat. And an ICG estimated waveform ICG _estimated (i) (hereinafter referred to as ICG _e (i)) by a linear combination model is created for each sample, and a transient response component ICG is generated. _Transient (i) (hereinafter referred to as ICG _t (i)) is extracted to create an ICG output waveform ICG (i). Further, the segment point recognition in the output waveform ICG (i) is performed to calculate the cardiac cycle parameter. Hereinafter, calculation method of original waveform ICG _m (i), creation method of estimated waveform ICG _e (i), extraction method of transient response component ICG _t (i), creation method of output waveform ICG (i), distinction point recognition method The calculation method of the cardiac cycle parameter will be described in order.
[0017]
(1) Calculation Method of Original Waveform ICG _m (i) Chest impedance (chest base resistance Z ₀ , dZ) is measured (step S1), and the measurement result is differentiated to obtain the original waveform ICG _m (i ) (= DZ / dt (i)) is calculated (step S2).
[0018]
▲ 2 ▼ creating estimated waveform ICG _e estimated waveform ICG _e (i) (i), when determined from the beat-to-beat ECG measurement RR interval of the electrocardiogram by (step S3) (see FIG. 2) (RR) (Step S4), expressed by the following equation (1) by a linear combination model, and ICG is estimated by this linear combination model (Step S5).
[0019]

Here, m represents the order of the harmonic having the heartbeat interval RR as a fundamental wave, and X _cos (m, i) and X _sin (m, i) are expressed by the following equations (2) and (3). And cosine wave component and sine wave component of the harmonic, respectively. X _cos (m, i) and X _sin (m, i) are written in the cos, sin table and updated according to RR each time a heartbeat is detected (step S6).
[0020]
X _cos (m, i) = cos {2πmi / RR} (2)
X _sin (m, i) = sin {2πmi / RR} (3)
W _cos (m, i) and W _sin (m, i) are weighting factors representing the magnitudes of the cosine wave component and sine wave component of each harmonic.
[0021]
The weighting factors W _cos (m, i) and W _sin (m, i) are obtained by subtracting the estimated waveform ICG _e (i) from the original waveform ICG _m (i) (step S7). Each sample is updated according to the following formulas (4) and (5) using the LMS method (the steepest descent method) so that the size of (= ICG _m (i) −ICG _e (i)) is reduced ( Step S8).
[0022]

Here, μ (m) is a correction coefficient that determines the correction amount of the weighting coefficients W _cos (m, i) and W _sin (m, i), and is an original waveform ICG for one heartbeat obtained by online or offline processing. It is expressed by the following formula (6) according to the amplitude spectrum of _m .
[0023]
μ (m) = μ0 ・ Amplitude {ICG _m (m)} (6)
(Μ0> 0)
Here, μ may be a positive fixed value μ0 itself, but it is converged by setting each harmonic according to the amplitude spectrum of the original waveform ICG _m for one heartbeat as described above. Speed and stability after convergence are improved.
[0024]
(2) Method for creating estimated waveform ICG _e (i) using heart rate interval and contraction time In creating the estimated waveform in item (2), a linear combination model based on the heart rate interval is used. It has the ability to follow fluctuations in heart rate intervals (heart rate) that occur over a long span. However, this model, which assumes that the waveform is shortened and lengthened in proportion to the shortening and lengthening of the heartbeat interval, cannot follow when a large fluctuation occurs in the heartbeat interval in a short span due to exercise or breathing. This is because in an actual biological signal, the heart contraction time and dilation time are not proportional to the heartbeat interval. Although it is easy to understand by taking an electrocardiogram as an example, the contraction time (in the case of an electrocardiogram, the time from the Q wave to the T wave) is shortened, the degree of extension is small, and the expansion time is shortened. In general, the extent of extension is large, and the same is true for ICG and pulse waves.
[0025]
In other words, since the actual waveform does not simply expand or contract in the time axis direction with respect to fluctuations in the heartbeat interval, even if the heartbeat interval for each heartbeat is used, it cannot be adapted well when the fluctuation is large. The steep part of the waveform cannot be estimated well.
[0026]
On the other hand, using not only the heartbeat interval but also the heart contraction time, the sine wave and cosine wave input signal of the linear combination model are divided into two parts, the first half and the second half. The shortening and extension of the accompanying waveform are made to conform to the actual biological signal, and the steep part of the waveform can be reproduced.
[0027]
This will be specifically described below.
[0028]
In this ICG measurement method, X _cos (m, i) = cos {2πmi / RR} consisting of the fundamental wave and its harmonics used as the input signal of the linear combination model.
X _sin (m, i) = sin {2πmi / RR}
In the phaser portion, one period is divided into two parts, a first half (systolic part) and a second half (diastolic part) as follows (see step S6 ′).
[0029]

Here, when the first half and the second half are separated using an electrocardiogram, RRf = QT × α (1 <α <2) and RRa = RR−RRf. As shown in FIG. 2, QT is the time from the Q waveform of the electrocardiogram to the end of the T waveform.
[0030]
Since other configurations are the same as those in the above-described embodiment, description thereof will be omitted.
[0031]
▲ 3 ▼ extraction methods transient response component ICG _t transient response component ICG _t (i) (i), as shown by the following formula (7), the original waveform ICG _m (i) from the estimated waveform ICG _e a (i) The subtracted difference e (i) is obtained by passing it through a band pass filter (step S9).
[0032]
ICG _t (i) = BPF {ICG _m (i) −ICG _e (i)} (7)
The linear combination model using the heart rate interval updates the fundamental and harmonics every time a heart rate is detected, so it has excellent response in the time direction, but it is assumed that the heart rate synchronization component is a repetitive waveform. Because of this model, the responsiveness to sudden waveform changes such as arrhythmia is not good. For this reason, it is conceivable to directly pass the signal of the original waveform ICG _m (i) through the band-pass filter, but if this is done, waveform distortion will occur. Accordingly, advance to separate the heartbeat-synchronous component from the original waveform ICG _m (i) by subtracting the estimated waveform ICG _e (i) from the original waveform ICG _m (i), the transient response component and a noise component which is a residual component Are passed through a band-pass filter to extract a transient response component.
[0033]
(4) Generation method of output waveform ICG (i) The output waveform ICG (i) is obtained by adding a transient response component ICG _t (i) to an estimated waveform ICG _e (i) which is a heartbeat synchronization component as shown in the following equation (8). ) Are added (step S10).
[0034]
ICG (i) = ICG _e (i) + ICG _t (i) (8)
(5) Classification point recognition method From the output waveform ICG (i), as shown in FIG. 3, the progenitor period PEP, left ventricular ejection time LVET, isovolumetric expansion period IVRP and filling time FT as cardiac cycle parameters are calculated. To determine, A point indicating mitral valve closure, B point indicating aortic valve opening, ICd (i) max negative value (dZ / dt) max point, X indicating aortic valve closure A point and an O point indicating the opening of the mitral valve are respectively obtained (step S11).
[0035]
The recognition method of each point is shown below.
[0036]
(i) (dZ / dt) max point First, with reference to the negative polarity maximum value (dZ / dt) max of the previous heartbeat, 70 to 80% is set as the trigger level, and this trigger level crosses from top to bottom. Search for points to detect heart rate. Next, the negative polarity maximum value is searched from the heartbeat detection point in the positive direction of the time axis, and the point having this negative polarity maximum value is defined as (dZ / dt) max point.
[0037]
(ii) Set the trigger level to 10% of the maximum negative polarity (dZ / dt) max detected at the B point, search the B point in the negative direction of the time axis from the (dZ / dt) max point, and set the trigger level. The point that crosses from the bottom to the top is the B point.
[0038]
(iii) Search for the maximum positive value from the X point (dZ / dt) max point in the positive direction of the time axis. When the maximum value is not updated for a certain period of time, the maximum value point held so far is X points. This fixed time is provided so that an inflection point (maximum value) due to noise is not erroneously recognized as the maximum value, and is set to 40 to 50 ms. Note that the maximum value here is not the maximum value during one heartbeat but the maximum value when proceeding from the (dZ / dt) max point in the positive direction of the time axis.
[0039]
(iv) Search the maximum slope in the negative direction of the time axis from the A point (dZ / dt) max point. When the maximum slope is not updated for a certain period of time, the maximum slope point value held so far is maximized. Inclination value. This fixed time is for preventing erroneous recognition due to noise, as in the case of X point recognition, and is set to 40 to 50 ms. The maximum inclination value here is not the maximum inclination value during one heartbeat, but the inclination value of the maximum inclination point when proceeding from the (dZ / dt) max point in the negative direction of the time axis. Say. Next, the A point is searched from the B point in the negative direction of the time axis, and the point where the slope value is less than 10% of the maximum slope value is set as the A point.
[0040]
(v) The mask period is set in the positive direction of the time axis from the O point X point, and when the minimum value is searched from there in the positive direction of the time axis and the update of the minimum value is not performed for a certain time, it is held until then. The point with the minimum value is the O point. This fixed time is for preventing erroneous recognition due to noise, as in the case of X point recognition and A point recognition, and is set to 40 to 50 ms. Further, the minimum value here is not the minimum value during one heartbeat but the minimum value when a mask period is set from the X point and then proceeds in the negative direction of the time axis. The mask period is not set to a fixed time, but is set according to the heart rate so that it is short when the heart rate is high and long when the heart rate is low. An example of determining this mask period Tmask is shown by the following equation (9).
[0041]
Tmask = 100−HR / β (9)
(Β = 3-4)
Here, HR represents the heart rate (bpm).
[0042]
(6) Calculation method of cardiac cycle parameters The cardiac cycle parameters are calculated by the following formulas (10), (11), (12), and (13) from the points obtained by the recognition of the dividing points (step S12).
[0043]
Precursor stage: PEP = T [B] -T [A] (10)
Left ventricular ejection time: LVET = T [X] -T [B] (11)
Isovolume expansion period: IVRP = T [O] -T [X] (12)
Charging time: FT = T [A ′] − T [O] (13)
Here, T [] represents the time of the point in [], and A ′ in equation (13) represents the A point of the next heartbeat.
[0044]
Furthermore, by using the derivation of Kubicek's cardiac output from LVET, (dZ / dt) max, as shown in the following formulas (14) and (15), the cardiac output SV (Stroke Volume) of one heartbeat (ml) One minute cardiac output CO (Cardiac Output) (ml) can be derived.
[0045]
SV = ρ (L / Z ₀ ) ² · LVET · (dZ / dt) max (14)
CO = SV · HR (15)
Here, ρ represents the blood viscosity resistance (Ω · cm), L represents the distance between measurement electrodes (cm), and Z ₀ represents the above-described chest base resistance (Ω). The unit of LVET is (s), and the unit of (dZ / dt) max is (Ω / s).
[0046]
Next, effects of the ICG measurement method of the present embodiment will be described below.
[0047]
(1) Application Example to Exercise Load Test FIG. 4 shows ICG measured while performing an exercise load test using an ergometer as a comparative example of the present invention. The waveform is recorded at a sampling frequency of 200 Hz. The measured waveform contains high-frequency noise in addition to noise caused by breathing and body movement.
[0048]
The results of passing this ICG signal through a band-pass filter are shown in FIGS. 5 to 7 show results of using a band-pass filter having a center frequency f ₀ of 3 Hz and a Q value of 0.2 to 0.5. As the Q value is increased and the band is narrowed, noise components are suppressed, but conversely, it is confirmed that the waveform attenuation and distortion (decrease in A wave level, increase in X wave level) increase. it can. As described above, in conventional filter processing, noise suppression and signal component holding are in a contradictory relationship, and it has been extremely difficult to achieve both.
[0049]
On the other hand, FIG. 8 shows the results of ICG estimation using the linear combination model according to the present embodiment. The order of the linear combination model was 16th order (M = 16), and estimation up to the 16th harmonic was performed. It can be seen that by performing such heart rate synchronization type estimation, a noise component that is asynchronous with the heart rate can be eliminated, and a stable heart rate synchronization component without waveform attenuation or distortion can be estimated and extracted.
[0050]
(2) Application example when the heart rate interval (heart rate) fluctuation is large FIG. 9 shows an ICG in which heart rate interval fluctuation has occurred in a short time, and FIG. 10 shows an estimated waveform when only the heart rate interval is used. FIG. 11 shows an estimated waveform when the interval and the contraction time are used.
[0051]
In the ICG shown in FIG. 9, the heartbeat interval is shortened (1) to (3) and extended (3) to (5) in a short time. In such a case, if the estimation is performed using only the heartbeat interval, the noise component can be suppressed. However, as shown in FIG. 10, some signal components are attenuated, and the entire waveform is dulled or widened. Reproducibility is worse. On the other hand, when the estimation is performed using the heartbeat interval and the contraction time, it can be understood that the portion attenuated in FIG. 10 can be retained as shown in FIG. 11, and the reproducibility of the entire waveform is very high.
[0052]
(3) Application example to waveform including arrhythmia FIG. 12 shows an electrocardiogram including arrhythmia. The early extrasystole occurs at the 4th beat from the left, and an electrocardiogram different from the rhythm appears. FIG. 13 shows the ICG measured simultaneously with the electrocardiogram. It can be confirmed that ICG has a completely different waveform at the time of arrhythmia as in the electrocardiogram.
[0053]
FIG. 14 shows ICG estimated by the linear combination model from the measured waveform of FIG. Since the linear combination model is a model that assumes repetition of the waveform, it cannot follow a rapid change in the waveform such as an arrhythmia, and the estimation result at the fourth beat is completely different from the measurement waveform.
[0054]
On the other hand, the output of the linear combination model is subtracted from the measured waveform of FIG. 13, a transient response component is extracted through a bandpass filter having a center frequency of 3 Hz and a Q value of 0.3, and the output of the linear combination model is added again. The results are shown in FIG. It can be seen that by extracting the transient response component, a stable ICG can be obtained while dealing with arrhythmia.
[0055]
(4) Automatic recognition of ICG, measurement example After ICG estimation from ICG measurement waveform in exercise load test, the result of automatically recognizing each base point (A point etc.) of the cardiac cycle by segment point recognition is shown in FIG. Shown in As shown in FIG. 16, it can be seen that the cardiac cycle can be fractionated even during exercise, and that the systolic and diastolic states of the heart can be observed.
[0056]
The left ventricular ejection time was automatically measured from each base point obtained by segment point recognition, and a trend graph of the left ventricular ejection time during exercise load was obtained as shown in FIG. Further, as a result of calculating the cardiac output from the left ventricular ejection time, (dZ / dt) max, and the heart rate using the Kubicek derivation formula, the result is as shown in FIG. The exercise load test was performed by applying a ramp load that linearly increases the load at a load increase rate of 20 W / min. From FIG. 17 and FIG. It can be confirmed, and cardiac function can be measured even during exercise.
[0057]
In the above-described embodiment, the heartbeat interval may be obtained from other than the electrocardiogram as long as it is a biological signal waveform resulting from the heartbeat. Further, when a biological signal other than ICG cannot be obtained, heartbeat detection can be performed by passing the original waveform ICG _m (i) itself through a bandpass filter to obtain a heartbeat interval. The band-pass filter used in this case is sufficient if it can detect the heartbeat, and may be a filter that causes waveform distortion.
[0058]
As described above, according to the above-described embodiment, the estimated waveform ICG _e (i) based on the linear combination model is created based on the heartbeat interval RR so as to approximate the original waveform ICG _m (i) of the biological signal (ICG). For this reason, for example, it is possible to eliminate a problem that a heartbeat component necessary for impedance measurement of the chest is buried in a noise component due to breathing or body movement.
[0059]
Further, since the estimated waveform ICG _e (i) is created based on the contraction time QT in addition to the heart rate interval RR, the estimated waveform ICG _e (i) accompanying the heart rate variability is more extended and shortened than the actual biological signal waveform. ICG _m (i) is extended and shortened, and the steep part of the waveform can be reproduced.
[0060]
Further, the estimated waveform ICG _e (i) is subtracted from the original waveform ICG _m (i), and the subtraction result is passed through a band-pass filter to extract the transient response component ICG _t (i). The noise component due to can be removed.
[0061]
Further, since the output waveform ICG (i) is created by adding the transient response component ICG _t (i) to the estimated waveform ICG _e (i), a stable output waveform ICG (i) is obtained while dealing with arrhythmia. be able to.
[0062]
Further, since the output waveform ICG (i) becomes stable, the base point (A point, etc.) of the cardiac cycle can be recognized from the output waveform ICG (i), and the cardiac cycle parameter ( Precursor period etc.) can be calculated.
[0063]
【The invention's effect】
According to the present invention, it is possible to eliminate a problem that a heartbeat component necessary for biological signal measurement, for example, chest impedance measurement, is buried in a noise component due to respiration or body movement, and to accurately measure the cardiac cycle.
[Brief description of the drawings]
FIG. 1 is a flowchart showing a procedure of an ICG measurement method according to an embodiment of the present invention.
FIG. 2 is an explanatory diagram of an RR interval and a QT interval in an electrocardiogram.
FIG. 3 is an explanatory diagram of each base point of a cardiac cycle and cardiac cycle parameters in ICG.
FIG. 4 represents ICG measured during exercise testing.
5 shows a signal waveform as a result of processing the ICG shown in FIG. 4 with a bandpass filter having a center frequency of 3 Hz and Q = 0.2.
6 shows a signal waveform as a result of processing the ICG shown in FIG. 4 with a band-pass filter having a center frequency of 3 Hz and Q = 0.3.
7 shows a signal waveform as a result of processing the ICG shown in FIG. 4 with a bandpass filter having a center frequency of 3 Hz and Q = 0.5.
8 represents an ICG estimated by a linear combination model from the ICG shown in FIG.
FIG. 9 shows an ICG in which heartbeat interval fluctuations occur greatly in a short time.
FIG. 10 shows an estimated waveform when only the heartbeat interval is used.
FIG. 11 shows an estimated waveform when a heartbeat interval and a contraction time are used.
FIG. 12 represents an electrocardiogram when an arrhythmia (early extrasystole) occurs.
FIG. 13 shows an ICG measurement waveform when an arrhythmia occurs.
FIG. 14 shows a signal waveform of an estimation result of ICG by a linear combination model when arrhythmia occurs.
FIG. 15 shows an output waveform that is a result of applying the transient response component extraction method.
FIG. 16 represents an automatic recognition result of each base point of ICG in the exercise load test.
FIG. 17 is a trend graph of left ventricular ejection time.
FIG. 18 is a trend graph of cardiac output.
[Explanation of symbols]
ICG Biological signal ICG _m (i) Original waveform RR Heartbeat interval ICG _e (i) Estimated waveform QT Contraction time ICG _t (i) Transient response component ICG (i) Output waveform A, B (dZ / dt) max, X, O Point Base point of heart cycle PEP, LVET, IVRP, FT Heart cycle parameters

Claims (2)

生体信号の原波形に近似させるよう心拍間隔及び心臓の収縮時間を基に線形結合モデルによる推定波形を作成し、
前記原波形から前記推定波形を減算し、この減算結果をバンドパスフィルタに通して過渡応答成分を抽出し、
前記推定波形に前記過渡応答成分を加算して出力波形を作成することを特徴とする生体信号処理方法。Create an estimated waveform based on the linear combination model based on the heartbeat interval and the contraction time of the heart to approximate the original waveform of the biological signal ,
Subtract the estimated waveform from the original waveform, pass the subtraction result through a bandpass filter to extract the transient response component,
A biological signal processing method , wherein an output waveform is created by adding the transient response component to the estimated waveform . 前記出力波形から心周期の基点を設定し、該心周期の基点から心周期パラメータを算出することを特徴とする請求項１記載の生体信号処理方法。The biological signal processing method according to claim 1 , wherein a base point of a cardiac cycle is set from the output waveform, and a cardiac cycle parameter is calculated from the base point of the cardiac cycle .

Images