特許 7106821 生体信号解析装置、生体信号解析方法及びプログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2022-07-19

(45)【発行日】2022-07-27

(54)【発明の名称】生体信号解析装置、生体信号解析方法及びプログラム

(51)【国際特許分類】

   A61B 5/02 20060101AFI20220720BHJP

【ＦＩ】

A61B5/02 310Z

【請求項の数】 8

(21)【出願番号】P 2017137302

(22)【出願日】2017-07-13

(65)【公開番号】P2019017555

(43)【公開日】2019-02-07

【審査請求日】2020-02-26

(73)【特許権者】

【識別番号】000006747

【氏名又は名称】株式会社リコー

(74)【代理人】

【識別番号】100107766

【弁理士】

【氏名又は名称】伊東 忠重

(74)【代理人】

【識別番号】100070150

【弁理士】

【氏名又は名称】伊東 忠彦

(72)【発明者】

【氏名】船橋 一樹

(72)【発明者】

【氏名】吉井 雅子

【審査官】山口 裕之

(56)【参考文献】

【文献】特開２０１７－１０４４８８（ＪＰ，Ａ）

【文献】米国特許出願公開第２０１７／００３５３０４（ＵＳ，Ａ１）

【文献】米国特許出願公開第２０１４／００８１０８８（ＵＳ，Ａ１）

【文献】特開２０１６－００２１１９（ＪＰ，Ａ）

【文献】特開２０１６－０４３１９１（ＪＰ，Ａ）

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ａ６１Ｂ ５／０２

(57)【特許請求の範囲】

【請求項1】

周期性を有し、かつ、周期に揺らぎのある脈波信号を解析する生体信号解析装置であって、
前記脈波信号を示す第１信号を生成する第１信号生成部と、
前記第１信号の周波数スペクトルを算出する算出部と、
算出した周波数スペクトルに基づいて、前記第１信号における生体起因の周期の揺らぎのうち、０．１５Ｈｚ以下の第１周波数帯における揺らぎ成分である第１揺らぎ成分を、前記第１信号のピーク周波数帯域とともに抽出する第１抽出部と、
前記第１揺らぎ成分が抽出された信号の脈拍間隔を該脈拍間隔の所定の固定値に合うように前記第１信号のサンプリング周期を補正することにより前記第１揺らぎ成分を平坦化することで０．１５Ｈｚより高く、０．４０Ｈｚ以下の第２周波数帯における生体起因の周期の揺らぎ成分である第２揺らぎ成分を示す信号値を抽出し、前記第２揺らぎ成分に起因するピークが強調された第２信号を生成する第２信号生成部と、
前記第２信号が有する振幅のノイズ成分を除去するノイズ除去部と、
前記ノイズ除去部によりノイズ成分が除去された信号に対して、第２信号生成部で平坦化された第１揺らぎ成分を付与する非平坦化処理部と、
前記非平坦化処理部により第１揺らぎ成分が付与された信号から、少なくとも前記第２揺らぎ成分を抽出する第２抽出部と、
を備える生体信号解析装置。

【請求項2】

生体を撮像した画像データに基づいて、前記第１信号生成部は、前記第１信号を生成することを特徴とする
請求項１に記載の生体信号解析装置。

【請求項3】

生体に対して非接触であることを特徴とする
請求項１又は２に記載の生体信号解析装置。

【請求項4】

前記ノイズ除去部は、
前記第２信号に含まれる周波数成分を解析して、設定周波数帯以外の周波数帯に含まれる成分を除去し、
前記設定周波数帯は、
前記第２信号の基本波及び前記第１信号における周期の揺らぎのうち、前記第２周波数帯における揺らぎ成分である第２揺らぎ成分に起因するピークを、
少なくとも含むことを特徴とする
請求項１乃至３のいずれか１項に記載の生体信号解析装置。

【請求項5】

前記設定周波数帯は、
前記基本波より０．１５Ｈｚ乃至０．４０Ｈｚ低い周波数帯域から最も振幅が大きい第２波と、前記基本波を中心に前記第２波と対称となる周波数帯域から第３波を含むことを特徴とする
請求項４に記載の生体信号解析装置。

【請求項6】

前記設定周波数帯は、
高調波に起因するピークを含むことを特徴とする
請求項４又は５に記載の生体信号解析装置。

【請求項7】

周期性を有し、かつ、周期に揺らぎのある脈波信号を解析する生体信号解析装置が行う生体信号解析方法であって、
生体信号解析装置が、前記脈波信号を示す第１信号を生成する第１信号生成手順と、
生体信号解析装置が、前記第１信号の周波数スペクトルを算出する算出手順と、
生体信号解析装置が、算出した周波数スペクトルに基づいて、前記第１信号における生体起因の周期の揺らぎのうち、０．１５Ｈｚ以下の第１周波数帯における揺らぎ成分である第１揺らぎ成分を、前記第１信号のピーク周波数帯域とともに抽出する第１抽出手順と、
生体信号解析装置が、前記第１揺らぎ成分が抽出された信号の脈拍間隔を該脈拍間隔の所定の固定値に合うように前記第１信号のサンプリング周期を補正することにより前記第１揺らぎ成分を平坦化することで０．１５Ｈｚより高く、０．４０Ｈｚ以下の第２周波数帯における生体起因の周期の揺らぎ成分である第２揺らぎ成分を示す信号値を抽出し、前記第２揺らぎ成分に起因するピークが強調された第２信号を生成する第２信号生成手順と、
生体信号解析装置が、前記第２信号が有する振幅のノイズ成分を除去するノイズ除去手順と、
生体信号解析装置が、ノイズ成分が除去された信号に対して、第２信号生成手順で平坦化された第１揺らぎ成分を付与する非平坦化処理手順と、
前記非平坦化処理手順により第１揺らぎ成分が付与された信号から、少なくとも前記第２揺らぎ成分を抽出する第２抽出手順と、
を含む生体信号解析方法。

【請求項8】

周期性を有し、かつ、周期に揺らぎのある脈波信号を解析するコンピュータに生体信号解析方法を実行させるプログラムであって、
コンピュータが、前記脈波信号を示す第１信号を生成する第１信号生成手順と、
コンピュータが、前記第１信号の周波数スペクトルを算出する算出手順と、
コンピュータが、算出した周波数スペクトルに基づいて、前記第１信号における生体起因の周期の揺らぎのうち、０．１５Ｈｚ以下の第１周波数帯における揺らぎ成分である第１揺らぎ成分を、前記第１信号のピーク周波数帯域とともに抽出する第１抽出手順と、
コンピュータが、前記第１揺らぎ成分が抽出された信号の脈拍間隔を該脈拍間隔の所定の固定値に合うように前記第１信号のサンプリング周期を補正することにより前記第１揺らぎ成分を平坦化することで０．１５Ｈｚより高く、０．４０Ｈｚ以下の第２周波数帯における生体起因の周期の揺らぎ成分である第２揺らぎ成分を示す信号値を抽出し、前記第２揺らぎ成分に起因するピークが強調された第２信号を生成する第２信号生成手順と、
コンピュータが、前記第２信号が有する振幅のノイズ成分を除去するノイズ除去手順と、
コンピュータが、ノイズ成分が除去された信号に対して、第２信号生成手順で平坦化された第１揺らぎ成分を付与する非平坦化処理手順と、
前記非平坦化処理手順により第１揺らぎ成分が付与された信号から、少なくとも前記第２揺らぎ成分を抽出する第２抽出手順と、
を実行させるプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、生体信号解析装置、生体信号解析方法及びプログラムに関するものである。

【背景技術】

【0002】

従来、生体から取得できる脈波等を解析して、生体の様々な状態を示す指標を評価する技術が知られている。例えば、指標は、脈拍数又は脈拍変動指標等である。

【0003】

また、脈拍変動指標は、脈拍のピーク間隔の揺らぎを評価する指標である。具体的には、脈拍変動指標には、０．０４Ｈｚ（ヘルツ）乃至０．１５Ｈｚ程度となる低周波成分の揺らぎを示すＬＦ（Ｌｏｗ Ｆｒｅｑｕｅｎｃｙ）値、０．１５Ｈｚ乃至０．４０Ｈｚ程度となる高周波成分の揺らぎを示すＨＦ（Ｈｉｇｈ Ｆｒｅｑｕｅｎｃｙ）値及びＬＦ値とＨＦ値の比であるＬＦ／ＨＦ値等の指標がある。そして、脈拍変動指標は、生体の自律神経の働きと関連するため、脈拍変動指標から、生体の自律神経の状態が評価できる。例えば、生体の自律神経の状態から、生体の疲れ又は病気等が把握できる。

