JP7460584B2

JP7460584B2 - 心拍データ解析装置及びプログラム

Info

Publication number: JP7460584B2
Application number: JP2021138627A
Authority: JP
Inventors: 耕司吉田; 一義和田
Original assignee: KDDI Corp; Tokyo Metropolitan Public University Corp
Current assignee: KDDI Corp; Tokyo Metropolitan Public University Corp
Priority date: 2021-08-27
Filing date: 2021-08-27
Publication date: 2024-04-02
Anticipated expiration: 2041-08-27
Also published as: JP2023032476A

Description

本発明は、心拍データ解析装置及びプログラムに関する。

心拍データを解析する手法として次のような枠組みに即した既存手法が存在し、特許文献１等でこの既存手法に即した解析が行われている。図１はこの既存手法の枠組みを説明するための模式図である。

この既存手法は、次のような知見に基づくものである。心拍の心室収縮の開始を表すR波の発生間隔（R-R間隔）は、周期性をもって変動し、その変動の周波数成分は、自律神経の働き具合を表すことが知られている。この関係を利用して自律神経の働き具合を評価でき、その一つの手法（図１）では、心拍データにおけるR-R間隔を順次計測し、そのR-R間隔を周波数解析（スペクトル分析）してその変動の周波数成分を求め、その周波数成分から自律神経の働き具合を評価する。

図１のデータD1は、心電センサなどの心拍センサにより取得された心拍データのグラフ例を示し、横軸方向を時間進行方向として、時系列上のR波の並びとして心拍データが構成される。ここで、R波は心室収縮の開始を表すことにより個々の脈動を表し、時間軸方向（ｘ軸方向）にR-R間隔RR1,RR2,RR3,…をもって出現し、その間隔は、周期性をもって変動する。

一般的に、時間軸方向のデータをそのまま周波数解析することは難しい。すなわち、図１のデータ例D1から直接、矢印Aで示すようにその周波数解析結果としてのデータ例D3を得ることは困難である。そこで、データ例D1に示されるようなR-R間隔変動を周波数解析する場合、データ例D2に示されるようにデータD1に対して中間関数を適用した中間データとしての時系列データを生成し、この中間データとしての時系列データD2を周波数解析することにより、データ例D3に示されるような周波数解析結果を得る。

このように周波数解析を可能とするための中間データD2を得るための中間関数としては、心拍データの解析目的等に応じて多種類の中間関数が提案されている。図１では、中間関数の例として、R-R間隔RR1,RR2,RR3,…のそれぞれ（時間軸方向に順次現れる幅）に対して定数の振幅値としてその逆数P1=1/RR1,P2=1/RR2,P3=1/RR3,…を対応付ける中間関数の例（非特許文献１）が示されている。なお、非特許文献１では上記の通り、隣り合うR波発生時刻の間で、そのR-R間隔の逆数（この逆数は瞬時心拍数を意味する）として中間関数を定義し、それをスペクトル解析して心拍揺らぎを実時間的にとらえることが記載されている。

図１にて、上記の中間データD2を周波数解析して得た結果データD3はパワースペクトル密度を表し、横軸は周波数(Hz)であり、縦軸はパワースペクトル密度(msec²/Hz)を表す。信号処理分野において既知の通り、パワースペクトル密度とは、時系列データのパワーを各周波数成分に分解して表したものである。ここで、データD3内にも同符号を用いて模式的に示されるように、所定の低周波側帯域として例えば0.04～0.15Hz成分の総和をLFとし、所定の高周波側帯域として例えば0.15～0.40Hz成分の総和をHFとすると、これらの商LF／HFは交感神経の働き具合を示し、HFは副交感神経の働き具合を示すことが知られているので、これらの値を自律神経指標とすることで、心拍データの解析結果を得ることができる。

特開2015-080624号公報

BME Vol．8，No．10,1994 p．13－16 「心拍ゆらぎの1／fスペクトル推定」

以上のような従来手法の心拍データ解析に関して、一般的に医学分野では、その自律神経指標の「（定常的）状態」を「正確」に評価することや、高い精度、再現性を得ることを目的としている。したがって、自律神経指標を得るために、比較的長い時間間隔(例えば90～180秒程度)の時間関数（心拍データ）を周波数解析の対象としている。

一方で、解析の用途によっては、比較的短い時間間隔（例えば10～80秒程度）の時間関数（心拍データ）を周波数解析の対象とすることで、自律神経指標の瞬間的状態あるいは状態変化を評価することが望まれる場合もある。この場合、可能な限り短時間間隔として例えば20秒の時間関数を周波数解析の対象とすることで、時間進行に沿って連続的に周波数解析を行うこととなる。

このような短時間間隔の時間関数を周波数解析の対象とする際に、顕著な問題となりうる事項として、脈動の検出漏れの発生がある。１つの脈動の検出漏れが発生した場合に、これに対処せずにそのまま周波数解析を行ったとすると、周波数解析結果が大きく乱れてしまう。

図２は、このように短時間間隔の時間関数を周波数解析対象とする場合の検出漏れ（データ欠落）の影響の例を、周波数解析結果として図１のデータD3の高周波側帯域HF及び低周波側帯域LFのうちの後者（低周波側帯域LF）のみに周波数成分が存在するデータの解析結果を示す図である。（すなわち、図２の当該データは、検出漏れの影響を容易に見て取ることを可能にするために、低周波側帯域LFのみに周波数成分が存在するデータとして、実験的なデモ用に恣意的に時間領域上で生成したものである。）

