ホーム > 特許ランキング > 国立大学法人東北大学
▶ 国立大学法人東北大学の特許一覧 ▶ 株式会社ハートビートサイエンスラボの特許一覧
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2023181010
(43)【公開日】2023-12-21
(54)【発明の名称】拍動解析装置、拍動解析方法、およびプログラム
(51)【国際特許分類】
A61B 5/02 20060101AFI20231214BHJP
【ＦＩ】
A61B5/02 310A
【審査請求】未請求
【請求項の数】6
【出願形態】ＯＬ
(21)【出願番号】P 2022094742
(22)【出願日】2022-06-10
(71)【出願人】
【識別番号】504157024
【氏名又は名称】国立大学法人東北大学
(71)【出願人】
【識別番号】522234297
【氏名又は名称】株式会社ハートビートサイエンスラボ
(74)【代理人】
【識別番号】100165179
【弁理士】
【氏名又は名称】田▲崎▼ 聡
(74)【代理人】
【識別番号】100188558
【弁理士】
【氏名又は名称】飯田 雅人
(74)【代理人】
【識別番号】100175824
【弁理士】
【氏名又は名称】小林 淳一
(74)【代理人】
【識別番号】100152272
【弁理士】
【氏名又は名称】川越 雄一郎
(74)【代理人】
【識別番号】100181722
【弁理士】
【氏名又は名称】春田 洋孝
(72)【発明者】
【氏名】湯田 恵美
(72)【発明者】
【氏名】早野 順一郎
【テーマコード（参考）】
4C017
【Ｆターム（参考）】
4C017AA02
4C017AA10
4C017AA19
4C017AC27
4C017BC14
4C017BC21
4C017FF05
(57)【要約】
【課題】心拍又は脈拍の時系列データからノイズを好適に除去可能な技術を提供する。
【解決手段】拍動解析装置は、心拍又は脈拍の時間間隔である拍間隔の時間に対する拍間隔時系列データを取得する拍間隔時系列データ取得部と、前記拍間隔時系列データ取得部により取得された前記拍間隔時系列データにもとづいて、前記拍間隔時系列データから拍間隔の変化速度、又は、拍間隔の加速度を導出する変動導出部と、前記変動導出部により導出された前記変化速度、又は、前記加速度を複数の区間に分類し、各区間のうち最大頻度を取得する最大頻度取得部と、前記最大頻度取得部により取得された前記最大頻度に応じて、前記間隔時系列データから前記ノイズデータを除去する除去部と、を備えた。
【選択図】図１
【解決手段】拍動解析装置は、心拍又は脈拍の時間間隔である拍間隔の時間に対する拍間隔時系列データを取得する拍間隔時系列データ取得部と、前記拍間隔時系列データ取得部により取得された前記拍間隔時系列データにもとづいて、前記拍間隔時系列データから拍間隔の変化速度、又は、拍間隔の加速度を導出する変動導出部と、前記変動導出部により導出された前記変化速度、又は、前記加速度を複数の区間に分類し、各区間のうち最大頻度を取得する最大頻度取得部と、前記最大頻度取得部により取得された前記最大頻度に応じて、前記間隔時系列データから前記ノイズデータを除去する除去部と、を備えた。
【選択図】図１
【特許請求の範囲】
【請求項1】
心拍又は脈拍の時間間隔である拍間隔の時間に対する拍間隔時系列データを取得する拍間隔時系列データ取得部と、
前記拍間隔時系列データ取得部により取得された前記拍間隔時系列データにもとづいて、前記拍間隔時系列データから拍間隔の変化速度、又は、拍間隔の加速度を導出する変動導出部と、
前記変動導出部により導出された前記変化速度、又は、前記加速度を複数の区間に分類し、各区間のうち最大頻度を取得する最大頻度取得部と、
前記最大頻度取得部により取得された前記最大頻度に応じて、前記拍間隔時系列データからノイズデータを除去する除去部と、
を備えた拍動解析装置。
前記拍間隔時系列データ取得部により取得された前記拍間隔時系列データにもとづいて、前記拍間隔時系列データから拍間隔の変化速度、又は、拍間隔の加速度を導出する変動導出部と、
前記変動導出部により導出された前記変化速度、又は、前記加速度を複数の区間に分類し、各区間のうち最大頻度を取得する最大頻度取得部と、
前記最大頻度取得部により取得された前記最大頻度に応じて、前記拍間隔時系列データからノイズデータを除去する除去部と、
を備えた拍動解析装置。
【請求項2】
前記変動導出部は、前記変化速度を、隣り合う脈拍間隔の差分から導出し、前記加速度を拍間隔の前後の拍間隔から導出する、請求項１に記載の拍動解析装置。
【請求項3】
前記最大頻度取得部は、前記変化速度、又は、前記加速度を複数の区間に分類したヒストグラムを作成する、請求項１に記載の拍動解析装置。
【請求項4】
前記除去部は、前記拍間隔時系列データのデータの総数、前記区間のビン幅、前記最大頻度を用いて除去するノイズを決定する請求項１に記載の拍動解析装置。
【請求項5】
心拍又は脈拍の時間間隔を示す拍間隔時系列データを取得する拍間隔時系列データ取得ステップと、
前記拍間隔時系列データ取得ステップにより取得された前記拍間隔時系列データにもとづいて、前記拍間隔時系列データから拍間隔の変化速度、又は、拍間隔の加速度を導出する変動導出を導出する導出ステップと、
前記導出ステップにより導出された拍間隔の変化速度、又は、拍間隔の加速度を複数の区間に分類し、各区間のうち最大頻度を取得する最大頻度取得ステップと、
前記最大頻度取得ステップにより取得された前記最大頻度に応じて、前記拍間隔時系列データからノイズデータを除去する除去ステップと、
を備えた拍動解析方法。
前記拍間隔時系列データ取得ステップにより取得された前記拍間隔時系列データにもとづいて、前記拍間隔時系列データから拍間隔の変化速度、又は、拍間隔の加速度を導出する変動導出を導出する導出ステップと、
前記導出ステップにより導出された拍間隔の変化速度、又は、拍間隔の加速度を複数の区間に分類し、各区間のうち最大頻度を取得する最大頻度取得ステップと、
前記最大頻度取得ステップにより取得された前記最大頻度に応じて、前記拍間隔時系列データからノイズデータを除去する除去ステップと、
を備えた拍動解析方法。
【請求項6】
コンピュータを、
心拍又は脈拍の時間間隔である拍間隔の時間に対する拍間隔時系列データを取得する拍間隔時系列データ取得部と、
前記拍間隔時系列データ取得部により取得された前記拍間隔時系列データにもとづいて、前記拍間隔時系列データから拍間隔の変化速度、又は、拍間隔の加速度を導出する変動導出部と、
前記変動導出部により導出された前記変化速度、又は、前記加速度を複数の区間に分類し、各区間のうち最大頻度を取得する最大頻度取得部と、
前記最大頻度取得部により取得された前記最大頻度に応じて、前記拍間隔時系列データからノイズデータを除去する除去部として機能させるプログラム。
心拍又は脈拍の時間間隔である拍間隔の時間に対する拍間隔時系列データを取得する拍間隔時系列データ取得部と、
前記拍間隔時系列データ取得部により取得された前記拍間隔時系列データにもとづいて、前記拍間隔時系列データから拍間隔の変化速度、又は、拍間隔の加速度を導出する変動導出部と、
前記変動導出部により導出された前記変化速度、又は、前記加速度を複数の区間に分類し、各区間のうち最大頻度を取得する最大頻度取得部と、
前記最大頻度取得部により取得された前記最大頻度に応じて、前記拍間隔時系列データからノイズデータを除去する除去部として機能させるプログラム。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、拍動解析装置、拍動解析方法、およびプログラムに関する。
【背景技術】
【0002】
光学式脈波センサ等の脈波信号から得られる脈拍間隔又は瞬時脈拍数の時系列データ（以下、拍間隔時系列データという）は、日常生活下における生体状態のモニタリングに用いられている。拍間隔時系列データには、環境光、呼吸、体動などが脈波信号測定に与える影響によって生ずる様々な生理的および物理的なノイズが混入する。特に活動中に記録された拍間隔時系列データは、体動に伴うパルス状の膨大なノイズを含むため、活動中に記録された拍間隔時系列データを用いた生体状態の評価の信頼性は大きく低下することとなる。
【0003】
具体的には、活動中に記録された拍間隔時系列データに含まれるノイズの多くは、１拍から数拍の脈波の検出エラーによる異常に長い脈拍間隔データと、脈波に含まれるノイズを脈波として誤検出することによる異常に短い脈波間隔データとによるノイズが多くを占めている。上述した異常に長い脈拍間隔データおよび異常に短い脈波間隔データは、いずれもパルス状のノイズを生ずる。
これらのノイズと脈波信号と識別する必要がある。
これらのノイズと脈波信号と識別する必要がある。
【0004】
脈波に含まれるノイズの除去には様々な手法がある。最も広く行われているのは、周波数解析により、脈波に関係する周波数帯域のみ通過するバンドパスフィルタ処理によって、脈波信号を抽出するものである。具体的には、例えば、０．７～２．０Ｈｚの周波数帯を通過する帯域フィルタを使用することで、脈拍成分以外の信号成分を除去することが行われている（例えば、特許文献１。）。しかしながら、上記手法は、脈波信号に多量に含まれるパルス状のノイズが広い周波数帯に影響を及ぼすために、真の脈波信号が適切に抽出できない。さらに、本願発明者らは、Ｒ-Ｒ間隔や、脈拍間隔の変動が自律神経機能、呼吸、およびその他の要因によって異なった周波数帯の変調をうける知見を得ており（非特許文献１。）、従来の帯域フィルタによるノイズの除去では効果的なノイズの除去ができないことを経験している。
【0005】
別の手法として、横軸を脈拍間隔、縦軸をその脈波間隔の出現頻度としたヒストグラムを作成し、そのヒストグラムの分布の形状から真の脈拍間隔データの分布を推定し、そこから外れるデータをノイズとして識別する手法が知られている（例えば、特許文献２。）。
【0006】
特許文献２等に記載の手法では、脈拍間隔データの分布を推定することでノイズとの識別を行っているが、非特許文献１に記載と同様に、脈拍間隔は個体間および個体の状態によって大きく変化するために、真の脈拍間隔データの分布そのものも大きな変動を示す。それによって、ノイズによるデータの分布との重なり生ずることも多く、両者の分布の違いから、真の脈拍間隔データとノイズを識別することは困難である。
【0007】
また、特許文献２等に記載の手法では、ノイズとみなして除去する際に、一定以下のサンプル数の場合としているが、実際には、ノイズを生ずる箇所のサンプル数が必ずしも少なくなるわけではないのでノイズでないものを削除する恐れがある。他の類似の手法として、一定間隔以下、又は、以上の脈拍間隔をノイズとして除去する手法もあるが、どの間隔以下、以上をノイズとするかは経験則に頼ることが多い。また、非特許文献１に記載のように、脈拍間隔には個体差および個体内での変動がみられ、固定された閾値を用いる手法では、適切なノイズ除去はできない。
【先行技術文献】
【特許文献】
【0008】
【特許文献1】特開２０１６－１９００２２号公報
【特許文献2】特開２００９－２９７１８４号公報
【非特許文献】
【0009】
【非特許文献1】E. Yuda et al., Journal of Physiological Anthropology(2020)39:21
【発明の概要】
【発明が解決しようとする課題】
【0010】
上記事情に鑑み、本発明は、心拍又は脈拍の拍間隔時系列データからノイズを好適に除去可能な技術を提供することを目的としている。
【課題を解決するための手段】
【0011】
本発明の一態様に係る、拍動解析装置、拍動解析方法又はプログラムは、心拍又は脈拍の時間間隔である拍間隔の時間に対する拍間隔時系列データを取得する拍間隔時系列データ取得部と、前記拍間隔時系列データ取得部により取得された前記拍間隔時系列データにもとづいて、前記拍間隔時系列データから拍間隔の変化速度、又は、拍間隔の加速度を導出する変動導出部と、前記変動導出部により導出された前記変化速度、又は、前記加速度を複数の区間に分類し、各区間のうち最大頻度を取得する最大頻度取得部と、前記最大頻度取得部により取得された前記最大頻度に応じて、前記間隔時系列データから前記ノイズデータを除去する除去部と、を備える。
【発明の効果】
【0012】
本発明により、心拍又は脈拍の時間間隔を示す拍間隔時系列データからノイズを好適に除去可能な技術を提供することが可能となる。
【図面の簡単な説明】
【0013】
【図1】拍動解析装置１００の機能構成を表す機能ブロック図である。
【図2】脈波データの一例を示す図である。
【図3】脈波データを拡大した図である。
【図4】脈波データの一例（基線動揺）を示す図である。
【図5】ｄ（ｉ）のヒストグラムを示す図である。
【図6】Ｒ-Ｒ間隔として測定した時系列データを示す図である。
【図7】ウェアラブルデバイスに設けられた光学式心拍センサで測定された脈波信号により得られた脈波間隔の時系列データを示す図である。
【図8】ノイズデータを除去した時系列データを示す図である。
【図9】拍動解析装置の処理の流れを示すフローチャートである。
【発明を実施するための形態】
【0014】
以下、本発明の具体的な構成例について、図面を参照しながら説明する。
【0015】
〔本発明の拍動解析装置〕
図１は、拍動解析装置１００の機能構成を表す機能ブロック図である。拍動解析装置１００は、バスで接続されたＣＰＵ（Central Processing Unit）やメモリや補助記憶装置などを備え、拍動解析プログラムを実行することによって通信部１１０、拍間隔時系列データ記憶部１４１、および制御部１２０を備える装置として機能する。なお、制御部１２０の各機能の全て又は一部は、ＡＳＩＣ（Application Specific Integrated Circuit）やＰＬＤ（Programmable Logic Device）やＦＰＧＡ（Field Programmable Gate Array）等のハードウェアを用いて実現されてもよい。拍動解析プログラムは、コンピュータ読み取り可能な記録媒体に記録されてもよい。コンピュータ読み取り可能な記録媒体とは、例えばフレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ－ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置である。拍動解析プログラムは、電気通信回線を介して送信されてもよい。
図１は、拍動解析装置１００の機能構成を表す機能ブロック図である。拍動解析装置１００は、バスで接続されたＣＰＵ（Central Processing Unit）やメモリや補助記憶装置などを備え、拍動解析プログラムを実行することによって通信部１１０、拍間隔時系列データ記憶部１４１、および制御部１２０を備える装置として機能する。なお、制御部１２０の各機能の全て又は一部は、ＡＳＩＣ（Application Specific Integrated Circuit）やＰＬＤ（Programmable Logic Device）やＦＰＧＡ（Field Programmable Gate Array）等のハードウェアを用いて実現されてもよい。拍動解析プログラムは、コンピュータ読み取り可能な記録媒体に記録されてもよい。コンピュータ読み取り可能な記録媒体とは、例えばフレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ－ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置である。拍動解析プログラムは、電気通信回線を介して送信されてもよい。
【0016】
拍間隔時系列データ記憶部１４１は、磁気ハードディスク装置や半導体記憶装置等の記憶装置を用いて構成される。拍間隔時系列データ記憶部１４１は、心拍又は脈拍の時間間隔を示す拍間隔時系列データを記憶する。
【0017】
通信部１１０は、他装置と通信を行うためのインタフェースである。通信部１１０は、例えば、ＮＩＣ（Network Interface Card）やＵＳＢ（Universal Serial Interface）インタフェースなどである。
【0018】
図１における制御部１２０は拍動解析装置１００の各部の動作を制御する。制御部１２０は、拍動解析プログラムを実行することによって、拍間隔時系列データ取得部１２１、変動導出部１２２、最大頻度取得部１２３、および除去部１２４として機能する。
【0019】
〔拍間隔時系列データ取得部〕
拍間隔時系列データ取得部１２１は、例えばウェアラブルデバイスなどの他装置から通信部１１０を介して心拍又は脈拍に関係する信号を受信し、受信した信号にもとづいて心拍又は脈拍間隔の時系列データを取得する。拍間隔時系列データ取得部１２１により取得された心拍又は脈拍間隔ｘ（ｉ）（以下では、「拍間隔」という場合もある。）の時系列データは、拍間隔時系列データ記憶部１４１に記憶される。拍間隔時系列データは、本実施形態において(t（ｉ）,ｘ（ｉ）)（t（ｉ）はｘ（ｉ）の時間的位置、ｉは整数）で示す。
拍間隔時系列データ取得部１２１は、例えばウェアラブルデバイスなどの他装置から通信部１１０を介して心拍又は脈拍に関係する信号を受信し、受信した信号にもとづいて心拍又は脈拍間隔の時系列データを取得する。拍間隔時系列データ取得部１２１により取得された心拍又は脈拍間隔ｘ（ｉ）（以下では、「拍間隔」という場合もある。）の時系列データは、拍間隔時系列データ記憶部１４１に記憶される。拍間隔時系列データは、本実施形態において(t（ｉ）,ｘ（ｉ）)（t（ｉ）はｘ（ｉ）の時間的位置、ｉは整数）で示す。
【0020】
以下の説明では、心拍又は脈拍間隔の拍間隔時系列データ(t（ｉ）,ｘ（ｉ）)について説明する。
図２は、脈波データの一例を示す図である。図２は、ノイズのない１分間の脈波を示している。図２において、横軸は時間（０～６０秒）を示し、縦軸は脈波強度（振幅強度）を示す。図３は、図２の脈波データを拡大した図（横軸を０～２０秒とした図）である。
図３に両矢印で示されるピーク間の時間間隔ｘ（ｉ）が脈拍間隔であり、ｘ（ｉ）の開始点の時刻がｔ（ｉ）となる。これを全ての脈拍に対して測定し、横軸を時間、縦軸を拍（脈拍、心拍）間隔とした時系列データが、拍間隔時系列データ（ｔ（ｉ）,ｘ（ｉ））となる。
図２は、脈波データの一例を示す図である。図２は、ノイズのない１分間の脈波を示している。図２において、横軸は時間（０～６０秒）を示し、縦軸は脈波強度（振幅強度）を示す。図３は、図２の脈波データを拡大した図（横軸を０～２０秒とした図）である。
図３に両矢印で示されるピーク間の時間間隔ｘ（ｉ）が脈拍間隔であり、ｘ（ｉ）の開始点の時刻がｔ（ｉ）となる。これを全ての脈拍に対して測定し、横軸を時間、縦軸を拍（脈拍、心拍）間隔とした時系列データが、拍間隔時系列データ（ｔ（ｉ）,ｘ（ｉ））となる。
【0021】
図２，図３では、ノイズがない場合を例に説明したが、実際のデータは、多くのノイズを含んでいる。
ノイズについて具体的に説明する。図４は、基線動揺ノイズを含む脈波データの一例を示す図である。図４において、横軸は時間（０～３２０秒）を示し、縦軸は脈波強度を示す。
ノイズについて具体的に説明する。図４は、基線動揺ノイズを含む脈波データの一例を示す図である。図４において、横軸は時間（０～３２０秒）を示し、縦軸は脈波強度を示す。
【0022】
図４において、脈波は２１０秒付近で上昇し、２８０秒付近で下降したのちに上昇に転じている。このように、基線動揺を含むノイズは、脈波強度の小さくなる部分や、脈拍の時間間隔にずれが生じていることがわかる。
他にも、こうしたノイズの原因としては、低振幅、光学式センサへの外光の侵入、および光学式センサと肌との遊離などが挙げられ、それぞれ、拍間隔データに様々なノイズを引き起こす。
他にも、こうしたノイズの原因としては、低振幅、光学式センサへの外光の侵入、および光学式センサと肌との遊離などが挙げられ、それぞれ、拍間隔データに様々なノイズを引き起こす。
【0023】
〔変動導出部〕
変動導出部１２２は、拍間隔ｘ（ｉ）の拍間隔時系列データ（ｔ（ｉ）,ｘ（ｉ））から拍間隔の変化速度、又は、拍間隔の加速度を導出する。
変動導出部１２２は、拍間隔ｘ（ｉ）の拍間隔時系列データ（ｔ（ｉ）,ｘ（ｉ））から拍間隔の変化速度、又は、拍間隔の加速度を導出する。
【0024】
変化速度の導出方法は、隣り合う間隔の差分を用いた方法である。個々の拍間隔ｘ（ｉ）の直前の脈拍間隔ｘ（ｉ―１）からの変化速度ｄ（ｉ）は、以下の（式１）で表される。
ｄ（ｉ）＝ｘ（ｉ）－ｘ（ｉ－１） （式１）
この方法では、ｘ（ｉ）＝ｘ（ｉ－１）となる場合、すなわち隣り合う間隔が変化しない場合にｄ（ｉ）＝０となる。
ｄ（ｉ）＝ｘ（ｉ）－ｘ（ｉ－１） （式１）
この方法では、ｘ（ｉ）＝ｘ（ｉ－１）となる場合、すなわち隣り合う間隔が変化しない場合にｄ（ｉ）＝０となる。
【0025】
加速度を導出する方法は、個々の拍間隔ｘ（ｉ）の前後の拍間隔の加速度ｄ２（ｉ）を用いた方法である。加速度、ｄ２（ｉ）は以下の（式２）で表される。
ｄ２（ｉ）＝ｘ（ｉ＋１）－２ｘ（ｉ）＋ｘ（ｉ－１） （式２）
この方法では、差分に変化がない場合、換言すると連続する３つの間隔のうちの２つめの間隔が、１つめの間隔と３つめの間隔と同じ場合にｄ（ｉ）＝０となる。
なお、加速度を導出する方法は、変化速度の導出方法で導いた（式１）を時間微分して求めてもよい。
ｄ２（ｉ）＝ｘ（ｉ＋１）－２ｘ（ｉ）＋ｘ（ｉ－１） （式２）
この方法では、差分に変化がない場合、換言すると連続する３つの間隔のうちの２つめの間隔が、１つめの間隔と３つめの間隔と同じ場合にｄ（ｉ）＝０となる。
なお、加速度を導出する方法は、変化速度の導出方法で導いた（式１）を時間微分して求めてもよい。
【0026】
本実施形態では、拍間隔の変化速度および拍間隔の加速度のいずれの方法を用いてもよい。
【0027】
本発明で、このような拍間隔の変化速度又は拍間隔の加速度を用いる理由について説明する。脈波信号のノイズの特徴として、拍間隔そのものは、生理的あるいは測定精度に伴う変動が大きく、その分布の違いから、ノイズによって生ずる異常な間隔と区別することは一般的に困難である。
【0028】
一方、自律神経による心拍数調整の伝達関数は、周波数特性に起因するため、複数の脈拍にわたる変動となって現れ、その変動値には、生理的な限度がある。上記拍間隔の変化速度ｄ（ｉ）又は拍間隔の加速度で得られるｄ２（ｉ）を求めると、年齢や体質、病態による個人差はあるものの、その変動値には生理的な限界があるため所定の範囲内に収まる。
【0029】
それに対して、脈波の検出エラーや誤検出によるノイズは、スパイク様又はパルス様の突発的な変動となるため、非生理的な変動速度や加速度を示す。そのため、変動速度や加速度を求めることにより、生理的な変動と、ノイズによる変動とを区別することが可能となる。
【0030】
本発明は、拍間隔のこの特性を利用し、個人毎に推定した、生理的な拍間隔の変化速度又は加速度の分布から逸脱した速度を、ノイズとして識別することで、個人に即した精度よいノイズの識別が可能になる。
【0031】
〔最大頻度取得部〕
次に、個人に即したノイズの識別手法について述べる。
前述で求めた、拍間隔の変化速度ｄ（ｉ）又は拍間隔の加速度ｄ２（ｉ）について、最大頻度取得部１２３で得られた値のヒストグラムを作成する。
次に、個人に即したノイズの識別手法について述べる。
前述で求めた、拍間隔の変化速度ｄ（ｉ）又は拍間隔の加速度ｄ２（ｉ）について、最大頻度取得部１２３で得られた値のヒストグラムを作成する。
【0032】
最大頻度取得部１２３は、拍間隔の変化速度ｄ（ｉ）又は拍間隔の加速度ｄ２（ｉ）の求まった値を、重なりのない区間ｂ（以下、「ビン（bins）幅ｂ」という。）に分類する。ビン幅の一例として８～２４msが挙げられるが、これに限るものではない。
【0033】
最大頻度取得部１２３は、拍間隔の変化速度ｄ（ｉ）又は拍間隔の加速度ｄ２（ｉ）をヒストグラムにする。確率的に、隣り合う拍の間の速度変動又は加速度変動がない場合の頻度が最も高いため、０付近が最も高い頻度になる。
【0034】
拍間隔の変化速度ｄ（ｉ）又は拍間隔の加速度ｄ２（ｉ）のデータ数（総数）をＮ、得られたヒストグラムの最も高いピークの高さ（最高度数）（以下「最高頻度」という。）をＦとする。なお、最高頻度Ｆになる箇所が複数の区間にあっても問題ない。
【0035】
このようにしてヒストグラム作成すると図５に示すような０付近を中心に左右に広がる分布を示すヒストグラムができる。
このヒストグラムの面積は、ビン幅ｂ×データ数をＮであり、以下の（式３）で表される。
Ｎ・ｂ （式３）
このヒストグラムの分布の中央の最頻度Ｆを含む部分は、生理的な変動による拍の変化（求めるべきピーク（パルス））であり、左右の周辺に広がるヒストグラムのいくつかの頻度は、脈波の検出エラーや誤検出によるスパイク様又はパルス様のノイズによるものと見なせる。
このヒストグラムの面積は、ビン幅ｂ×データ数をＮであり、以下の（式３）で表される。
Ｎ・ｂ （式３）
このヒストグラムの分布の中央の最頻度Ｆを含む部分は、生理的な変動による拍の変化（求めるべきピーク（パルス））であり、左右の周辺に広がるヒストグラムのいくつかの頻度は、脈波の検出エラーや誤検出によるスパイク様又はパルス様のノイズによるものと見なせる。
【0036】
また、このヒストグラムの最高頻度Ｆは、個人の生理的な変動による拍の変化の最も特徴的なものであり個人差が大きい。例えば、若年者は最高頻度Ｆが低くなり、ヒストグラムの頻度分布は広い範囲になる傾向がある。一方で、高齢者や糖尿病の患者は、最高頻度Ｆが高くなり、ヒストグラムの頻度分布は狭い範囲になる傾向がある。
【0037】
つまり、ヒストグラムからは、最高頻度Ｆを中心とする周辺の頻度の情報が、生理的な情報を示し、中心から離れるほど、生理的な限度を超えた値、つまり、ノイズに伴う情報であると言える。一方で、最高頻度Ｆの高さは個人の脈波の傾向を表し、個人により、頻度がどの程度、集中しているのかの傾向を示していると言える。
【0038】
〔除去部〕
そこで、除去部１２４は、最大頻度取得部１２３で作成したヒストグラムから、個人の生理的な変動の分布の推定範囲を決定する。
除去部１２４は、個人の生理的な特徴を示す、最高頻度Ｆの高さに応じて、生理的な変動の分布の範囲（逆に言えば、ノイズを除去する範囲）を決定する。
そこで、除去部１２４は、最大頻度取得部１２３で作成したヒストグラムから、個人の生理的な変動の分布の推定範囲を決定する。
除去部１２４は、個人の生理的な特徴を示す、最高頻度Ｆの高さに応じて、生理的な変動の分布の範囲（逆に言えば、ノイズを除去する範囲）を決定する。
【0039】
一例としては、最高頻度Ｆ周辺の高さを頂点Ａとする三角形とみなし、その範囲とする底辺の長さＢ－Ｃ間を２ＴＩとする。この時、三角形ＡＢＣの面積、以下の（式４）で表される。
ＴＩ・Ｆ （式４）
このとき、本来、得られたデータは、三角形の範囲内に所定の係数αで収まるものと仮定すると、三角形の形成する面積（式４）と、ヒストグラムの面積（式３）が等しくなるから、底辺の長さＴＩは、以下の（式５）で表される。
αＴＩ・＝Ｎ・ｂ／Ｆ （式５）
ＴＩ・Ｆ （式４）
このとき、本来、得られたデータは、三角形の範囲内に所定の係数αで収まるものと仮定すると、三角形の形成する面積（式４）と、ヒストグラムの面積（式３）が等しくなるから、底辺の長さＴＩは、以下の（式５）で表される。
αＴＩ・＝Ｎ・ｂ／Ｆ （式５）
【0040】
図５にはヒストグラムにα＝１のときの、三角形ＡＢＣを重ねて示している。横軸が０における高さＦで、面積がＮ・ｂの三角形で、２ＴＩは２Ｎｂ／Ｆで表される底辺ＢＣの長さである。Ｂの座標は（Ｎｂ／Ｆ，０）であり、Ｃの座標は（－Ｎｂ／Ｆ，０）である。
【0041】
つまり、脈拍の速度変動ｄ（ｉ）の場合は、（式６）の、脈拍間隔の加速度ｄ２（ｉ）の場合は、（式７）の範囲に対応するｘ（ｉ）が、生理的な変動によるものと推定できる。
－αＴＩ＜ｄ（ｉ）＜αＴＩ （式６）
－αＴＩ＜ｄ２（ｉ）＜αＴＩ （式７）
逆に言えば、脈拍の速度変動ｄ（ｉ）の場合は、（式８）の、脈拍間隔の加速度ｄ２（ｉ）の場合は、（式９）の範囲に対応するｘ（ｉ）が、パルス状のノイズによるものと推定できる。
ｄ（ｉ）＜－αＴＩ、αＴＩ＜ｄ（ｉ） （式８）
ｄ２（ｉ）＜－αＴＩ、αＴＩ＜ｄ２（ｉ） （式９）
ここでαは係数であるが、発明者の研究では最適値は１．０前後である。
－αＴＩ＜ｄ（ｉ）＜αＴＩ （式６）
－αＴＩ＜ｄ２（ｉ）＜αＴＩ （式７）
逆に言えば、脈拍の速度変動ｄ（ｉ）の場合は、（式８）の、脈拍間隔の加速度ｄ２（ｉ）の場合は、（式９）の範囲に対応するｘ（ｉ）が、パルス状のノイズによるものと推定できる。
ｄ（ｉ）＜－αＴＩ、αＴＩ＜ｄ（ｉ） （式８）
ｄ２（ｉ）＜－αＴＩ、αＴＩ＜ｄ２（ｉ） （式９）
ここでαは係数であるが、発明者の研究では最適値は１．０前後である。
【0042】
なお、本発明では、生理学的な変動の分布の推定範囲を、脈拍の速度変動ｄ（ｉ）又は脈拍間隔の加速度ｄ２（ｉ）標準偏差を用いていない。
これは、パルス状のノイズのｄ（ｉ）又はｄ２（ｉ）は絶対値の非常に大きい正や負の値を示すことが多いためである。標準偏差は、絶対値の非常に大きい正や負の値に強い影響を受け、精度が悪くなる場合があるためである。例えば、数分にわたり、脈がとれない期間が含まれると、標準偏差が非常に大きくなる。一方、最高頻度Ｆ周辺の高さを頂点とする三角形とみなした本発明のＴＩは、そのような大きいノイズのｄ（ｉ）又はｄ２（ｉ）の値に影響を受けることがなく精度よく、生理的な変動に伴うピーク（心拍又は脈拍）と、ノイズとを分離できる。
これは、パルス状のノイズのｄ（ｉ）又はｄ２（ｉ）は絶対値の非常に大きい正や負の値を示すことが多いためである。標準偏差は、絶対値の非常に大きい正や負の値に強い影響を受け、精度が悪くなる場合があるためである。例えば、数分にわたり、脈がとれない期間が含まれると、標準偏差が非常に大きくなる。一方、最高頻度Ｆ周辺の高さを頂点とする三角形とみなした本発明のＴＩは、そのような大きいノイズのｄ（ｉ）又はｄ２（ｉ）の値に影響を受けることがなく精度よく、生理的な変動に伴うピーク（心拍又は脈拍）と、ノイズとを分離できる。
【0043】
上記例では、最高頻度Ｆ周辺の高さを頂点とする三角形としてＴＩを求めたがこれに限るものではない。
例えば、三角形は図５では、最高頻度Ｆを中心とし、底辺±ＴＩの二等辺三角形としたが、最高頻度Ｆや、頻度分布の状態によっては、頻度分布の割合に応じて、三角形の形状を変更させてもよい。つまり、最高頻度Ｆは頂点でありつつも、必ずしも底辺の中心である必要はない。頻度分布の状態に応じて、三角形の形状は変化可能である。
例えば、三角形は図５では、最高頻度Ｆを中心とし、底辺±ＴＩの二等辺三角形としたが、最高頻度Ｆや、頻度分布の状態によっては、頻度分布の割合に応じて、三角形の形状を変更させてもよい。つまり、最高頻度Ｆは頂点でありつつも、必ずしも底辺の中心である必要はない。頻度分布の状態に応じて、三角形の形状は変化可能である。
【0044】
また、三角形である必要はない。最高頻度Ｆを頂点とし、ヒストグラムのデータが内包できる別の関数ｇ（Ｆ）（＞０）を用いてもよい。その場合、ｇの導関数ｇ'（Ｆ）＜０（Ｆ＞０）となる関数ｇを用いて閾値（±ｇ（Ｆ））を定める。
また、最高頻度Ｆが複数の区間に生じた場合は、頻度分布の状態に応じて、例えば、前後の区間の頻度を合わせたデータ数の平均が来る区間を最高頻度Ｆの区間としてもよい。
なお、ノイズ除去後のデータは、不図示の記憶部に保存、或いは、通信部１１０を介して他の装置に供される。
また、最高頻度Ｆが複数の区間に生じた場合は、頻度分布の状態に応じて、例えば、前後の区間の頻度を合わせたデータ数の平均が来る区間を最高頻度Ｆの区間としてもよい。
なお、ノイズ除去後のデータは、不図示の記憶部に保存、或いは、通信部１１０を介して他の装置に供される。
【0045】
〔本発明の実施例〕
本発明のノイズデータの除去方法により、実際にノイズが除去された場合のデータ例について説明する。
図６は、心電図のＲ波を検出し、その１拍ごとの間隔をＲ-Ｒ間隔として測定した時系列データを示す図である。図６では、横軸は分（min）、縦軸は、Ｒ－Ｒ波間隔（ミリ秒）を示す。図６に示される時系列データは、ノイズのない脈拍間隔の真値の規準を示すＲ－Ｒ波間隔の時系列データである。なお、生理学的にＲ－Ｒ波の間隔と脈拍間隔は、１分毎の平均値は一致するが、周波数の高い変動は一致しないため、図６の時系列データは、完全に一致することはなく、あくまで目安である。
本発明のノイズデータの除去方法により、実際にノイズが除去された場合のデータ例について説明する。
図６は、心電図のＲ波を検出し、その１拍ごとの間隔をＲ-Ｒ間隔として測定した時系列データを示す図である。図６では、横軸は分（min）、縦軸は、Ｒ－Ｒ波間隔（ミリ秒）を示す。図６に示される時系列データは、ノイズのない脈拍間隔の真値の規準を示すＲ－Ｒ波間隔の時系列データである。なお、生理学的にＲ－Ｒ波の間隔と脈拍間隔は、１分毎の平均値は一致するが、周波数の高い変動は一致しないため、図６の時系列データは、完全に一致することはなく、あくまで目安である。
【0046】
図７は、図６において測定されたときと同じ時間帯にウェアラブルデバイス（Silmee W22（登録商標）（TDK株式会社製））に設けられた光学式心拍センサで測定された脈波信号により得られた脈波間隔の時系列データを示す図である。この図７に示される時系列データは、脈拍間隔時系列データの一例である。図７には、複数の脈拍とノイズが含まれている。図８は、図７に示される時系列データから本発明のノイズ除去方法を用いて、ノイズデータを除去した時系列データを示す図である。それぞれの図の横軸は分（min）で、縦軸は、拍間隔（ミリ秒）を示す。
【0047】
図７に示される間隔時系列データは、図６に示される時系列データと同様に、８００ミリ秒が最頻値にみえる。しかし、図６にはみられない１０００ミリ秒に到達しているデータが非常に多いことがわかる。すなわち、ノイズデータが非常に多いと考えられる。
【0048】
一方、図８に示される時系列データは、１０００ミリ秒に到達しているデータは除去され、図６に示されるＲ－Ｒ波間隔と比較して、若干振幅の大きいデータがあるものの、概ねＲ－Ｒ波間隔と同様のデータとなっている。すなわち、ノイズが好適に除去された時系列データとなっている。なお、本実施形態によるノイズデータの除去方法によれば、期外収縮やブロックによる脈波間隔もノイズデータとして除去されるが、呼吸により生じる変動は、その変動によるｄ（ｉ）が区間［－ＴＩ，ＴＩ］に含まれるため、ノイズデータとして除去されない。
【0049】
〔本発明の処理ステップ〕
次に、上述した拍動解析装置１００の処理の流れを、フローチャートを用いて説明する。図９は、拍動解析装置１００の処理の流れを示すフローチャートである。なお、以下は、変動速度ｄ（ｉ）で説明するが、加速度ｄ２（ｉ）も同様に実行できる。図９において、拍間隔時系列データ取得部１２１は、拍間隔時系列データ（ｔ（ｉ）,ｘ（ｉ））を取得する（ステップＳ１０１）。変動導出部１２２は、拍間隔時系列データ（ｔ（ｉ）,ｘ（ｉ））にもとづいて、ｄ（ｉ）を導出する（ステップＳ１０２）。図９に示される処理において、ｄ（ｉ）の総数は、Ｎ個（ｄ（１）～ｄ（Ｎ））とする。
次に、上述した拍動解析装置１００の処理の流れを、フローチャートを用いて説明する。図９は、拍動解析装置１００の処理の流れを示すフローチャートである。なお、以下は、変動速度ｄ（ｉ）で説明するが、加速度ｄ２（ｉ）も同様に実行できる。図９において、拍間隔時系列データ取得部１２１は、拍間隔時系列データ（ｔ（ｉ）,ｘ（ｉ））を取得する（ステップＳ１０１）。変動導出部１２２は、拍間隔時系列データ（ｔ（ｉ）,ｘ（ｉ））にもとづいて、ｄ（ｉ）を導出する（ステップＳ１０２）。図９に示される処理において、ｄ（ｉ）の総数は、Ｎ個（ｄ（１）～ｄ（Ｎ））とする。
【0050】
最大頻度取得部１２３は、ビン幅ｂを取得して（ステップＳ１０３）、ヒストグラムを作成し、最大頻度Ｆを取得する（ステップＳ１０４）。ビン幅ｂは、デフォルトで設定され、記憶装置に予め記憶されているか、拍動解析装置１００の操作者により設定可能としてもよい。最大頻度取得部１２３は、作成したヒストグラムに対し、最大頻度Ｆを頂点とする、最も合理的に内包する形状（関数）、係数αを選択して、ＴＩを決定する。形状で三角形を選択し、αを１とした場合、ＴI＝Ｎｂ／Ｆを算出する。除去部１２４は、ＴI＝Ｎｂ／Ｆを取得する（ステップＳ１０５）。
【0051】
除去部１２４は、ループカウンタｋを１で初期化する（ステップＳ１０６）。除去部１２４は、ｄ（ｋ）＜－ＴＩ又はＴＩ＜ｄ（ｋ）か否かを判定する（ステップＳ１０７）。－ＴＩ≦ｄ（ｋ）≦ＴＩの場合には（ステップＳ１０７：ＮＯ）、除去部１２４は、ステップＳ１０９に進む。ｄ（ｋ）＜－ＴＩ又はＴＩ＜ｄ（ｋ）の場合には（ステップＳ１０７：ＹＥＳ）、除去部１２４は、ｘ（ｋ）をノイズデータとして除去する（ステップＳ１０８）。除去部１２４は、ループカウンタｋを１だけ増分し（ステップＳ１０９）、ｋ＞Ｎか否かを判定する（ステップＳ１１０）。ｋ＞Ｎの場合には（ステップＳ１１０：ＹＥＳ）、除去部１２４は、処理を終了する。ｋ≦Ｎの場合には（ステップＳ１１０：ＮＯ）、除去部１２４は、ステップＳ１０７に戻る。
【0052】
〔その他〕
以上説明した実施形態において、拍間隔の測定法はどのようなものであってもよい。また、本実施形態は、心拍間隔、脈波間隔に限らず、以下の時系列データ（（１）～（６））にも適用できる。
（１）１拍毎の心拍間隔又は脈拍間隔の時系列データ
（２）１拍毎の瞬時心拍数又は瞬時脈拍数の時系列データ
（３）一定間隔毎の瞬時又は平均心拍間隔の時系列データ
（４）一定間隔毎の瞬時又は平均脈拍間隔の時系列データ
（５）一定間隔毎の瞬時又は平均心拍数の時系列データ
（６）一定間隔毎の瞬時又は平均脈拍数の時系列データ
上記（１）～（６）に示される時系列データを本実施形態における（ｔ（ｉ）,ｘ（ｉ））とし、第１の方法又は第２の方法によってｄ（ｉ）又はｄ２（ｉ）を導出することで、本実施形態を適用することができる。
以上説明した実施形態において、拍間隔の測定法はどのようなものであってもよい。また、本実施形態は、心拍間隔、脈波間隔に限らず、以下の時系列データ（（１）～（６））にも適用できる。
（１）１拍毎の心拍間隔又は脈拍間隔の時系列データ
（２）１拍毎の瞬時心拍数又は瞬時脈拍数の時系列データ
（３）一定間隔毎の瞬時又は平均心拍間隔の時系列データ
（４）一定間隔毎の瞬時又は平均脈拍間隔の時系列データ
（５）一定間隔毎の瞬時又は平均心拍数の時系列データ
（６）一定間隔毎の瞬時又は平均脈拍数の時系列データ
上記（１）～（６）に示される時系列データを本実施形態における（ｔ（ｉ）,ｘ（ｉ））とし、第１の方法又は第２の方法によってｄ（ｉ）又はｄ２（ｉ）を導出することで、本実施形態を適用することができる。
【0053】
また、本発明では、ヒストグラムを作成しているが、ヒストグラムを作成する必要はなく、単に複数の区間に分類した頻度を求め、最大頻度を求め、計算処理により、ノイズを除去する範囲を決定してもよい。
【0054】
上述した本実施形態に係る拍動解析装置１００の機能をウェアラブルデバイスで実現するようにしてもよい。この場合のウェアラブルデバイスは、脈波を検出するセンサと拍動解析装置１００の機能を備える。また、拍動解析装置１００の機能を、ウェアラブルデバイスから間隔時系列データを取得するスマートフォンやサーバで実現してもよい。
【0055】
以上、説明したように、本実施形態によれば、心拍又は脈拍の時間間隔を示す間隔時系列データからノイズを好適に除去可能な技術を提供することができる。
【0056】
上述した実施形態における拍動解析装置１００の機能をコンピュータで実現するようにしても良い。その場合、この機能を実現するためのプログラムをコンピュータ読み取り可能な記録媒体に記録して、この記録媒体に記録されたプログラムをコンピュータシステムに読み込ませ、実行することによって実現しても良い。なお、ここでいう「コンピュータシステム」とは、ＯＳや周辺機器等のハードウェアを含むものとする。また、「コンピュータ読み取り可能な記録媒体」とは、フレキシブルディスク、光磁気ディスク、ＲＯＭ、ＣＤ－ＲＯＭ等の可搬媒体、コンピュータシステムに内蔵されるハードディスク等の記憶装置のことをいう。さらに「コンピュータ読み取り可能な記録媒体」とは、インターネット等のネットワークや電話回線等の通信回線を介してプログラムを送信する場合の通信線のように、短時間の間、動的にプログラムを保持するもの、その場合のサーバやクライアントとなるコンピュータシステム内部の揮発性メモリのように、一定時間プログラムを保持しているものも含んでも良い。また上記プログラムは、前述した機能の一部を実現するためのものであっても良く、さらに前述した機能をコンピュータシステムにすでに記録されているプログラムとの組み合わせで実現できるものであっても良い。
【0057】
本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。
【符号の説明】
【0058】
１００ 拍動解析装置、１２０ 制御部、１２１ 拍間隔時系列データ取得部、１２２ 変動導出部、１２３ 最大頻度取得部、１２４ 除去部、１４１ 拍間隔時系列データ記憶部
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】