JP6687645B2 - 生体信号処理方法および生体信号処理装置 - Google Patents

Description

本発明は、心電図波形から得られる生体信号に混入したノイズ成分を適切に除去して、生体信号の解析精度を向上させる生体信号処理方法および生体信号処理装置に関するものである。

心臓の拍動リズムは、自律神経すなわち交感神経・迷走神経の影響を受けて変動することが知られている。例えば、安静かつ寛いだ状態では、迷走神経が亢進し、呼吸に伴う心拍変動（呼吸性洞性不整脈）が顕著にみられるようになる。このときの呼吸数は、心電図波形から抽出される心拍（Ｒ波）時刻の時間間隔、すなわちＲ−Ｒ間隔を分析することで、推定することができる。

一方、心電図波形の振幅にも呼吸の影響が現れる。この心電図波形振幅への影響は、呼吸に伴って肺や胸郭が伸縮することで、心電図測定系から見たインピーダンスが変動することによると考えられる。

図１７の（ａ）は、被験者が４−７−８呼吸法（吸気４秒−息止め７秒−呼気８秒）と呼ばれる呼吸方法を実践した際の心電図波形の例を示す図である。図１７の（ａ）によると、心電図波形の高さが変動していることが容易に見て取れる。この心電図波形から、Ｒ波と１つ前のＲ波の時間間隔であるＲ−Ｒ間隔のデータを抽出して時系列にプロットした図を図１７の（ｂ）に示す。Ｒ−Ｒ間隔は、呼気時に長くなり、吸気時に短くなる（心拍は、呼気時に遅くなり、吸気時に速くなる）。図１７の（ｂ）によると、Ｒ−Ｒ間隔が、心電図波形と同期して変動していることが分かる。

また、図１７の（ｃ）は、Ｒ波のピーク値からＳ波のピーク値までの振幅、すなわちＲＳ振幅のデータを心拍毎に抽出して時系列にプロットした図である。ＲＳ振幅は、例えば、抽出したＲ波に対して、前後一定の区間の心電図波形を探索し、その区間の最大値と最小値の差をとる、といった方法で取得できる。図１７の（ｃ）では、上述の吸気−息止め−呼気の時間幅がＲＳ振幅の変動パターンに現れていることが分かる。

図１８の（ａ）は、図１７の（ｂ）に示したＲ−Ｒ間隔のデータを最大エントロピー法（ＭＥＭ：Maximum Entropy Method）で解析したスペクトルを示す図、図１８の（ｂ）は、図１７の（ｃ）に示したＲＳ振幅のデータをＭＥＭで解析したスペクトルを示す図である。ＭＥＭを用いることにより、６０秒という短いデータからでも滑らかなスペクトル曲線が得られる。図１８の（ａ）、図１８の（ｂ）の例では、ＲＳ振幅から得られたスペクトルは、Ｒ−Ｒ間隔から得られたスペクトルよりも明確に、１９秒周期≒０．０５２Ｈｚのピークを示している。このように、Ｒ−Ｒ間隔と併せて、ＲＳ振幅からも呼吸についての情報を得ることができる。

ところで、心電図波形を計測する際、波形にノイズが加わることがある。特に、携帯型の装置やウエアラブルな装置を用いて日常生活の中での心電図波形を取得する場合には、体動などによるノイズが入りやすい。そうしたノイズは、ＲＳ振幅にも混入する可能性がある。また、Ｒ波抽出の誤りを誘発し、その結果、Ｒ−Ｒ間隔などに不適切なデータを入り込ませる要因になる。

特許文献１には、心電図波形のＴ波のピーク値を基に呼吸推定を行う方法において、Ｔ波が発生しない場合もしくはＴ波のピーク値が非常に小さい場合を所定閾値により検出し、所定閾値より小さいピークの値を校正する構成が開示されている。しかしながら、特許文献１に開示された技術では、心電図波形に大きなノイズが重畳した場合、ＲＳ振幅やＲ−Ｒ間隔などの生体信号を補正することは不可能であった。

ＲＳ振幅やＲ−Ｒ間隔などの生体信号を基に、被験者の呼吸等に関する分析を行う場合、それら生体信号の、呼吸等に同期した変動成分に着目することになる。そこで、ＭＥＭや高速フーリエ変換（ＦＦＴ：Fast Fourier Transform）といったスペクトル解析などの手法を生体信号に適用し、呼吸等に関する情報を抽出するといった方法が採られる。ところが、生体信号のデータ列に、ノイズ等に由来する不適切な値が含まれていると、解析の結果が本来のものとは乖離してしまうという問題点があった。

特許第５６３２５７０号公報

本発明は、上記のような点に鑑みてなされたものであり、ＲＳ振幅やＲ−Ｒ間隔などの生体信号に混入したノイズ成分を適切に除去することができ、生体信号の解析の精度を向上させることができる生体信号処理方法および生体信号処理装置を提供することを目的とする。

本発明の生体信号処理方法は、生体の心電図波形から生体信号を抽出する第１のステップと、この第１のステップで抽出した生体信号の時系列データを用いて平均化データを算出する第２のステップと、前記第１のステップで抽出した生体信号のデータが適切かどうかを、このデータよりも過去の時刻の生体信号のデータを用いて算出した前記平均化データに基づいてデータ毎に判定する第３のステップと、この第３のステップで不適切と判定した生体信号のデータの削除および補間のいずれかを行う第４のステップと、前記第４のステップで処理した後の生体信号の１階微分値と２階微分値とを算出する第５のステップと、前記１階微分値と前記２階微分値とが共に０を中心とする所定の範囲内の状態が所定時間以上継続した場合に、生体信号の変動が低いことを警告する第６のステップとを含むことを特徴とするものである。