図２にて、グラフGAはデータ欠落がない時間間隔の周波数解析結果であり、5つのパワースペクトル密度のグラフ例が、上側のグラフから順に、12,18,24,48,96パルス（脈動数）に相当する、それぞれ約10,14,19,38,77秒の間隔のウィンドウサイズの時間関数を用いた場合の例として示されている。（すなわち、解析対象となるある共通の時間関数（心拍データ）を、5種類のそれぞれ異なるウィンドウサイズで切り出して周波数解析した結果を列挙したのがグラフGAである。）グラフGBは、グラフGAと同様の解析（同様の5種類の間隔のウィンドウサイズによる）を、データ欠落がある時間関数（心拍データ）に対して適用した結果の例である。グラフGAとグラフGBとではほぼ同様の心拍データを解析対象としているのにもかかわらず、グラフGBの心拍データにおいて１つの偶発的なデータ欠落が存在しているという違いだけで、グラフGA,GBの周波数解析結果が大きく変化してしまっている様子を見て取ることができる。

このようなデータ欠落が存在する時間関数（心拍データ）を扱う場合であっても、グラフGCとして示すように、データ欠落が存在しないグラフGAと概ね同様の周波数解析結果となるように補正する手法が望まれるが、従来技術ではこのように適切に補正することができなかった。なお、図２では、グラフGA,GB,GCに関して、ウィンドウサイズが共通となるグラフ同士の対応関係を、グラフGA,GB間とグラフGA,GC間とにおいて点線で結ぶことにより示している。

なお、データ欠落（脈動の検出漏れ）の原因としては、心拍データの計測装置における計測の際のノイズ等による検出漏れ（心拍データには脈動が存在しているが計測できていない場合）のほか、心拍データにおける不整脈の存在が挙げられる。不整脈については、健康な場合であってもストレス等の一時的な要因等により偶発的に発生することがあることが知られている。このようなデータ欠落の原因の如何を問わず、データ欠落の影響を適切に補正することが望まれる。

なお、図２の例のように周波数解析結果が大きく変化してしまう原因として、データ欠落が存在する場合に中間関数に変換すると、特異なデータが発生してしまうことが挙げられる。図３は、データ欠落がある場合の中間関数の例を、データ欠落がない場合と対比で示す図である。データD4は正常な心拍データ（欠落のない時間関数）の例であり、これより得られる中間関数がデータD5として示される。このデータD4,D5は、図１のデータD1,D2と同様の例であり、中間関数としてR-R間隔の逆数を用いる例となっている。

一方で、図３ではデータD4,D5の対比例としてさらに、データD6は不整脈のある心拍データ（欠落がある時間関数）の例であり、これより得られる中間関数がデータD7として示される。データD6の時間関数では箇所R6においてデータ欠落が発生していることから、R-R間隔が突発的に2倍程度に広がってしまうこととなり、この結果としてデータD7の時間関数では同様に、対応する箇所R7においてR-R間隔が突発的に2倍程度に広がり、且つ同時に、その際の振幅値が1/2倍程度に小さくなってしまう。

このように、正常な場合の中間関数のデータD5とデータ欠落がある場合の中間関数のデータD7とを対比すると、データD7はデータ欠落の影響を受けた箇所R7が特異な値となってしまうことにより、図２に示したように、このような特異な値を有する中間関数を周波数解析すると、その結果が大きく乱れてしまうこととなる。

このようなデータ欠落への対処を検討した場合、例えば特許文献１では、時間関数上のノイズ箇所について、前後のデータ列の位相を乱さないデータ列の補間パターンを挿入する（ノイズ箇所を補間パターンで置換する）ことで、短期間の補間で生体変動の復元性を高め、スペクトル構造の復元性を高めるとされるが、このような補間が適切に行われるためには長めの時間関数を用いる必要があり、短時間間隔では適切に補正することができなかった。（すなわち、短時間間隔の場合において前後のデータ列で十分に補完するためには、当該短時間間隔の範囲を逸脱して補間する必要があるため、自律神経指標の瞬間的状態あるいは状態変化を必要とする場合，その周波数成分の変化量に大きく影響してしまう。また、短時間間隔の範囲内のみで補完を試みようとしても、そもそも補間に十分なデータが存在せず、適切に補間できないことから、周波数解析結果を適切に補正できない。）また、補間処理は、図１のデータD1に例示されるような時間関数（一般に複雑な変化挙動を示す関数）を対象とするものであり、必ずしも簡素に補間処理を行うことができないという問題もあった。

上記従来技術の課題に鑑み、本発明は、比較的短い時間関数を解析対象とする場合において、データ欠落が存在する場合であっても適切な補正を施してデータ欠落の影響を低減することのできる心拍データ解析装置及びプログラムを提供することを目的とする。

上記目的を達成するため、本発明は、時系列の心拍データからウィンドウで切り出した心拍データ区間より抽出される脈動間隔の列であるＲ－Ｒ間隔列を読み込み、当該Ｒ－Ｒ間隔列を変換して中間関数を生成する中間関数生成部と、前記中間関数を周波数解析した結果を得る周波数解析部と、を備える心拍データ解析装置であって、前記心拍データ区間における脈動の欠落箇所の有無を判定する欠落判定部を備え、前記欠落判定部において脈動の欠落があると判定された場合、前記中間関数生成部は、前記Ｒ－Ｒ間隔列を変換して中間関数を生成したうえでさらに、当該中間関数における脈動欠落に相当する箇所を、脈動欠落がない状態へと補正することにより、前記周波数解析部で周波数解析の対象となる中間関数を生成することを特徴とする。また、コンピュータを前記心拍データ解析装置として機能させるプログラムであることを特徴とする。

本発明によれば、欠落影響のある箇所を中間関数上において脈動欠落がない状態へと補正することにより、データ欠落が存在する場合であっても適切な補正を施してデータ欠落の影響を低減することができる。

