(57)【要約】【課題】本発明は、人の拍動データの呼吸性洞性不整脈（ＲＳＡ）成分の強度を評定するための方法及びシステムに向けられている。【解決手段】方法は、光電脈波計（ＰＰＧ）の様なトランスデューサ機器を通して心拍数をデジタル式にサンプリングして、拍動から拍動までの間隔の時系列を取得する段階と、ウィンドウ化されたデータセグメントに自己回帰（ＡＲ）分析を行う段階と、強度をそれぞれのデータウィンドウの現在のＲＳＡ成分と相関させるフィルタ極を識別する段階を含んでいる。方法は、ＲＳＡ成分の堅牢な測定値を、式、2013111467.tif 000018によって計算することを規定している。ここにＲは成分出力測定値、Ｅwは総ウィンドウエネルギー、ｋは比例定数、ＰiはＲＳＡ成分との相関が最も強いフィルタ極である。【選択図】図２
2013111467.tif 000002
【特許請求の範囲】
【請求項１】
拍動センサーデータからＲＳＡを計算するためのコンピュータ実装型の方法において、
（ａ）時間のウィンドウ期間中の被験者の拍動をセンサーによってサンプリングして、拍動センサーデータを取得する段階と、
（ｂ）前記センサーによって、前記ウィンドウ期間中のそれぞれのサンプルに対応する拍動センサーデータを演算機器に送信する段階と、
（ｃ）前記演算機器によって、前記拍動センサーデータから、拍動間間隔を計算する段階と、
（ｄ）前記演算機器によって、前記拍動間間隔データに最良適合する全極型フィルタモデルを求める段階と、
（ｅ）前記演算機器によって、前記フィルタモデル極の位置を計算する段階と、
（ｆ）前記演算機器によって、単位円の事前に規定されているセクタ内で最も大きな絶対値を有する極を識別する段階と、
（ｇ）前記演算機器によって、前記識別された極に基づいてＲＳＡの値を特定する段階と、から成る方法。
【請求項２】
前記特定する段階は、ＳRSAと表される成分強度を、方程式、
2013111467.tif 000016
に従って特定し、ここに、ＰRSAは識別された極、Ｅwは総ウィンドウエネルギー、ｋは比例定数である、請求項１に記載の方法。
【請求項３】
前記センサーは、心電計、光電脈波計、圧力センサー、及びマイクロフォンから成る群より選択されたトランスデューサ機器である、請求項１に記載の方法。
【請求項４】
前記計算された拍動間間隔は、均一レートでサンプリングし直される、請求項１に記載の方法。
【請求項５】
前記演算機器によって、ウィンドウ中の前記サンプルを平均値除去及びトレンド除去する段階を更に含んでいる、請求項１に記載の方法。
【請求項６】
前記最良適合全極型モデルフィルタ係数を計算するのに、ブルグ法が使用されている、請求項１に記載の方法。
【請求項７】
前記単位円の前記事前に規定されているセクタは、
2013111467.tif 000017
によって定義されており、ここにｆsはサンプリング周波数である、請求項１０に記載の方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、心拍数変動（ＨＲＶ）の呼吸性洞性不整脈（ＲＳＡ）成分を測定するための方法及び本方法を利用するシステムに関する。本発明は、更に、心拍数変動のＲＳＡ成分を測定するための方法及びシステムに基づくバイオフィードバック若しくは不安管理のアプリケーションに関する。
【背景技術】
【０００２】
ＲＳＡに基づくバイオフィードバックは、不安やストレスの制御に有用であることが示されており、臨床試験は高血圧の軽減におけるその有効性を示している。瞬間ＲＳＡは、被験者又は使用者への視聴覚的表示を変調して、或る歩調の呼吸運動が当人の心臓のパターンにどの様に影響しているかを指し示すためのバイオフィードバックパラメータとして使用することができる。これの単純な例として、使用者に表示される画面色を変調することや被験者に向けて発せられる発信音の周波数を変調することが挙げられる。
【０００３】
一般的な測定値であるＲＳＡ強度Ｓ′(ｔ)は、以下の方程式（１）、即ち、
2013111467.tif 000003
にある様に、同期した時間ウィンドウに亘って計算された心拍数信号Ｈ(ｔ)と呼吸数Ｂ(ｔ)の関数の相互相関関数の最大値と定義することができ、ここに、*は、同じ時間間隔に亘って２つの信号に適用された相互相関関数を表し、ｍａｘは最大値関数を表している。
【０００４】
呼吸数は心拍数よりも測定するのが簡便でないことから、測度Ｓ′(ｔ)は、Ｂ(ｔ)を振幅１．０の正弦関数として演繹的に定義し、使用者が呼吸するように指導される目標呼吸周波数を大抵０．０５ヘルツから０．２ヘルツの範囲と定義した上で、推定される場合が多い。これによりＳ′(ｔ)を、呼気流又は肺の体積を直接測定することなく推定できる。
【０００５】
心臓の活動を、処理及びＲＳＡ特徴分析に利用できる電気信号に変換するには、幾つかのセンサー技術を利用することができる。これらのセンサーには、マイクロフォン（音声心臓信号）、圧力センサー（脈圧）、心電計（ＥＣＧ）、光電脈波計（ＰＰＧ）、並びにＲＦ技術を利用している非接触センサーが含まれる。
【０００６】
ＲＳＡを求めるための現行方法は、相互相関又はパワースペクトル分析に基づいており、それらは、大抵、多くのモバイル機器上でリアルタイムで実行させるには演算処理上の要求が多すぎる。
【発明の概要】
【発明が解決しようとする課題】
【０００７】
以上より、市場では、心拍数信号のＲＳＡ成分を抽出する方法において、資源が制約されているモバイル機器を含め、広範な演算機器で使用することのできる方法が必要とされている。
【課題を解決するための手段】
【０００８】
本発明は、心臓コヒーレンスの状態が実現されるように歩調取りされた呼吸をする技法について使用者を案内するＲＳＡバイオフィードバックアプリケーションを開示している。本発明は、更に、目標が満たされたときに点数を割り当てることによって被験者に見返りを与える点数制システムを開示している。
【０００９】
本発明は、ヒト被験者の心拍数のＲＳＡ成分のメトリックを生の拍動間隔のデータストリームから効率的に推定する方法を規定している。
【００１０】
ＲＳＡ強度推定方法は、デジタル心拍数信号を重なり合ったウィンドウ式のデータセグメントに分割し、拍動間隔のそれぞれのウィンドウについて最良適合全極型線形フィルタ（自己回帰モデル）を計算する。複素平面の目標セクタ内の単一極がＲＳＡ極として識別され、ＲＳＡ強度の推定値は、以下の方程式、即ち、
2013111467.tif 000004
に基づいて計算され、ここに、|ＰRSA|は識別されたＲＳＡ極の絶対値、Ｅwは総ウィンドウエネルギー、kは比例定数である。
【００１１】
この方法の主要な利点は、フィルタモデルの度数をスケール化することによるものであって、そのおかげで本方法はその処理効率をスケールすることができる。このことは、複雑な高速フーリエ変換（ＦＦＴ）処理をリアルタイムで行うには不十分な演算能力しか持ち合わせていないモバイル機器の様な演算機器に処理をオフロードしようとする場合には特に有用である。自己回帰方法論は、更に、より短いデータセグメントを使用してより高い分解能のＲＳＡ識別を可能にする。より短いデータウィンドウなら、被験者の生理状態の更新をより高速化できることから、或る特定のバイオフィードバックアプリケーションにとっては好適であろう。
【００１２】
本発明のこれらの実施形態及び他の態様は、以下の詳細な説明及び添付図面から容易に明らかになるものであり、詳細な説明及び添付図面は本発明を例示的に説明しようとするものであって本発明を限定しようとするものではない。
【図面の簡単な説明】
【００１３】
【図１】本発明の或る例示としての実施形態のシステム図である。
【図２】生の心拍数センサーデータからＲＳＡメトリックを抽出するための或る好適な方法のフロー図である。
【図３】良適合全極型線形フィルタの極位置を示すアルガン図を描いており、識別されたＲＳＡ極Ａとこの極の位相角度Ｂが表示されている。
【発明を実施するための形態】
【００１４】
本発明は、次に続く詳細な説明を添付図面と関連付けて読んで頂ければ、より完全に理解されることであろう。ここには本発明の詳細な実施形態が開示されているが、但し、開示されている実施形態は、本発明を単に例示しているものであり、本発明は様々な形態に具現化することができるものと理解頂きたい。従って、ここに開示されている特定の機能関連の詳細事項は、制限を課すものと解釈されてはならず、単に、特許請求の範囲のための根本事項として、また本発明を様々な形で採用することを実質的且つ適切に詳述されている実施形態の中で当業者に教示するための代表的な根本事項として解釈されるべきである。
【００１５】
図１は、本発明の例示としての実施形態１００のシステム図である。センサー１０４は、ヒト被験者１０２の心拍数を測っている。最も好適には、センサー１０４は、生の心拍数センサーデータを通信リンク１０６を介して演算機器１０８に連絡し、演算機器１０８は心拍数データを分析し、ＲＳＡメトリックを抽出する。或る好適な実施形態では、通信リンク１０６は、例えばブルートゥース通信の使用による、近距離通信リンクである。しかしながら、通信リンク１０６は、近距離通信、例えば、インターネットの様な広域ネットワーク、ワイヤレス通信、ローカルエリア通信ネットワークに限定されない。概して、通信リンク１０６は、データ通信をサポートしている如何なる型式の通信リンク又はネットワークによって提供されていてもよい。或る好適な実施形態では、ＲＳＡメトリックは、演算機器１０８で実行されているアプリケーションであって、処理されたセンサーデータに基づいて被験者１０２にバイオフィードバックを提供するバイオメトリックアプリケーションによって使用されている。センサー１０４は、市販のオキシメーターとして描かれているが、本発明はその様に限定されず、他の心拍数センサーが使用されていてもよい。加えて、演算機器１０８はスマートフォンとして描かれているが、とりわけ、携帯電話、ラップトップコンピュータ、パッドコンピュータ、ゲームコンソール、ケーブルＴＶシステムで使用されている様なセットトップボックス、及びパソコンを含め、如何なる型式の演算機器が使用されていてもよい。演算機器１０８は、最低でも、プロセッサと、プログラムコードやデータを記憶するためのデータストレージと、センサー１０４と通信するための通信性能と、バイオメトリックアプリケーションによるグラフィック式の結果を表示するためのディスプレイを保有している。
【００１６】
ＲＳＡ強度を心拍数信号から計算するための方法
図２は、ＲＳＡメトリックを生の心拍数センサーデータから抽出するための或る好適な方法２００のフロー図である。方法２００では、生のセンサーデータフィード２０２が、心拍数測定装置、即ちセンサーの或る実施形態１０４から入って来る。その様な心拍数測定値測定機器の例には、光学的に取得される脈波図である光電脈波図（ＰＰＧ）機器、器官の体積測定機器、及び心電図（ＥＣＧ）機器が含まれる。概して、方法２００は、ヒトの拍動を測定するか又はヒトの拍動の態様に関する電気信号を提供するセンサー又はトランスデューサからの何らかのデータ信号に働きかける。成功している或る試作品では、センサー１０４は、ブルートゥース近距離通信リンクを使用してNokia N97モバイル機器に接続しているNONIN9560オキシメーターとして実装されている。NONINブルートゥースオキシメータのデータストリームは、NOKIA N97側で利用できるJSR82 JAVAブルートゥースソフトウェアモジュールを使用してアクセスされ、１６ビット符号なしデータとして表現されるＰＰＧ心拍数センサーデータがNOKIA N97に毎秒７５サンプルのレートで送信される。
【００１７】
生のサンプリングされた拍動センサーデータストリームは、ステップ２０４でノイズを減らすための第１ローパスフィルタリングによって、ここではＩＢＩデータと呼ばれる拍動間時間間隔の時系列へ変換され、次にステップ２０６で、ピーク検出アルゴリズムが、それぞれの拍動について局所的な最大ピークを割り出し、ピーク間の時間間隔を出力する。入って来る生データは、ローパスフィルタを使用してフィルタリングされ、移動平均フィルタを適用することによってトレンドが取り除かれる。１つの実施形態では、この目的に次数４のチェビシェフフィルタがカットオフ周波数４．０Ｈｚで使用されている。
【００１８】
これにより、ノイズが減ったゼロ平均信号をステップ２０６のピーク検出アルゴリズムに基づく閾の中へ通せるようになる。１階微分がゼロに等しければ、それはピークである可能性を示唆することになるが、その様な１階微分がゼロに等しい全ての点を比較することによってピークが検出され、それぞれの拍動サイクルの正の範囲内の最大値が見つけられる。これは、生理信号に対して堅牢なピーク検出を行うためのデジタル信号処理（ＤＳＰ）の技術分野では知られている方法である。
【００１９】
ステップ２０８で、単位を秒とするＩＢＩデータは、次に、異常値を識別し置き換えるために分析され、ここに、異常値とは、ＩＢＩデータストリーム内の現在の正規変動及びトレンドから外れている値のことである。異常値は、測定ノイズやピーク検出アルゴリズム内の誤差、更には非常に短いか非常に長いかのどちらかの拍動を発生させる異所性拍動の様な生理的原因によっても起こり得る。異所性拍動は、正常の健康な被験者にも起こり、これらの拍動の比率が非常に高ければ心臓疾患が疑われる。異常値除去は、ここでは、それぞれの拍動の期待値を計算することによって実現されている。異常値を処理するための第１パスは、ＩＢＩデータの単純な移動平均フィルタを使用して平均期待値を計算し、次に、以下、即ち、
2013111467.tif 000005
に示されている様に、この期待値を中心とする範囲より外の値を異常値として定義することによって行われている。異常値は、次に、ＩＢＩストリームの中で期待値に置き換えられ、後の処理に向けて誤差のない信号が作成される。
【００２０】
ステップ２１２では、ステップ２１８で求められたＡＲフィルタモデルを使用して、期待ＩＢＩ値が計算される。ステップ２１８は、音声信号処理では線形予測符号化（ＬＰＣ）として知られている方法を使用しているものであって、その方法は、安定及び準安定の時系列における期待値の堅牢な測定値を与える。ＡＲモデル計算がステップ２１８から戻されてくる。nをモデルの次数として、前のn次ＩＢＩ値にＡＲモデルを適用して、次のＩＢＩ値の推定値を計算する。４から１５の間のＡＲモデルフィルタ次数が使用される。ステップ２１８は、ステップ２１８がステップ２１２の後に起こっているので、値、ＡＲモデル値は、前のＩＢＩデータウィンドウに対応することになる。これは、ＩＢＩ信号が準定常であり且つ新しいデータウィンドウが比較的に十分短い間隔で取られているとの仮定に基づくと有効である。この仮定は、ウィンドウの重なりが５秒又はそれ以下である場合に有効である。
【００２１】
異常値が除去されたら、ステップ２１０で、ＩＢＩデータストリームが均一レートでサンプリングし直される。これは、生データのサンプリングが生データストリーム中のそれぞれのピークで行われており、これらのピークは時間的に均一に間隔が空いているわけではないことから、必要になる。線形システム分析を利用するためには、データは均一レートで提供されなくてはならず、ひいてはデータは後の処理が適用できるようになる前にこれを実現するべく処理されなくてはならない。均一サンプリングは、均等に間隔が空けられた点にＩＢＩデータの補完点を取ることによって実現される。非均一ＩＢＩデータに三次スプライン補完が適用され、サンプリングし直し点は０．２５ｓ（４Ｈｚ）間隔で確定される。
【００２２】
ステップ２１４で、均一ＩＢＩデータストリームは、長さｎのデータのセグメントへウィンドウ化され、それぞれのウィンドウは前のウィンドウとｍサンプル分重なり合っている。ｎの好適な値は１２８であり、ｍの好適な値は３２である。ウィンドウ化とは、特定の間隔の外では値がゼロとなる関数を適用することを指す。ウィンドウ化処理が原因のＡＲパラメータの歪みを小さくするために、ウィンドウデータには更にハミングウィンドウ関数が適用されている。
【００２３】
ステップ２１６で、ウィンドウ化されたデータには、次いで、データのＤＣオフセットを除去する平均値除去が行われ、そしてデータウィンドウの長さに亘ってデータから内在する何らかのトレンドを除去するトレンド除去が行われる。ウィンドウ中のデータの平均値除去及びトレンド除去は、ウィンドウ中のそれぞれのデータ点を通る最良適合線を、最小二乗平均アルゴリズムを使用して計算し、この中心線をウィンドウ中のそれぞれのデータ点から差し引くことによって行われる。
【００２４】
ステップ２１８で、データウィンドウは、次に、自己回帰（ＡＲ）法により処理されて、データの最良適合全極型フィルタモデルが計算される。ＡＲ法は、タイムドメインで、以下の方程式、即ち、
2013111467.tif 000006
にある様に、時系列におけるそれまでの値の線形加重和に基づいて次の値の最良推定値を見つけるために使用される線形予測法として定式化することができる。ここに
2013111467.tif 000007
は期待値であり、
2013111467.tif 000008
は予測係数である。これは、以下の方程式、即ち、
2013111467.tif 000009
にある様に、入力信号を出力信号に関係付けるｚ空間の伝達関数と見ることもできる。
【００２５】
方程式５は、全極型フィルタを記述しており、それぞれのフィルタ極は、以下に方程式（６）、即ち、
2013111467.tif 000010
に示されている大きさと位相角度を有するものとして極座標に記述することができる。
【００２６】
極Ｐiの大きさｒが大きいほど、それがシステムの力学全体に及ぼす影響は強くなる。それぞれの極について、位相角度θは、この極がその効果を働かせる周波数に関係付けられる。それぞれの極は、モデル化されるデータストリームの成分にマップさせることができ、この成分の強度と周波数は、それが整合している極の大きさと移相によって追跡することができる。
【００２７】
全極型モデルのパラメータａ1からａpまでを入力データのウィンドウから計算するのに方法が幾つかある。或る好適な実施形態では、最大エントロピー法又はＭＥＭとしても知られているブルグ法が使用されている。ブルグ法は、当技術ではよく知られており、入力データのウィンドウに基づき、次数ｐの全極型フィルタモデルのパラメータを推定するものである。ブルグ法の演算効率は、特に低次数システムについては、他の方法に勝って有利である。ブルグ法は、更に、正確な極推定がより難しい短いデータセグメント（例えば３０秒など）の分析に特に有利である。１つの実施形態では、フィルタ次数は、方法を実行する演算機器の演算能力に適合されている。前述のＮ９７機器と共に使用される１つの実施形態では、フィルタ次数９が成功裏に使用された。
【００２８】
ブルグアルゴリズムの出力は、ユール・ウォーカの多項式として知られている、伝達関数に整合する多項式を記述するパラメータセットａ1，ａ2，…ａpである。ステップ２２０で、伝達関数の極が、この多項式の根を解くことによって計算される。このステップでは、根を解く標準的な如何なるアルゴリズムが使用されてもよい。１つの実施形態では、堅牢な結果を与える効率のよいアルゴリズムであるジェンキンス・トラウブ法が使用されている。
【００２９】
一旦、伝達関数の極が計算されると、ステップ２２２で、拍動データのＲＳＡ成分に最も密接に関係している極が識別される。この極は、毎秒４呼吸から３０呼吸の正常な呼吸数の制約付きのセクタに在る。この範囲には、更に、呼吸が２０呼吸を上回れば交感神経系がＲＳＡの効果を減少させ、その結果極の大きさが下がるという事実の制約が付けられている。これは、０．０７Ｈｚから０．３３Ｈｚの周波数範囲に関係している。（毎分５−７呼吸程度の遅い制御された呼吸数はＲＳＡ成分強度にピークを生じさせることが示されている。）サンプリングレートfsを所与とすると、これは、以下の方程式（７）、即ち、
2013111467.tif 000011
で与えられるセクタ角度の範囲に相当する。このセクタ内の極が、大きさの観点から比較される。それぞれの極Ｐiの影響のメトリックは極エネルギーＥiとして定義することができ、ここに、
2013111467.tif 000012
である。
【００３０】
目標セクタ内の最も大きいＥiを有する極が識別される。より低い位相角度のセクタに、０．７×Ｅi  ｍａｘより大きい極エネルギーを有する極が他にあれば、これらの極の角度の最も低いものがＲＳＡ極として識別される。
【００３１】
平均値除去及びトレンド除去されたＩＢＩデータウィンドウの総ウィンドウエネルギーＥwは、以下の方程式、即ち、
2013111467.tif 000013
によって計算される。
【００３２】
ステップ２２４で、ＲＳＡ成分の総強度が、以下の方程式（１０）、即ち、
2013111467.tif 000014
を使用して計算される。ここに、ＥRSAは、選択されたＲＳＡ極ＰRSAと関連付けられる極エネルギー、Ｅwは、総ウィンドウエネルギー、ｋは、目標範囲への出力を正規化するのに使用される比例定数である。１つの実施形態では、ｋの値は、ＩＢＩデータウィンドウ中のサンプルの数ｎの逆数と定義されており、ｎ＝１２８の場合、ｋ＝０．００７８の値を意味する。以上の方程式（９）と方程式（１０）を組み合わせると、以下の式、即ち、
2013111467.tif 000015
が与えられる。強度ＳRSAは、心拍数データの現在のＲＳＡ成分の堅牢な測定値である。現在の実施形態では、ＲＳＡが強いと、ＳRSA測定値のこの値は１．０を超えた。０．２より下のＳRSAは、ＩＢＩデータウィンドウ中にＲＳＡ成分が殆ど又は全くないことの表れであった。
【００３３】
ステップ２２６で、ＳRSA、即ちＲＳＡ成分の総強度は、演算機器３０８に送信される。他の実施形態では、ＳRSA値は、記憶され、要求があり次第演算機器１０８に送られる。
【００３４】
図３は、最良適合全極型フィルタの極位置を示すアルガン図を描いており、識別されたＲＳＡ極Ａとこの極の位相角度Ｂが表示されている。図３に示されている実施形態は、心拍数データの１つのウィンドウ化されたセグメントについて、９次全極型フィルタモデルをサンプリングレート２Ｈｚで使用して計算された極位置を示している。関心セクタ内で、１つの極（Ａと標示）は、最大の大きさを有し、閾強度値より上であり、ＲＳＡ極として識別される。この極位置から計算された強度値ＳRSAは、図２のステップ４２６として送信される出力ＲＳＡ測定値であり、当該出力ＲＳＡ測定値が演算機器側のバイオフィードバックアプリケーションで利用される。
【符号の説明】
【００３５】
１００  システム
１０２  ヒト被験者
１０４  センサー
１０６  通信リンク
１０８  演算機器
【図２】
2013111467.tif 000019
【図３】
2013111467.tif 000020
【図１】
2013111467.tif 000021
【外国語明細書】
2013111467.tif 000022
2013111467.tif 000023
2013111467.tif 000024
2013111467.tif 000025
2013111467.tif 000026
2013111467.tif 000027
2013111467.tif 000028
2013111467.tif 000029
2013111467.tif 000030
2013111467.tif 000031
2013111467.tif 000032
2013111467.tif 000033
2013111467.tif 000034
2013111467.tif 000035
2013111467.tif 000036
2013111467.tif 000037
2013111467.tif 000038
2013111467.tif 000039
2013111467.tif 000040
2013111467.tif 000041