また、本発明の生体信号処理装置は、生体の心電図波形から生体信号を抽出するように構成された生体信号抽出部と、この生体信号抽出部で抽出された生体信号の時系列データを用いて平均化データを算出するように構成された平均化処理部と、前記生体信号抽出部で抽出された生体信号のデータが適切かどうかを、このデータよりも過去の時刻の生体信号のデータを用いて算出された前記平均化データに基づいてデータ毎に判定するように構成された異常値判定部と、この異常値判定部で不適切と判定された生体信号のデータの削除および補間のいずれかを行うように構成された異常値処理部と、前記異常値処理部によって処理された後の生体信号の１階微分値と２階微分値とを算出するように構成された微分部と、前記１階微分値と前記２階微分値とが共に０を中心とする所定の範囲内の状態が所定時間以上継続した場合に、生体信号の変動が低いことを警告するように構成された通知手段とを備えることを特徴とするものである。

本発明によれば、第１のステップで抽出した生体信号のデータが適切かどうかを、このデータよりも過去の時刻の生体信号のデータを用いて算出した平均化データに基づいてデータ毎に判定し、不適切と判定した生体信号のデータの削除および補間のいずれかを行うことにより、生体信号に混入したノイズ成分を適切に除去することができ、生体信号の解析の精度を向上させることができる。生体信号においては本来その値は変動するものであり、不適切なデータが混入しうる場合にそれを判別するためには、固定的な値ではなく、それまでの値の平均化した値に基づいて処理することが好ましい。

図１は、心電図波形およびＲ−Ｒ間隔の時系列データの例を示す図である。 図２は、ＲＳ振幅の時系列データの例を示す図である。 図３は、本発明の第１の実施例に係る生体信号処理装置の構成を示すブロック図である。 図４は、本発明の第１の実施例に係る生体信号処理装置の動作を説明するフローチャートである。 図５は、ＲＳ振幅の時系列データ、ＲＳ振幅の平均化データ、および補間後のＲＳ振幅のデータの例を示す図である。 図６は、ＲＳ振幅の時系列データ、およびＲＳ振幅の時系列データを再サンプリングしたデータの例を示す図である。 図７は、ＲＳ振幅のスペクトルの例を示す図である。 図８は、本発明の第２の実施例に係る生体信号処理装置の構成を示すブロック図である。 図９は、本発明の第２の実施例に係る生体信号処理装置の動作を説明するフローチャートである。 図１０は、心電図波形、Ｒ−Ｒ間隔の時系列データ、Ｒ−Ｒ間隔の平均化データ、および補間後のＲ−Ｒ間隔のデータの例を示す図である。 図１１は、本発明の第２の実施例によるＲ−Ｒ間隔の平均化データ、および補間後のＲ−Ｒ間隔のデータの例を示す図である。 図１２は、本発明の第３の実施例に係る生体信号処理装置の動作を説明するフローチャートである。 図１３は、ＲＳ振幅の時系列データの１階微分値および２階微分値と、変化量低下判定部の判定結果の例を示す図である。 図１４は、本発明の第４の実施例においてＲＳ振幅を高速フーリエ変換により周波数解析したスペクトルの例を示す図である。 図１５は、本発明の第４の実施例においてＲＳ振幅を最大エントロピー法により周波数解析したスペクトルの例を示す図である。 図１６は、本発明の第１〜第５の実施例に係る生体信号処理装置を実現するコンピュータの構成例を示すブロック図である。 図１７は、心電図波形、Ｒ−Ｒ間隔の時系列データおよびＲＳ振幅の時系列データの例を示す図である。 図１８は、Ｒ−Ｒ間隔のスペクトルおよびＲＳ振幅のスペクトルを示す図である。

［第１の実施例］
図１ないし図７を参照して本発明の第１の実施例について説明する。図１の（ａ）は心電図波形、および心電図波形から抽出したＲ−Ｒ間隔の時系列データの例を示す図であり、図１の（ｂ）は図１の（ａ）の一部を拡大した図である。図１の（ａ）、図１の（ｂ）において、１００は心電図波形（単位［μＶ］）を示し、１０１はＲ−Ｒ間隔（単位［ｍｓ］）を示している。図１の（ａ）、図１の（ｂ）の横軸は時間である。

心電図波形にはノイズが重畳しているものの、Ｒ波に相当する棘波が確認でき、Ｒ−Ｒ間隔は正しく抽出されている。図１の（ａ）を拡大した図１の（ｂ）によると、１０２の部分に見られる心拍では、Ｒ波がノイズのピークと重なっており、Ｓ波はノイズのスロープにかかっていることが分かる。

図２は、図１の（ａ）の心電図波形からＲＳ振幅のデータを心拍毎に抽出して時系列にプロットした図である。ここでは、ＲＳ振幅は、心拍時刻の±２５ｍｓの区間での心電図波形の最大値と最小値の差から算出している。１０３の部分の点はノイズの影響を受けたものであり、不適切な値になっている。

図３は本実施例に係る生体信号処理装置の構成を示すブロック図である。生体信号処理装置は、心電計１が測定した心電図波形から生体信号を抽出する生体信号抽出部２と、生体信号抽出部２で抽出された生体信号の時系列データに対して平均化処理を行う平均化処理部３と、生体信号抽出部２で抽出された生体信号のデータとこのデータよりも前の時刻のデータを用いて算出された平均化データとを比較して、生体信号抽出部２で抽出された生体信号のデータが適切かどうかを判定する異常値判定部４と、異常値判定部４で不適切と判定された生体信号のデータを削除するか、もしくは適切なデータに置き換えて補間する異常値処理部５と、異常値処理部５によって処理された後の生体信号の１階微分値と２階微分値を算出する微分部６と、１階微分値と２階微分値を基に生体信号の変動が低いかどうかを判定する変化量低下判定部７と、異常値処理部５によって処理された生体信号の時系列データをサンプリングするリサンプリング部８と、リサンプリング部８が取得した生体信号の時系列データを周波数解析して生体信号のスペクトルを求める周波数解析部９と、異常値処理部５によって処理された後の生体信号、変化量低下判定部７の判定結果および周波数解析部９の周波数解析結果を表示する表示部１０とを備えている。