心拍データを周波数解析する既存手法の枠組みを説明するための模式図である。短時間間隔の時間関数を周波数解析対象とする場合の検出漏れ（データ欠落）の影響(と補正の効果)の例を示す図である。データ欠落がある場合の中間関数の例を、データ欠落がない場合の中間関数との対比で示す図である。一実施形態に係る心拍データ解析装置の機能ブロック図である。一実施形態に係る心拍データ解析装置の動作のフローチャートである。リアルタイムで処理する際のタイミング及びウィンドウ幅の設定例を示す図である。 R-R間隔の逆数として中間関数を生成する例を示す図である。補正部による中間関数の補正処理を説明するための例を示す図である。解析部によって得られるリアルタイムの各時刻でのパワースペクトル密度のグラフ例を示す図である。中間関数としてビート関数を用いる場合の補正を説明するための例を示す図である。一般的なコンピュータ装置におけるハードウェア構成を示す図である。

図４は、一実施形態に係る心拍データ解析装置10の機能ブロック図である。心拍データ解析装置10は、R-R間隔抽出部1、欠落判定部2、生成部31及び補正部32を含む中間関数生成部3並びに解析部4を備える。心拍データ解析装置10は、その全体的な動作として、心拍データを計測する計測装置20より心拍データを入力として読み込み、この心拍データ（時間関数）を中間関数に変換したうえでその周波数解析結果を出力する。なお、図４では計測装置20を心拍データ解析装置10の外部構成としているが、計測装置20は心拍データ解析装置10に含まれるものであってもよい。

図５は、一実施形態に係る心拍データ解析装置10の動作のフローチャートである。以下、図５の各ステップを説明しながら、図４の心拍データ解析装置10の各機能部の処理について説明する。

図５のフローが開始されると、ステップS1では、計測装置20より現時刻までの最新の心拍データ（時間関数）をR-R間隔抽出部1において読み込んでから、ステップS2へと進む。当該読み込む最新の心拍データは、ステップS6の更新間隔に渡るものである。すなわち、この更新回数をi(i=1,2,3,…)としてステップS1の時刻をT(i)とすると、ステップS1では１つ前の時刻T(i-1)から現在時刻T(i)までの最新の心拍データを（まだ読み込んでいなかったものとして追加で）読み込めばよい。（i=1の初回では過去時刻T(0)から現時刻T(1)まで待機して読み込めばよい。）この読み込む最新の現時刻T(i)がウィンドウ端に一致するように設定しておいてもよい。

すなわち、図５のフローは、計測装置20において心拍データの取得及び解析の対象となるユーザに対してリアルタイムで心拍データの計測を行い、この心拍データ（時間関数）から現在時刻の側の直近のウィンドウサイズ分を切り出して、このウィンドウサイズ内の時間関数をリアルタイムで周波数解析することを、リアルタイムの各時刻（各処理タイミング）において繰り返す実施形態に関するものである。（なお、後述する図５のステップS6が、このリアルタイムでの繰り返し処理の各時刻の更新管理を行うステップである。）一方で、リアルタイムの計測及び解析ではなく、計測装置20において一定期間に渡って既に計測済みの心拍データを順次、心拍データ解析装置10においてスライディングウィンドウを用いて解析することも、図５のフローと同様にして可能である。

ここで、図５のフローに関してリアルタイムで処理する際のタイミング及びウィンドウ幅は、図６の模式例に示されるように、種々の所定の設定を用いてよい。すなわち、図６の例EX1に示すように、リアルタイムの順次の処理タイミングである時刻t1,t2,t3,…において当該各時刻までのウィンドウW1,W2,W3,…を読み込み、各ウィンドウが接している（あるウィンドウの最後の時刻が次のウィンドウの最初の時刻と一致している）ようにしてもよい。また、例EX2に示すように、リアルタイムの順次の処理タイミングである時刻t4,t5,t6,…において当該各時刻までのウィンドウW4,W5,W6,…を読み込み、各ウィンドウが離れている（あるウィンドウの最後の時刻と次のウィンドウの最初の時刻との間に間隔が存在する）ようにしてもよい。また、例EX3に示すように、リアルタイムの順次の処理タイミングである時刻t7,t8,t9,…において当該各時刻までのウィンドウW7,W8,W9,…を読み込み、各ウィンドウが部分的に重複している（あるウィンドウの最後側の部分と次のウィンドウの最初側の部分とが重複する）ようにしてもよい。また、例EX3では隣接ウィンドウ間の重複割合がウィンドウ幅の2割程度であるが、この重複割合も任意のものを設定できる。（例えばウィンドウ幅の9割が重複して、各ウィンドウが、両隣の2つの隣接ウィンドウよりもさらに遠方のウィンドウとも重複していてもよい。）

計測装置20は、心拍データを取得する対象とされるユーザにウェアラブル装置等の形式で装着され、任意の既存手法により任意の既存形式の心拍データを取得することができる。例えば、ECG（心電図）センサとして皮膚上に装着した電極を介して心電図の形式で心拍データを得てもよいし、PPG（光学式心拍）センサとして血管中の血流体積を光学的に検出して脈波の形式で心拍データを得てもよい。その他の種類の心拍センサを利用してもよい。

ステップS2では、ステップS1で読み込んだ最新の時間関数（心拍データ）を含む現時刻までの時間関数を対象として以下の（１）～（３）の各処理を行ってからステップS3へと進む。

（１） R-R間隔抽出部1が、現時刻までの時間関数よりR-R間隔を抽出し、このR-R間隔を生成部31（及び図４中に線L1で示されるように実施形態によっては欠落判定部2）へと出力する。前掲の図１のデータ例D1に関して説明した通り、R-R間隔は、時間関数の振幅値がピークとなり、且つ、当該ピークの値が閾値を超える箇所をR波の発生時刻と判定することで、当該R波の発生時刻の間隔がなす時系列として抽出することができる。