【0004】

そして、まず、装置が、生体信号におけるピーク間隔の検出サンプル数を区間ごとにリサンプリングする。次に、装置が、リサンプリングされた信号を直交変換する。さらに、装置が、基本周波数成分及び高調波成分を抽出することで、ノイズとなる信号を確実に除去する方法が知られている（例えば、特許文献１等）。

【発明の概要】

【発明が解決しようとする課題】

【0005】

しかしながら、従来の技術では、高周波成分の揺らぎが精度良く抽出できない課題がある。

【0006】

本発明の１つの側面は、高周波成分の揺らぎを精度良く抽出することを目的としている。

【課題を解決するための手段】

【0007】

上述した課題を解決するために、本発明の一態様である、周期性を有し、かつ、周期に揺らぎのある脈波信号を解析する生体信号解析装置は、
前記脈波信号を示す第１信号を生成する第１信号生成部と、
前記第１信号の周波数スペクトルを算出する算出部と、
算出した周波数スペクトルに基づいて、前記第１信号における生体起因の周期の揺らぎのうち、０．１５Ｈｚ以下の第１周波数帯における揺らぎ成分である第１の揺らぎ成分を、前記第１信号のピーク周波数帯域とともに抽出する第１抽出部と、
前記第１揺らぎ成分が抽出された信号の脈拍間隔を該脈拍間隔の所定の固定値に合うように前記第１信号のサンプリング周期を補正することにより前記第１の揺らぎ成分を平坦化することで０．１５Ｈｚより高く、０．４０Ｈｚ以下の第２周波数帯における生体起因の周期の揺らぎ成分である第２の周期の揺らぎ成分を示す信号値を抽出し、前記第２揺らぎ成分に起因するピークが強調された第２信号を生成する第２信号生成部と、
前記第２信号が有する振幅のノイズ成分を除去するノイズ除去部と、
前記ノイズ除去部によりノイズ成分が除去された信号に対して、第２信号生成部で平坦化された第１揺らぎ成分を付与する非平坦化処理部と、
前記非平坦化処理部により第１揺らぎ成分が付与された信号から、少なくとも前記第２揺らぎ成分を抽出する第２抽出部と、
を備えることを特徴とする。

【発明の効果】

【0008】

高周波成分の揺らぎが精度良く抽出できる。

【図面の簡単な説明】

【0009】

【図1】本発明の一実施形態に係る生体信号解析装置の一例を示す概略図である。

【図2】本発明の一実施形態に係る生体信号解析装置のハードウェア構成の一例を示すブロック図である。

【図3】本発明の一実施形態に係る生体信号解析装置の機能構成の一例を示す機能ブロック図である。

【図4】本発明の一実施形態に係る生体信号解析装置による全体処理の一例を示すフローチャートである。

【図5】本発明の一実施形態に係る生体信号解析装置による撮像例を示す図である。

【図6】本発明の一実施形態に係る生体信号解析装置による脈波信号の抽出例を示すフローチャートである。

【図7】本発明の一実施形態に係る脈波信号の一例を示す図である。

【図8】本発明の一実施形態に係るＳＮ比を計算する処理の一例を示すフローチャートである。

【図9】本発明の一実施形態に係るＳＮ比の計算において算出されるパワースペクトルの一例を示す図である。

【図10】本発明の一実施形態に係る閾値の設定例を示す図である。

【図11】本発明の一実施形態に係る生体信号解析装置による制御の一例を示すタイミングチャートである。

【図12】本発明の一実施形態に係る生体信号解析装置による脈波解析の一例を示すフローチャートである。

【図13】本発明の一実施形態に係る生体信号解析装置による脈波変動指標の算出例を示すフローチャートである。

【図14】本発明の一実施形態に係る生体信号解析装置による第１信号の周波数解析結果例を示す図である。

【図15】本発明の一実施形態における高周波成分の揺らぎ及び低周波成分の揺らぎの例を示す図である。

【図16】本発明の一実施形態に係る生体信号解析装置による脈拍間隔の検出例を示す図である。

【図17】本発明の一実施形態に係る生体信号解析装置によるピーク時間の補正例を示す図である。

【図18】本発明の一実施形態に係る平坦化の例を示す図である。

【図19】本発明の一実施形態に係る第２信号の一例を示す図である。

【図20】本発明の一実施形態に係る出力信号の一例を示す図である。

【図21】本発明の一実施形態に係る効果の一例を示す図である。

【発明を実施するための形態】

【0010】

以下、本発明の実施形態について添付の図面を参照しながら説明する。なお、本明細書及び図面において、実質的に同一の機能構成を有する構成要素については、同一の符号を付し、重複する説明を省略する。

【0011】

＜生体信号解析装置例＞
図１は、本発明の一実施形態に係る生体信号解析装置の一例を示す概略図である。

【0012】

本実施形態では、まず、生体信号解析装置の例である端末装置１００は、カメラ等の撮像装置１００Ｈ１を有する。そして、端末装置１００は、撮像装置１００Ｈ１によって、脈波等を測定する対象の生体１を撮像する。次に、撮像装置１００Ｈ１が撮像した画像に基づいて、端末装置１００は、生体１の脈波等を示す生体信号を取得する。

【0013】

＜ハードウェア構成例＞
図２は、本発明の一実施形態に係る生体信号解析装置のハードウェア構成の一例を示すブロック図である。図示するように、例えば、端末装置１００は、撮像装置１００Ｈ１と、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）１００Ｈ２と、記憶装置１００Ｈ３と、入力装置１００Ｈ４とを有する。さらに、端末装置１００は、出力装置１００Ｈ５と、Ｉ／Ｆ（ｉｎｔｅｒｆａｃｅ）１００Ｈ６とを有する。これらのハードウェアは、バス１００Ｈ７によって接続される。

【0014】

撮像装置１００Ｈ１は、例えば、カメラ、光センサ又はこれらの組み合わせである。以下、撮像装置１００Ｈ１がカメラである例で説明する。また、この例では、カメラの色フィルタ構成は、例えば、Ｒ（Ｒｅｄ）、Ｇ（Ｇｒｅｅｎ）及びＢ（Ｂｌｕｅ）（以下「Ｒ」、「Ｇ」及び「Ｂ」で各色を示す場合がある。）の３チャンネルの色信号を出力できるベイヤ（Ｂａｙｅｒ）構成等である。なお、カメラの色フィルタ構成は、１チャンネル以上の色信号を出力できる構成であればよい。

【0015】

また、カメラは、脈拍に基づく輝度変化を取得できるチャンネルを有するのが望ましい。脈拍に基づく輝度変化は、例えば、「Ｇ」又は近赤外（ＮＩＲ、Ｎｅａｒ－ｉｎｆｒａｒｅｄ）の光から取得しやすい。したがって、カメラは、近赤外のチャンネルを有してもよい。近赤外の光は、主に７５０ｎｍ（ナノメートル）乃至１．４μｍ（マイクロメートル）程度の波長である。

【0016】

なお、撮像装置１００Ｈ１が有する撮像素子は、例えば、ＣＣＤ（Ｃｈａｒｇｅ Ｃｏｕｐｌｅｄ Ｄｅｖｉｃｅ）又はＣＭＯＳ（Ｃｏｍｐｌｅｍｅｎｔａｒｙ Ｍｅｔａｌ Ｏｘｉｄｅ Ｓｅｍｉｃｏｎｄｕｃｔｏｒ）センサ等の光センサである。

【0017】

さらに、撮像装置１００Ｈ１は、カメラ等に限られない。例えば、撮像装置は、光を測定する光センサ等でもよい。すなわち、光センサは、生体から反射する光を測定できるセンサ等である。

【0018】

以下、撮像装置１００Ｈ１がカメラであって、ＲＧＢ画像を出力する例で説明する。

【0019】

ＣＰＵ１００Ｈ２は、中央処理装置である。すなわち、ＣＰＵ１００Ｈ２は、処理を実現するための演算及びデータの加工等を行う演算装置並びにハードウェアを制御する制御装置である。また、ＣＰＵ１００Ｈ２は、記憶装置１００Ｈ３等が記憶するプログラム等に基づいて、処理を実行する。

【0020】

記憶装置１００Ｈ３は、例えば、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）、ＲＯＭ（Ｒｅａｄ－Ｏｎｌｙ Ｍｅｍｏｒｙ）、ハードディスク又はこれらの組み合わせである。すなわち、記憶装置１００Ｈ３は、主記憶装置及び補助記憶装置等である。

【0021】

入力装置１００Ｈ４は、ユーザによる操作を入力する装置である。例えば、入力装置１００Ｈ４は、キーボード、マウス又はこれらの組み合わせ等である。