次に、図４を参照して本実施例の生体信号処理装置の動作を説明する。本実施例では、図２に示したＲＳ振幅の時系列データから不適切なデータを取り除く例について説明する。なお、微分部６と変化量低下判定部７の動作については他の実施例で説明する。

心電計１は、図示しない被験者（生体）の心電図波形を測定する。心電図波形の具体的な測定方法は周知の技術であるので、詳細な説明は省略する。
生体信号抽出部２は、心電計１が測定した心電図波形から生体信号（本実施例ではＲＳ振幅）を抽出する（図４ステップＳ１）。本実施例の生体信号抽出部２は、上記のとおり、心拍時刻を中心とする所定の区間（ここでは心拍時刻の±２５ｍｓの区間）での心電図波形の最大値（Ｒ波の最大値に相当）と心電図波形の最小値（Ｓ波の最小値に相当）との差をＲＳ振幅として心拍毎に抽出する。

なお、心拍時刻を求める方法としては、例えば特開２０１５−１５６９３６号公報に開示された技術を用いることができる。この文献に開示された技術では、Ｒ波の代表点とこの点の後に存在するＳ波の代表点との間で、所定の閾値を跨ぐ心電図波形の２点のサンプリングデータを検出して、この２点のサンプリングデータを結ぶ直線が前記閾値と交わる時刻を心拍時刻として算出する。また、この心拍時刻を、ＲＳ振幅のデータの時刻とする。

異常値判定部４は、生体信号抽出部２が抽出した生体信号（ＲＳ振幅）のデータと平均化処理部３が直前の時刻までの生体信号（ＲＳ振幅）のデータを用いて算出した平均化データとを比較し、生体信号抽出部２が抽出した生体信号のデータが適切かどうかをデータ毎に判定する（図４ステップＳ２）。

具体的には、異常値判定部４は、ある時刻での生体信号のデータの値Ｘ（ｉ）が、直前の時刻までの生体信号の平均化データＸ´（ｉ−１）を中心とする所定の正常値の範囲内の場合、生体信号のデータＸ（ｉ）は適切であると判定し、データＸ（ｉ）が前記正常値の範囲から外れている場合、データＸ（ｉ）は不適切であると判定する。本実施例では、平均化データＸ´（ｉ−１）の±３０％の範囲を正常値の範囲とする。なお、異常値判定部４は、過去の時刻の生体信号の平均化データＸ´（ｉ−１）を用いて、生体信号のデータＸ（ｉ）が適切かどうかを判定するので、生体信号抽出部２が抽出した１個目の生体信号のデータについては判定処理を実施しない。異常値判定部４が判定処理を実施するのは、生体信号抽出部２が抽出した２個目以降の生体信号のデータについてである。

次に、平均化処理部３は、生体信号抽出部２が抽出した生体信号（ＲＳ振幅）の時系列データに対して平均化処理を行う（図４ステップＳ３）。平均化処理部３は、平均化処理前のｉ番目の生体信号の値をＸ（ｉ）、ｉ番目までの生体信号を平均化した値をＸ´（ｉ）としたとき、次式のようにして生体信号の平均化処理を行う。
Ｘ´（ｉ）＝ｒ×Ｘ（ｉ）＋（１−ｒ）×Ｘ´（ｉ−１） ・・・（１）

式（１）におけるｒは予め定められた係数である。係数ｒの値を小さくするほど、生体信号のデータ列の値の細かな変動が抑えられるものの、一方で生体信号の大まかな変化に追随し難くなるので、このような点を鑑みて、例えばｒ＝０．２などとすると、生体信号の瞬時の変動が抑えられ、適度に平均化された生体信号のデータ列が得られる。

ただし、平均化処理部３は、誤った値が平均化処理に混入することを防ぐために、異常値判定部４が不適切と判定した生体信号（ＲＳ振幅）のデータを平均化処理に取り込まないようにする。例えば生体信号のデータＸ（ｉ）が不適切と判定された場合、データＸ（ｉ）を使用せずに、直前の時刻までの生体信号の平均化データＸ´（ｉ−１）をそのまま平均化データＸ´（ｉ）とすればよい。これにより、平均化データの値の推移をより安定させることができる。

異常値処理部５は、異常値判定部４が不適切と判定した生体信号（ＲＳ振幅）のデータを適切なデータに置き換えて補間する（図４ステップＳ４）。このときの補間の方法としては、不適切と判定された生体信号のデータを、その前後の適切なデータを用いて補間する線形補間がある。また、線形補間に代わってスプライン補間などの補間方法を用いてもよい。

生体信号抽出部２と異常値判定部４と平均化処理部３と異常値処理部５とは、例えば被験者から測定終了の指示があるまで（図４ステップＳ５においてＹＥＳ）、ステップＳ１〜Ｓ４の処理を所定の周期毎（例えば心電計１のサンプリング毎）に行う。

表示部１０は、異常値処理部５によって処理された生体信号（ＲＳ振幅）の時系列データを表示する（図４ステップＳ６）。
図５はＲＳ振幅の時系列データ、このＲＳ振幅の時系列データを平均化処理部３が平均化したデータ、および異常値処理部５が補間したＲＳ振幅のデータの例を示す図である。図５の○印５０は図２と同じＲＳ振幅の時系列データを示し、破線５１はこの時系列データを平均化したデータを示し、＊印５２は補間後のＲＳ振幅のデータを示している。図５によれば、ノイズの影響を受けた不適切なＲＳ振幅のデータを適切に補間できていることが分かる。