（２）生成部31が、R-R間隔抽出部1で算出したR-R間隔を変換して中間データとしての中間関数を生成し、この中間関数を解析部4又は補正部32へと出力する。ここで、生成部31で生成した中間関数は、後述するステップS3で肯定判定となりステップS4に進む場合に補正部32へと出力され、後述するステップS3で否定判定となりステップS4がスキップされる場合に、解析部4へと出力される。また、生成部31で生成した中間関数は、実施形態によっては図４中に線L2で示されるように欠落判定部2へも出力される。

前述の通り中間データへ変換するための中間関数には様々な種類が存在するが、以下では説明例として生成部31が図７に例示されるように、R-R間隔の逆数として中間関数を生成する場合について説明することとする。図７（前掲の図１，３の例でも同様）に示されるように、読み込んだウィンドウ内の時間関数（心拍データ）に対してR-R間隔抽出部1の抽出結果として、閾値超えのピーク時刻t_k(k=0,1,2,3,4,…)が順次に判定され、これらピーク時刻の隣接間隔としてR-R間隔I_k+1=t_k+1-t_k(k=0,1,2,3,4…)が抽出されているとする。このR-R間隔を、その逆数を振幅に持つ中間関数に変換すると、関数f(t)=1/I_k+1(t_k≦t<t_k+1)の形の中間関数が得られる。

（３）欠落判定部2が、ステップS1で読み込んだ時間関数（心拍データ）における欠落の有無を判定し、判定結果（欠落ありの判定結果の場合は、欠落発生箇所の情報も含む）を補正部32へと出力する。この欠落有無の判定は、図４中に線L1で示すように、R-R間隔抽出部1で得たR-R間隔（複数のR-R間隔が並んだ時系列データとしてのR-R間隔列）のうち、所定の閾値を超えて長いと判定されるR-R間隔が存在する場合に、当該長いと判定されたR-R間隔の箇所に欠落が発生しているものとして判定すればよい。（なお、単一の「R-R間隔」と、これらの複数が時間順で並んだ「R-R間隔列」とに関して、区別が文脈上明らかである場合にいずれも「R-R間隔」と呼ぶこととする。）

ここで、一般的にヒトの心拍間隔は、3割程度の変化に留まると言われている。したがって、閾値の設定例として、平均R-R間隔の3割を超過した場合（R-R間隔が平均R-R間隔の1.3倍を超過した場合）、または、前後のR-R間隔（前後２つのR-R間隔のうち少なくとも片方、あるいは両方、あるいは両方の平均）の1.3倍を超過したした場合にデータ欠落があるものと判定してもよい。この平均R-R間隔は、予め固定値として求めておいてもよいし、R-R間隔抽出部1において直近の一定期間（例えば現在から３分前の過去までの期間）におけるR-R間隔の履歴を記録しておき、この記録から各時刻において動的に平均R-R間隔を算出し、これに応じた閾値判定を欠落有無判定部2において行うようにしてもよい。

あるいは、欠落判定部2では同様の判定を、線L0に示されるようにステップS1で読み込んで更新して得られた現時刻までの時間関数（心拍データ）を直接解析して行うようにしてもよいし、線L2で示されるように生成部31で得た中間関数を解析して行うようにしてもよい。

ステップS3では、上記ステップS2の（３）の処理における欠落判定部2の判定結果が「欠落あり」であったか否かによる場合分けが行われ、肯定の場合（欠落が存在すると判定された場合）はステップS4へと進み、否定の場合（欠落が存在しないと判定された場合）は図５中に示されるようにステップS4はスキップして、ステップS5へと進む。

ステップS4では、補正部32が生成部31から得られた中間関数を補正して、補正された中間関数を解析部4へと出力してから、ステップS5へと進む。上記ステップS3の判定結果により、この中間関数は当初の時間関数（心拍データ）に欠落が発生しているものであり、前掲の図２，３で説明したようにそのまま周波数解析するには不適切であるため、補正部32においてこれを補正し、周波数解析するのに適した補正された中間関数とする。

図８は、補正部32による中間関数の補正処理を説明するための例を示す図である。補正対象となる中間関数の例としてデータD7が示され、このデータ例D7は図３で示したものと同一であって、R-R間隔の逆数を振幅とすることによって中間関数を得ており、箇所R7がデータ欠落の存在によって特異な値となっている。説明のため、図８のデータD7にも描かれているように、この欠落発生箇所R7のR-R間隔（時間幅）をI_nとすると、その逆数を振幅に持つ中間関数であることからこの欠落発生箇所R7の振幅はf(t)=1/I_n（この時刻tは時間幅I_nの箇所R7に属する範囲内にある）となる。なお、図８のデータ例D7に示されるような欠落発生箇所R7（ウィンドウ内の時間範囲）の情報は、欠落判定部2において欠落が存在する旨が判定された際に取得されており、補正部32へとこの情報が出力されている。

図８にさらに示されるように、補正部32は、データD7における欠落発生箇所R7を、データD7Cのような形の中間関数に加工することによって補正結果を得る。データD7Cに関して、図８では上段側に全体的な様子が、下段側にこの一部分としての補正対象となる欠落発生箇所R7及びその近傍の拡大図が示されている。この補正処理においては、データD7Cに示されるように、欠落発生箇所R7のみが補正の対象となり、当該ウィンドウサイズ内において、欠落発生箇所R7以外の箇所は補正されない。また、補正の際に参照する中間関数の情報は、欠落発生箇所R7の最近傍にあるR-R間隔及びその逆数としての振幅値となる。

すなわち、欠落発生箇所R7（時間幅I_n及び一定振幅1/I_nで構成される）は、前掲の図３のデータ例D6,D7を参照して説明した通り、本来であれば発生しているべき１回の脈動が欠落している状態にあるため、欠落状態として不自然に長い1回分のR-R間隔（及び不自然に小さい振幅）が反映された時間幅I_n及び一定振幅1/I_nで構成される状態から、本来発生しているべきであった自然な短い２回分のR-R間隔（及び自然な大きさの振幅）に相当する、「時間幅I_m1及び一定振幅1/I_m1」及び「時間幅I_m2及び一定振幅1/I_m2」で構成される状態へと置き換えることにより補正する。