【0022】

出力装置１００Ｈ５は、ユーザに対して画像、各種処理結果又はこれらの組み合わせ等を表示する装置である。例えば、出力装置１００Ｈ５は、液晶ディスプレイ等である。

【0023】

なお、入力装置１００Ｈ４及び出力装置１００Ｈ５は、一体となっているタッチパネル等でもよい。

【0024】

Ｉ／Ｆ１００Ｈ６は、外部装置と接続するためのインタフェースである。例えば、Ｉ／Ｆ１００Ｈ６は、ＵＳＢ（Ｕｎｉｖｅｒｓａｌ Ｓｅｒｉａｌ Ｂｕｓ）等である。また、Ｉ／Ｆ１００Ｈ６は、画像データ等を外部装置と入出力する。さらに、Ｉ／Ｆ１００Ｈ６は、ネットワーク等を介して外部装置と通信を行って、データを送受信してもよい。

【0025】

また、Ｉ／Ｆ１００Ｈ６は、記録媒体２００等からプログラム等を入力する。

【0026】

このように、端末装置１００は、スマートフォン又はＰＣ（Ｐｅｒｓｏｎａｌ Ｃｏｍｐｕｔｅｒ）等の情報処理装置であり、コンピュータである。なお、端末装置１００は、情報処理装置と、情報処理装置に接続される撮像装置との組み合わせ等でもよい。

【0027】

＜機能構成例＞
図３は、本発明の一実施形態に係る生体信号解析装置の機能構成の一例を示す機能ブロック図である。図示するように、端末装置１００は、第１信号生成部１００Ｆ１と、抽出部１００Ｆ２と、第２信号生成部１００Ｆ３と、ノイズ除去部１００Ｆ４とを備える。

【0028】

第１信号生成部１００Ｆ１は、撮像装置等の測定装置により取得される生体１の画像データ等のデータに基づいて、生体信号を示す第１信号を生成する第１信号生成手順を行う。

【0029】

抽出部１００Ｆ２は、第１信号生成部１００Ｆ１が生成する第１信号が有する低周波揺らぎ成分を抽出する抽出手順を行う。

【0030】

第２信号生成部１００Ｆ３は、抽出部１００Ｆ２が抽出する低周波揺らぎ成分を平坦化させて第２信号を生成する第２信号生成手順を行う。

【0031】

ノイズ除去部１００Ｆ４は、第２信号生成部１００Ｆ３が生成する第２信号が有するノイズ成分を除去するノイズ除去手順を行う。

【0032】

例えば、第１信号生成部１００Ｆ１と、抽出部１００Ｆ２と、第２信号生成部１００Ｆ３と、ノイズ除去部１００Ｆ４とは、ＣＰＵ１００Ｈ２等によって実現される。

【0033】

なお、端末装置１００は、更に指標を計算する機能等を備えてもよい。

【0034】

＜全体処理例＞
図４は、本発明の一実施形態に係る生体信号解析装置による全体処理の一例を示すフローチャートである。図示するように、端末装置は、第１信号を生成するのに、ステップＳ０３乃至ステップＳ０５のように、撮像条件を調整して生体信号等を精度良く抽出できるようにするのが望ましい。以下、図示するように、撮像条件の変更等を行う場合を例に説明する。なお、図示するステップＳ０１乃至ステップＳ０２は、第１信号の生成例であり、第１信号は、図示する以外の処理によって生成されてもよい。

【0035】

＜生体の撮像例＞（ステップＳ０１）
ステップＳ０１では、端末装置は、生体を撮像する。例えば、ステップＳ０１では、端末装置は、生体を撮像し、生体を示す画像データを生成する。例えば、端末装置は、３０ｆｐｓ（フレーム毎秒）程度で撮像を行い、動画データを生成する。

【0036】

図５は、本発明の一実施形態に係る生体信号解析装置による撮像例を示す図である。端末装置は、生体１の部位のうち、図示する領域１１が中心となるように撮像するのが望ましい。具体的には、領域１１は、生体１の鼻及び頬を含む領域である。

【0037】

この領域１１は、脈拍によって画素値が変化しやすい領域である。そのため、領域１１を含む画像に基づいて指標が計算されると、脈波変動指標等が精度良く計算できる。

【0038】

さらに、領域１１は、髪又は衣服等によって肌が隠れる場合が少ない領域である。したがって、領域１１を含む画像に基づいて指標が計算されると、脈波変動指標等が精度良く計算できる。

【0039】

なお、端末装置は、脈拍による画素値の変化が観察できる領域であれば、領域１１以外の領域を撮像してもよい。例えば、端末装置は、額又は指先等の部位を撮像してもよい。

【0040】

また、端末装置は、生体に対して非接触であるのが望ましい。

【0041】

＜脈波信号の抽出例＞（ステップＳ０２）
図４に戻り、ステップＳ０２では、端末装置は、脈波信号を抽出する。

【0042】

図６は、本発明の一実施形態に係る生体信号解析装置による脈波信号の抽出例を示すフローチャートである。図示する処理は、例えば、動画では、フレームごとに行われる。

【0043】

＜顔における特徴点座標の算出例＞（ステップＳ２１）
ステップＳ２１では、端末装置は、顔における特徴点座標を算出する。具体的には、まず、端末装置は、撮像された画像から、目、口及び鼻等の特徴点の座標を検出する。なお、各部位の検出は、例えば、公知の顔認証技術等によって実現できる。

【0044】

＜脈波信号の抽出に用いる画素の領域の設定例＞（ステップＳ２２）
ステップＳ２２では、端末装置は、脈波信号の抽出に用いる画素の領域を設定する。すなわち、端末装置は、領域１１から脈波信号が抽出できるように、設定を行う。具体的には、端末装置は、ステップＳ２１による算出結果に基づいて、生体の鼻及び頬を含む領域を設定する。また、この場合において、設定される領域は、目及び口が領域内に入らない程度に設定される。なお、領域の設定は、顔認証等に基づいて行われるに限られない。例えば、ユーザによる操作によって、領域の始点位置、幅及び高さ等が入力され、設定が行われてもよい。また、領域は、複数に分割されて設定されてもよい。例えば、領域は、鼻を含む領域と、左頬を含む領域と、右頬を含む領域とに分割されて設定されてもよい。

【0045】

＜領域内の画素値の平均化例＞（ステップＳ２３）
ステップＳ２３では、端末装置は、領域内の画素値を平均化する。具体的には、ステップＳ２３では、端末装置は、ステップＳ２２で設定される領域から生成される画像が有するＲ、Ｇ及びＢ等の画素値を平均化して、それぞれの平均値を計算する。

【0046】

生体の脈拍に起因する画素値の変化は、微小な変化である。そのため、１画素単位では、ノイズの影響が大きい。そこで、複数の画素が示すそれぞれの画素値を平均化すると、脈波信号に対するノイズの影響が低減できる。

【0047】

＜脈波信号の生成例＞（ステップＳ２４）
ステップＳ２４では、端末装置は、脈波信号を生成する。例えば、端末装置は、ステップＳ２３で計算される平均値に基づいて、下記（１）式を計算して、脈波信号の値を生成する。

ｐ０（ｎ）＝ａ_ｒ×ｒ（ｎ）＋ａ_ｇ×ｇ（ｎ）＋ａ_ｂ×ｂ（ｎ） （１）

上記（１）式では、「ｎ」は、フレーム番号を示す値である。また、「ｒ（ｎ）」は、「ｎ」フレーム目の画像が示すＲの画素値である。同様に、「ｇ（ｎ）」は、「ｎ」フレーム目の画像が示すＧの画素値である。さらに、「ｂ（ｎ）」は、「ｎ」フレーム目の画像が示すＢの画素値である。

【0048】

また、上記（１）式では、「ａ_ｒ」は、Ｒに対する重みとなる係数である。同様に、「ａ_ｇ」は、Ｇに対する重みとなる係数である。さらに、「ａ_ｂ」は、Ｂに対する重みとなる係数である。

【0049】

例えば、「ａ_ｒ」、「ａ_ｇ」及び「ａ_ｂ」は、あらかじめ「ａ_ｒ＝０」、「ａ_ｇ＝１」、「ａ_ｂ＝０」と設定されると、端末装置は、Ｇの成分だけを抽出した脈波信号を生成できる。生体の脈拍に起因する画素値の変化は、Ｇの成分から観察できる。したがって、上記のような設定とすると、端末装置は、生体の脈拍に起因する画素値の変化を観察しやすい脈波信号を生成できる。

【0050】