なお、平均化処理部３が平均化処理の対象とするのは、生体信号抽出部２が抽出した生体信号のデータであって、異常値処理部５が補間したデータではなく、補間したデータが、その後の平均化処理に組み込まれることはない。補間後のデータを平均化処理に使用しない理由は、補間後のデータは、平均化されたデータに基づいて推定されたものであり、その補間後のデータを、このデータ自身を導くのに用いた値に再び組み込むことは適切ではないからである。

次に、リサンプリング部８は、異常値処理部５によって処理された生体信号（ＲＳ振幅）の時系列データを、心電計１のサンプリング周波数よりも遅いサンプリング周波数（例えば１秒間隔）でサンプリングする（図４ステップＳ７）。
図６はＲＳ振幅の時系列データ、およびＲＳ振幅の時系列データを再サンプリングしたデータの例を示す図である。図６の○印６０はＲＳ振幅の時系列データを示し、＊印６１は異常値処理部５による補間後にリサンプリング部８が再サンプリングしたデータを示している。

周波数解析部９は、リサンプリング部８が取得した生体信号（ＲＳ振幅）の時系列データを高速フーリエ変換または最大エントロピー法（ＭＥＭ）により周波数解析し、生体信号のスペクトルを求める（図４ステップＳ８）。
表示部１０は、周波数解析部９による周波数解析結果のスペクトルを表示する（図４ステップＳ９）。

図７はＭＥＭにより周波数解析したスペクトルの例を示す図である。図７の実線７０はリサンプリング部８が取得したＲＳ振幅の時系列データ（図６の６１）をＭＥＭにより周波数解析したスペクトルを示している。ここでは、比較のため、補間処理前のＲＳ振幅の時系列データ（図６の６０）をＭＥＭにより周波数解析したスペクトルを破線７１で示している。

補間処理前のＲＳ振幅のデータから得られたスペクトルと、不適切なデータを補間した後のＲＳ振幅のデータから得られたスペクトルとでは様相が異なっており、補間処理前のＲＳ振幅のデータから得られたスペクトルには、不適切なデータに起因するスペクトルの成分ないしは分布（図７の７２，７３，７４）が含まれていることが見て取れる。不適切なデータを補間した後のＲＳ振幅のデータから得られたスペクトルの方が、生体の状態を正しく反映した結果を示していると考えられる。

以上のように、本実施例では、ＲＳ振幅などの生体信号に混入したノイズ成分を適切に除去することができ、生体信号の解析の精度を向上させることができる。
なお、本実施例では、異常値判定部４が不適切と判定した生体信号のデータをもっともらしいデータに置き換えて補間しているが、これに限るものではなく、異常値処理部５は、異常値判定部４が不適切と判定した生体信号のデータを削除（データ欠損）してもよい。

［第２の実施例］
次に、本発明の第２の実施例について説明する。図８は本実施例に係る生体信号処理装置の構成を示すブロック図である。本実施例の生体信号処理装置は、生体信号抽出部２と、平均化処理部３ａと、異常値判定部４と、異常値処理部５と、微分部６と、変化量低下判定部７と、リサンプリング部８と、周波数解析部９と、表示部１０とを備えている。本実施例では、Ｒ−Ｒ間隔の時系列データから不適切なデータを取り除く例について説明する。

次に、図９を参照して本実施例の生体信号処理装置の動作を説明する。生体信号抽出部２は、心電計１が測定した心電図波形から生体信号（本実施例ではＲ−Ｒ間隔）を抽出する（図９ステップＳ１０）。Ｒ−Ｒ間隔を求める方法としては、例えば文献「“ECG Implementation on the TMS320C5515 DSP Medical Development Kit (MDK) with the ADS1298 ECG-FE”，Texas Instruments Incorporated，＜http://www.ti.com/lit/an/sprabj1/sprabj1.pdf＞，2011」に開示された技術を用いることができる。この文献に開示された技術では、心電図波形を時間差分した値の変化を基にＲ−Ｒ間隔を求めている。あるいは特開２０１５−１５６９３６号公報に開示された技術等を用いて心拍時刻を求め、心拍時刻の間隔をＲ−Ｒ間隔としてもよい。

異常値判定部４は、生体信号抽出部２が抽出した生体信号（Ｒ−Ｒ間隔）のデータと平均化処理部３ａが直前の時刻までの生体信号（Ｒ−Ｒ間隔）のデータを用いて算出した平均化データとを比較し、生体信号抽出部２が抽出した生体信号のデータが適切かどうかをデータ毎に判定する（図９ステップＳ１１）。

本実施例の異常値判定部４は、ある時刻での生体信号のデータの値Ｘ（ｉ）が、直前の時刻までの生体信号の平均化データＸ´（ｉ−１）の所定数倍（本実施例では１．３５倍）の値を超える場合、データＸ（ｉ）は不適切であると判定し、データＸ（ｉ）が平均化データＸ´（ｉ−１）の所定数倍以下の値である場合、データＸ（ｉ）は適切であると判定する。すなわち、本実施例では、平均化データＸ´（ｉ−１）の所定数倍以下の範囲が正常値の範囲となる。第１の実施例で説明したとおり、異常値判定部４が判定処理を実施するのは、生体信号抽出部２が抽出した２個目以降の生体信号のデータについてである。