具体的に、欠損発生箇所R7の最近傍にある中間関数の情報のみを補正のための元データとして参照することによって自然な補正結果を得るべく、以下の式(1),(2)を満たすように、２回分のR-R間隔に相当する２つの時間幅I_m1及びI_m2（及びこれらの逆数としての振幅）を求めることで、補正された中間関数を得るようにすればよい。

式(1)は、時間幅I_nを有する欠落発生箇所R7を、合計２回のR-R間隔の和として等しい時間幅となるような２つのR-R間隔I_m1及びI_m2（I_m1が前半側で、I_m2が後半側）に置き換えることを意味している。この置き換えを行う場合は、逆数を振幅に持つ中間関数として、前半側のR-R間隔I_m1の振幅が逆数1/I_m1となり、後半側のR-R間隔I_m2の振幅が逆数1/I_m2となるが、これら振幅が補正結果として自然な値となるようにするための制約が式(2)である。２つの変数としての２つのR-R間隔I_m1及びI_m2に対して、式(1)及び(2)の２つの制約を課すことで、この２つの変数の値を決定し、補正された中間関数（データ例D7C）を得ることができる。

式(2)の制約は幾何的には、図８のデータ例D7C（下段側の拡大図）に描いて示される通りである。すなわち、補正対象である欠落発生箇所R7の直前に位置するR-R間隔R7B（時間幅I_n-1及び一定振幅1/I_n-1で構成される）の中点m_n-1と、補正対象である欠落発生箇所R7の直後に位置するR-R間隔R7A（時間幅I_n+1及び一定振幅1/I_n+1で構成される）の中点m_n+1と、を結ぶ直線Lを中間関数のグラフ上に引いた際に、この直線L上に、補正された前半側のR-R間隔I_m1の中点m_m1と補正された後半側のR-R間隔I_m2の中点m_m2とが乗るようにする、という制約である。

なお、上記の中点m_n-1, m_m1, m_m2, m_n+1とは、図８からも明らかなように、逆数の振幅を有する中間関数の定義から、そのグラフ形状が「短冊関数」（R-R間隔に応じて変動する横幅及びこの横幅の逆数を高さとする矩形（短冊）が時間軸上で並ぶ形状の関数）となる中間関数における、各「短冊」（矩形）の上辺（振幅値がなす辺）の中点を意味するものである。（すなわち、短冊関数は、ステップ関数f(x)=0(x<0),f(x)=1(x≧1)のステップ位置x=0とステップ高さ1とを様々に変えた（ステップ高さは負の値でもよい）、複数のステップ関数の和として表現できるものである。）

以上、図８で説明した補正部32の補正手法を第１実施形態として、さらにこの変形例である第２，第３実施形態を説明する。

第２実施形態は、第１実施形態の式(1)及び式(2)の制約を用いることに変えて、式(1)及び以下の式(3)を用いることで補正された前半側のR-R間隔I_m1及び補正された後半側のR-R間隔I_m2を定めることができる。
I_m1=I_m2 …(3)

すなわち、第２実施形態は、図８等で示される欠落発生箇所R7の時間幅I_n（R-R間隔I_n）の時間軸上の中点に、本来発生しているべき脈動があったと推定して補正結果を得るものである。式(3),(1)より、補正された２つの等しいR-R間隔及び対応する振幅は以下の通りとなり、データ例D7の欠落発生箇所R7を中間関数として自然な値に補正することが可能となる。
I_m1=I_m2=I_n/2
→R-R間隔I_m1,I_m2は当初の半分の2回分となる。
1/I_m1=1/I_m2=2*(1/I_n)
→振幅値1/I_m1,1/I_m2は当初の不自然な値1/I_nの2倍となり自然な値となる。

なお、第２実施形態は、第１実施形態（及び第３実施形態）とは異なり、中間関数を補正するに際して、欠落発生箇所R7以外の中間関数のデータを補正元となるデータとして参照することが不要となる。

第３実施形態は、上記の第２実施形態が時間軸上の中点（すなわち、中間関数の横軸方向の中点）として本来の脈動を推定して補正結果を得たのに対して、中間関数の縦軸方向の中点として、本来の脈動に対応する中間関数の補正結果を得るものである。

すなわち、図８で示したのと同様に、欠落発生箇所R7の直前に位置するR-R間隔R7B（時間幅I_n-1及び一定振幅1/I_n-1（＝y_n-1とする）で構成される）と、欠落発生箇所R7の直後に位置するR-R間隔R7A（時間幅I_n+1及び一定振幅1/I_n+1（＝y_n+1とする）で構成される）と、を補正のための参照データとして利用し、欠落発生箇所R7の補正された振幅y_nを以下のようにこれら前後の振幅y_n-1及びy_n+1の平均値として算出することで、補正結果を得る。
y_n=(y_n-1+y_n+1)/2

この第３実施形態では、補正された振幅を上記平均値y_nとし、時間軸方向の長さについては補正前のI_nと同様の長さであるものとすればよい。（すなわち、中間関数の欠落発生箇所R7（時間軸での長さI_n）について、補正前の振幅1/I_nを振幅y_nへと変更したものを補正結果とすればよい。この補正後の箇所については、「振幅×時間軸長さ=y_n×I_n≠1」となるため、仮に１回のR-R間隔がこの補正後の箇所に対応しているものとすると、時間間隔の逆数を振幅として有する中間関数の定義に反した振幅値となっているが、振幅値自体は自然な値に補正されているため、自然な周波数解析結果を得ることが可能となる。（あるいは、補正後の箇所では、時間軸長さI_nの範囲内に、それぞれが「振幅y_n×時間軸長さ=1」となるような時間軸長さのR-R間隔が複数(整数でなくてもよい)だけ存在するものと解釈してもよい。）すなわち、第３実施形態では、時間軸方向の長さがI_nとなる欠落発生箇所R7のどの位置に本来の脈動が存在したのかという点で第１実施形態や第２実施形態のように明確な推定結果を得ることはできないが、周波数解析する対象である中間関数の形において、自然な補正結果を得ることが可能となる。）