他にも、「ａ_ｒ」、「ａ_ｇ」及び「ａ_ｂ」は、あらかじめ「ａ_ｒ＝－ｋ」、「ａ_ｇ＝１」、「ａ_ｂ＝０」等と設定されてもよい。このような設定であると、脈波信号は、Ｇの成分から、「ｋ」で補正されたＲの成分を引いた値で生成される。このようにすると、端末装置は、Ｇの成分に含まれる体動等を起因とするノイズを低減させることができる。なお、ノイズは、例えば、周辺の光量の変化又は光源のちらつき等の周辺環境の変化が起因する場合もある。

【0051】

したがって、「ｋ」は、ノイズとなる成分が少なくなるように設定される。なお、「ｋ」は、正の値である。また、「ｋ」は、例えば、フレームごとに設定されてもよい。

【0052】

また、領域が複数設定される場合がある。このような場合には、端末装置は、まず、領域ごとに、それぞれの脈波信号を生成する。そして、端末装置は、複数の脈波信号を合成して１つの脈波信号としてもよい。具体的には、端末装置は、加算平均等によって、複数の脈波信号を合成する。すなわち、端末装置は、領域ごとの脈波信号を平均して、合成してもよい。他にも、端末装置は、領域ごとの脈波信号に重み付けをして合成してもよい。

【0053】

図７は、本発明の一実施形態に係る脈波信号の一例を示す図である。図では、横軸は、経過時間である。一方で、縦軸は、信号値である。

【0054】

ステップＳ０１、すなわち、図６に示す処理が行われると、端末装置は、図示するような脈波信号を抽出できる。図示するように、脈波信号は、脈拍に応じて信号値が周期的に変化する時系列データである。

【0055】

なお、１フレームが「Δｔ」で一定、すなわち、「Δｔ」ごとに撮像が行われるとすると、時間は、「Δｔ×ｎ」となる。「ｎ」は、フレーム番号を示す値である。

【0056】

＜ＳＮ比の計算例＞（ステップＳ０３）
図４に戻り、ステップＳ０３では、端末装置は、ＳＮ比を計算する。

【0057】

例えば、端末装置は、脈波信号のパワースペクトルに基づいて、まず、信号成分となる「Ｖｓ」と、ノイズ成分となる「Ｖｎ」とを計算する。そして、「Ｖｓ」と、「Ｖｎ」との比を計算すると、端末装置は、ＳＮ比を計算できる。

【0058】

また、「Ｖｓ」は、例えば、信号成分となる周波数帯域を決定し、脈波信号のパワースペクトルにおいて、信号成分となる周波数帯域に出力されるパワーの平均値を計算すると特定できる。同様に、「Ｖｎ」は、例えば、ノイズ成分となる周波数帯域を決定し、脈波信号のパワースペクトルにおいて、ノイズ成分となる周波数帯域に出力されるパワーの平均値を計算すると特定できる。なお、ＳＮ比は、上記以外の方法によって計算されてもよい。

【0059】

図８は、本発明の一実施形態に係るＳＮ比を計算する処理の一例を示すフローチャートである。

【0060】

＜パワースペクトルの算出例＞（ステップＳ３１）
ステップＳ３１では、端末装置は、パワースペクトルを算出する。具体的には、ステップＳ３１では、まず、端末装置は、データ数が２の累乗となるように調整する。次に、端末装置は、調整されたデータに対して、高速フーリエ変換（ＦＦＴ、Ｆａｓｔ Ｆｏｕｒｉｅｒ Ｔｒａｎｓｆｏｒｍ）を行う。続いて、端末装置は、各周波数におけるパワーをそれぞれ算出する。

【0061】

図９は、本発明の一実施形態に係るＳＮ比の計算において算出されるパワースペクトルの一例を示す図である。例えば、図示するグラフがステップＳ３１によって生成される。なお、図９では、縦軸は、パワーを示し、一方で、横軸は、周波数を示す。

【0062】

＜脈波周波数の検出例＞（ステップＳ３２）
図８に戻り、ステップＳ３２では、端末装置は、脈波周波数を検出する。以下、図９に示すパワースペクトルがステップＳ３１によって算出された例で説明する。この例では、脈波周波数Ｆｐは、図示するように検出される。

【0063】

生体が安静状態であると、脈拍数は、１分当たり３０乃至１２０回程度であることが多い。したがって、パワースペクトルでは、この脈拍数に応じた周期で、強いピークの基本波が検出される。すなわち、周波数が０．５乃至２．０Ｈｚとなる帯域で、端末装置は、パワースペクトルにおいて、最もパワーが大きくなる周波数を検出する。このようにすると、端末装置は、脈波周波数Ｆｐを図示するように検出できる。具体的には、図９に示す例では、周波数が１．０Ｈｚ程度となる周波数で、強いピークが検出される。したがって、この例では、周波数が１．０Ｈｚ程度となる周波数が脈波周波数Ｆｐと検出される。

【0064】

＜信号及びノイズのそれぞれの周波数帯域の決定例＞（ステップＳ３３）
図８に戻り、ステップＳ３３では、端末装置は、信号及びノイズのそれぞれの周波数帯域を決定する。例えば、信号成分の周波数帯域は、ステップＳ３２で検出される脈波周波数を中心として決定される。具体的には、図９では、信号成分の周波数帯域Ｆｓは、「Ｆｐ－ΔＦ」乃至「Ｆｐ＋ΔＦ」の帯域と決定される例である。なお、図９では、「ΔＦ」は、０．２Ｈｚと設定された例である。また、「ΔＦ」は、あらかじめ設定される値である。例えば、「ΔＦ」は、パワースペクトルのピーク形状等を考慮して設定される。また、信号成分の周波数帯域Ｆｓには、脈波周波数Ｆｐの整数倍の周波数に現れる高調波のピークを含む周波数帯域が含まれてもよい。

【0065】

一方で、ステップＳ３３では、端末装置は、信号成分の周波数帯域Ｆｓ以外の周波数帯域をノイズ成分の周波数帯域Ｆｎと決定する。ただし、脈波周波数Ｆｐより低周波のパワーは、計測される際の生体の体動等によって影響を受ける。したがって、端末装置は、図９に示すように、信号成分の周波数帯域Ｆｓより高周波となる周波数帯域をノイズ成分の周波数帯域Ｆｎと決定してよい。なお、ノイズ成分の周波数帯域Ｆｎには、信号成分の周波数帯域Ｆｓより低周波の周波数帯域が含まれてもよい。また、ノイズ成分の周波数帯域Ｆｎの上限となる周波数は、（サンプリング周波数／２）である。すなわち、例えば、カメラのフレームレートが「３０ｆｐｓ（フレーム／秒）」とすると、ノイズ成分の周波数帯域Ｆｎの上限となる周波数は、「３０Ｈｚ／２＝１５Ｈｚ」である。

【0066】

＜ＳＮ比の計算例＞（ステップＳ３４）
ステップＳ３４では、端末装置は、ＳＮ比を計算する。具体的には、ステップＳ３４では、まず、端末装置は、ステップＳ３３で決定される信号成分の周波数帯域と、ノイズ成分の周波数帯域とに基づいて、「Ｖｓ」及び「Ｖｎ」のパワーのそれぞれの平均値を計算する。次に、端末装置は、計算される「Ｖｓ」及び「Ｖｎ」の比を計算して、ＳＮ比を計算する。

【0067】

＜ＳＮ比が閾値以上であるか否かの判断例＞（ステップＳ０４）
図４に戻り、ステップＳ０４では、端末装置は、ＳＮ比が閾値以上であるか否かを判断する。すなわち、端末装置は、ＳＮ比と、あらかじめ設定される閾値とを比較して、ＳＮ比が閾値以上であるか否かを判断する。

【0068】

ＳＮ比が閾値以上である脈波信号であると、端末装置は、脈波信号に基づいて精度良く指標を計算できる。したがって、ＳＮ比が閾値未満となる小さい値である場合には、端末装置は、撮像条件の変更等の調整を行う。

【0069】

次に、ＳＮ比が閾値以上であると判断されると（ステップＳ０４でＹＥＳ）、端末装置は、ステップＳ０６に進む。一方で、ＳＮ比が閾値以上ではないと判断されると（ステップＳ０４でＮＯ）、端末装置は、ステップＳ０５に進む。なお、ＳＮ比が閾値以上でないと判断された場合には、端末装置は、アラーム等を表示して、処理を中止してからステップＳ０５に進んでもよい。

【0070】

図１０は、本発明の一実施形態に係る閾値の設定例を示す図である。なお、図１０は、ＳＮ比と、ピーク検出誤差との関係の一例を示す図である。

