(57)【要約】【課題】心電図波形にノイズが加わっている場合でも、適切に呼吸情報を抽出する。【解決手段】呼吸推定装置は、心電図波形のサンプリングデータ列からサンプリングデータの時間差分値を算出する時間差分値算出部３と、時間差分値の最大値と最小値との差である差分ＲＳ振幅を検出する差分ＲＳ振幅検出部４と、差分ＲＳ振幅の時系列データに対して再サンプリング処理を行う再サンプリング処理部５と、再サンプリング処理後の時系列データを周波数解析して差分ＲＳ振幅の周波数スペクトルを求める周波数解析部７と、差分ＲＳ振幅の周波数スペクトルから被験者の呼吸周波数を抽出する呼吸周波数推定部８とを備えている。【選択図】  図１
2018011753.tif 000002
【特許請求の範囲】
【請求項１】
生体の心電図波形のサンプリングデータ列からサンプリングデータの時間差分値をサンプリング時刻ごとに算出する時間差分値算出ステップと、
前記時間差分値の最大値と最小値との差である振幅を検出する振幅検出ステップと、
この振幅検出ステップで検出した振幅の時系列データに基づいて生体の呼吸情報を抽出する呼吸情報抽出ステップとを含むことを特徴とする呼吸推定方法。
【請求項２】
請求項１記載の呼吸推定方法において、
前記呼吸情報は、呼吸周波数であることを特徴とする呼吸推定方法。
【請求項３】
請求項２記載の呼吸推定方法において、
前記呼吸情報抽出ステップは、
前記振幅検出ステップで検出した振幅の時系列データを補間して等間隔のデータに再構成する再サンプリング処理を行う再サンプリング処理ステップと、
再サンプリング処理後のデータを周波数解析し、前記振幅の周波数スペクトルを求める周波数解析ステップと、
この周波数解析ステップで得られた周波数スペクトルから生体の呼吸周波数を推定する呼吸周波数推定ステップとを含むことを特徴とする呼吸推定方法。
【請求項４】
請求項２記載の呼吸推定方法において、
前記呼吸情報抽出ステップは、
前記振幅検出ステップで検出した振幅の時系列データを補間して等間隔のデータに再構成する再サンプリング処理を行う再サンプリング処理ステップと、
再サンプリング処理後の時系列データから位相情報を抽出する特徴量抽出ステップと、
この特徴量抽出ステップで得られた位相情報についてノイズを濾した位相情報を推定するカルマンフィルタ処理ステップと、
このカルマンフィルタ処理ステップで得られた推定位相値を周波数に変換した結果を、生体の呼吸周波数とする呼吸周波数変換ステップとを含むことを特徴とする呼吸推定方法。
【請求項５】
生体の心電図波形のサンプリングデータ列からサンプリングデータの時間差分値をサンプリング時刻ごとに算出する時間差分値算出手段と、
前記時間差分値の最大値と最小値との差である振幅を検出する振幅検出手段と、
この振幅検出手段で検出された振幅の時系列データに基づいて生体の呼吸情報を抽出する呼吸情報抽出手段とを備えることを特徴とする呼吸推定装置。
【請求項６】
請求項５記載の呼吸推定装置において、
前記呼吸情報は、呼吸周波数であることを特徴とする呼吸推定装置。
【請求項７】
請求項６記載の呼吸推定装置において、
前記呼吸情報抽出手段は、
前記振幅検出手段で検出された振幅の時系列データを補間して等間隔のデータに再構成する再サンプリング処理を行う再サンプリング処理手段と、
再サンプリング処理後のデータを周波数解析し、前記振幅の周波数スペクトルを求める周波数解析手段と、
この周波数解析手段で得られた周波数スペクトルから生体の呼吸周波数を推定する呼吸周波数推定手段とを備えることを特徴とする呼吸推定装置。
【請求項８】
請求項６記載の呼吸推定装置において、
前記呼吸情報抽出手段は、
前記振幅検出手段で検出された振幅の時系列データを補間して等間隔のデータに再構成する再サンプリング処理を行う再サンプリング処理手段と、
再サンプリング処理後の時系列データから位相情報を抽出する特徴量抽出手段と、
この特徴量抽出手段で得られた位相情報についてノイズを濾した位相情報を推定するカルマンフィルタと、
このカルマンフィルタで得られた推定位相値を周波数に変換した結果を、生体の呼吸周波数とする呼吸周波数変換手段とを備えることを特徴とする呼吸推定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、心電図波形から呼吸に関する情報を抽出する際に、心電図波形に混入したノイズ等の影響を除去して抽出精度を向上させる呼吸推定方法および装置に関するものである。
【背景技術】
【０００２】
心電図波形の振幅には、呼吸の影響が現れる。心電図波形に呼吸の影響が現れるのは、呼吸に伴って胸郭が伸縮したり、肺（空気を含む）の組成が変化したりすることで、心電図測定系から見たインピーダンスが変動することによると考えられる。この現象を利用して、心電図波形のＲ波のピーク値からＳ波のピーク値までの振幅であるＲＳ振幅の時系列データから、呼吸情報を抽出することができる。ＲＳ振幅は、抽出したＲ波に対して、前後一定の区間の心電図波形を探索し、その区間の最大値と最小値の差をとることで得られる。
【０００３】
しかし、心電図波形を計測する際、波形にノイズが加わることがある。特に、携帯型の装置や、人体に装着するウェアラブルデバイスを用いて日常生活の中での心電図波形を取得する場合には、体動などによるノイズが入りやすい。
【０００４】
特許文献１には、心電図波形のＴ波のピーク値を基に呼吸推定を行う方法において、Ｔ波が発生しない場合もしくはＴ波のピーク値が非常に小さい場合を所定閾値により検出し、所定閾値より小さいピークの値を校正する構成が開示されている。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特許第５６３２５７０号公報
【発明の概要】
【発明が解決しようとする課題】
【０００６】
特許文献１に開示された技術では、心電図波形にノイズが重畳した場合、ＲＳ振幅などの生体信号を補正することは不可能であった。このため、心電図波形にノイズが加わっている状況では、ＲＳ振幅から呼吸情報を抽出することは困難であった。
【０００７】
本発明は、心電図波形にノイズが加わっている場合でも、適切に呼吸情報を抽出することができる呼吸推定方法および装置を提供することを目的とする。
【課題を解決するための手段】
【０００８】
本発明の呼吸推定方法は、生体の心電図波形のサンプリングデータ列からサンプリングデータの時間差分値をサンプリング時刻ごとに算出する時間差分値算出ステップと、前記時間差分値の最大値と最小値との差である振幅を検出する振幅検出ステップと、この振幅検出ステップで検出した振幅の時系列データに基づいて生体の呼吸情報を抽出する呼吸情報抽出ステップとを含むことを特徴とするものである。
また、本発明の呼吸推定方法の１構成例において、前記呼吸情報は呼吸周波数である。
【０００９】
また、本発明の呼吸推定方法の１構成例において、前記呼吸情報抽出ステップは、前記振幅検出ステップで検出した振幅の時系列データを補間して等間隔のデータに再構成する再サンプリング処理を行う再サンプリング処理ステップと、再サンプリング処理後のデータを周波数解析し、前記振幅の周波数スペクトルを求める周波数解析ステップと、この周波数解析ステップで得られた周波数スペクトルから生体の呼吸周波数を推定する呼吸周波数推定ステップとを含むことを特徴とするものである。
また、本発明の呼吸推定方法の１構成例において、前記呼吸情報抽出ステップは、前記振幅検出ステップで検出した振幅の時系列データを補間して等間隔のデータに再構成する再サンプリング処理を行う再サンプリング処理ステップと、再サンプリング処理後の時系列データから位相情報を抽出する特徴量抽出ステップと、この特徴量抽出ステップで得られた位相情報についてノイズを濾した位相情報を推定するカルマンフィルタ処理ステップと、このカルマンフィルタ処理ステップで得られた推定位相値を周波数に変換した結果を、生体の呼吸周波数とする呼吸周波数変換ステップとを含むことを特徴とするものである。
【００１０】
また、本発明の呼吸推定装置は、生体の心電図波形のサンプリングデータ列からサンプリングデータの時間差分値をサンプリング時刻ごとに算出する時間差分値算出手段と、前記時間差分値の最大値と最小値との差である振幅を検出する振幅検出手段と、この振幅検出手段で検出された振幅の時系列データに基づいて生体の呼吸情報を抽出する呼吸情報抽出手段とを備えることを特徴とするものである。
【発明の効果】
【００１１】
本発明によれば、生体の心電図波形のサンプリングデータ列から時間差分値を算出し、時間差分値の最大値と最小値との差である振幅を検出し、この振幅の時系列データに基づいて生体の呼吸情報を抽出することにより、生体の体動などにより心電図波形に重畳したノイズの影響を相殺することができ、生体の呼吸情報の適切な抽出を実現することができる。
【図面の簡単な説明】
【００１２】
【図１】本発明の第１の実施の形態に係る呼吸推定装置の構成を示すブロック図である。
【図２】本発明の第１の実施の形態に係る呼吸推定装置の動作を説明するフローチャートである。
【図３】心電図波形の１例を示す図である。
【図４】図３の心電図波形の時間差分値を示す図である。
【図５】図３の心電図波形から求めたＲＳ振幅と図４の時間差分値波形から求めた差分ＲＳ振幅とを示す図である。
【図６】心電図波形の別の例を示す図である。
【図７】図６の心電図波形の時間差分値を示す図である。
【図８】図６の心電図波形と同時に測定した呼吸計の出力波形を示す図である。
【図９】図６の心電図波形から求めたＲＳ振幅と図７の時間差分値波形から求めた差分ＲＳ振幅とを示す図である。
【図１０】心電図波形の別の例を示す図である。
【図１１】図１０の心電図波形の時間差分値を示す図である。
【図１２】図１０の心電図波形と同時に測定した呼吸計の出力波形を示す図である。
【図１３】図１０の心電図波形から求めたＲＳ振幅と図１１の時間差分値波形から求めた差分ＲＳ振幅とを示す図である。
【図１４】本発明の第２の実施の形態に係る呼吸推定装置の構成を示すブロック図である。
【図１５】本発明の第２の実施の形態に係る呼吸推定装置の動作を説明するフローチャートである。
【図１６】本発明の第２の実施の形態に係る呼吸推定装置の特徴量抽出部の構成を示すブロック図、および特徴量抽出部の各部の信号波形を示す図である。
【図１７】本発明の第２の実施の形態に係る呼吸推定装置のカルマンフィルタのブロック線図である。
【発明を実施するための形態】
【００１３】
［第１の実施の形態］
以下、本発明の実施の形態について図面を参照して説明する。図１は本発明の第１の実施の形態に係る呼吸推定装置の構成を示すブロック図である。呼吸推定装置は、心電図波形のサンプリングデータ列を出力する心電計１と、心電図波形のサンプリングデータ列を記憶する記憶部２と、心電図波形のサンプリングデータ列からサンプリングデータの時間差分値をサンプリング時刻ごとに算出する時間差分値算出部３と、時間差分値算出部３が算出した時間差分値の最大値と最小値との差である差分ＲＳ振幅を検出する差分ＲＳ振幅検出部４（振幅検出手段）と、差分ＲＳ振幅の時系列データに対して再サンプリング処理を行う再サンプリング処理部５と、再サンプリング処理後の時系列データを帯域制限するバンドパスフィルタ６と、再サンプリング処理後の時系列データを周波数解析して差分ＲＳ振幅の周波数スペクトルを求める周波数解析部７と、差分ＲＳ振幅の周波数スペクトルから被験者の呼吸周波数を抽出する呼吸周波数推定部８とを備えている。再サンプリング処理部５とバンドパスフィルタ６と周波数解析部７と呼吸周波数推定部８とは、呼吸情報抽出手段を構成している。
【００１４】
次に、図２を参照して本実施の形態の呼吸推定装置の動作を説明する。心電計１は、図示しない被験者の心電図波形を測定し、心電図波形のサンプリングデータ列Ｘ（ｉ）を出力する。心電図波形の具体的な測定方法は周知の技術であるので、詳細な説明は省略する。記憶部２は、心電計１から出力された心電図波形のサンプリングデータ列Ｘ（ｉ）を記憶する。
【００１５】
時間差分値算出部３は、心電図波形のサンプリングデータＸ（ｉ）の時間差分値Ｙ（ｉ）を算出するため、サンプリングデータＸ（ｉ）の１サンプリング後のデータＸ（ｉ＋１）と１サンプリング前のデータＸ（ｉ−１）とを記憶部２から取得して、サンプリングデータＸ（ｉ）の時間差分値Ｙ（ｉ）を次式のようにサンプリング時刻ごとに算出する（図２ステップＳ１）。
Ｙ（ｉ）＝Ｘ（ｉ＋１）−Ｘ（ｉ−１）                            ・・・（１）
【００１６】
次に、差分ＲＳ振幅検出部４は、時間差分値算出部３が算出した時間差分値の最大値と最小値との差である差分ＲＳ振幅を検出する（図２ステップＳ２）。時間差分値の最大値は、心電図波形のＲ波の急峻な立ち上がりに由来するものであり、時間差分値の最小値は、Ｒ波からＳ波への急峻な心電位の低下に由来するものである。差分ＲＳ振幅検出部４は、時間差分値の最大値とこれに続く最小値の組ごとに振幅検出を行う。
【００１７】
続いて、再サンプリング処理部５は、差分ＲＳ振幅検出部４が算出した差分ＲＳ振幅の時系列データを補間して等間隔のデータに再構成する再サンプリング処理を行う（図２ステップＳ３）。このときの補間方法としては、線形補間、スプライン補間などがある。
【００１８】
バンドパスフィルタ６は、再サンプリング処理後のデータの帯域制限を行う（図２ステップＳ４）。バンドパスフィルタ６を用いる理由は、人の呼吸周波数が低周波のみに限られるためである。このバンドパスフィルタ６の通過帯域は、例えば０．１５〜０．４Ｈｚである。
【００１９】
次に、周波数解析部７は、バンドパスフィルタ６によって帯域制限された再サンプリング処理後のデータにハニング窓関数を掛け合わせた後に、再サンプリング処理後のデータを高速フーリエ変換して、差分ＲＳ振幅の周波数スペクトルを求める（図２ステップＳ５）。周知のとおり、ハニング窓関数は、所望の区間のデータを切り出すために使用される。周波数解析部７は、以上のような周波数解析を、再サンプリング処理後のデータ周期ごとに行う。
【００２０】
呼吸周波数推定部８は、周波数解析部７で得られた差分ＲＳ振幅の周波数スペクトルのピークの周波数を被験者の呼吸周波数として検出する（図２ステップＳ６）。このとき、周波数スペクトルに複数のピークが現れた場合には、強度が最大のピークの周波数を呼吸周波数とする。呼吸周波数推定部８は、以上のような呼吸周波数推定を、再サンプリング処理後のデータ周期ごとに行う。こうして、被験者の呼吸周波数の時系列データを得ることができる。
【００２１】
図３は、被験者が緩徐な呼吸を行っている際の心電図波形の例を示す図である。図３の１００は心電図波形、○印１０１はＲ−Ｒ間隔を示している。Ｒ−Ｒ間隔は、自律神経の働きにより、呼吸に伴って変動する。図３に示す心電図波形の振幅にも、Ｒ−Ｒ間隔と同様の変動、すなわち呼吸に同期した変動がみてとれる。
【００２２】
図４は、図３の心電図波形の時間差分をとった波形を示す図である。○印１０１はＲ−Ｒ間隔を示し、１０２は時間差分値を示している。図３と同様に、時間差分値の波形にも呼吸に伴う振幅の変動が現れている。
【００２３】
図５は、図３の心電図波形から求めたＲＳ振幅と図４の時間差分値から求めた差分ＲＳ振幅とを示す図である。図５の１０３はＲＳ振幅を示し、１０４は差分ＲＳ振幅を示している。ＲＳ振幅と差分ＲＳ振幅の様子はよく一致している。このように、心電図波形の時間差分値から得られる差分ＲＳ振幅は、心電図波形そのものから得られるＲＳ振幅と同様に、被験者の呼吸状態を反映して変動する。したがって、差分ＲＳ振幅から被験者の呼吸周波数を推定可能なことが分かる。
【００２４】
図６は、心電図波形の別の例を示す図である。図６の例では、大きなノイズの重畳はないものの、基線がやや搖動している様子が見うけられる。
図７は、図６の心電図波形の時間差分をとった波形を示す図である。図７によると、心電図波形で観測された基線の搖動がキャンセルされていることが分かる。
【００２５】
図８は、図６の心電図波形と同時に測定した呼吸計の出力波形を示す図である。呼吸計の出力波形は、被験者にマスクを装着させた状態で、呼吸気の流路に設けたサーミスタの抵抗値（任意単位）を示しており、被験者の呼吸リズムを示す標準的な情報である。図８によると、被験者が、一定のリズムで呼吸していることが分かる。
【００２６】
図９は、図６の心電図波形から求めたＲＳ振幅と図７の時間差分値から求めた差分ＲＳ振幅とを示す図である。図５と同様に、１０３はＲＳ振幅を示し、１０４は差分ＲＳ振幅を示している。ＲＳ振幅は、心電図波形の基線の搖動の影響により、呼吸に伴う変動が分かりにくくなっている。それに対して、差分ＲＳ振幅は、ＲＳ振幅よりもはっきりとした変動リズムを呈し、この変動リズムは、図８の呼吸計の出力波形が示す呼吸リズムとも一致している。したがって、差分ＲＳ振幅から被験者の呼吸周波数を推定すれば、適切な推定が可能なことが分かる。
【００２７】
図１０は、心電図波形のさらに別の例を示す図である。この図１０の例の場合、上向きのＲ波のピークが小さく、下向きのＳ波のピークが大きいという特徴がある。また、一部に急峻な基線の搖動が生じている。
図１１は、図１０の心電図波形の時間差分をとった波形を示す図である。図１０によると、心電図波形で観測された基線の搖動はおおむねキャンセルされている。
【００２８】
図１２は、図１０の心電図波形と同時に測定した呼吸計の出力波形を示す図である。図１２によると、被験者の呼吸そのものは安定している様子がうかがえる。
図１３は、図１０の心電図波形から求めたＲＳ振幅と図１１の時間差分値から求めた差分ＲＳ振幅とを示す図である。図５と同様に、１０３はＲＳ振幅を示し、１０４は差分ＲＳ振幅を示している。この図１３の例においても、ＲＳ振幅は、心電図波形の基線搖動の影響を受けて乱れている。一方、差分ＲＳ振幅は、被験者の呼吸に近い変動を示している。
【００２９】
以上のように、本実施の形態によれば、心電図波形の時間差分値を算出し、差分ＲＳ振幅を検出して、この差分ＲＳ振幅から被験者の呼吸周波数を推定することにより、被験者の体動などにより心電図波形に重畳したノイズの影響を相殺することができ、呼吸周波数の適切な推定を実現することができる。
【００３０】
なお、本実施の形態で説明した記憶部２と時間差分値算出部３と差分ＲＳ振幅検出部４と再サンプリング処理部５とバンドパスフィルタ６と周波数解析部７と呼吸周波数推定部８とは、ＣＰＵ（Central Processing Unit）、記憶装置及びインタフェースを備えたコンピュータと、これらのハードウェア資源を制御するプログラムによって実現することができる。ＣＰＵは、記憶装置に格納されたプログラムに従って本実施の形態で説明した処理を実行する。
【００３１】
［第２の実施の形態］
次に、本発明の第２の実施の形態について説明する。図１４は本発明の第２の実施の形態に係る呼吸推定装置の構成を示すブロック図であり、図１と同一の構成には同一の符号を付してある。本実施の形態の呼吸推定装置は、心電計１と、記憶部２と、時間差分値算出部３と、差分ＲＳ振幅検出部４と、再サンプリング処理部５と、バンドパスフィルタ６と、バンドパスフィルタ６で得られた差分ＲＳ振幅の時系列データから位相情報を抽出する特徴量抽出部９と、特徴量抽出部９で得られた位相情報についてノイズを濾した位相情報を推定するカルマンフィルタ１０と、カルマンフィルタ１０で得られた推定位相値を周波数に変換した結果を、被験者の呼吸周波数とする呼吸周波数変換部１１とを備えている。再サンプリング処理部５とバンドパスフィルタ６と特徴量抽出部９とカルマンフィルタ１０と呼吸周波数変換部１１とは、呼吸情報抽出手段を構成している。
【００３２】
図１５は本実施の形態の呼吸推定装置の動作を説明するフローチャートである。心電計１と記憶部２と時間差分値算出部３と差分ＲＳ振幅検出部４と再サンプリング処理部５とバンドパスフィルタ６の動作（ステップＳ１〜Ｓ４）は、第１の実施の形態で説明したとおりである。
【００３３】
特徴量抽出部９は、バンドパスフィルタ６によって帯域制限された再サンプリング処理後の時系列データから瞬時的な線形化された位相情報を抽出する（図１５ステップＳ７）。本実施の形態では、人の呼吸の動きを波形で表したときに、呼吸曲線の理想モデルにおいては、短時間での呼吸数の変化が無視できるほど小さいことと呼吸の動きが正弦波になると仮定し、ヒルベルト変換により瞬時的な位相情報を抽出する。
【００３４】
図１６（Ａ）は特徴量抽出部９の構成を示すブロック図である。特徴量抽出部９は、ヒルベルト変換部９０と、角度算出部９１と、アンラップ処理部９２とから構成される。
【００３５】
まず、ヒルベルト変換部９０は、バンドパスフィルタ６によって帯域制限された再サンプリング処理後の時系列データ（差分ＲＳ振幅の時系列データ）をヒルベルト変換して、位相がπ／２異なる２つの信号（実部及び虚部）を生成する。ヒルベルト変換部９０に入力される信号をＡｅｘｐ（−ｉθ）の正弦波で表現すると（図１６（Ｂ））、生成された実部の信号はＡｃｏｓθで表現される正弦波であり、虚部の信号はｉＡｓｉｎθで表現される正弦波の信号である。Ａは振幅、θは角度、ｉは虚数単位である。
【００３６】
続いて、角度算出部９１は、ヒルベルト変換部９０の出力信号の実部Ａｃｏｓθと虚部ｉＡｓｉｎθとから角度θ（−πから＋π）を算出する。
【００３７】
最後に、アンラップ処理部９２は、角度算出部９１が算出した角度θを連続した位相値に線形化する位相アンラップを行う。角度算出部９１が算出する角度θは−πから＋πの間の値となるので、図１６（Ｃ）に示すように隣り合う点に２πの位相飛びが生じる場合がある。そこで、アンラップ処理部９２は、例えば２πを足したり引いたりすることで、位相を繋ぎ合わせるようにする。これにより、図１６（Ｄ）に示すような連続した位相が得られる。
【００３８】
次に、カルマンフィルタ１０は、差分ＲＳ振幅の位相情報についてノイズを濾した位相情報を推定する（図１５ステップＳ８）。図１７はカルマンフィルタ１０のブロック線図である。カルマンフィルタは、ベイジアン推定と自動再帰推定に基づくものである。システムへの入力はガウシアン雑音を有するという前提の下、カルマン利得Ｋは最小自乗近似により求める。生体システムを記述する物理量ｘ（ｋ）は再帰的に決定されており、次式で表される。
【００３９】
【数１】
2018011753.tif 000003
【００４０】
ここで、ｕ（ｋ）はシステム入力、ｘ（ｋ）は雑音無しの物理量である。また、ｍは測定数、ｎは生体システムからの信号数、Ａはシステムモデルを示すｎ×ｍ行列、Ｈは測定系モデルを示すｍ×ｎ行列である。式（２）に示すように、物理量ｘ（ｋ＋１）には、前時刻における物理量Ａｘ（ｋ）と、ｗs（ｋ）の生体システム雑音とが含まれる。式（３）に示すように、システム入力ｕ（ｋ）には、測定値Ｈｘ（ｋ）と、ｗm（ｋ）の測定系雑音とが含まれる。本実施の形態の場合、測定値Ｈｘ（ｋ）は特徴量抽出部９から入力される位相値である。
【００４１】
【数２】
2018011753.tif 000004
【００４２】
ここで、Ｋ（ｋ）はカルマン利得を示すｎ×ｍ行列、ハットｘ（ｋ）は物理量ｘ（ｋ）の推定値（本実施の形態では位相の推定値）である。以下、同様に文字上に付した「∧」をハットと呼ぶ。カルマン利得Ｋ（ｋ）は次の式（５）〜式（７）により求めることができる。
【００４３】
【数３】
2018011753.tif 000005
【００４４】
ここで、Ｒはセンサ雑音に関する共分散マトリクス、Ｑは生体システム雑音に関する共分散マトリクス、Ｐは推定値に関する共分散マトリクスである。ＨT，ＡTはそれぞれ行列Ｈ，Ａの転置行列である。カルマン利得Ｋ（ｋ）は、測定系雑音ｗm（ｋ）が最小となるように（すなわち、推定誤差が最小となるように）再帰的に決定される。共分散マトリクスＰ（ｋ｜ｋ）の対角成分はフィルタ処理により最小化される自乗推定誤差である。
【００４５】
センサ雑音に関する共分散マトリクスＲを予め正確に設定するために、被験者の呼吸レートが安定している時刻域を利用し、特徴量抽出部９で得られた差分ＲＳ振幅の位相情報について、当該時刻域における標準偏差σを計算する。差分ＲＳ振幅の位相情報について得られた標準偏差σの平均値を、カルマンフィルタ１０の共分散マトリクスＲの対角成分として予め設定する。センサ雑音に関する共分散マトリクスＲは、測定環境に依存するものの、個人差による違いはさほど大きくないと考えられる。
【００４６】
一方、生体システム雑音に関する共分散マトリクスＱは、生体システム雑音であり、個人差が反映される。この共分散マトリクスＱは対角行列である。共分散マトリクスＱの対角成分を変更可能なパラメータとして、被験者毎に共分散マトリクスＱの対角成分の最適値を数値的にテストし、被験者に応じた対角成分の値を予め設定しておくようにすればよい。以上のように、予め被験者に呼吸をさせつつデータを計測し、計測したデータに基づいて共分散マトリクスＲ，Ｑを決定しておくことで、位相の推定精度を高めることができる。ここで、各パラメータの例示を行う。ｘ（ｋ）は位相θｋとドットθｋのベクトルとして表現できる。ドットθｋは位相θｋの微分である。ここでは、文字上に付した「・」をドットと呼ぶ。共分散マトリクスＱ、Ｒの定数ｑ１およびｒ１は設定可能なパラメータであり、上述の通りに設定を行う。例えば、ｑ１＝１×１０-3、ｒ１＝２．４である。
【００４７】
【数４】
2018011753.tif 000006
【００４８】
カルマンフィルタ１０は、フィルタ処理を再サンプリング処理後のデータ周期ごとに行う。
次に、呼吸周波数変換部１１は、カルマンフィルタ１０によって処理された位相値を周波数に変換して呼吸周波数の時系列データを出力する（図１５ステップＳ９）。カルマンフィルタ１０から出力される位相値を時間微分すれば、瞬時角周波数が得られるので、この瞬時角周波数を２πで割ることにより、呼吸周波数を求めることができる。
こうして、本実施の形態では、被験者の呼吸周波数の時系列データを生成して出力することができ、第１の実施の形態と同様の効果を得ることができる。
【００４９】
本実施の形態で説明した記憶部２と時間差分値算出部３と差分ＲＳ振幅検出部４と再サンプリング処理部５とバンドパスフィルタ６と特徴量抽出部９とカルマンフィルタ１０と呼吸周波数変換部１１とは、ＣＰＵ、記憶装置及びインタフェースを備えたコンピュータと、これらのハードウェア資源を制御するプログラムによって実現することができる。ＣＰＵは、記憶装置に格納されたプログラムに従って本実施の形態で説明した処理を実行する。
【産業上の利用可能性】
【００５０】
本発明は、人の呼吸を観察する呼吸モニタリングに適用することができる。
【符号の説明】
【００５１】
１…心電計、２…記憶部、３…時間差分値算出部、４…差分ＲＳ振幅検出部、５…再サンプリング処理部、６…バンドパスフィルタ、７…周波数解析部、８…呼吸周波数推定部、９…特徴量抽出部、１０…カルマンフィルタ、１１…呼吸周波数変換部、９０…ヒルベルト変換部、９１…角度算出部、９２…アンラップ処理部。
【図１】
2018011753.tif 000007
【図２】
2018011753.tif 000008
【図３】
2018011753.tif 000009
【図４】
2018011753.tif 000010
【図５】
2018011753.tif 000011
【図６】
2018011753.tif 000012
【図７】
2018011753.tif 000013
【図８】
2018011753.tif 000014
【図９】
2018011753.tif 000015
【図１０】
2018011753.tif 000016
【図１１】
2018011753.tif 000017
【図１２】
2018011753.tif 000018
【図１３】
2018011753.tif 000019
【図１４】
2018011753.tif 000020
【図１５】
2018011753.tif 000021
【図１６】
2018011753.tif 000022
【図１７】
2018011753.tif 000023