ステップS5では、解析部4が、中間関数生成部3で生成された中間関数に対して周波数解析を行い、その結果を出力してからステップS6へと進む。解析部4での解析の対象となる中間関数は、ステップS4からステップS5へ進んだ場合（データ欠落が存在する場合）は生成部31で生成され、さらに、補正部32でデータ欠落の影響を低減するように補正された中間関数であり、ステップS4をスキップしてステップS3からステップS5へ進んだ場合（データ欠落が存在しない場合）は、生成部31で生成された中間関数である。また、ステップS5で周波数解析となる中間関数の時間軸上での範囲は、図６で説明したような所定ウィンドウ範囲であり、現在時刻まで得られている中間関数を包含する最新時刻側のウィンドウ範囲（過去のステップS5では未処理であったウィンドウ範囲）を、新たな解析対象とすればよい。

データ欠落が存在する場合、補正部32では以上の説明のようにしてデータ欠落が存在しない状態に相当するものとなるように中間関数を補正するので、解析部4における周波数解析結果をデータ欠落の影響を低減させた適切な結果として得ることが可能となる。すなわち、前掲の図２で説明したような、データ欠落が存在する場合であってもグラフGBのような不自然な周波数解析結果ではなく、補正が反映されたグラフGCのように自然な周波数解析結果（本来のグラフGAに近い周波数結果）を得ることが可能となる。

解析部4では前掲の図１の例におけるデータD2からデータD3を得るのと同様の処理で周波数解析結果を得ることができる。すなわち、中間関数を周波数解析してパワースペクトル密度F(f)を得て、以下の式のように所定の低周波側帯域として例えば0.04～0.15Hz成分の総和をLFとし、所定の高周波側帯域として例えば0.15～0.40Hz成分の総和をHFとする（以下の式）として算出する。

商LF／HFは交感神経の働き具合を示し、HFは、副交感神経の働き具合を示すことが知られているので、これらの値（LF／HF及びHF）を自律神経指標とすることで、周波数解析に基づく心拍データの解析結果を得ることができる。

図９に、解析部4によって得られるリアルタイムの各時刻でのパワースペクトル密度のグラフ例を示す。パワースペクトル密度がリアルタイムで得られることから、自律神経指標である値（LF／HF及びHF）も同様に、リアルタイムで得ることが可能となる。

ステップS6では、時刻を次の時刻（リアルタイムの次の処理タイミング）に更新してからステップS1へと戻ることで、以上説明した処理が心拍データ解析装置10においてリアルタイムの各時刻について繰り返し実行される。

以上、本発明の実施形態によれば、自律神経指標を求める際に、不整脈あるいは計測誤差等に起因したデータ欠落（脈動の検出漏れ）の影響を少なくすることができる。この際、周波数成分に影響を及ぼさないように、心拍間隔（R-R間隔）ではなく、中間関数データ上で補正することで、適切かつ簡素な補正が可能である。この補正にあたり、他の区間（ウィンドウ）のデータを参照しないことで、比較的短いウィンドウを対象としてリアルタイム解析に即した補正及び自律神経指標の瞬間的状態あるいは状態変化を導く要求に対する補正が可能である。

なお、以下（４）において後述するように、データ欠落がウィンドウ端で発生している場合は、例外的にウィンドウ外の中間関数も参照して補正することとなるが、参照箇所はウィンドウ外部に範囲を広げる際の最小範囲としてR-R間隔１個分だけでよく、ウィンドウ外部の広範囲な中間関数を参照することは不要であるため、実質的に単一ウィンドウ内のみで補正しているのと同様である。

以下、種々の追加例、補足例、代替例などに関して説明する。

（１）本発明の実施形態によれば、自律神経指標のリアルタイムでの高精度な把握が可能となることから、例えば精神的なストレスを原因とする疾病等の予防等へと利用することで、国連が主導する持続可能な開発目標（ＳＤＧｓ）の目標３「あらゆる年齢のすべての人々の健康的な生活を確保し、福祉を推進する」に貢献することが可能となる。（なお、データ欠落の一因となる不整脈が偶発的なものではない場合には、別途の医療的な対処が望まれる。）

（２）以上の説明では、生成部31においてR-R間隔の逆数を振幅値とする中間関数を用いる場合の補正に関して説明したが、その他の中間関数でも同様の補正処理が可能である。例えば逆数（R-R間隔の-1乗）以外で振幅値を定める中間関数であっても、逆数の場合と同様にグラフ形状が短冊関数となるような中間関数であれば、当該中間関数の個別制約のもとで、図８で説明した第１実施形態や、その変形例としての第２，第３実施形態による補正が可能である。

例えば第１実施形態であれば、図８で説明したのと同様に、短冊関数となる中間関数におけるデータ欠落発生箇所R7を補正するのに、その直前の１回のR-R間隔に相当する中間関数部分R7Bと、その直後の１回のR-R間隔に相当する中間関数部分R7Aと、のみを補正元データとして参照すればよい。変形例として、直前及び直後のそれぞれにおいて、２回分以上のR-R間隔に相当する中間関数部分を補正元データとして参照するようにしてもよい。（ただし、データ欠落発生箇所R7からできる限り近い範囲内のみを補正元データとすることが望ましい。）