次に、平均化処理部３ａは、生体信号抽出部２が抽出した生体信号（Ｒ−Ｒ間隔）の時系列データに対して平均化処理を行う（図９ステップＳ１２，Ｓ１３）。本実施例の平均化処理については後述する。

異常値処理部５は、異常値判定部４が不適切と判定した生体信号（Ｒ−Ｒ間隔）のデータを適切なデータに置き換えて補間する（図９ステップＳ１４）。まず、異常値処理部５は、異常値判定部４が生体信号（Ｒ−Ｒ間隔）に不適切なデータが生じていると判定した場合、不適切なデータが生じた時刻とその直前の時刻との間に挿入すべきデータの個数（心拍数に相当）を推定する。具体的には、異常値処理部５は、生体信号（Ｒ−Ｒ間隔）に不適切なデータが生じた時刻ｔ₂とその直前のデータの時刻ｔ₁との時間間隔（ｔ₂―ｔ₁）を、時刻ｔ₁までのＲ−Ｒ間隔の平均化データの値で除することによって、時刻ｔ₂と時刻ｔ₁との間に挿入するデータの個数Ｎを決める。

そして、異常値処理部５は、決定したデータの個数Ｎで時刻ｔ₂と時刻ｔ₁との時間間隔（ｔ₂―ｔ₁）を等分した値を、時刻ｔ₂と時刻ｔ₁との間に挿入するＲ−Ｒ間隔のもっともらしい値として算出する。こうして、異常値処理部５は、Ｒ−Ｒ間隔に不適切なデータが生じた時刻ｔ₂とその直前の時刻ｔ₁との間に、個数ＮだけＲ−Ｒ間隔のもっともらしい値を挿入することによりＲ−Ｒ間隔を補間することができる。

図１０の（ａ）は心電図波形、および心電図波形から抽出したＲ−Ｒ間隔の時系列データの例を示す図である。図１０の（ａ）において、８０は心電図波形（単位［μＶ］）を示し、８１はＲ−Ｒ間隔（単位［ｍｓ］）を示している。この図１０の（ａ）に示した例では、１６時３１分４０秒付近の時刻で、生体の呼吸運動に由来してＲ−Ｒ間隔が大きく揺らいでいる。

図１０の（ｂ）は図１０の（ａ）に示したＲ−Ｒ間隔の時系列データを平均化したデータ、および異常値処理部５が補間したデータの例を示す図である。図１０の（ｂ）の破線８２はＲ−Ｒ間隔の時系列データを第１の実施例で説明した手法で平均化したデータを示し、＊印８３は補間後のデータを示している。

図１０の（ｂ）の例では、Ｒ−Ｒ間隔を平均化した値を求める際に、Ｒ−Ｒ間隔そのものの値を用い平均化している。また、誤った値が平均化に混入することを防ぐために、それまでの平均化した値に比べ一定の割合（この例では±３０％超）離れた値の場合には平均化に取り込まないこととしている。しかしそのため、本来の生体の状態に由来した揺らぎを示しているＲ−Ｒ間隔の値を不適切な値と見なして平均化に取り込まず、過去の値に固執する形となってしまい、それに基づいて無用な補間を続けるといった、誤った結果を招いていることが分かる。

そこで、本実施例の平均化処理部３ａは、生体信号抽出部２が抽出した生体信号（Ｒ−Ｒ間隔）の平均化処理に際して、Ｒ−Ｒ間隔の逆数に基づく値に対して平均化処理を行い、この平均化処理で求めた値の逆数からＲ−Ｒ間隔の平均化データを算出することを、個々のデータ毎に行う。具体的には、Ｒ−Ｒ間隔の逆数に基づく値（以下、処理対象データＣと呼ぶ）として、Ｃ＝６００００／Ｒ−Ｒ間隔＝心拍数を用いる。ｉ番目の処理対象データをＣ（ｉ）、ｉ−１番目までの処理対象データを平均化した値をＣ´（ｉ−１）、所定の係数をｒとすれば、ｉ番目までの処理対象データを平均化した値Ｃ´（ｉ）は、式（１）と同様の次式によって求めることができる。
Ｃ´（ｉ）＝ｒ×Ｃ（ｉ）＋（１−ｒ）×Ｃ´（ｉ−１） ・・・（２）

平均化処理部３ａの逆数平均化処理部３０は、式（２）によりＣ´（ｉ）を算出する（図９ステップＳ１２）。そして、平均化処理部３ａの平均化データ算出部３１は、Ｒ−Ｒ間隔＝６００００／Ｃ´（ｉ）により、Ｒ−Ｒ間隔の平均化データを算出すればよい（図９ステップＳ１３）。
ただし、平均化処理部３ａは、誤った値が平均化処理に混入することを防ぐために、ある時刻での処理対象データＣ（ｉ）が、直前の時刻までの処理対象データを平均化した値Ｃ´（ｉ−１）を中心とする所定の正常値範囲から外れている場合、処理対象データＣ（ｉ）を不適切であると判定して、平均化処理に取り込まないようにする。例えば処理対象データＣ（ｉ）が不適切と判定した場合、直前の時刻までの処理対象データを平均化した値Ｃ´（ｉ−１）をそのままＣ´（ｉ）とすればよい。ここでは、平均化データＣ´（ｉ−１）の±３０％の範囲を正常値範囲とする。

図９のステップＳ１４，Ｓ１５，Ｓ１６，Ｓ１７，Ｓ１８，Ｓ１９の処理は、それぞれ図４のステップＳ４，Ｓ５，Ｓ６，Ｓ７，Ｓ８，Ｓ９と同様であるので、説明は省略する。

