JP2024145525A - 心拍データ解析装置及びプログラム - Google Patents
心拍データ解析装置及びプログラム Download PDFInfo
- Publication number
- JP2024145525A JP2024145525A JP2023057915A JP2023057915A JP2024145525A JP 2024145525 A JP2024145525 A JP 2024145525A JP 2023057915 A JP2023057915 A JP 2023057915A JP 2023057915 A JP2023057915 A JP 2023057915A JP 2024145525 A JP2024145525 A JP 2024145525A
- Authority
- JP
- Japan
- Prior art keywords
- window
- analysis device
- data analysis
- intermediate function
- unit
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
Images
Landscapes
- Measuring Pulse, Heart Rate, Blood Pressure Or Blood Flow (AREA)
- Measurement Of The Respiration, Hearing Ability, Form, And Blood Characteristics Of Living Organisms (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
【課題】心拍データを解析して効果的に帯域評価を行うことのできる心拍データ解析装置を提供する。
【解決手段】被験者の心拍データから算出されるRR間隔を順次読み込んで、心臓の拍動1拍分ごとに値が定まる中間関数を生成(11)し、直流成分を除去(12)して交流中間関数とする第1処理と、前記交流中間関数のサンプリング信号(13)に処理対象のウィンドウを設定(14)する第2処理と、前記ウィンドウ内のサンプリング信号に周波数解析を適用(16)してパワースペクトル密度を算出(17)する第3処理と、前記パワースペクトル密度のLF帯域とHF帯域のそれぞれでの加算値を自律神経指標として算出(18)する第4処理と、を実行する。
【選択図】図1
【解決手段】被験者の心拍データから算出されるRR間隔を順次読み込んで、心臓の拍動1拍分ごとに値が定まる中間関数を生成(11)し、直流成分を除去(12)して交流中間関数とする第1処理と、前記交流中間関数のサンプリング信号(13)に処理対象のウィンドウを設定(14)する第2処理と、前記ウィンドウ内のサンプリング信号に周波数解析を適用(16)してパワースペクトル密度を算出(17)する第3処理と、前記パワースペクトル密度のLF帯域とHF帯域のそれぞれでの加算値を自律神経指標として算出(18)する第4処理と、を実行する。
【選択図】図1
Description
本発明は、心拍データ解析装置及びプログラムに関する。
自律神経系を評価する従来技術の例として、特許文献1では疲労を客観的に評価するための指標として自律神経活動量(ccvTP)を、LF値+HF値の総和を時間中の心拍数で補正したものとして利用している。具体的には、トータルパワー(TP)の平方根を,時間中(Window)の平均心拍数「平均(RR)」で除算することで式(PR1)のようにccvTPを算出する。TPは低周波成分LFと高周波成分HFの和「TP=LF+HF」として求める。この低周波/高周波成分LF,HFには式(PR2)~(PR4)を用いる。式(PR4)のC(t)はRR間隔の自己相関関数である。
このLF,HFに関して、上記(PR2),(PR3)では自律神経と相関がある周波数帯についてLF:0.04-0.15Hz,HF:0.15-0.4Hzとの定義を利用しており、この定義は多くの論文で引用されている非特許文献1でも提唱されおり、一般的に広く利用されている。
Heart rate variability-Standards of measurement, physiological interpretation, and clinical use-European Heart Journal 1996
座位状態での心拍測定を用いたリアルタイムなストレス緩和システム、マルチメディア,分散,協調とモバイルシンポジウム 平成5年,7月
しかしながら、従来技術ではこのように自律神経と相関があることから自律神経活動の推定指標として使える周波数帯LF,HFをリアルタイムで効果的に周波数解析したい際に発生する要望について考慮されていなかった。
● 単位時間当たりの自律神経系のLF,HFの各周波数成分を定量的に比較したい。
すなわち、周波数解析を行うに際してFFT(高速フーリエ変換)の精度(分解能)の時間分解能と周波数分解能との間のトレードオフが制約として存在するが、当該制約された分解能のもとで、比較的狭い周波数帯LF,HFのそれぞれの帯域積分値を求めようとしても、ポイント数が不足して積分値の精度が確保できない(実際の周波数特性曲線の面積と、これの近似値としてのポイント毎の長方形面積の和と、の間の数値積分誤差が大きくなる)という問題がある。例えば非特許文献2では、周波数解析を行っているものの、このように比較的狭い周波数帯LF,HFに限定した帯域積分はそもそも行われておらず、帯域積分の精度を確保することが想定されていない。
すなわち、周波数解析を行うに際してFFT(高速フーリエ変換)の精度(分解能)の時間分解能と周波数分解能との間のトレードオフが制約として存在するが、当該制約された分解能のもとで、比較的狭い周波数帯LF,HFのそれぞれの帯域積分値を求めようとしても、ポイント数が不足して積分値の精度が確保できない(実際の周波数特性曲線の面積と、これの近似値としてのポイント毎の長方形面積の和と、の間の数値積分誤差が大きくなる)という問題がある。例えば非特許文献2では、周波数解析を行っているものの、このように比較的狭い周波数帯LF,HFに限定した帯域積分はそもそも行われておらず、帯域積分の精度を確保することが想定されていない。
● 解析対象時間を極力短くしたい。
また、リアルタイムでの周波数解析の精度を確保するにはウィンドウサイズを小さくすることが望まれるが、従来技術はウィンドウサイズを比較的大きく設定しており、小さなウィンドウを扱うことが想定されていない。例えば非特許文献2では180秒の長さのウィンドウサイズを用いており、リアルタイムでの急激な自律神経活動の変化を捉えようとしても、比較的長い180秒のウィンドウでは急激な変化となるべき波形が鈍ってしまう。
また、リアルタイムでの周波数解析の精度を確保するにはウィンドウサイズを小さくすることが望まれるが、従来技術はウィンドウサイズを比較的大きく設定しており、小さなウィンドウを扱うことが想定されていない。例えば非特許文献2では180秒の長さのウィンドウサイズを用いており、リアルタイムでの急激な自律神経活動の変化を捉えようとしても、比較的長い180秒のウィンドウでは急激な変化となるべき波形が鈍ってしまう。
上記従来技術の課題に鑑み、本発明は、心拍データを解析して効果的に帯域評価を連続的に行うことのできる心拍データ解析装置及びプログラムを提供することを目的とする。
上記目的を達成するため、本発明は心拍データ解析装置またはプログラムであって、被験者の心拍データから算出されるRR間隔を順次読み込んで、心臓の拍動1拍分ごとに値が定まる中間関数を生成し、直流成分を除去して交流中間関数とする第1処理と、前記交流中間関数のサンプリング信号に処理対象のウィンドウを設定する第2処理と、前記ウィンドウ内のサンプリング信号に周波数解析を適用してパワースペクトル密度を算出する第3処理と、前記パワースペクトル密度のLF帯域とHF帯域のそれぞれでの加算値を自律神経指標として算出する第4処理と、を実行し、前記第2処理では、サンプリング信号の元となるRR間隔が等しい個数だけ連続した範囲のそれぞれを、拍動基底ウィンドウとして設定する、または、互いに等しい長さの時間基底ウィンドウを設定することを特徴とする。
本発明によれば、心拍データを解析して効果的に自律神経指標に相当する帯域の変化を評価することができる。
図1は、一実施形態に係る心拍データ解析装置10の機能ブロック図であり、図示するように心拍データ解析装置10は中間関数生成部11、直流成分除去部12、標本化部13、ウィンドウ定義部14、解像度向上部15、周波数解析部16、振幅特性部17、等値化処理部18及び帯域積分部19を備え、この順番で処理を行う。
心拍データ解析装置10の全体的な処理としては図1にも示される通り、被験者の心拍データであるRR間隔をリアルタイムで読み込み、ウィンドウを設定して、当該ウィンドウ内でのLF積分値とHF積分値を、自律神経活動指標として出力する。当該分野において既知のように、心拍データの取得は、被験者に心電センサを取り付ける等により心電データ時系列を取得し、そのR波の間隔としてRR間隔を、心電データ時系列の順にRR間隔を並べた数列データ(i=1,2,3,…番目のRR間隔の値RRIiからなる数列{RRIi|i=1,2,3,…,})の形で取得すればよい。
以下、図1の各部の処理を順番に説明する。なお、以下の説明で「後段」とは、当該説明している機能ブロックより図1で後段側に配置されている機能ブロックを意味するものとする。(例えば、中間関数生成部11に関して「後段」は各部12~19を意味するものとする。)同様に、「前段」は「後段」の逆側に配置されている機能ブロックを意味するものとする。
<中間関数生成部11>
中間関数生成部11は、リアルタイムで入力される(単なる数列データとしての)RR間隔より、周波数解析可能な時間関数である中間関数を生成し、後段に出力する。図2は中間関数生成部11で生成する中間関数の模式図である。RR間隔はそのままでは単なる数列であるため、これを周波数解析可能とするために時系列データである中間関数に変換することが当該技術分野では行われている。中間関数としては、心拍データの解析目的等に応じて多種類の中間関数が提案されており、例えば、中間関数の例として、RR間隔RRI1,RRI2,RRI3,…のそれぞれ(時間軸方向に順次現れる幅)に対して定数の値としてその逆数P1=1/RRI1,P2=1/RRI2,P3=1/RRI3,…を対応付ける中間関数などが存在する。
中間関数生成部11は、リアルタイムで入力される(単なる数列データとしての)RR間隔より、周波数解析可能な時間関数である中間関数を生成し、後段に出力する。図2は中間関数生成部11で生成する中間関数の模式図である。RR間隔はそのままでは単なる数列であるため、これを周波数解析可能とするために時系列データである中間関数に変換することが当該技術分野では行われている。中間関数としては、心拍データの解析目的等に応じて多種類の中間関数が提案されており、例えば、中間関数の例として、RR間隔RRI1,RRI2,RRI3,…のそれぞれ(時間軸方向に順次現れる幅)に対して定数の値としてその逆数P1=1/RRI1,P2=1/RRI2,P3=1/RRI3,…を対応付ける中間関数などが存在する。
本実施形態では、i番目のRR間隔の値をRRIiとして、時刻tがRRIi内にある際の中間関数の値Pi>0を以下の式(21a)ように定義することで、図2にも模式的に示される通り、時刻tの値P(t)が以下の式(21b)となるように、中間関数を生成する。(すなわち、時刻tがi番目のRR間隔RRIiの時間範囲内にある場合にP(t)=Piの値となるものとして中間関数P(t)を定義する。なお、式(21b)でi=1の場合は「0≦t≦RRI1についてP(t)=P1を意味する。)
なお、式(21a)でCoは、一拍あたりのエネルギーが変化する前提で使う係数である。この係数は、心臓の一回拍出量や血管抵抗で定まる。例えば一回拍出量や血管抵抗の変化がないと仮定した場合はCo=1とすればよい。
なお、本実施形態において式(21a),(21b)で定義される、RR間隔の平方根の逆数による中間関数は、本出願人による特許第7221195号の手法と同様であり、同文献の段落[0020]に記載される通りの以下の効果を奏することができる。
「時系列データ全体にわたって1回の拍動あたりのパワースペクトル密度の総和が等しくなるように、時系列データの時間領域表現を最適化しているので、時系列データの時々の所定拍動回数区間のパワースペクトル密度の総和を等しくすることができ、その下で、それぞれのパワースペクトル密度を求めることができる。したがって、それらのパワースペクトル密度を比較するだけで、自律神経の働き具合の経時的変化を正しく評価することができるようになる。」
「時系列データ全体にわたって1回の拍動あたりのパワースペクトル密度の総和が等しくなるように、時系列データの時間領域表現を最適化しているので、時系列データの時々の所定拍動回数区間のパワースペクトル密度の総和を等しくすることができ、その下で、それぞれのパワースペクトル密度を求めることができる。したがって、それらのパワースペクトル密度を比較するだけで、自律神経の働き具合の経時的変化を正しく評価することができるようになる。」
上記の効果は、1回の拍動あたりのエネルギーである式(21b)の値(振幅)の2乗の積分値が以下の式(21c)の通り、一定値Coとなることによる。
<直流成分除去部12>
直流成分除去部12は、前段で得た中間関数から直流成分を除去し、直流成分が除去され、交流成分のみで構成される中間関数(以下、「交流中間関数」と呼ぶ)を後段に出力する。
直流成分除去部12は、前段で得た中間関数から直流成分を除去し、直流成分が除去され、交流成分のみで構成される中間関数(以下、「交流中間関数」と呼ぶ)を後段に出力する。
図3は、直流成分除去部12の処理の意義を説明するためのデータ例を示す図であり、データD31にそのグラフ例が示されるような中間関数は、多くの直流成分を含んでいる。解像度が高い場合問題ないが、交流成分が小さいことから、解像度が低くなると,周波数解析をしたい交流成分が直流成分のサイドローブに紛れてしまう。これはデータ例D32として直流成分を除去せずに中間関数のデータD31をそのままで周波数解析した例に示される通りである。データ例D33,D34に示されるように、直流成分を除去して交流中間関数としてから周波数解析をすることで、当該紛れる問題を回避し、解析すべき波形W33を明瞭に取得できる。
現在時刻tの時点においてn個のRR間隔{RRIi|i=1,2,…,n}が読み込まれてその中間関数P(t)が生成されているものとすると、直流成分除去部12では、以下の式(12a),(12b)のように中間関数P(t)の単位時間あたりの平均値avePnを直流成分として求め、中間関数P(t)から当該直流成分を減算することで交流中間関数P'(t)を得ることができる。図4はこの直流成分除去部12の処理の模式図である。
なお、上記の式(12a)の平均値は、後述するウィンドウ定義部14で拍動基底のウィンドウを設定する場合のものであり、ウィンドウ幅Wを有する時間基底のウィンドウを設定する場合は、例えばn+1個のRR間隔{RRIi|i=1,2,…,n}がこの幅Wのウィンドウ内に読み込まれて最後のn+1番目のRR間隔が端数となっている(n+1番目のRR間隔の途中にウィンドウ後端が位置している)場合に、同様にして以下の式(12aW),(12bW)で交流中間関数P'(t)を得ることができる。
<標本化部13>
標本化部13は、ディジタル信号処理を可能とすべく、交流中間関数P'(t)を時間軸上で等間隔にサンプリングして、サンプリング信号x(i)(i=1,2,…)を後段へ出力する。例えば、サンプリングレートを10msecとし、読み込まれているRR間隔{RRIk|k=1,2,…,n}の時間長さ単位をsecとすると、時間長さRRIk内において100×RRIk個のサンプリング信号x(i)を得ることができるので、図5にも模式的に示される通り、以下の式(131),(132),(133)に例示されるような形でサンプリング信号を得ることができる。
標本化部13は、ディジタル信号処理を可能とすべく、交流中間関数P'(t)を時間軸上で等間隔にサンプリングして、サンプリング信号x(i)(i=1,2,…)を後段へ出力する。例えば、サンプリングレートを10msecとし、読み込まれているRR間隔{RRIk|k=1,2,…,n}の時間長さ単位をsecとすると、時間長さRRIk内において100×RRIk個のサンプリング信号x(i)を得ることができるので、図5にも模式的に示される通り、以下の式(131),(132),(133)に例示されるような形でサンプリング信号を得ることができる。
<ウィンドウ定義部14>
ウィンドウ定義部14は、前段で得たサンプリング信号x(i)を周波数解析するためにウィンドウWを設定し、当該設定されたウィンドウ内のサンプリング信号x(i)∈Wを後段に出力する。すなわち、周波数成分の変化を求める場合には,時間長が異なる複数のデータを定量的に比較できるようにする必要があることからウィンドウを設定する。
ウィンドウ定義部14は、前段で得たサンプリング信号x(i)を周波数解析するためにウィンドウWを設定し、当該設定されたウィンドウ内のサンプリング信号x(i)∈Wを後段に出力する。すなわち、周波数成分の変化を求める場合には,時間長が異なる複数のデータを定量的に比較できるようにする必要があることからウィンドウを設定する。
本実施形態では、(1)ウィンドウサイズが等しい時間基底、または、(2)拍動数が等しい拍動基底、のいずれかとしてウィンドウを設定する。いずれのウィンドウを用いるかは予め設定しておく。図6は、ウィンドウ定義部14で設定する2通りのウィンドウの例を示す模式図であり、上段部分の交流中間関数P'(t)に対して設定するウィンドウとして、中段部分に時間基底のウィンドウ例Wt1,Wt2,Wt3が示され、下段部分に拍動基底のウィンドウ例Wp1,Wp2,Wp3が示されている。
時間基底ウィンドウは、
Wp1={i|0≦i≦N-1}
Wp2={i|Δn1≦i≦Δn1+N-1}
Wp3={i|Δn1+Δn2≦i≦Δn1+Δn2+N-1}
のように、全て等しいデータ長(サンプリング信号個数)Nを有するが、そのウィンドウ端に、RR間隔の途中で短く途切れてしまう箇所が存在する。例えばウィンドウWp1は、その最後部分が6番目のRR間隔であるRRI6の途中にあり、RRI6の途中で途切れている。(図6の例では示されていないが、時間基底ウィンドウの最後部分と同様に、最初部分でもRR間隔の途中で途切れることが発生しうる。)なお、このように短く途切れた箇所も含めてそのまま周波数解析するとホワイトノイズが生じるが、主に高周波成分に影響するため、本実施形態ではウィンドウ幅をそれなりに大きいものとして設定するため、LF,HF帯への影響は小さいことが想定される。
Wp1={i|0≦i≦N-1}
Wp2={i|Δn1≦i≦Δn1+N-1}
Wp3={i|Δn1+Δn2≦i≦Δn1+Δn2+N-1}
のように、全て等しいデータ長(サンプリング信号個数)Nを有するが、そのウィンドウ端に、RR間隔の途中で短く途切れてしまう箇所が存在する。例えばウィンドウWp1は、その最後部分が6番目のRR間隔であるRRI6の途中にあり、RRI6の途中で途切れている。(図6の例では示されていないが、時間基底ウィンドウの最後部分と同様に、最初部分でもRR間隔の途中で途切れることが発生しうる。)なお、このように短く途切れた箇所も含めてそのまま周波数解析するとホワイトノイズが生じるが、主に高周波成分に影響するため、本実施形態ではウィンドウ幅をそれなりに大きいものとして設定するため、LF,HF帯への影響は小さいことが想定される。
一方で、拍動基底ウィンドウの場合は、当該短く途切れる事態が発生しないように、等しい拍動数で構成されるものとしてウィンドウを生成する、すなわち、連続する所定数のRR間隔の全体に渡る時間範囲をウィンドウとして定める。例えばWt1はRRI1からRRI5までであり、Wt2はRRI2からRRI6までであり、Wt3はRRI3からRRI7までであり、いずれも5拍分(RR間隔5個分)の範囲のウィンドウとしている。時間基底の場合の短い途切れが発生しない代わりに、ウィンドウのデータ長(ウィンドウ内のサンプリング信号の個数)が異なるウィンドウ間で(偶然で一致する場合を除いて一般に)不等長となる。本実施形態において拍動基底ウィンドウを利用する場合は、後述する手法(等値化処理部18)により当該不等長となることに対処する措置を講ずることができる。
<解像度向上部15>
解像度向上部15は、ウィンドウ定義部14で設定されて処理対象となるウィンドウが時間基底ウィンドウまたは拍動基底ウィンドウのいずれの場合であっても共通の処理として、FFTの解像度を向上させるべく、ウィンドウに値が全てゼロで構成されるダミーのサンプリング信号を追加したものを解像度向上ウィンドウとして後段に出力する。
解像度向上部15は、ウィンドウ定義部14で設定されて処理対象となるウィンドウが時間基底ウィンドウまたは拍動基底ウィンドウのいずれの場合であっても共通の処理として、FFTの解像度を向上させるべく、ウィンドウに値が全てゼロで構成されるダミーのサンプリング信号を追加したものを解像度向上ウィンドウとして後段に出力する。
すなわち、現時刻tで処理対象となるウィンドウW(t)がM(t)個のサンプリング信号で構成されているものとして、これに対して所定数D個の値が全てゼロとなるダミーサンプリング信号DSを追加したものを解像度向上ウィンドウWres(t)=W(t)∪DSとする。なお、ダミーサンプリング信号DSを追加する位置は、ウィンドウW(t)の外部(過去側の外部及び/又は未来側の外部のいずれでもよい)に連続して追加すればよい。図7に解像度向上部15による処理の模式例として、ウィンドウW(t)の未来側の外部にダミーサンプリング信号DSの全部を追加した例を示す。(なお、所定数D個は、当該追加した合計数M(t)+Dが一定値となるように定めればよい。図7ではW(t)がM(t)=N=90000個でD=110000個として一定の合計200000個となるようにしており、例えば次のW(t+1)がM(t+1)=80000個であったとすると、D=120000個として、一定の合計200000個となるようにすればよい。)
以下の式(15a)~(15c)は解像度向上部15の当該処理によってFFT解像度が向上することの説明例であり、FFTの解像度Δfは式(15a)の通り、サンプリング速度Fs(サンプリング周波数Fs)をサンプリング点数Nで割った値となる。従って、例えば生データ(すなわち、解像度向上部15の処理を適用する前の、ウィンドウ定義部14から得られたウィンドウ内のサンプリング信号)がFs=1000Hz, 長さ90秒でN=90×1000である場合、その解像度は式(15b)の通り、0.011[Hz]となるが、当該生データに解像度向上部15で例えば110000個のダミーデータを追加することにより、その解像度を式(15c)の通り、0.005[Hz]として約2倍に向上させることができる。
<周波数解析部16及び振幅特性部17>
周波数解析部16は解像度向上ウィンドウ(当該ウィンドウ内のサンプリング信号)に対して、DFT(離散フーリエ変換)を以下の式(16a)のように施して、式(16b)のように、その振幅スペクトルと位相スペクトルによる表現を求め、後段に出力する。振幅特性部17は、この振幅スペクトルと位相スペクトルによる表現から、以下の式(17a)による絶対値の二乗である式(17b)を、当該解像度向上ウィンドウのPSD(f)(パワースペクトル密度)として後段に出力する。
周波数解析部16は解像度向上ウィンドウ(当該ウィンドウ内のサンプリング信号)に対して、DFT(離散フーリエ変換)を以下の式(16a)のように施して、式(16b)のように、その振幅スペクトルと位相スペクトルによる表現を求め、後段に出力する。振幅特性部17は、この振幅スペクトルと位相スペクトルによる表現から、以下の式(17a)による絶対値の二乗である式(17b)を、当該解像度向上ウィンドウのPSD(f)(パワースペクトル密度)として後段に出力する。
<等値化処理部18>
等値化処理部18は、ウィンドウ定義部14で設定されたウィンドウが拍動基底の場合に、前述の通り拍動基底では出力データ長が異なるため、その違いが周波数成分に現れてしまうため、これを等価化するために周波数成分PSD(f)をデータ長Tnで除算して、以下の式(18a)~(18c)のように、等値化された周波数成分PSD'(f)=PSD(f)/Tnとして後段に出力する。当該除算することで、データ長Tn(図7のW(t)のデータ長)が異なるものでも所定の帯域積分値を等しくすること(等値化)が実現できる。なお、当該実現できることの定性的な理由については、図11を参照して後述する。(なお、当該データ長Tnによる除算は、後述する帯域積分部19における帯域積分の式(19a),(19b)や式(19a2),(式19b2)においては、除算項のDatalength(=Tn)として表現されており、式(18a)~(18c)は、拍動基底の場合に1回だけデータ長による除算が必要となることを表すものである。すなわち、等値化処理部18でデータ長Tnにより1回除算して、さらに帯域積分部19でデータ長Datalength(=Tn)による除算をもう1回行って、合計2回のデータ長による除算が行われるわけではない。)
等値化処理部18は、ウィンドウ定義部14で設定されたウィンドウが拍動基底の場合に、前述の通り拍動基底では出力データ長が異なるため、その違いが周波数成分に現れてしまうため、これを等価化するために周波数成分PSD(f)をデータ長Tnで除算して、以下の式(18a)~(18c)のように、等値化された周波数成分PSD'(f)=PSD(f)/Tnとして後段に出力する。当該除算することで、データ長Tn(図7のW(t)のデータ長)が異なるものでも所定の帯域積分値を等しくすること(等値化)が実現できる。なお、当該実現できることの定性的な理由については、図11を参照して後述する。(なお、当該データ長Tnによる除算は、後述する帯域積分部19における帯域積分の式(19a),(19b)や式(19a2),(式19b2)においては、除算項のDatalength(=Tn)として表現されており、式(18a)~(18c)は、拍動基底の場合に1回だけデータ長による除算が必要となることを表すものである。すなわち、等値化処理部18でデータ長Tnにより1回除算して、さらに帯域積分部19でデータ長Datalength(=Tn)による除算をもう1回行って、合計2回のデータ長による除算が行われるわけではない。)
図8に、拍動基底に関して、例示されるようなウィンドウ長T1,T2,T3が互いに異なることにより等値化処理部18の処理を適用することの模式例を示す。
なお、等値化処理部18は、ウィンドウ定義部14で設定されたウィンドウが時間基底の場合には、前述の通りウィンドウサイズは全て同じであって当該等値化処理は不要であるため、等値化処理を施すことをスキップして周波数成分PSD(f)をそのまま後段に出力する。
<帯域積分部19>
帯域積分部19は自律神経指標の元となるLF,HFを、それぞれの周波数帯(前述の通り、LFが0.04Hz≦f≦0.15Hz, HFが0.15Hz≦f≦0.4Hz)でのPSD(f)(等値化されている場合はPSD'(f))の合算値(数値積分)として算出する。図9Aは帯域積分部19によるLF,HFの算出例であり、図中にも示される通り、LF,HFの各帯域での積分値を数値積分による合算値としてPSD(f)(またはPSD'(f))の和で以下の式(19a),(19b)のように計算することができる。なお、Datalengthによる除算は、拍動基底の場合に必要であり、図7のW(t)として説明した時間領域でのデータ長の値を用いればよい。(時間基底の場合はデータ長が同じなので、この除算を省略してよい。)なお、式(19a),(19b)等におけるAN[k]は、式(17a)のA'N[k]と同一である。
帯域積分部19は自律神経指標の元となるLF,HFを、それぞれの周波数帯(前述の通り、LFが0.04Hz≦f≦0.15Hz, HFが0.15Hz≦f≦0.4Hz)でのPSD(f)(等値化されている場合はPSD'(f))の合算値(数値積分)として算出する。図9Aは帯域積分部19によるLF,HFの算出例であり、図中にも示される通り、LF,HFの各帯域での積分値を数値積分による合算値としてPSD(f)(またはPSD'(f))の和で以下の式(19a),(19b)のように計算することができる。なお、Datalengthによる除算は、拍動基底の場合に必要であり、図7のW(t)として説明した時間領域でのデータ長の値を用いればよい。(時間基底の場合はデータ長が同じなので、この除算を省略してよい。)なお、式(19a),(19b)等におけるAN[k]は、式(17a)のA'N[k]と同一である。
本実施形態に係る図9Aや式(19a),(19b)の例では、前記の式(15c)のパラメータ設定に即して周波数解像度を0.005Hzに向上させており、図8の右側に示すように、帯域積分においてLFとHFの境界である0.15Hzを適切に区切り、且つ、LFの下限側境界の0.04HzやHFの上限側境界の0.15Hzも、適切に区切ることが可能となっている。
図9Aとは別の計算例として、例えばW(t)が極端に長く、サイドローブが小さい場合であっても図9B及び式(19a2),(19b2)のように数値積分を算出できる。wkは、数値積分で加算するAN[k](式(17a)のA'N[k]と同一)が対応する積分範囲に入っている割合を表し、全部入っている場合には1であり、半分のみ入っている場合は1/2となる。図9Bにも示される通り、k=40,150,400についてw40=1/2, w150=1/2, w400=1/2であり、当該積分範囲でのその他の全てのkについて、wk =1である。この例も同様に、各境界箇所で1/2となる形で、適切に区切ることが可能となっている。
図10は本実施形態を適用せずにLF,HFの各帯域での積分値を算出する対比例であり、サンプリング周波数Fs=100Hz,データ長24secの設定で周波数解像度Δfを、
Δf=100/(100×60)=0.04167Hz
とするもとで合算値を求めている。この場合、LF帯(0.04-0.15Hz),HF帯(0.15-0.4Hz)に含まれる周波数は、それぞれ3ポイント、4ポイントとなり、解像度が不足して曲線面積の合算値による近似精度が低く、また、各帯域も適切にカバーできておらず、境界箇所で過不足、逸脱、重複が発生してしまっている。
Δf=100/(100×60)=0.04167Hz
とするもとで合算値を求めている。この場合、LF帯(0.04-0.15Hz),HF帯(0.15-0.4Hz)に含まれる周波数は、それぞれ3ポイント、4ポイントとなり、解像度が不足して曲線面積の合算値による近似精度が低く、また、各帯域も適切にカバーできておらず、境界箇所で過不足、逸脱、重複が発生してしまっている。
本発明の実施形態では図9のように、高い解像度である周波数解像度0.005HzのもとでLF帯(0.04-0.15Hz),HF帯(0.15-0.4Hz)に含まれる周波数は、それぞれ21ポイント、39ポイントとなり、合算値による曲線面積の近似精度も高く、且つ、各帯域も適切にカバーできている。
以上、本発明の実施形態によれば、自律神経指標であるLF,HFのリアルタイムな「変化」を導くために必要な(1)単位時間当たりの自律神経系の各周波数成分の定量的な比較をできる。これにより、自律神経指標の変化を表すことができる。また、(2)ウィンドウサイズが短い場合の周波数解像度を向上させることができる。これにより、積分誤差を小さくできる。さらに、(3)定義済みのLF,HFの周波数帯の帯域積分を正確に実施することができる。
以下、種々の補足例、追加例、代替例などについて説明する。
(1) 本発明の実施形態によれば、自律神経指標のリアルタイムでの高精度な把握が可能となることから、例えば精神的なストレスを原因とする疾病等の予防等へと利用することで、国連が主導する持続可能な開発目標(SDGs)の目標3「あらゆる年齢のすべての人々の健康的な生活を確保し、福祉を推進する」に貢献することが可能となる。(なお、データ欠落の一因となる不整脈が偶発的なものではない場合には、別途の医療的な対処が望まれる。)
(2) 等値化処理部18でウィンドウサイズが変化しうる拍動基底の場合に、サイズ長で除算すればよい旨の理由について、図11の評価例で簡単に説明する。データ例D11aに示されるように、周波数成分(単一の交流成分とする)が同一で時間長が異なる中間関数(3個が例示される)についてPSDを比較すると、そのピーク値は時間長の2乗に概ね比例することが数学的解析からも導出できる。そこで、データ例D11bに示すように、PSDを時間長の2乗で除算すれば、ピーク値を互いにに概ね近い値とすることができる。しかしながら、データ例D11cに示すように、ピーク値を揃えたとしても、メインローブの幅が異なると、面積が異なるため、LFやHFの積分値も互いに異なるものとなってしまう。これに対して、データ例D11dに示すように、PSDを時間長で割ることで、面積を互いに等しいものとし、LFやHFの積分値も互いに概ね等しいものとすることができる。等値化処理部18ではこのような数学的性質を利用している。
(3) 中間関数は、式(21a),(21b)で定義されるものを用いる場合について説明したが、心拍データのRR間隔に一致する(横軸の)時間間隔及び/または(縦軸)値をもって心臓の拍動1拍分ごとに値が定まる任意のその他の中間関数を用いるようにしてよい。図12に、その他に利用可能な中間関数の例を示す。(図中Ikはk番目のRR間隔RRIkを表す。)例(a)は単位インパルスによる中間関数であり、間隔Ikをもって現れる心臓の拍動1拍分ごとに値(単位値として1の値)が定まるものである。例(b)はビート関数による中間関数であり、時間間隔はΔkとして等しくし、当該間隔Δkごとに値Ikのインパルスを与える。例(c)はスプライン補間による中間関数であり、間隔Ikをもって現れる心臓の拍動1拍分ごとに値Ikの点を与えたうえで、この点をスプライン補間したものである。例(d)は直線補間による中間関数であり、例(c)でスプライン補間に代えて直線補間を用いたものである。例(e)は1/fゆらぎによる中間関数であり、間隔Ikをもって現れる心臓の拍動1拍分ごとに値1/Ikが連続して現れるようにするものである。例(f)は、式(21a),(21b)で定義した通りのものであり、他の例と対比を可能とすべく図12内にも再掲したものである。
(4) 図13は、一般的なコンピュータ装置70におけるハードウェア構成の例を示す図である。心拍データ解析装置10は、このような構成を有する1台以上のコンピュータ装置70として実現可能である。なお、2台以上のコンピュータ装置70で心拍データ解析装置10を実現する場合、ネットワーク経由で処理に必要な情報の送受を行うようにしてよい。コンピュータ装置70は、所定命令を実行するCPU(中央演算装置)71、CPU71の実行命令の一部又は全部をCPU71に代わって又はCPU71と連携して実行する専用プロセッサとしてのGPU(グラフィックス演算装置)72、CPU71(及びGPU72)にワークエリアを提供する主記憶装置としてのRAM73、補助記憶装置としてのROM74(フラッシュROMとして読み書き可能なSSD等であり、HDD等でもよい)、通信インタフェース75、ディスプレイ76、マウス、キーボード、タッチパネル等によりユーザ入力を受け付ける入力インタフェース77、心拍センサ78と、これらの間でデータを授受するためのバスBSと、を備える。
心拍データ解析装置10の各機能部は、各部の機能に対応する所定のプログラムをROM74から読み込んで実行するCPU71及び/又はGPU72によって実現することができる。プログラム実行に際して必要な参照データは、補助記憶装置としてのROM74でデータベースを構築して、当該データベースから読み込むようにしてよい。なお、CPU71及びGPU72は共に、演算装置(プロセッサ)の一種である。ここで、表示関連の処理が行われる場合にはさらに、ディスプレイ76が連動して動作し、データ送受信に関する通信関連の処理が行われる場合にはさらに通信インタフェース75が連動して動作する。
心拍データ(データ処理されてRR間隔まで加工された状態となる場合も含む)を外部で計測されたものとして取得するのではなく、心拍データ解析装置10自身において取得する場合には、心拍計測機能を提供するセンサとして、任意種類のものを用いてよく、心電信号を取得するECGセンサや脈波信号を取得するPPGセンサや心弾道図信号を取得するBCGセンサ等として構成される心拍センサ78を用いることができる。
10…心拍データ解析装置、11…中間関数生成部、12…直流成分除去部、13…標本化部、14…ウィンドウ定義部、15…解像度向上部、16…周波数解析部、17…振幅特性部、18…等値化処理部、19…帯域積分部
Claims (6)
- 被験者の心拍データから算出されるRR間隔を順次読み込んで、心臓の拍動1拍分ごとに値が定まる中間関数を生成し、直流成分を除去して交流中間関数とする第1処理と、
前記交流中間関数のサンプリング信号に処理対象のウィンドウを設定する第2処理と、
前記ウィンドウ内のサンプリング信号に周波数解析を適用してパワースペクトル密度を算出する第3処理と、
前記パワースペクトル密度のLF帯域とHF帯域のそれぞれでの加算値を自律神経指標として算出する第4処理と、を実行し、
前記第2処理では、サンプリング信号の元となるRR間隔が等しい個数だけ連続した範囲のそれぞれを、拍動基底ウィンドウとして設定する、または、互いに等しい長さの時間基底ウィンドウを設定することを特徴とする心拍データ解析装置。 - 前記第2処理では、サンプリング信号の元となるRR間隔が等しい個数だけ連続した範囲のそれぞれを、拍動基底ウィンドウとして設定することを特徴とする請求項1に記載の心拍データ解析装置。
- 前記第3処理では、パワースペクトル密度を算出したうえで、対応するデータ長さで除算する等値化処理をさらに適用することを特徴とする請求項2に記載の心拍データ解析装置。
- 前記第2処理では、互いに等しい長さの時間基底ウィンドウを設定することを特徴とする請求項1に記載の心拍データ解析装置。
- 前記第2処理では、前記交流中間関数のサンプリング信号に加えて値がゼロとなるダミー信号も追加したものとして、処理対象のウィンドウを設定することを特徴とする請求項1に記載の心拍データ解析装置。
- コンピュータを請求項1ないし5のいずれかに記載の心拍データ解析装置として機能させることを特徴とするプログラム。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2023057915A JP2024145525A (ja) | 2023-03-31 | 2023-03-31 | 心拍データ解析装置及びプログラム |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| JP2023057915A JP2024145525A (ja) | 2023-03-31 | 2023-03-31 | 心拍データ解析装置及びプログラム |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| JP2024145525A true JP2024145525A (ja) | 2024-10-15 |
Family
ID=93057455
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| JP2023057915A Pending JP2024145525A (ja) | 2023-03-31 | 2023-03-31 | 心拍データ解析装置及びプログラム |
Country Status (1)
| Country | Link |
|---|---|
| JP (1) | JP2024145525A (ja) |
Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH0654815A (ja) * | 1992-08-07 | 1994-03-01 | Fukuda Denshi Co Ltd | Rr間隔スペクトル分析方法及びその装置 |
| JPH0788094A (ja) * | 1993-07-30 | 1995-04-04 | Isuzu Motors Ltd | 心拍変動波形解析方法及び装置 |
| JPH10511177A (ja) * | 1994-12-09 | 1998-10-27 | エクソン・ケミカル・パテンツ・インク | パワースペクトル密度監視によるプラントパラメータの検出 |
| JP2006125976A (ja) * | 2004-10-28 | 2006-05-18 | Nsk Ltd | 機械設備の異常診断システム |
| JP2007083065A (ja) * | 1999-03-02 | 2007-04-05 | Quantum Intech Inc | 生理的コヒーレンスおよび自律神経バランスを促進するための方法および装置 |
| JP2016000149A (ja) * | 2014-06-12 | 2016-01-07 | オムロンヘルスケア株式会社 | 脈拍測定装置、脈拍測定方法、および、脈拍測定プログラム |
| JP2019051011A (ja) * | 2017-09-14 | 2019-04-04 | 日本電信電話株式会社 | 瞬時心拍の時系列データの補完装置、補完方法及びそのプログラム |
| JP2021094203A (ja) * | 2019-12-17 | 2021-06-24 | Kddi株式会社 | 心拍変動解析装置、方法およびプログラム |
| JP2023032476A (ja) * | 2021-08-27 | 2023-03-09 | Kddi株式会社 | 心拍データ解析装置及びプログラム |
-
2023
- 2023-03-31 JP JP2023057915A patent/JP2024145525A/ja active Pending
Patent Citations (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH0654815A (ja) * | 1992-08-07 | 1994-03-01 | Fukuda Denshi Co Ltd | Rr間隔スペクトル分析方法及びその装置 |
| JPH0788094A (ja) * | 1993-07-30 | 1995-04-04 | Isuzu Motors Ltd | 心拍変動波形解析方法及び装置 |
| JPH10511177A (ja) * | 1994-12-09 | 1998-10-27 | エクソン・ケミカル・パテンツ・インク | パワースペクトル密度監視によるプラントパラメータの検出 |
| JP2007083065A (ja) * | 1999-03-02 | 2007-04-05 | Quantum Intech Inc | 生理的コヒーレンスおよび自律神経バランスを促進するための方法および装置 |
| JP2006125976A (ja) * | 2004-10-28 | 2006-05-18 | Nsk Ltd | 機械設備の異常診断システム |
| JP2016000149A (ja) * | 2014-06-12 | 2016-01-07 | オムロンヘルスケア株式会社 | 脈拍測定装置、脈拍測定方法、および、脈拍測定プログラム |
| JP2019051011A (ja) * | 2017-09-14 | 2019-04-04 | 日本電信電話株式会社 | 瞬時心拍の時系列データの補完装置、補完方法及びそのプログラム |
| JP2021094203A (ja) * | 2019-12-17 | 2021-06-24 | Kddi株式会社 | 心拍変動解析装置、方法およびプログラム |
| JP2023032476A (ja) * | 2021-08-27 | 2023-03-09 | Kddi株式会社 | 心拍データ解析装置及びプログラム |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110236573B (zh) | 心理压力状态的检测方法及相关装置 | |
| Clifford et al. | Quantifying errors in spectral estimates of HRV due to beat replacement and resampling | |
| EP3033991B1 (en) | System and method for blood pressure estimation | |
| Sahoo et al. | Wavelet based pulse rate and Blood pressure estimation system from ECG and PPG signals | |
| US20150105666A1 (en) | Narrow band feature extraction from cardiac signals | |
| JP2016137038A (ja) | R−r間隔補間方法および心拍変動計測装置 | |
| Wachowiak et al. | Assessing heart rate variability through wavelet-based statistical measures | |
| WO2025066208A1 (zh) | 一种hrv数据预处理方法、装置及电子设备 | |
| US11298065B2 (en) | Fetal heart rate extraction within a processor constrained environment | |
| Cesarelli et al. | PSD modifications of FHRV due to interpolation and CTG storage rate | |
| WO2018130897A1 (en) | Method and system for determining heart rate variability | |
| JP2024145525A (ja) | 心拍データ解析装置及びプログラム | |
| Escrivá Muñoz et al. | Novel characterization method of impedance cardiography signals using time-frequency distributions | |
| Daoud et al. | Use of cardiorespiratory coherence to separate spectral bands of the heart rate variability | |
| JP7221195B2 (ja) | 心拍変動解析装置、方法およびプログラム | |
| JP7460584B2 (ja) | 心拍データ解析装置及びプログラム | |
| Kudrynski et al. | Real-time estimation of the spectral parameters of Heart Rate Variability | |
| EP2938247B1 (en) | Method and apparatus for reducing motion artifacts in ecg signals | |
| Singh et al. | A new baroreflex sensitivity index based on improved Hilbert–Huang transform for assessment of baroreflex in supine and standing postures | |
| Al Osman et al. | A pattern-based windowed impulse rejection filter for nonpathological HRV artifacts correction | |
| JP2000296118A (ja) | 生体信号解析方法および装置 | |
| US10993676B2 (en) | Signal processing method and apparatus | |
| Warmerdam et al. | Reliability of spectral analysis of fetal heart rate variability | |
| Ye et al. | PPG based Respiration Signal Estimation using VMDPCA | |
| JP2024145524A (ja) | 心拍データ解析装置及びプログラム |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| A621 | Written request for application examination |
Free format text: JAPANESE INTERMEDIATE CODE: A621 Effective date: 20250110 |
|
| A977 | Report on retrieval |
Free format text: JAPANESE INTERMEDIATE CODE: A971007 Effective date: 20250822 |
|
| A131 | Notification of reasons for refusal |
Free format text: JAPANESE INTERMEDIATE CODE: A131 Effective date: 20250827 |