（３）生成部31における中間関数としてビート関数を用いる場合の補正部32における補正を、図１０を用いて説明する。ビート関数とは、R-R間隔がI_k(時間軸上で順番にk=1,2,3,…)である場合に、この時間軸上での順番kを離散時間kに置き換えることによって、中間関数をそのまま離散時間kでの等間隔のパルス列I_k(k=1,2,3,…)として得るものである。（すなわち、逆数の中間関数では時間軸は当初の心拍データと共通していたのと異なり、ビート関数の中間関数では時間軸はパルス番号に変換され、パルス高さに時間幅の情報が反映されることとなる。）

図１０ではデータ例D10としてこの離散時間パルス列が示され、3番目のパルスI₃がその他のパルスと比べて顕著に高いことから、欠落判定部2において欠落発生箇所として判定されているものとする。この場合、データ例D10Cとして示すように、欠落起因のパルスI₃を２つの補正後のパルスI_3-1,I_3-2に置換することで、以下に示すように補正した中間関数を得ることができる。
補正前の中間関数（データD10） I₁,I₂,I₃,I₄の順に並ぶ４個のパルス列
補正後の中間関数（データD10C） I₁,I₂,I_3-1,I_3-2,I₄の順に並ぶ５個のパルス列

２つの補正後のパルスI_3-1,I_3-2の高さについては、欠落発生箇所であった3番目のパルスI₃の直前のパルスI₂と直後のパルスI₄とから高さを内挿することにより、以下のように算出すればよい。当該内挿による算出は、図１０の補正後データ例D10Cにおいて線L10として模式的に示されている。
I_3-1=(2*I₂+I₄)/3
I_3-2=(I₂+2*I₄)/3

なお、上記ビート関数を中間関数とする際の補正例は、以下の読み替えを適用することで、図８で説明した逆数を中間関数とする際の補正と同様の方針で補正している。
3番目のパルスI₃＝データ欠落発生箇所R7
2番目のパルスI₂＝データ欠落発生箇所R7の直前のR-R間隔箇所R7B
4番目のパルスI₄＝データ欠落発生箇所R7の直後のR-R間隔箇所R7A

上記では中間関数としてビート関数を用いる場合を説明したが、スプライン補間と直線補間を用いる場合も概ね同様とすればよい。スプライン補間ではR-R間隔がI_k(k=1,2,…)であり、この間隔I_kの時刻tの範囲がt_k≦t≦t_k+1であるとすると、（手順１）時刻t_k(k=1,2,…)の値が間隔I_kである（これ以外の時刻tの値は0である）パルス列関数を求め、（手順２）このパルス列関数を３次スプライン補間し、（手順３）この３次スプライン補間から平均R-R間隔で（再）サンプリングしたパルス列関数を中間関数とする。

上記スプライン補間の中間関数については、（手順１）のパルス列関数を補正の対象とすればよく、欠落発生箇所について1個のR-R間隔で構成される状態から2個のR-R間隔で構成される状態とすればよい。すなわち、欠落発生箇所として1個の間違ったパルスを2個の自然なパルスに補正すればよく、その2個のパルスの時間軸上での位置を、前後のパルス位置を３等分した位置とし、当該2個のパルスの高さを前後のパルスの高さで内挿した高さとすることで、図１０のビート関数の場合と概ね同様に補正できる。

中間関数として直線補間を用いる場合も同様に補正できる。すなわち、直線補間の中間関数はスプライン補間における手順１～３のうち手順２，３で用いるものをスプライン補間から直線補間に置き換えて得られるため、スプライン補間の中間関数の場合に手順１で３等分して内挿した手法と全く同様の手法により、直線補間の中間関数の場合についても同様に補正することが可能となる。

（４）図８の例で、データ欠落発生箇所R7を補正するに際して、補正元データとして直前のR-R間隔箇所R7B及び直後のR-R間隔箇所R7Aの中間関数を参照するものとした。頻繁に発生する可能性は低いが、例外的な場合としてデータ欠落発生箇所R7がステップS5で解析に用いるウィンドウの端部に位置する場合、直前のR-R間隔箇所R7Bが当該ウィンドウの１つ前のウィンドウ内にある、または、直後のR-R間隔箇所R7Aが当該ウィンドウの１つ後のウィンドウ内にある、という状況が発生しうるが、全く同様にして補正するようにすればよい。補正においてはデータ欠落発生箇所R7の最近傍のみを参照するため、ウィンドウ外であっても周波数解析結果に大きく影響することはない。

（５）ステップS1で新たに読み込むデータに関しては、一定時間間隔すなわち時間関数（心拍データ）上の所定時間範囲としてもよいし、R-R間隔上での所定個数（所定数のR-R間隔）としてもよい。データ欠落があるものと判定された場合は、データ欠落を補正して1個から2個へと増やしたR-R間隔においてR-R間隔の個数を定めるようにしてもよい。

（６）図１１は、一般的なコンピュータ装置70におけるハードウェア構成の例を示す図である。心拍データ解析装置10及び計測装置20の各々（前述の通りこれらは一体で構成されていてもよい）は、このような構成を有する１台以上のコンピュータ装置70として実現可能である。なお、２台以上のコンピュータ装置70で心拍データ解析装置10及び計測装置20の各々を実現する場合、ネットワーク経由で処理に必要な情報の送受を行うようにしてよい。コンピュータ装置70は、所定命令を実行するCPU（中央演算装置）71、CPU71の実行命令の一部又は全部をCPU71に代わって又はCPU71と連携して実行する専用プロセッサとしてのGPU（グラフィックス演算装置）72、CPU71（及びGPU72）にワークエリアを提供する主記憶装置としてのRAM73、補助記憶装置としてのROM74、通信インタフェース75、ディスプレイ76、マウス、キーボード、タッチパネル等によりユーザ入力を受け付ける入力インタフェース77、心拍センサ78と、これらの間でデータを授受するためのバスBSと、を備える。