図１１は図１０の（ａ）に示したＲ−Ｒ間隔の時系列データを本実施例の手法で平均化したデータ、および異常値処理部５が補間したデータの例を示す図である。図１１の○印および＊印９１はＲ−Ｒ間隔の時系列データおよび補間後のＲ−Ｒ間隔のデータを示し、破線９２はＲ−Ｒ間隔の平均化データを示している。

最も変動が急な部分のＲ−Ｒ間隔の変動域がおよそ６５０ｍｓ→９００ｍｓ（変動量４２％）であるのに対し、Ｒ−Ｒ間隔の逆数に比例する心拍数の変動域はおよそ９２ｂｐｍ→６７ｂｐｍ（変動量２７％）となる。したがって、本実施例の手法によれば、Ｒ−Ｒ間隔が揺らいでいる箇所でも、それらの揺らぎは平均化に取り込まれており、補間処理したＲ−Ｒ間隔のデータが補間処理前のＲ−Ｒ間隔のデータから逸脱することがない。つまり、心拍数を用いることにより、変動のスケールが抑えられるので、平均化処理が安定化され、生体信号のデータ列に対する処理が正しく行なわれている。

生体信号のデータの値が変動する場合、ある数値域においては、データの値そのものよりも、データの値の逆数に基づく値に対して平均化処理を行うことで、平均化データの変動幅を抑えることができ、平均化処理を安定化することができる。本実施例によれば、生体信号の時系列データにおいて、ノイズ等に起因する不適切なデータを取り除いて、もっともらしく修復することができ、それにより生体の状態のより正確な分析に結び付けることができる。

［第３の実施例］
次に、本発明の第３の実施例について説明する。本実施例においても、生体信号処理装置の構成は第１の実施例と同様であるので、図３の符号を用いて説明する。本実施例は、図３の生体信号処理装置の微分部６と変化量低下判定部７の動作を説明するものである。

図１２を参照して本実施例の生体信号処理装置の動作を説明する。図１２のステップＳ１〜Ｓ６の処理は第１の実施例で説明したとおりである。
次に、微分部６は、異常値処理部５によって処理された生体信号（ＲＳ振幅）の１階微分値と２階微分値を算出する（図１２ステップＳ２０）。

ある時刻ｔ_kにおける生体信号（ＲＳ振幅）の補間された値をｆ（ｔ_k）とすると、１階微分値ｆ’（ｔ_k）は式（３）のように表すことができる。
ｆ’（ｔ_k）＝｛ｆ（ｔ_k+1）−ｆ（ｔ_k）｝／（ｔ_k+1−ｔ_k）・・・（３）

また、２階微分値ｆ”（ｔ_k）は式（４）のように表すことができる。
ｆ”（ｔ_k）＝｛ｆ（ｔ_k+1）−２ｆ（ｔ_k）＋ｆ（ｔ_k-1）｝／（ｔ_k+1−ｔ_k）²
・・・（４）

ＲＳ振幅の変化は呼吸運動に起因するため、息を止めているときなど、呼吸運動が停止している際は、１階微分値と２階微分値が共に０近傍の値をとる。そこで、変化量低下判定部７は、微分部６が算出した１階微分値と２階微分値が共に０を中心とする所定の範囲内の状態が所定時間以上継続した場合に、生体信号の変動が低い（呼吸運動が停止している）と判定し、１階微分値と２階微分値のうち少なくとも一方が０を中心とする所定の範囲から外れているか、あるいは１階微分値と２階微分値が共に所定の範囲内の状態の継続時間が所定時間未満の場合に、生体信号の変動が正常（呼吸運動が経時変化している）と判定する（図１２ステップＳ２１）。

表示部１０は、変化量低下判定部７の判定結果を表示する（図１２ステップＳ２２）。これにより、生体信号の変動が低いことを警告することが可能になる。警告の方法としては、通知手段となる表示部１０に警告メッセージを表示させてもよいし、ＬＥＤ等の発光器からなる通知手段を別に設けて点滅させてもよいし、スピーカ等の通知手段を別に設けて音声で警告を発するようにしてもよい。

図１３は、図５に示したＲＳ振幅の時系列データの１階微分値および２階微分値と、変化量低下判定部７の判定結果の例を示す図である。図１３の□印１１０は１階微分値を示し、△印１１１は２階微分値を示し、実線１１２は変化量低下判定部７の判定結果を示している。

この図１３の例では、１階微分値と２階微分値が共に±８０［μＶ］の範囲内の状態が３秒以上継続した場合に、生体信号の変動が低いとして判定結果を“１”とし、１階微分値と２階微分値のうち少なくとも一方が±８０［μＶ］の範囲から外れているか、あるいは１階微分値と２階微分値が共に±８０［μＶ］の範囲内の状態の継続時間が３秒未満の場合に、生体信号の変動が正常として判定結果を“０”としている。
こうして、本実施例では、呼吸運動が経時変化しているのか、停止しているのかを判定することができ、生体の呼吸運動を監視することができる。

なお、本実施例では、図３に示した生体信号処理装置の微分部６と変化量低下判定部７の動作について説明しているが、図８に示した第２の実施例においても同様の動作を実現することができる。

［第４の実施例］
次に、本発明の第４の実施例について説明する。本実施例においても、生体信号処理装置の構成および処理の流れは第１の実施例と同様であるので、図３、図４の符号を用いて説明する。本実施例では、ＲＳ振幅の時系列データを周波数解析する例について説明する。

第１の実施例で説明したとおり、生体信号処理装置の周波数解析部９は、リサンプリング部８が取得した生体信号（ＲＳ振幅）の時系列データを周波数解析して生体信号のスペクトルを求め（図４ステップＳ８）、表示部１０は、このスペクトルを表示する（図４ステップＳ９）。