【0071】

図では、横軸は、ＳＮ比を示す。一方で、縦軸は、端末装置による解析によって検出されるピーク間隔と、接触式等の脈波を計測する装置によって検出されるピーク間隔、すなわち、真のピーク間隔との誤差を示す。

【0072】

図示するように、ＳＮ比が小さい値であると、ピーク検出誤差は、増大する関係となる。したがって、まず、端末装置は、ピーク検出誤差と、ＳＮ比との関係を示す関係式ＦＲを算出する。なお、関係式ＦＲは、例えば、各点に基づいて最小二乗法等によって算出される。

【0073】

次に、端末装置は、後段に行われる解析の目的等に合わせて、許容できるピーク検出誤差、すなわち、縦軸の値を決定する。続いて、端末装置は、関係式ＦＲに基づいて、決定される縦軸の値に対する値を閾値ＴＨと決定する。このようにして、端末装置は、閾値ＴＨを設定する。また、この決定によって定まる閾値ＴＨがあらかじめ端末装置に設定され、ステップＳ０４が行われる。

【0074】

＜撮像条件の変更例＞（ステップＳ０５）
図４に戻り、ステップＳ０５では、端末装置は、撮像条件を変更する。すなわち、端末装置は、ＳＮ比が改善するように、撮像条件を変更する。

【0075】

撮像条件は、例えば、照明の明るさ、撮像装置と生体の距離、露光時間等のカメラパラメータ又はこれらの組み合わせである。したがって、端末装置は、例えば、照明の明るさが不足していることを知らせるメッセージを表示する等によって、ユーザに撮像条件を変更させる操作を行わせる。他にも、例えば、端末装置は、露光時間を変更する制御を行う。

【0076】

図１１は、本発明の一実施形態に係る生体信号解析装置による制御の一例を示すタイミングチャートである。以下、シャッタスピードの初期設定等によって定まる露光時間が図１１（Ａ）に示す露光時間（以下「第１露光時間ＥＴ１」という。）である例で説明する。

【0077】

まず、初期設定等では、フレームレートが設定される。そして、フレームレートが定まると、サンプリング周期が定まる。以下、図１１（Ａ）に示すサンプリング周期（以下「第１サンプリング周期ＳＴ１」という。）であるとする。なお、露光時間は、カメラのシャッタが切られ、光が取り込まれている時間、すなわち、露光が行われている時間である。

【0078】

次に、ＳＮ比が所定の閾値以上でないと判断されると、端末装置は、例えば、露光時間を図１１（Ｂ）に示す露光時間（以下「第２露光時間ＥＴ２」という。）に変更するように制御する。具体的には、端末装置は、シャッタスピードの設定値等を変更して露光時間が第２露光時間ＥＴ２となるように制御する。

【0079】

図示するように、第２露光時間ＥＴ２は、第１露光時間ＥＴ１と比較すると、長い時間である。このように、露光時間が長くなると、１フレームにおける受光量が増加するため、パワースペクトルにおいて信号成分が強くなる場合が多い。すなわち、露光時間が長くなると、ＳＮ比を大きい値にすることができる。また、ＳＮ比と、ピーク検出誤差との関係は、図１０に示す関係にあることが多いため、ＳＮ比を大きい値にすると、ピーク検出誤差を減少させることができる。したがって、端末装置は、露光時間を変更してＳＮ比を大きい値にし、ピーク検出誤差を減少させることができる。

【0080】

また、端末装置は、脈波信号をサンプリングする周期の例であるフレームレートによって定まる周期を遅くしてもよい。図１１に示す例では、端末装置は、第１サンプリング周期ＳＴ１を図１１（Ｂ）に示すサンプリング周期（以下「第２サンプリング周期ＳＴ２」という。）に変更する制御を行う。露光時間は、１周期において、フレームレートによって定まる周期より、長く設定できない場合が多い。そこで、端末装置は、図１１（Ｂ）に示す第２サンプリング周期ＳＴ２のように、サンプリングする周期が長くなるようにフレームレート等を変更する制御を行う。例えば、変更前、すなわち、図１１（Ａ）では、フレームレートが「６０ｆｐｓ」であるとすると、端末装置は、フレームレートを「３０ｆｐｓ」とする。このように、フレームレートを遅くすると、１周期において、第１露光時間ＥＴ１を第２露光時間ＥＴ２とするように、端末装置は、露光時間を長くできる。

【0081】

＜脈波解析例＞（ステップＳ０６）
図４に戻り、ステップＳ０６では、端末装置は、脈波を解析する。

【0082】

図１２は、本発明の一実施形態に係る生体信号解析装置による脈波解析の一例を示すフローチャートである。

【0083】

＜周波数スペクトルの算出例＞（ステップＳ６１）
ステップＳ６１では、端末装置は、第１信号の周波数スペクトルを算出する。例えば、端末装置は、第１信号の例である脈波信号の時系列データに対して、種々のスペクトル解析手法を適用すると、周波数スペクトルを算出できる。

【0084】

具体的には、ＦＦＴを用いると、端末装置は、振幅の周波数分布を示す振幅スペクトル、位相の周波数分布を示す位相スペクトル又はパワーの周波数分布を示すパワースペクトル等の周波数スペクトルを算出できる。

【0085】

＜脈拍間隔の低周波揺らぎ成分の抽出例＞（ステップＳ６２）
ステップＳ６２では、端末装置は、ステップＳ６１による周波数スペクトルの算出結果に基づいて、第１周波数帯における揺らぎ成分の例である脈拍間隔の低周波揺らぎ成分（以下「ＬＦ揺らぎ成分」という。）を抽出する。

【0086】

ＬＦ揺らぎ成分とは、脈拍間隔の時系列データからパワースペクトルを算出すると、所定の周波数以下の周波数帯に含まれる成分である。所定の周波数は、例えば、０．１５Ｈｚである。

【0087】

また、ＬＦ揺らぎ成分は、例えば、振幅スペクトル又はパワースペクトルに対して所定の周波数帯のみを通すバンドパスフィルタを適用すると、抽出できる。なお、バンドパスフィルタが通す周波数帯は、脈波ピークの周波数帯を中心として幅が「±Δｆ」となる範囲等に設定される。具体的には、「±Δｆ」は、ＬＦ揺らぎ成分を確実に抽出するため、所定の周波数（例えば、０．１５Ｈｚである。）以上の値が設定されるのが望ましい。ＬＦ揺らぎ成分による影響は、脈波ピークの幅の広がりとして現れる。したがって、脈波ピークの周波数帯を中心として幅が「±Δｆ」となる範囲に設定されたバンドパスフィルタを振幅スペクトル又はパワースペクトルに対して適用すると、端末装置は、ＬＦ揺らぎ成分を抽出できる。

【0088】

また、周波数スペクトルに白色ノイズによるオフセット成分等が含まれる場合には、端末装置は、スペクトルサブトラクション法（ｓｐｅｃｔｒａｌ ｓｕｂｔｒａｃｔｉｏｎ ｍｅｔｈｏｄ）等のノイズ除去法を組み合わせてもよい。

【0089】

＜時系列データに逆変換する例＞（ステップＳ６３）
ステップＳ６３では、端末装置は、ステップＳ６２で抽出したＬＦ揺らぎ成分の振幅スペクトルを時系列データに逆変換する。例えば、時系列データは、逆フーリエ変換（ＩＦＦＴ（Ｉｎｖｅｒｓｅ ＦＦＴ））等によって算出できる。この場合には、位相には、ステップＳ６１におけるＦＦＴで算出した位相が用いられるとよい。

【0090】

＜脈拍間隔の検出例＞（ステップＳ６４）
ステップＳ６４では、端末装置は、ステップＳ６３で算出される脈波信号の時系列データから、脈拍間隔の時系列データを生成する。脈拍信号から検出した「ｍ」番目の脈拍のピーク時間を「Ｔ_ｍ」とし、ピーク時間「Ｔ_ｍ」における脈拍間隔を「Ｉ（Ｔ_ｍ）」とすると、脈拍間隔「Ｉ（Ｔ_ｍ）」は、下記（２）式のように計算できる。

Ｉ（Ｔ_ｍ）＝Ｔ_ｍ－Ｔ_ｍ－１ （２）

また、上記（２）式におけるピーク時間「Ｔ_ｍ」は、例えば、極大点ＬＭが発生する時間である。例えば、端末装置は、信号値が極小となる時点と、次に信号値が極小となる時点との間において、信号値が最大となる時点を特定すると、極大点ＬＭを特定できる。

【0091】