心拍データ解析装置10及び計測装置20の各機能部は、各部の機能に対応する所定のプログラムをROM74から読み込んで実行するCPU71及び／又はGPU72によって実現することができる。なお、CPU71及びGPU72は共に、演算装置（プロセッサ）の一種である。ここで、表示関連の処理が行われる場合にはさらに、ディスプレイ76が連動して動作し、データ送受信に関する通信関連の処理が行われる場合にはさらに通信インタフェース75が連動して動作する。計測装置20における心拍計測機能を提供するセンサとして、ECGセンサやPPGセンサ等として構成される心拍センサ78を用いることができる。心拍センサ78を有する計測装置20に関して、CPU71及び／又はGPU72に代えて、心拍データを計測して通信インタフェース75において送信可能なデータ形式として得るための専用集積回路（ASIC（特定用途向けIC）やFPGA（Field Programmable Gate Array）等）を有していてもよい。

10…心拍データ解析装置、20…計測装置、1…R-R間隔抽出部、2…欠落判定部、3…中間関数生成部、4…解析部、31…生成部、32…補正部

Claims

時系列の心拍データからウィンドウで切り出した心拍データ区間より抽出される脈動間隔の列であるＲ－Ｒ間隔列を読み込み、当該Ｒ－Ｒ間隔列を変換して中間関数を生成する中間関数生成部と、
前記中間関数を周波数解析した結果を得る周波数解析部と、を備える心拍データ解析装置であって、
前記心拍データ区間における脈動の欠落箇所の有無を判定する欠落判定部を備え、
前記欠落判定部において脈動の欠落があると判定された場合、前記中間関数生成部は、前記Ｒ－Ｒ間隔列を変換して中間関数を生成したうえでさらに、当該中間関数における脈動欠落に相当する箇所を、脈動欠落がない状態へと補正することにより、前記周波数解析部で周波数解析の対象となる中間関数を生成し、
前記中間関数は、ｋを整数としてk番目のＲ－Ｒ間隔ごとに定数の振幅を有する短冊関数であり、
前記中間関数生成部は、k番目のＲ－Ｒ間隔をＩ _ｋとし、nを整数として、ｎ番目のＲ－Ｒ間隔Ｉ _n が前記脈動欠落に相当する箇所であると判定された場合に、以下の式(1)及び(2)を満たすように当該ｎ番目のＲ－Ｒ間隔Ｉ _n を２つのＲ－Ｒ間隔Ｉ _m1 及びＩ _m2 へと置き換えて、前記補正することを特徴とする心拍データ解析装置。
時系列の心拍データからウィンドウで切り出した心拍データ区間より抽出される脈動間隔の列であるＲ－Ｒ間隔列を読み込み、当該Ｒ－Ｒ間隔列を変換して中間関数を生成する中間関数生成部と、
前記中間関数を周波数解析した結果を得る周波数解析部と、を備える心拍データ解析装置であって、
前記心拍データ区間における脈動の欠落箇所の有無を判定する欠落判定部を備え、
前記欠落判定部において脈動の欠落があると判定された場合、前記中間関数生成部は、前記Ｒ－Ｒ間隔列を変換して中間関数を生成したうえでさらに、当該中間関数における脈動欠落に相当する箇所を、脈動欠落がない状態へと補正することにより、前記周波数解析部で周波数解析の対象となる中間関数を生成し、
前記中間関数は、ｋを整数としてk番目のＲ－Ｒ間隔ごとに定数の振幅を有する短冊関数であり、
前記中間関数生成部は、k番目のＲ－Ｒ間隔をＩ _ｋとし、nを整数として、ｎ番目のＲ－Ｒ間隔Ｉ _n が前記脈動欠落に相当する箇所であると判定された場合に、以下の式(1)及び(3)を満たすように当該ｎ番目のＲ－Ｒ間隔Ｉ _n を２つのＲ－Ｒ間隔Ｉ _m1 及びＩ _m2 へと置き換えて、前記補正することを特徴とする心拍データ解析装置。
時系列の心拍データからウィンドウで切り出した心拍データ区間より抽出される脈動間隔の列であるＲ－Ｒ間隔列を読み込み、当該Ｒ－Ｒ間隔列を変換して中間関数を生成する中間関数生成部と、
前記中間関数を周波数解析した結果を得る周波数解析部と、を備える心拍データ解析装置であって、
前記心拍データ区間における脈動の欠落箇所の有無を判定する欠落判定部を備え、
前記欠落判定部において脈動の欠落があると判定された場合、前記中間関数生成部は、前記Ｒ－Ｒ間隔列を変換して中間関数を生成したうえでさらに、当該中間関数における脈動欠落に相当する箇所を、脈動欠落がない状態へと補正することにより、前記周波数解析部で周波数解析の対象となる中間関数を生成し、
前記中間関数は、ｋを整数としてk番目のＲ－Ｒ間隔ごとに定数の振幅を有する短冊関数であり、
前記中間関数生成部は、k番目のＲ－Ｒ間隔をＩ _ｋとし、その短冊関数における定数の振幅をy _ｋとし、nを整数として、ｎ番目のＲ－Ｒ間隔Ｉ _n が前記脈動欠落に相当する箇所であると判定された場合に、その短冊関数における定数の振幅y _n を以下の式(4)によって定めることで前記補正することを特徴とする心拍データ解析装置。
前記欠落判定部は、前記Ｒ－Ｒ間隔列に含まれる複数のＲ－Ｒ間隔に、閾値判定で長いと判定されるＲ－Ｒ間隔が存在する場合に、前記心拍データ区間において当該長いと判定されたＲ－Ｒ間隔の箇所に脈動の欠落があるものとして判定することを特徴とする請求項１ないし３のいずれかに記載の心拍データ解析装置。
コンピュータを請求項１ないし４のいずれかに記載の心拍データ解析装置として機能させることを特徴とするプログラム。