図１４はＲＳ振幅を高速フーリエ変換により周波数解析したスペクトルの例を示す図である。横軸は周波数、縦軸は最大ピーク値が１となるように規格化されたスペクトルエネルギー密度である。このように、周波数分布を得たことにより、計測されたＲＳ振幅は、どの周波数帯域の成分により構成されているのか把握可能となり、呼吸運動のより特徴的な知見を得ることができる。

図１４の例では、０．２秒おきにサンプリングされたＲＳ振幅の時系列データを３２個すなわち、計６．２秒間のデータからフーリエ変換の計算を行っている。高速フーリエ変換における周波数分解能は、次の式（５）で規定されるため、横軸のプロット間隔はおよそ０．１６Ｈｚである（文献「三上直樹，“はじめて学ぶディジタル・フィルタと高速フーリエ変換”，ＣＱ出版，ｐ．１３５−１３７，２００５年」参照）。
ｄｆ＝１／Ｔ（＝１／６．４＝１／（０．２×３２）≒０．１６） ・・・（５）

式（５）が意味するところは、周波数分解能は計測時間によって一意に定まるということである。周波数分解能を０．１６Ｈｚよりも高くしたい場合には、リサンプリング部８のサンプリングレートが固定されている状況においては、計測時間を長くするより他ない。

図１５は、図１４の場合と同じＲＳ振幅のデータを用いて最大エントロピー法により周波数解析したスペクトルの例を示す図である。最大エントロピー法は高速フーリエ変換とは異なり、周波数分解能は計測時間に無関係であるため、計測時間を長くせずとも高い周波数分解能が得られる。

平時の呼吸の周期は３〜４秒であり、高速フーリエ変換で周波数分解能を高めようとして計測時間を長くすると複数周期分のデータを含むこととなり、得られる周波数分布はそれらの統計量となってしまい、単回の呼吸活動の情報が埋もれてしまう。しかし、最大エントロピー法を用いれば高い分解能を短い計測時間で得られるため、統計量ではない単回の呼吸についての周波数情報を得ることができる。

図１５においては、図１４に比べて周波数分解能が高まったために、呼吸周波数のピークが０．３５Ｈｚにあることが明瞭に観測でき、かつ、図１４では確認できない０．２３Ｈｚのピークも検出できている。

最大エントロピー法による周波数の計算は、以下の方法により行う（文献「南茂夫，“科学計測のための波形データ処理”，ＣＱ出版，ｐ．１７３−１７４，１９８６年」参照）。最大エントロピー法には、Ｂｕｒｇ法とＹｕｌｅ−Ｗａｌｋｅｒ方があるが、ここではＢｕｒｇ法を例にとる。
求めたいスペクトルエネルギー密度Ｓ（ω）は次の式で与えられる。

Δｔはサンプリングレートであり、図１５の例では０．２秒である。ωは角周波数であり、周波数ｆとω＝２πｆの関係をもつ。
スペクトルエネルギー密度Ｓ（ω）を求めるためには、自己回帰モデルの係数ａ_mi、予測誤差の分散Ｐ_m、および用いるモデル次数ｍが分かればよい。モデル次数ｍは任意の整数をとるが、ここでは最大モデル次数を１６以上の値（文献「井上博，“循環器疾患と自律神経機能”，医学書院，ｐ．８５−８６，２０１０年」参照）から２０を選択し、次数ｍは１から２０をとるものとする。

自己回帰モデルの係数ａ_miを求めるためには、まず式（７）〜式（９）からａ_mmを求めておく必要がある。

式（７）におけるＮはＲＳ振幅のデータの個数であり、ここでは３２である。また、式（７）〜式（９）中の係数ｂ_mi，ｂ’_miの初期値は次式のようになる。ｘ_iはＮ個のデータのうちのｉ番目のデータを現す。

求めたａ_mmから式（１３）、式（１４）の漸化関係式を用いて、自己回帰モデルの係数ａ_miと予測誤差の分散Ｐ_mを求める。

式（１４）で用いるＰ₀は、ＲＳ振幅のデータの平均値をｘ_aveとすると式（１５）で求められる。

次に、式（６）に代入すべき係数ａ_miを決定するための統計量Ｑ_mを計算する。

式（１６）を用いて統計量Ｑ_mをＱ₁からＱ₂₀まで計算し、それらのうち最小のＱ_mを与えるｍが、式（６）で使用するモデル次数ｍである。計算はｍ＝１から１ずつ増加させてＱ₂₀まで計算する途中において統計量Ｑ_mの極小値が最初に出現した段階で計算を打ち切り、Ｑ_mが極小値となったときのｍを式（６）のモデル次数ｍとして使用してもよい。仮にＱ₂₀まで計算しても統計量の極小値が現れなかった場合には、次数の候補の最大値（この場合は２０）を使用する。

以上により、自己回帰モデルの係数ａ_mi、予測誤差の分散Ｐ_mおよびモデル次数ｍが入手できたので、式（６）により、周波数分布を求めることができる。図１５では、周波数分解能をΔｆ＝０．００１とした。すなわち、角周波数ωを０から０．００６２８ずつ増加させて、３．１４まで計算した。高速フーリエ変換の際のΔｆ＝０．１６と比べて高い分解能を得ることができた。

［第５の実施例］
第１〜第４の実施例では、判定の対象となる生体信号のデータＸ（ｉ）が、このデータＸ（ｉ）よりも過去の時刻の生体信号のデータを用いて算出された平均化データＸ´（ｉ−１）に基づく所定の正常値の範囲外の場合に、判定の対象となる生体信号のデータＸ（ｉ）を不適切と判定していたが、判定処理は上記の例に限るものではない。例えば、平均化データＸ´（ｉ−１）±α（αは規定値）の範囲を正常値の範囲としてもよい。