また、前の極大点ＬＭが発生する時間である「Ｔ_ｍ－１」からピーク時間「Ｔ_ｍ」までが脈拍の１周期に相当する。さらに、極大点ＬＭとなる時点は、光電式容積脈波計（ＰＰＧ）で計測すると観察できる容積脈波の立ち上がり開始時点に相当する。

【0092】

なお、ピーク時間「Ｔ_ｍ」は、極大点ＬＭに基づいて計算されるに限られず、例えば、極小点に基づいて計算されてもよい。脈波信号が極小となる時点は、光電式容積脈波計で計測できる容積脈波の振幅が最大となる時点に相当する。

【0093】

また、脈拍のピーク時間は、補正されてもよい。補正例は、後述する。

【0094】

＜脈拍間隔の低周波ゆらぎ成分を平坦化した第２信号の生成例＞（ステップＳ６５）
ステップＳ６５では、端末装置は、抽出された脈波信号「ｐ_０（ｎ）」と、ステップＳ６４で算出した脈拍間隔の時系列データに基づいて、脈拍間隔の低周波ゆらぎ成分を平坦化して第２信号を生成する。第２信号を「ｐ_１（ｔ_ｎ）」とし、第１信号を「ｐ_０（ｎ）＝ｐ_０（ｎΔｔ）」とすると、第２信号は、例えば、下記（３）式のように生成される。

ｐ_１（ｔ_ｎ）＝ｐ_０（ｎΔｔ）
ｔ_ｎ＝ｔ_ｎ－１＋Δｔ_ｎ＝ｔ_ｎ－１＋Δｔ×（Ｉ_ａｖｅ／Ｉ（ｎΔｔ）） （３）

上記（３）式において、「Ｉ_ａｖｅ」は、脈拍間隔の平均値である。なお、「Ｉ_ａｖｅ」は、脈拍間隔の平均値に限られず、別の固定値でもよい。また、上記（３）式において、「Ｉ（ｎΔｔ）」は、経過時間「ｔ＝ｎΔｔ」における脈拍間隔である。例えば、「ｔ＝ｎΔｔ」が区間［Ｔ_ｍ－１，Ｔ_ｍ］に含まれる場合には、「Ｉ（Ｔ_ｍ）」に「Ｉ（ｎΔｔ）」の値が用いられる。

【0095】

また、上記（３）式において、「ｐ_１（ｔ_ｎ）」は、不等間隔「Δｔ_ｎ」でサンプリングされた信号である。したがって、「Δｔ」でリサンプリングすると、端末装置は、「ｐ_１（ｎ）＝ｐ_１（ｎΔｔ）」が得られる。なお、リサンプリングを実現する方法は、公知の方法等が用いられる。

【0096】

なお、ＬＦ揺らぎ成分を平坦化する方法は、上記の方法に限られず、他の方法でもよい。

【0097】

＜周波数スペクトルの算出例＞（ステップＳ６６） ステップＳ６６では、端末装置は、ステップＳ６５で生成される第２信号「ｐ_１（ｔ）」の周波数スペクトルを算出する。例えば、算出方法は、ステップＳ６１と同様である。

【0098】

＜ノイズ成分の除去例＞（ステップＳ６７）
ステップＳ６７では、端末装置は、ステップＳ６６で算出される振幅スペクトルから、脈波成分及びノイズ成分の周波数成分の違いを利用してノイズ成分を除去する。

【0099】

ＬＦ揺らぎ成分が平坦化されて生成される第２信号の振幅スペクトルを用いるため、脈波成分に起因するピークがはっきりと観測できる。そのため、従来ではノイズ成分と判別が難しい、第２周波数帯における揺らぎ成分の例である脈拍間隔の高周波揺らぎ成分（以下「ＨＦ揺らぎ成分」という。）のピークが、精度良く抽出できる。なお、ＨＦ揺らぎ成分とは、脈拍間隔の時系列データからパワースペクトルを算出すると、所定の周波数帯に含まれる成分である。所定の周波数帯は、例えば、０．１５Ｈｚ乃至０．４０Ｈｚである。

【0100】

ノイズ成分を除去するには、例えば、所定の周波数帯の成分のみを通す周波数フィルタを適用する。また、周波数フィルタは、脈波信号の基本波及び高調波に起因するピークと、ＨＦ揺らぎ成分に起因するピークとを抽出できるように設計される。このような設計の周波数フィルタを振幅スペクトルに対して適用すると、ＨＦ揺らぎ成分を含む脈波成分を抽出し、ノイズ成分を除去できる。

【0101】

各ピークの周波数は、例えば、以下のような方法で計算される。

【0102】

脈波信号の基本波Ｌ１に起因するピークの周波数をｆ（Ｌ１）とすると、ｆ（Ｌ１）は、所定の周波数帯において振幅が最大となる周波数である。この場合では、所定の周波数帯は、生体が取り得る脈拍の平均周波数に基づいて設定されるのが望ましい。例えば、生体が安静状態であると、所定の周波数帯は、０．５Ｈｚ乃至２．０Ｈｚとなる。

【0103】

脈波信号の高調波は、Ｌ２、Ｌ３、・・・・とあり、各高調波のピークの周波数ｆ（Ｌ２）、ｆ（Ｌ３）、・・・・は、ｆ（Ｌ１）を整数倍した周波数の近傍の周波数帯で振幅が最大となる周波数から特定できる。例えば、２次高調波Ｌ２に起因するピークの周波数ｆ（Ｌ２）は、ｆ（Ｌ１）×２Ｈｚの近傍の周波数帯で振幅が最大となる周波数から特定できる。

【0104】

ＨＦ揺らぎ成分に起因するピークＨ１ａの周波数ｆ（Ｈ１ａ）は、（ｆ（Ｌ１）－０．４０）Ｈｚ乃至（ｆ（Ｌ１）－０．１５）Ｈｚの周波数帯において振幅が最大となる周波数から特定できる。ＨＦ揺らぎ成分は、呼吸に起因する揺らぎ成分である。そして、揺らぎ成分の周波数Ｆｈは、一般的に０．１５Ｈｚ乃至０．４０Ｈｚの周波数帯に含まれる。さらに、ＨＦ揺らぎ成分に起因するピークは、基本波Ｌ１に起因するピークの周波数ｆ（Ｌ１）を中心に低周波側と、高周波側とに、それぞれ揺らぎ成分の周波数Ｆｈ分離れた周波数付近に現れる。また、基本波Ｌ１に起因するピークの周波数ｆ（Ｌ１）と、ＨＦ揺らぎ成分に起因するピークＨ１ａの周波数ｆ（Ｈ１ａ）との周波数差Δｆｒを特定すると、ＨＦ揺らぎ成分に起因するピークＨ１ｂの周波数ｆ（Ｈ１ｂ）は、ｆ（Ｌ１）＋Δｆｒと特定できる。

【0105】

各ピークを抽出できる周波数フィルタは、例えば、下記（４）式のようになる。

【0106】

【数1】

上記（４）式において、「ｓ」は、次数を示す値である。また、「Δｆ'」は、ＨＦ揺らぎ成分に起因するピークの幅を考慮して、「Δｆ」を補正する値である。「Δｆ'」及び「Δｆ」は、「Δｆ'≧Δｆ」という関係になる。上記のようなフィルタを用いると、端末装置は、第２信号に含まれるノイズ成分を除去できる。なお、端末装置は、上記（４）式に示す以外のフィルタによってノイズを除去してもよい。例えば、周波数に応じて、上記（４）式における「１」及び「０」は、他の値でもよい。さらに、「Δｆ'」は、「ｓ」の値、すなわち、高調波の次数によって変更されてもよい。

【0107】

＜時系列データに逆変換する例＞（ステップＳ６８）
ステップＳ６８では、端末装置は、ステップＳ６７で抽出される脈波成分の振幅スペクトルを時系列データに逆変換する。例えば、変換方法は、ステップＳ６３と同様である。

【0108】

＜脈拍間隔の低周波ゆらぎ成分を非平坦化した脈波信号の生成例＞（ステップＳ６９）
ステップＳ６９では、端末装置は、ステップＳ６８で算出された脈波成分の時系列データ「ｐ_１'（ｎ）」に対して、ステップＳ６５で平坦化によって除去したＬＦ揺らぎ成分を付与する非平坦化処理を行う。

【0109】

非平坦化した脈波信号を「Δｔ」の間隔でサンプリングした時系列データを「ｐ_２（ｎ）」とすると、「ｐ_２（ｎ）」は、下記（５）式のようになる。

ｐ_２（ｎ）＝ｐ_２（ｎΔｔ）＝ｐ_１'（ｔ_ｎ） （５）