また、第１〜第４の実施例の異常値判定部４は、判定の対象となる生体信号のデータＸ（ｉ）とこのデータよりも過去の時刻の生体信号のデータ（直前の時刻までの生体信号のデータ）とを用いて算出した平均化データから得られる分散σ²を算出してもよい。そして、異常値判定部４は、この分散σ²が、過去の時刻の生体信号のデータを用いて算出した平均化データから得られる分散σ_p ²に基づく所定の正常値の範囲外の場合に、生体信号のデータＸ（ｉ）を不適切と判定するようにしてもよい。例えば２σ_p ²以下の範囲を正常値の範囲とし、分散σ²が２σ_p ²を上回る場合に、生体信号のデータＸ（ｉ）を不適切と判定すればよい。

［第６の実施例］
第１〜第５の実施例で説明した生体信号処理装置は、ＣＰＵ（Central Processing Unit）、記憶装置及びインターフェースを備えたコンピュータと、これらのハードウェア資源を制御するプログラムによって実現することができる。このコンピュータの構成例を図１６に示す。コンピュータは、ＣＰＵ４０と、記憶装置４１と、インターフェース装置（以下、Ｉ／Ｆと略する）４２とを備えている。Ｉ／Ｆ４２には、心電計１などが接続される。このようなコンピュータにおいて、本発明の生体信号処理方法を実現させるためのプログラムは、フレキシブルディスク、ＣＤ−ＲＯＭ、ＤＶＤ−ＲＯＭ、メモリカードなどの記録媒体に記録された状態で提供され、記憶装置４１に格納される。ＣＰＵ４０は、記憶装置４１に記憶されたプログラムに従って第１〜第５の実施例で説明した処理を実行する。

本発明は、心電図波形から得られる生体信号を解析する技術に適用することができる。

１…心電計、２…生体信号抽出部、３，３ａ…平均化処理部、４…異常値判定部、５…異常値処理部、６…微分部、７…変化量低下判定部、８…リサンプリング部、９…周波数解析部、１０…表示部、３０…逆数平均化処理部、３１…平均化データ算出部。

Claims (7)

1. 生体の心電図波形から生体信号を抽出する第１のステップと、
   この第１のステップで抽出した生体信号の時系列データを用いて平均化データを算出する第２のステップと、
   前記第１のステップで抽出した生体信号のデータが適切かどうかを、このデータよりも過去の時刻の生体信号のデータを用いて算出した前記平均化データに基づいてデータ毎に判定する第３のステップと、
   この第３のステップで不適切と判定した生体信号のデータの削除および補間のいずれかを行う第４のステップと、
   前記第４のステップで処理した後の生体信号の１階微分値と２階微分値とを算出する第５のステップと、
   前記１階微分値と前記２階微分値とが共に０を中心とする所定の範囲内の状態が所定時間以上継続した場合に、生体信号の変動が低いことを警告する第６のステップとを含むことを特徴とする生体信号処理方法。
2. 請求項１記載の生体信号処理方法において、
   前記第２のステップは、
   前記第１のステップで抽出した生体信号の時系列データのそれぞれの値の逆数に基づく値に対して平均化処理を行う第７のステップと、
   この第７のステップの平均化処理で求めた値の逆数から前記平均化データを算出する第８のステップとを含むことを特徴とする生体信号処理方法。
3. 請求項１または２記載の生体信号処理方法において、
   前記第２のステップは、前記第１のステップで抽出した生体信号のデータの平均化処理に際して、このデータよりも過去の時刻のデータを用いて算出した前記平均化データの値を中心とする所定の範囲から外れている生体信号のデータ、および前記第４のステップで補間した生体信号のデータを、前記平均化データの算出に使用しないことを特徴とする生体信号処理方法。
4. 請求項１乃至３のいずれか１項に記載の生体信号処理方法において、
   前記生体信号は、Ｒ波と１つ前のＲ波の時間間隔であるＲ−Ｒ間隔、およびＲ波のピーク値からＳ波のピーク値までのＲＳ振幅のいずれかであることを特徴とする生体信号処理方法。
5. 請求項１乃至４のいずれか１項に記載の生体信号処理方法において、
   前記第４のステップで処理した後の生体信号を高速フーリエ変換および最大エントロピー法のいずれかにより周波数解析して、生体信号のスペクトルを求める第９のステップをさらに含むことを特徴とする生体信号処理方法。
6. 生体の心電図波形から生体信号を抽出するように構成された生体信号抽出部と、
   この生体信号抽出部で抽出された生体信号の時系列データを用いて平均化データを算出するように構成された平均化処理部と、
   前記生体信号抽出部で抽出された生体信号のデータが適切かどうかを、このデータよりも過去の時刻の生体信号のデータを用いて算出された前記平均化データに基づいてデータ毎に判定するように構成された異常値判定部と、
   この異常値判定部で不適切と判定された生体信号のデータの削除および補間のいずれかを行うように構成された異常値処理部と、
   前記異常値処理部によって処理された後の生体信号の１階微分値と２階微分値とを算出するように構成された微分部と、
   前記１階微分値と前記２階微分値とが共に０を中心とする所定の範囲内の状態が所定時間以上継続した場合に、生体信号の変動が低いことを警告するように構成された通知手段とを備えることを特徴とする生体信号処理装置。
7. 請求項６記載の生体信号処理装置において、
   前記平均化処理部は、
   前記生体信号抽出部で抽出された生体信号の時系列データのそれぞれの値の逆数に基づく値に対して平均化処理を行うように構成された逆数平均化処理部と、
   この逆数平均化処理部の平均化処理で求めた値の逆数から前記平均化データを算出するように構成された平均化データ算出部とから構成されることを特徴とする生体信号処理装置。

Images