なお、経過時間「ｔ＝ｔ_ｎ」における信号値は、例えば、補間によって求まる。補間の方法は、例えば、線形補間又はスプライン補間等である。

【0110】

＜指標の計算例＞（ステップＳ０７）
図４に戻り、ステップＳ０７では、端末装置は、ステップＳ０６でノイズ成分が除去された脈波信号を用いて脈波指標を計算する。

【0111】

図１３は、本発明の一実施形態に係る生体信号解析装置による脈波変動指標の算出例を示すフローチャートである。

【0112】

＜脈拍間隔の検出例＞（ステップＳ７１）
ステップＳ７１では、端末装置は、ステップＳ０６でノイズ成分が除去された脈波信号から脈拍間隔を検出する。例えば、検出方法は、ステップＳ６４と同様である。

【0113】

＜リサンプリング例＞（ステップＳ７２）
ステップＳ７２では、端末装置は、脈拍間隔の時系列データにおいて、時間間隔が等間隔となるようにリサンプリングを行う。リサンプリングの時間間隔は、例えば、０．２５秒とする。また、リサンプリングする経過時間における信号値は、例えば、補間によって求められる。補間方法は、例えば、線形補間又はスプライン補間等である。

【0114】

＜パワースペクトルの算出例＞（ステップＳ７３）
ステップＳ７３では、端末装置は、ステップＳ７２で算出した脈拍間隔の時系列データからパワースペクトルを算出する。例えば、パワースペクトルの算出方法は、既知の周波数解析手法等である。具体的には、最大エントロピー法によって、指定の周波数におけるパワーが算出できる。なお、ラグ等のパラメータは、あらかじめ設定される。

【0115】

＜脈拍変動指標の算出例＞（ステップＳ７４）
ステップＳ７４では、端末装置は、ステップＳ７３で算出したパワースペクトルに基づいて、脈拍変動指標を算出する。例えば、脈拍変動指標は、既知の方法で算出される。具体的には、端末装置は、０．０４Ｈｚ乃至０．１５Ｈｚの周波数帯におけるパワー積分値をＬＦ値、０．１５Ｈｚ乃至０．４０Ｈｚの周波数帯におけるパワー積分値をＨＦ値とし、これらの比を計算してＬＦ／ＨＦ値を算出する。

【0116】

＜処理結果例＞
まず、図７に示すような脈波信号が、脈波解析の対象となる。そして、周波数スペクトルの算出が行われると（ステップＳ６１）、例えば、以下のような結果が得られる。

【0117】

図１４は、本発明の一実施形態に係る生体信号解析装置による第１信号の周波数解析結果例を示す図である。図は、脈波信号に対してＦＦＴを行った結果である。図では、横軸は、周波数であり、縦軸は、各周波数の振幅スペクトルである。

【0118】

脈波信号は、生体が取り得る脈拍の周波数帯に、ピークが発生する。例えば、生体が安静状態であると、脈拍の平均周波数は、０．５Ｈｚ乃至２．０Ｈｚであるため、この範囲にピークが発生しやすい。図示する例は、１．０Ｈｚ付近に最も強いピークが発生した例である。なお、ピークは、図示するように、１．０Ｈｚ付近に限られず、例えば、生体の状態等によって、ピークが発生する周波数帯は異なる。

【0119】

以下、脈波によるピークが含まれる成分を「脈波ピーク成分ＰＫ１」という。一方で、脈波信号には、ノイズ成分ＮＺが含まれる場合が多い。

【0120】

図１５は、本発明の一実施形態における脈拍間隔の高周波の揺らぎ成分及び低周波の揺らぎ成分の例を示す図である。例えば、図１５（Ａ）に示すような脈波信号には、図１５（Ｂ）に示すようなＬＦ揺らぎ成分と、図１５（Ｃ）に示すようなＨＦ揺らぎ成分とが含まれる。

【0121】

図示するように、脈波信号の脈拍間隔には、図１５（Ｂ）に示すようなＬＦ揺らぎ成分ＦＫＬが発生する場合がある。したがって、ＬＦ揺らぎ成分ＦＫＬによって、ピーク間隔は、時間の経過において、長くなったり、短くなったりする。

【0122】

同様に、脈波信号の脈拍間隔には、図１５（Ｃ）に示すようなＨＦ揺らぎ成分ＦＫＨが発生する場合がある。したがって、ＨＦ揺らぎ成分ＦＫＨによって、ピーク間隔は、時間の経過において、長くなったり、短くなったりする。

【0123】

ＬＦ揺らぎ成分ＦＫＬは、脈拍間隔の時系列データのパワースペクトルを算出すると、所定の周波数以下の周波数帯域に含まれる成分である。

【0124】

ＬＦ揺らぎ成分ＦＫＬは、脈波ピーク成分ＰＫ１の近傍の周波数に現れる。すなわち、揺らぎがあると、ピーク間隔が広くなったり、狭くなったりするため、ピークが現れる周波数がばらついて幅が生じる。

【0125】

図１６は、本発明の一実施形態に係る生体信号解析装置による脈拍間隔の検出例を示す図である。まず、上記の通り、逆フーリエ変換等が行われると（ステップＳ６３）、図示するような時系列データが得られる。

【0126】

また、脈拍間隔は、揺らぎがあると、図示するように、一定でなく、変動する。具体的には、揺らぎによって、例えば、脈拍間隔「Ｉ（Ｔ_ｍ＋１）」のように、間隔が長くなったり、脈拍間隔「Ｉ（Ｔ_ｍ）」のように、間隔が短くなったりする。

【0127】

ＬＦ揺らぎ成分は、ピーク間隔の揺らぎの周期が例えば１０秒程度の周期になる。そして、脈拍間隔は、ＬＦ揺らぎ成分の影響により変動する。また、ＨＦ揺らぎ成分は、ピーク間隔の揺らぎの周期が例えば３乃至４秒程度の周期になる。そして、脈拍間隔は、ＨＦ揺らぎ成分の影響により変動する。なお、これらの値は、個人差がある。

【0128】

なお、端末装置は、ステップＳ６４及びステップＳ７１において脈拍間隔を検出するとき、脈拍信号のピーク時間を補正してもよい。

【0129】

図１７は、本発明の一実施形態に係る生体信号解析装置によるピーク時間の補正例を示す図である。図示する補正は、いわゆる折れ線近似である。以下、補正によるピーク間隔の補正量をピーク補正量「ΔＴ」とする。また、計測される信号値のうち、ピークとなる時点の信号値を第１信号値「Ｐ_ｉ」とする。さらに、第１信号値Ｐ_ｉより１つ前のピークとなる信号値を第２信号値「Ｐ_ｉ－１」とする。同様に、第１信号値Ｐ_ｉより１つ後のピークとなる信号値を第３信号値「Ｐ_ｉ＋１」とする。この場合には、ピーク補正量「ΔＴ」は、図示する式等によって計算される。

【0130】

図示する式によって、ピーク補正量「ΔＴ」が計算されると、ピーク補正量「ΔＴ」は、信号値が計測される間隔「Δｔ」より短い値にできるため、端末装置は、信号の時間分解能を高くできる。すなわち、ピーク時間を補正して、補正ピーク時間「Ｔ'_ｍ」を用いると、端末装置は、ピーク時間を精度良く検出できる。そのため、端末装置は、時間分解能が低い信号からでも、補正ピーク時間「Ｔ'_ｍ」を用いることで、脈拍間隔を精度良く算出できる。

【0131】

なお、算出した脈拍間隔に高周波成分が含まれる場合には、端末装置は、フィルタリングによって高周波成分を除去するのが望ましい。

【0132】

図１８は、本発明の一実施形態に係る平坦化の例を示す図である。図は、低周波成分のみを示す図である。上記のように、平坦化が行われると、例えば、図示するようにＬＦ揺らぎ成分ＦＫＬがほぼない状態が得られる。このような平坦化が行われると、図１５（Ｂ）に示すようなＬＦ揺らぎ成分ＦＫＬが平坦化され、一方で、図１５（Ｃ）に示すようなＨＦ揺らぎ成分ＦＫＨは、そのままとなる。

【0133】

図１９は、本発明の一実施形態に係る第２信号の一例を示す図である。まず、第２信号は、例えば、図示するような周波数特性となる信号である。

【0134】

図示するように、まず、所定の周波数帯において最も振幅が大きい周波数成分を特定すると、端末装置は、基本波Ｌ１を特定できる。なお、所定の周波数帯は、生体が取り得る脈拍の平均周波数に基づき設定すればよい。例えば、生体が安静状態の場合は、０．５Ｈｚ乃至２．０Ｈｚを所定の周波数帯として用いることができる。

【0135】

一方で、第２波の例であるＨＦ揺らぎ成分に起因するピークＨ１ａ及び第３波の例であるＨＦ揺らぎ成分に起因するピークＨ１ｂ、すなわち、ＨＦ揺らぎ成分ＦＫＨによる周波数成分は、基本波Ｌ１より振幅が小さい周波数成分である。特に、ＨＦ揺らぎ成分に起因するピークＨ１ｂは、図示するように、ノイズ成分等と切り分けが難しい程度の振幅であることが多い。

【0136】

そこで、まず、図示するように、端末装置は、ＨＦ揺らぎ成分に起因するピークＨ１ａを特定する。具体的には、端末装置は、例えば、（基本波の周波数ｆ（Ｌ１）－０．４０Ｈｚ）乃至（基本波の周波数ｆ（Ｌ１）－０．１５Ｈｚ）の周波数帯域において、最も振幅が大きい周波数成分を特定すると、端末装置は、ＨＦ揺らぎ成分に起因するピークＨ１ａを特定できる。

【0137】

そして、ＨＦ揺らぎ成分に起因するピークＨ１ｂは、基本波Ｌ１を中心にＨＦ揺らぎ成分に起因するピークＨ１ａと対称となる周波数帯域に発生する性質がある。そこで、あらかじめ特定されるＨＦ揺らぎ成分に起因するピークＨ１ｂと基本波Ｌ１を中心に対称となる周波数帯域を特定すると、端末装置は、ＨＦ揺らぎ成分に起因するピークＨ１ｂを特定しやすい。

【0138】

具体的には、基本波Ｌ１と、ＨＦ揺らぎ成分に起因するピークＨ１ａとの周波数差を「Δｆｒ」とすると、端末装置は、基本波Ｌ１より高い周波数帯域において、「Δｆｒ」分の周波数差がある周波数帯域、すなわち、「基本波の周波数ｆ（Ｌ１）＋Δｆｒ」となる周波数帯域を特定すると、ＨＦ揺らぎ成分に起因するピークＨ１ｂを特定できる。

【0139】

また、基本波Ｌ１の整数倍となる周波数成分を含む周波数帯域にも、同様の周波数成分が発生することが多い。例えば、整数倍となる周波数成分を含む周波数帯域は、２次高調波Ｌ２及び３次高調波Ｌ３等である。

【0140】

したがって、端末装置は、基本波Ｌ１、ＨＦ揺らぎ成分に起因するピークＨ１ａ、ＨＦ揺らぎ成分に起因するピークＨ１ｂ、２次高調波Ｌ２及び３次高調波Ｌ３の周波数帯域以外の周波数帯域を少なくともフィルタリングすると、ノイズ成分を除去できる。

【0141】

なお、２次高調波Ｌ２及び３次高調波Ｌ３を特定し、「Δｆｒ」に基づいて、図示するように、２次高調波Ｌ２及び３次高調波Ｌ３等でも、ＨＦ揺らぎ成分に起因するピークが特定されてもよい。図示する例では、２次高調波Ｌ２及び３次高調波Ｌ３の周波数から、高周波側及び低周波側に、「Δｆｒ」分の周波数差がある周波数帯域に所定値以上の振幅がある場合には、２次高調波Ｌ２及び３次高調波Ｌ３等にも、ＨＦ揺らぎ成分に起因するピークがあるとする。なお、所定値は、あらかじめ設定される値である。このような場合には、２次高調波Ｌ２及び３次高調波Ｌ３に対するＨＦ揺らぎ成分に起因するピークが含まれる周波数帯は、フィルタリングの対象外とする。

【0142】

したがって、図示する例では、端末装置は、設定周波数帯以外の周波数帯の例であるフィルタリング領域ＦＬＴをフィルタリングしてノイズ成分を減衰させる。

【0143】

図２０は、本発明の一実施形態に係る出力信号の一例を示す図である。例えば、図７に示すような第１信号に対して、図１２に示すような処理が行われると、端末装置は、図示するような出力信号を生成できる。すなわち、端末装置は、第２信号に含まれるノイズ成分を除去すると、図示するような出力信号を生成できる。

【0144】

＜効果例＞
図２１は、本発明の一実施形態に係る効果の一例を示す図である。まず、仮に、ＬＦ揺らぎ成分ＦＫＬがなく、かつ、ＨＦ揺らぎ成分ＦＫＨがないとすると、脈波信号は、図２１（Ａ）のような周波数特性となる。図示するように、脈波信号は、脈拍の一定な動作による周波数（以下「基本周波数ＢＦ」という。）が主な成分となる信号となる。

【0145】

一方で、実際には、脈波信号には、ＬＦ揺らぎ成分ＦＫＬがある。したがって、ＬＦ揺らぎ成分ＦＫＬによって、図２１（Ｂ）のような周波数特性となる場合が多い。図示するように、脈拍間隔の低周波成分は、基本周波数ＢＦのピークの広がりとして現れる。

【0146】

さらに、ＨＦ揺らぎ成分ＦＫＨがある。したがって、脈波信号は、例えば、図２１（Ｃ）のような周波数特性となる信号である。図２１（Ｂ）と比較すると、ＨＦ揺らぎ成分ＦＫＨの成分が加わる点が異なる。図示するように、ＨＦ揺らぎ成分ＦＫＨが加わる場合は、基本周波数ＢＦの前後となる周波数帯域にピークが発生しやすく、基本周波数ＢＦにおける振幅と比較すると、ピークにおける振幅が弱い場合が多い。このときの脈波信号の振幅スペクトルの例が図１４である。図示するように、脈波信号にノイズ成分があると、ＨＦ揺らぎ成分ＦＫＨの成分は、ノイズ成分との判別が困難になる場合が多い。

【0147】

そこで、端末装置は、ＨＦ揺らぎ成分ＦＫＨを正確に検出するために、ＬＦ揺らぎ成分ＦＫＬを平坦化し、ＨＦ揺らぎ成分ＦＫＨが保存された第２信号を生成する。第２信号は、例えば、図１９に示すような周波数特性となる。図示するように、ＬＦ揺らぎ成分ＦＫＬの平坦化によって、ＨＦ揺らぎ成分ＦＫＨに起因するピークが強調されて現れる。したがって、出力信号では、ＨＦ揺らぎ成分ＦＫＨの成分がノイズ成分に埋もれにくい。ゆえに、端末装置は、出力信号から、ＨＦ揺らぎ成分ＦＫＨの成分を精度良く抽出することができる。

【0148】

＜その他の実施形態＞
また、本発明に係る実施形態は、生体を撮像した画像データ以外のデータに基づいて第１信号が生成されてもよい。例えば、生体の一部に対して光を照射して、透過した光又は反射した光によって脈波を測定する装置等のような非接触式の測定装置が生成するデータが用いられてもよい。さらに、パッド等を生体に取り付けることで脈波を測定する装置等のような接触式の測定装置が生成するデータが用いられてもよい。

【0149】

また、本発明に係る実施形態は、周期性を有する信号であればどのような生体信号に対しても適用可能である。例えば、接触式の測定装置によって取得された心電信号に対して適用することができる。

【0150】

また、本発明に係る実施形態は、本発明の一実施形態に係る処理が、プログラムに基づいて情報処理装置又は情報処理システムによって実行されることで実現されてもよい。すなわち、本発明に係る実施形態は、生体信号解析方法をコンピュータに実行させるためのプログラム等によって実現されてもよい。なお、プログラムは、コンピュータが読み取り可能な記録媒体等に記憶されてコンピュータにインストールされる。

【0151】

さらに、本発明に係る実施形態は、１以上の情報処理装置を有する生体信号解析システムによって実現されてもよい。

【0152】

なお、生体信号解析システムは、処理を冗長、分散又は並列して行ってもよい。

【0153】

以上、本発明の好ましい実施例について詳述したが、本発明は係る特定の実施形態に限定されない。すなわち、特許請求の範囲に記載された本発明の要旨の範囲内において、種々の変形又は変更が可能である。

【符号の説明】

【0154】

１ 生体
１００ 端末装置
１００Ｈ１ 撮像装置
ＮＺ ノイズ成分
ＦＫＬ ＬＦ揺らぎ成分
ＦＫＨ ＨＦ揺らぎ成分
Ｌ１ 基本波
Ｌ２ ２次高調波
Ｌ３ ３次高調波

【先行技術文献】

【特許文献】

【0155】

【文献】国際公開第２０１４－０８４１６２号

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】

【図13】

【図14】

【図15】

【図16】

【図17】

【図18】

【図19】

【図20】

【図21】