信号処理装置、信号処理方法および生体情報測定装置
【課題】信号処理量を低減し消費電力をより抑えた、時系列信号からノイズ成分を除去する信号処理装置を提供することを目的とする。
【解決手段】信号処理装置は、周期性を有する第１信号成分と第１ノイズ成分とを含む第１の時系列信号と、前記第１信号成分と所定の関係を有する第２信号成分および前記第１ノイズ成分と所定の関係を有する第２ノイズ成分を含む第２の時系列信号と、前記第１ノイズ成分と前記第２ノイズ成分との所定の関係とに基づいて、基準となる１つの前記信号を生成し、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を前記所定の差分ずつ異ならせた他の複数の信号を、前記第１時系列信号と前記所定の差分とから算出したデータとを用いて漸化式によって生成し、第１の所定時間範囲での該生成信号の周期性を用いて、前記第１の所定時間範囲での前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を推定する。
【解決手段】信号処理装置は、周期性を有する第１信号成分と第１ノイズ成分とを含む第１の時系列信号と、前記第１信号成分と所定の関係を有する第２信号成分および前記第１ノイズ成分と所定の関係を有する第２ノイズ成分を含む第２の時系列信号と、前記第１ノイズ成分と前記第２ノイズ成分との所定の関係とに基づいて、基準となる１つの前記信号を生成し、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を前記所定の差分ずつ異ならせた他の複数の信号を、前記第１時系列信号と前記所定の差分とから算出したデータとを用いて漸化式によって生成し、第１の所定時間範囲での該生成信号の周期性を用いて、前記第１の所定時間範囲での前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を推定する。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、時系列信号からノイズ成分を除去する信号処理装置、信号処理方法および生体情報測定装置に関する。
【背景技術】
【０００２】
従来、ノイズ成分が重畳した時系列データからノイズ成分を除去する信号処理に関する技術が、様々な信号処理装置に応用されてきた。特に時系列データが生体情報に関する情報を含んでいる場合には、上記の信号処理装置は、生体情報測定装置と呼ばれている。生体情報測定装置は、生体組織から生体情報を非侵襲で検出する装置であり、具体的には光電脈波計と呼ばれる生体の脈波波形および脈拍数を測定する測定装置や、パルスオキシメータと呼ばれる動脈血中酸素飽和濃度を測定する測定装置等である。これらの測定装置の原理は、生体組織を透過または反射した光を受光することによって得られる、生体組織の脈動による変動分に対応した信号成分に基づいて、血中における吸光物質の濃度等の生体情報を求めるものである。
【０００３】
一般に、生体組織を透過または反射した光を受光することによって得られる、生体情報の検出に必要なデータには様々なノイズ成分が重畳されている。ノイズ成分は、主に、生体情報測定装置を使用している際に、生体が体を動かす等の体動を行うことによるものである。ノイズ成分が信号成分に重畳すると生体情報の算出において誤差要因となるため、ノイズ成分を除去することが望まれている。
【０００４】
互いに波長の異なる複数の光を生体にそれぞれ照射した場合に、生体組織を透過または反射した光の強度における直流交流比に基づいて、生体情報を算出する技術が提案されてきた。特に、信号成分にノイズ成分が重畳している場合には、各波長についての直流交流比は、信号成分とノイズ成分とで表される。このように表されたノイズ成分を除去する技術として、測定されたデータに基づいて生成された周期性を有する信号成分を含むデータから、周期性を用いることによって信号成分を抽出し、データからノイズ成分を除去する技術が提案されている（特許文献１等参照）。
【０００５】
この技術によれば、信号成分とノイズ成分との相互相関を用いた従来の手法に比べ、ノイズ成分除去のための演算処理量を低減することができる。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】国際公開第２０１０／０７３９０８号明細書
【発明の概要】
【発明が解決しようとする課題】
【０００７】
しかしながら、特許文献１の技術によって演算処理量が低減したとはいえ、乗算回数が比較的多いため、消費電力が大きくなってしまうことは否めない。特に、携帯用の生体情報測定装置では、通常、電池で駆動されるため、消費電力は重大な問題である。
【０００８】
また、生体情報測定装置は、手術室、集中治療室等の病棟のみならず、呼吸不全患者、在宅酸素療法患者の日常生活中の呼吸状態のデータ収集や管理、睡眠時無呼吸症候群のスクリーニング、登山等のスポーツ分野等常時身に着ける用途にまで用いられつつある。このような状況に鑑みても、生体情報測定装置には、小型化、軽量化、省電力化および低価格化等も求められている。
【０００９】
例として取りあげた生体情報測定装置に限らず、一般に、時系列データからノイズ成分を除去する信号処理を行う信号処理装置においても、ノイズ成分を除去する信号処理は消費電力の点で重大な問題といえる。
【００１０】
本発明は、上述の事情に鑑みて為された発明であり、その目的は、信号処理量を低減し消費電力を抑えた信号処理装置、信号処理方法および生体情報測定装置を提供することである。
【課題を解決するための手段】
【００１１】
本発明者は、種々検討した結果、上記目的は、以下の本発明により達成されることを見出した。
【００１２】
本発明の一態様に係る信号処理装置は、周期性を有する第１信号成分と、第１ノイズ成分とを含む第１の時系列信号と、前記第１信号成分と所定の関係を有する第２信号成分および前記第１ノイズ成分と所定の関係を有する第２ノイズ成分を含む第２の時系列信号と、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係とに基づいて、前記第１の時系列信号と前記第２の時系列信号とから前記第１信号成分を含む信号であって、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を所定の差分ずつ異ならせた複数の信号を生成する信号生成部と、第１の所定時間範囲での前記信号生成部で生成した各信号の周期性を用いて、前記第１の所定時間範囲での前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を推定する推定部とを備え、前記信号生成部は、基準となる１つの信号を生成し、当該信号と、前記第１時系列信号と前記所定の差分とから算出したデータとを用いた漸化式によって他の信号を生成することを特徴とする。
【００１３】
また、上述の信号処理装置において、前記信号生成部が、生成した信号に２値化処理を行うことで２値化信号をさらに生成し、前記推定部が、前記２値化信号の周期性を用いて前記信号生成部で生成した信号の周期を推定することを特徴とする。
【００１４】
そして、本発明の他の一態様に係る信号処理方法は、信号処理方法であって、第１の所定時間範囲での、周期性を有する第１信号成分と第１ノイズ成分とを含む第１の時系列信号と、前記第１信号成分と所定の関係を有する第２信号成分および前記第１ノイズ成分と所定の関係を有する第２ノイズ成分を含む第２の時系列信号と、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係とに基づいて、前記第１の時系列信号と前記第２の時系列信号とから前記第１信号成分を含む信号であって、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を所定の差分ずつ異ならせた複数の信号を生成する信号生成ステップと、第１の所定時間範囲での前記生成ステップで生成した各信号の周期性を用いて、前記第１の所定時間範囲での前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を推定する推定ステップとを備え、前記信号生成ステップは、基準となる１つの信号を生成し、当該信号と、前記第１時系列信号と前記所定の差分とから算出したデータとを用いた漸化式によって他の信号を生成することを特徴とする。
【００１５】
このような構成の信号処理装置では、前記周期性を用いることによってノイズ成分を除去する際に必要な信号処理において、漸化式を用いることで、乗除算よりも加減算を多用するので、結果的に従来技術より信号処理量を低減することができる。このため、信号処理量の低減により、必要とされるＣＰＵ性能を低く抑えることが可能となり、また、ノイズ成分除去の演算に伴う消費電力も抑制することが可能となる。
【００１６】
また、上述の信号処理装置において、前記信号生成部は、基準となる信号の閾値を算出し、当該閾値と、前記第１の時系列信号又は前記第２の時系列信号と前記所定の差分とを用いて算出したデータとを用いて、漸化式によって他の信号の閾値を算出することを特徴とする。
【００１７】
このような構成の信号処理装置によれば、２値化処理を行う際の閾値を、生成信号毎に設定でき、且つ、閾値を求める際の信号処理量が従来よりも低減できるので、より正確で迅速な２値化処理が可能となる。
【００１８】
また、上述の信号処理装置において、前記信号生成部が生成した信号は、時系列信号であり、前記信号生成部は、生成した信号から所定数毎に間引かれた結果のデータを用いて２値化処理を行うことを特徴とする。
【００１９】
このような構成の信号処理装置によれば、測定データのうちの一部のデータを用いて信号処理を行うので、より信号処理量を低減することができる。
【００２０】
また、上述の信号処理装置において、前記閾値を挟む２つの前記間引かれた結果のデータの間の間引いたデータに、前記２値化処理を行うことを特徴とする。
【００２１】
このような構成の信号処理装置によれば、閾値近辺では、測定した間隔のデータを用いて２値化するので、迅速であり、且つ、より正確な２値化処理が可能となる。
【００２２】
また、本発明の一態様に係る生体情報測定装置は、互いに波長の異なる複数の光を生体へそれぞれ照射して前記生体を透過または反射した各光をそれぞれ受光することによって得られた少なくとも第１測定データまたは第２測定データに基づいて、前記生体の生体情報を測定する生体情報測定装置であって、周期性を有する第１信号成分と、第１ノイズ成分からなる第１測定データを測定する第１測定部と、前記第１信号成分と所定の関係を有する第２信号成分および第１ノイズ成分と所定の関係を有する第２ノイズ成分を含む第２測定データを測定する第２測定部と、前記第１測定データと、前記第２測定データと、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係とに基づいて、前記第１測定データと前記第２測定データとから前記第１信号成分を含む信号であって、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を所定の差分ずつ異ならせた複数の信号を生成する信号生成部と、前記信号生成部が生成した信号の周期性を用いて、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を推定する推定部とを備え、前記信号生成部は、基準となる１つの信号を生成し、当該信号と、前記第１測定データと前記所定の差分とから算出したデータとを用いた漸化式によって他の信号を生成することを特徴とする。
【００２３】
このような構成の生体情報測定装置によれば、演算処理量を低減し消費電力を抑制した生体の血中酸素飽和度を測定する生体情報測定装置を提供することができる。また、演算処理量が低減されていることから、高性能な処理機能を備える必要が少なくなり、低価格化、小型化、軽量化の実現が可能となる。
【００２４】
また、上述の生体情報測定装置において、前記推定部が推定した所定の関係に基づいて、所定の周期で血中酸素飽和度を算出する血中酸素飽和度算出部を備え、前記生成信号部は、前記所定の周期よりも長い周期で、前記複数の信号を生成して、２値化信号を生成し、前記推定部は、前記２値化信号の周期性に基づいて、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を推定することを特徴とする。
【００２５】
測定データから血中酸素飽和度を求める周期よりも長い周期で、前記第１ノイズ成分と前記第２ノイズ成分との比である所定の関係k_vを算出する。k_vは、急激には変化しないからである。従って、複数の信号を生成し２値化する等の処理、データを測定するごとk_vを求めるよりも信号処理量が減ることになり、より迅速な信号処理量が可能となる。
【００２６】
また、上述の生体情報測定装置において、前記推定部が推定した所定の関係を有する信号に２値化処理を行うことで２値化信号を生成し、当該２値化信号の周期性を用いて脈拍数を算出する脈拍数算出部を備えることを特徴とする。
【００２７】
このような構成の信号処理装置では、前記周期性を用いることによってノイズ成分を除去し、心拍数を推定する際に必要な信号処理において、漸化式を用いることで、乗除算よりも加減算を多用するので、結果的に従来技術より信号処理量を低減することができる。
【発明の効果】
【００２８】
本発明に係る信号処理装置、信号処理方法および生体情報測定装置は、ノイズ成分除去のために必要な信号処理量を低減することができ、これにより、ノイズ成分除去の信号処理に伴う消費電力を抑制することが可能となる。
【図面の簡単な説明】
【００２９】
【図１】実施形態における生体情報測定装置の構成を示すブロック図である。
【図２】実施形態の生体情報測定装置における生体情報測定処理を示すフローチャートである。
【図３】実施形態の生体情報測定装置における初期処理を示すフローチャートである。
【図４】実施形態の生体情報測定装置におけるkv_min算出処理を示すフローチャートである。
【図５】測定データによるR-kv*IRの波形を示す図である。
【図６】実施形態の生体情報測定装置における２値化処理を示すフローチャートである。
【図７】図６の２値化処理を説明するための図である。
【図８】実施形態の生体情報測定装置におけるバラツキ指数算出処理を示すフローチャートである。
【図９】実施形態の生体情報測定装置における脈拍数PR算出処理を示すフローチャートである。
【図１０】変形例の生体情報測定装置における２値化処理を示すフローチャートである。
【図１１】変形例の生体情報測定装置における脈拍数PR算出処理を示すフローチャートである。
【発明を実施するための形態】
【００３０】
以下に、本発明の前提となる従来技術および実施形態について述べる。なお、便宜上、信号処理装置のうち、特に生体情報測定装置を例として取りあげるが、生体情報測定装置のみならず、ノイズ成分を除去する信号処理装置にも本発明は適用可能である。
＜本発明の前提となる従来技術＞
まず、本発明の前提となる従来技術について説明する。この説明をするにあたり、一例として、互いに波長IR、Rの異なる複数の光を生体へそれぞれ照射して前記生体を透過または反射した各光をそれぞれ受光することによって得られた各測定データに基づいて前記生体の生体情報として血中酸素飽和度を測定する場合について説明する。
【００３１】
ランバート・ベールの法則によって、生体組織を透過または反射したある波長の光の強度における交流成分と直流成分との比は、その波長での生体組織の吸光度の変化分に等しいと近似される。
【００３２】
上記ランバート・ベールの法則による近似を用いることによって、赤外波長IRについての、透過光または反射光の強度の直流成分と交流成分との比である赤外直交比IR_signalは、波長IRについての生体組織の吸光度の変化分と等しいと見なすことができる。同様に、赤色波長Rについての、透過光または反射光の強度の直流成分と交流成分との比である赤色直交比R_signalも、波長Rについての生体組織の吸光度の変化分と等しいと見なすことができる。
【００３３】
前記赤外直交比IR_signalは、式（１）のように表される。
【００３４】
【数１】
【００３５】
ここで、sは、吸光度の変化分の信号成分で、nは、信号成分に重畳しているノイズ成分である。
【００３６】
そして、前記赤色直交比R_signalは、式（２）で表される。
【００３７】
【数２】
【００３８】
ここで、k_aは、波長IRにおける吸光度の変化分の信号成分sと波長Rの吸光度の変化分の信号成分との比であり、k_vは、波長IRの信号成分に重畳したノイズ成分nと波長Rの信号成分に重畳したノイズ成分との比である。
【００３９】
この式（２）のk_aと動脈血中酸素飽和度とは、一対一に対応することが知られており、k_aを求めることによって、動脈血中酸素飽和度を求めることができる。
【００４０】
また、式（１）にk_vを乗算することによって式（３）が得られ、式（３）から式（２）を減算することによって、以下の式（４）が得られる。
【００４１】
【数３】
【００４２】
【数４】
【００４３】
同様に、上の式（１）にk_aを乗算することによって式（５）が得られ、式（５）から（式２）を減算することによって、以下の式（６）が得られる。
【００４４】
【数５】
【００４５】
【数６】
【００４６】
ここで、信号成分sとノイズ成分nは、独立であるという関係、すなわち以下の関係式（７）を用い、かつ、短い時間内では、k_aおよびk_vが一定という条件の下で、上の（４）式と（６）式の相関をとることによって、式（８）が得られる。
【００４７】
【数７】
【００４８】
【数８】
【００４９】
ここで、Σは、k_aおよびk_vが一定であるような短い時間に関しての総和である。iは、光の強度の変化分の時系列データIR_signal、R_signalのデータ番号であり、データの測定時間間隔をΔt、測定開始時刻をt0として、t＝Δt＊i+ｔ0という関係で時間tと結ばれている。
式（８）にはk_vとk_aという未知数が２つ含まれているため、式（８）のみからk_aとk_vとを求めることができない。
【００５０】
ここで、k_vを求めるために式（４）の右辺は、ほぼ周期的であることに着目し、式（４）の左辺が周期性をもつようなk_vを求める。
【００５１】
このように求められたk_vを式（８）に代入して、式（８）を満たす場合のk_aを求める。このk_aに基づいて、ノイズ成分を低減した動脈血中酸素飽和度を求めることができる。そして、上述の演算方法では、例えば、演算処理量が比較的大きなフーリエ変換等を用いる必要がないので、その分演算処理量を低減することが可能となっている。
＜実施形態＞
本実施形態は、式（４）の左辺が周期性をもつようなk_vを求める際の演算処理量を上記式をそのまま演算してk_vを求める場合の演算処理量に比べて大幅に削減するものである。
【００５２】
ここで、k_vを求める際の演算処理の削減方法について説明する。
【００５３】
式（４）の左辺に相当するq（以下の式参照）が最も周期性を持つようなk_vを求める。
【００５４】
q(i)=IR_signal(i)*k_v-R_signal(i)
ここで、iは、１からN個まである時系列データのi番目であることを表し、iと時間との対応は、上記＜本発明の前提となる従来技術＞で述べた通りである。Nは、動脈酸素飽和度の算出に必要なデータ数が選ばれ、例えば２００等の値が用いられる。
【００５５】
動脈血中酸素飽和度を求めるには、まず、qを求め、qを２値化してそのパルス周期のバラツキを表わす指標が最小になるk_vを決定する。そして、決定したk_vをkv_minとして（以下、パルス周期のバラツキを表わす指標が最小になるk_vを「kv_min」という。）、
ΣR_signal(i)2-kv_min*ΣR_signal(i)*IR_signal(i)
-k_a*{ΣR_signal(i)*IR_signal(i)-kv_min*ΣIR_signal(i)2}=0
を解いて、
k_a={ΣR_signal(i)2-kv_min*ΣR_signal(i)*IR_signal(i)}
/{ΣR_signal(i)*IR_signal(i)-kv_min*ΣIR_signal(i)2}
よりk_aを決定する。そして、k_aに基づいて血中酸素飽和度を算出する。Σはiに関して１からNまで取られる総和である（以下の式でも同様である。）。また、ΣR_signal(i)2、ΣIR_signal(i)2、ΣR_signal(i)*IR_signal(i)を、それぞれR_signal2、ΣIR_signal2、ΣR_signal*IR_signalと記載することがある。
【００５６】
k_a={ΣR_signal(i)2-kv_min*ΣR_signal(i)*IR_signal(i)}
/{ΣR_signal(i)*IR_signal(i)-kv_min*ΣIR_signal(i)2}
において、R_signal、IR_signalをそれぞれ時間差分に置き換えたものでも良い。
【００５７】
kv_minを求める方法は、次の通りである。
【００５８】
まず、k_vをΔk_vずつ変化させて、それぞれのk_vについてqを求める。
【００５９】
例えば、k_vをある値からΔk_vずつM回プラスし、また、M回マイナスして、(2*M+1)通りのqを算出する。この場合、各k_vについての演算では、N回の乗算を行うことになるので、すべてのk_vについてのqを求めるには、N*(2*M+1)回の乗算を行うことになる。
【００６０】
次に、それぞれのqにおける周期性を判断するために、それぞれのqを２値化する。２値化信号を生成する場合に用いる閾値は、各k_v応じて閾値を設定することが好ましい。qの振幅の大きさは各k_vによって異なるからである。
【００６１】
閾値は例えば、以下の式
閾値=f*√(1/N)Σq(i)2
より求めることができる。fは、例えば０．２である。
【００６２】
(1/N)Σq(i)2をQと表わすと
Q=(1/N)*Σq2
=k_v*k_v*(1/N)*ΣIR_signal2
-2*k_v*(1/N)*ΣR_signal*IR_signal+(1/N)*ΣR_signal2
である。
【００６３】
従って、(1/N)*ΣIR_signal2、(1/N)*ΣR_signal*IR_signal、(1/N)*ΣR_signal2は一度求めておけばよい。しかし、各k_vのQの算出には４回の乗算が必要であるので、すべてのk_vについては、4*(2*M+1)回の乗算が必要であることになり、計算量が非常に多くなる。
【００６４】
さらに、より精度を求めて、k_vを細かく変化させてq(i)のパルス幅のバラツキを調べる場合には、計算量が大幅に増えてしまうことになる。
【００６５】
そこで、本実施形態では計算量を削減するために、以下に示す手順でqおよびQを求める。
【００６６】
m番目(mは-M〜Mの範囲の整数)のk_vに対応するqおよびQは、q(m,i)、Q(m)と表し、m=0のときのq(0,i)、Q(0)、及び、mが0でない場合のQ(m)を求めるために必要なT(0)は、以下の式（Ａ１）〜（Ａ３）で求める。
【００６７】
q(0,i)=R_signal(i)-k_v(0)*IR_signal(i) ・・・（Ａ１）
Q(0)=k_v(0)2*SUM_sq_IR-2*k_v(0)*SUM_prdct_R_IR+SUM_sq_R ・・・（Ａ２）
T(0)=2*k_v(0)*dkv_SUM_sq_IR-2*dkv_SUM_prdct_R_IR ・・・（Ａ３）
なお、以下の値は事前に計算しておく。
【００６８】
SUM_sq_R=(1/N)*ΣR_signal(i)2 ・・・（Ａ４）
SUM_sq_IR=(1/N)*ΣIR_signal(i)2 ・・・（Ａ５）
SUM_prdct_R_IR=(1/N)*ΣR_signal(i)*IR_signal(i) ・・・（Ａ６）
次に、1≦m≦Mのとき(m=0から増加させる)のq及びQは、以下の式（Ｂ１）〜（Ｂ３）で求める。
【００６９】
q(m,i)=q(m-1,i)+dkv_IR(i) ・・・（Ｂ１）
Q(m)=Q(m-1)+T(m-1) ・・・（Ｂ２）
T(m)=T(m-1)+U ・・・（Ｂ３）
また、-M≦m≦-1のとき(m=0から減少させる)のq及びQは、以下の式（Ｂ４）〜（Ｂ６）で求める。
【００７０】
q(m,i)=q(m+1,i)-dkv_IR(i) ・・・（Ｂ４）
Q(m)=Q(m+1)-T(m+1)+U ・・・（Ｂ５）
T(m)=T(m+1)-U ・・・（Ｂ６）
なお、以下の値は事前に計算しておく。
【００７１】
dkv_IR(i)=Δk_v*IR_signal(i) ・・・（Ａ７）
U=2*sq_dkvSUM_sq_IR ・・・（Ａ８）
dkv_SUM_prdct_R_IR=Δk_v*SUM_prdct_R_IR ・・・（Ａ９）
dkv_SUM_sq_IR=Δk_v*SUM_sq_IR ・・・（Ａ１０）
sq_dkvSUM_sq_IR=Δk_v*dkv_SUM_sq_IR ・・・（Ａ１１）
なお、dkv_IR(i)はIR_signal(i)を算出するごとに算出しておいてもよい。
【００７２】
このように、式（Ａ１）を用いてq(0,i)、式（Ａ７）を用いてdkv_IR(i)を算出しておけば、q(m,i)の算出は式（Ｂ１）および式（Ｂ４）で示すように加減算のみで実行できる。
【００７３】
また、同様に式（Ａ２）を用いてQ(0)、式（Ａ３）を用いてT(0)、式（Ａ８）を用いてUを算出しておけば、Q(m)の算出は式（Ｂ２）および式（Ｂ３）、並びに、式（Ｂ５）および式（Ｂ６）で示すように加減算のみで実行できる。従って、計算量を大幅に少なくすることができる。
【００７４】
ここで、式（Ｂ２）および式（Ｂ３）、式（Ｂ５）および式（Ｂ６）について、説明する。
mを増加させるときのQ(m)(式Ｂ(２))は以下のように求めることができる。
Q(m)={k_v'(m)}2*SUM_sq_IR
-2*k_v'(m)*SUM_prdct_R_IR+SUM_sq_IR
={k_v'(m-1)+Δk_v}2*SUM_sq_IR
-2*{k_v'(m-1)+Δk_v}*SUM_prdct_R_IR+SUM_sq_IR
={k_v'(m-1)}2*SUM_sq_IR+2*k_v'(m-1)*dkv_SUM_sq_IR+sq_dkv_SUM_sq_IR
-2*k_v'(m-1)*SUM_prdct_R_IR-2*dkv_SUM_prdct_R_IR+SUM_sq_IR
ここで
{k_v'(m-1)}2*SUM_sq_IR-2*k_v'(m-1)*SUM_prdct_R_IR+SUM_sq_IR=Q(m-1)だから、
Q(m)=Q(m-1)+2*k_v'(m-1)*dkv_SUM_sq_IR+sq_dkv_SUM_sq_IR-2*dkv_SUM_prdct_R_IR
2*k_v'(m)*dkv_SUM_sq_IR+sq_dkv_SUM_sq_IR-2*dkv_SUM_prdct_R_IR=T(m)と定義すると、
2*k_v'(m-1)*dkv_SUM_sq_IR+sq_dkv_SUM_sq_IR-2*dkv_SUM_prdct_R_IR=T(m-1)だから
Q(m)=Q(m-1)+T(m-1)
また、
T(m)=2*k_v'(m)*dkv_SUM_sq_IR+sq_dkv_SUM_sq_IR-2*dkv_SUM_prdct_R_IR
=2*k_v'(m-1)*dkv_SUM_sq_IR+sq_dkv_SUM_sq_IR-2*dkv_SUM_prdct_R_IR
+2*sq_dkv_SUM_sq_IR
=T(m-1)+2*sq_dkv_SUM_sq_IR
である。
【００７５】
2*sq_dkv_SUM_sq_IR=Uと表わすとT(m)は次の漸化式で求めることができる。
T(m)=T(m-1)+U
したがって、以下の２つの漸化式によりQ(m)を順次求めることができる。
Q(m)=Q(m-1)+T(m-1)
T(m)=T(m-1)+U
次に、mを減少させるときのQ(m)(式Ｂ(５))は以下のように求めることができる。
mを増加させるときのT(m)=T(m-1)+Uより
T(m-1)=T(m)-U
だからこれをQ(m)=Q(m-1)+T(m-1)に代入すると
Q(m)=Q(m-1)+T(m)-U
Q(m-1)=Q(m)-T(m)+U
mをm+1に置き換えると、以下の漸化式でmを-1から順次減少させてQ(m)を求めることができる。
Q(m)=Q(m+1)-T(m+1)+U
T(m)=T(m+1)-U
以下、本発明に係る実施の一形態を図面に基づいて説明する。なお、各図において同一の符号を付した構成は、同一の構成であることを示し、適宜、その説明を省略する。
【００７６】
＜構成＞
まず、本発明の実施形態の構成について説明する。図１は、実施形態における生体情報測定装置の構成を示すブロック図である。なお、図中の実線は、後述する脈派の時系列データに相当する電気信号成分の各ブロック間での流れを表す。
【００７７】
この生体情報測定装置４０は、例えば体動等によるノイズ成分の影響を除去した生体情報を測定するために、測定対象（被験体）の生体を透過または反射した光の強度の時系列データを用いて測定対象の体動等による生体情報測定値への影響を除去し、ノイズ成分が除去された時系列データに基づいて、例えば血中酸素飽和度等の生体情報を測定する測定装置である。
【００７８】
本発明の一態様に即してより具体的に言えば、生体情報測定装置４０の生体情報測定装置本体は、互いに波長の異なる複数の光を生体へそれぞれ照射して生体を透過または反射した各光をそれぞれ受光することによって得られた各測定データに基づいて生体の生体情報を測定する。該生体情報測定装置本体は、測定データに基づいて生成され周期性を有する信号成分を含むデータから、信号成分の周期性を用いることによってデータからノイズ成分を除去し、ノイズ成分を除去したデータに基づいて生体の生体情報を測定する測定部１０を備える。
【００７９】
すなわち、生体情報測定装置４０は、互いに波長の異なる複数の光を生体へそれぞれ照射して生体を透過または反射した各光をそれぞれ受光することによって各測定データを取得するデータ取得部１と、データ取得部１で得られた各測定データに基づいて生体の生体情報を測定する該生体情報測定装置本体とを備える。
【００８０】
以下、より詳細に、この生体情報測定装置４０の構成について述べる。生体情報測定装置４０は、例えば、図１に示すように、データ取得部１と、測定部１０と、出力部３０とを備えて構成される。
【００８１】
データ取得部１は、測定部１０での生体情報の測定に必要な、生体の脈動に関する時系列データを、所定の時間間隔（測定周期）で測定して取得するための装置である。ここで、脈動を検出する方法としては、各種の方法が採用可能であるが、例えば生体組織のヘモグロビンの吸光特性を利用する方法を好適に採用することができる。周知の通り、酸素は、ヘモグロビンによって生体の各細胞に運ばれるが、ヘモグロビンは、肺で酸素と結合して酸化ヘモグロビンとなり、生体の細胞で酸素が消費されるとヘモグロビンに戻る。酸素飽和度は、血中の酸化ヘモグロビンの割合をいう。これらヘモグロビンの吸光度および酸化ヘモグロビンの吸光度は、波長依存性を有しており、例えば、ヘモグロビンは、赤色領域の波長Rの赤色光に対し酸化ヘモグロビンよりも光を多く吸収するが、赤外線領域の波長IRの赤外光に対しては酸化ヘモグロビンよりも光の吸収が少ない。この方法は、このようなヘモグロビンと酸化ヘモグロビンとの赤色光と赤外光とに対する吸光特性の違いを利用して例えば血中酸素飽和度や脈拍数等の生体情報を求めるものである。データ取得部１は、例えば、図１に示すように、赤色光（以下、R）を所定の生体組織に照射する発光素子（R）、および、前記発光素子（R）で照射され測定対象の生体組織を透過または反射した光を受光する受光素子（R）を備えたセンサ（R）部２と、赤外光（以下、IR）を前記所定の生体組織に照射する発光素子（IR）、および、前記発光素子（IR）で照射され測定対象の生体組織を透過または反射した各光を受光する受光素子（IR）を備えたセンサ（IR）部３とを含む反射型若しくは透過型センサである。このような構成のデータ取得部１は、所定の生体組織にセットされ、前記受光素子（R、IR）によって各受光量をそれぞれモニタして、これら受光された各光を光強度に従って電気信号へそれぞれ光電変換することによって脈波に関する前記各時系列データをそれぞれ取得する。データ取得部１は、測定部１０に接続され、これら各時系列データを後述するAC/DC(R)部１１へ出力する。
【００８２】
出力部３０は、測定部１０に接続され、測定部１０で測定された生体に関する生体情報の測定値等を出力するための装置である。出力部３０は、例えば、液晶表示装置（ＬＣＤ；Liquid Crystal Display）、７セグメントＬＥＤ、有機フォトルミネセンス表示装置、ＣＲＴ（Cathode Ray Tube）表示装置およびプラズマ表示装置等の表示装置や、マイクやスピーカー等の音出力装置や、プリンタ等の印刷装置等である。例えば、脈波データのデータ解析結果等の各種測定情報は、光点灯（点消灯、点滅を含む）、文字、画像、音声あるいは印刷等の適宜に任意の形態で出力される。
【００８３】
この出力部３０は、例えば、図１に示すように、SpO2出力部３１と、脈拍数出力部３２と、信頼度出力部３３とを備えて構成される。SpO2出力部３１は、後述するSvO2推定SpO2決定部１９で算出された血中酸素飽和度を出力するための装置である。脈拍数出力部３２は、後述する脈拍数算出部２０で算出された脈拍数を出力するための装置である。そして、信頼度出力部３３は、後述する信頼度算出部２１で算出された信頼度を出力するための装置である。
【００８４】
この測定部１０は、例えば、図１に示すように、機能的に、AC/DC(R)変換部１１、AC/DC(IR)変換部１２、BPF(R)部１３、BPF(IR)部１４、ΣR_signal*IR_signal算出部１５、ΣR2算出部１６、ΣIR2算出部１７、初期値算出部１８、SvO2推定SpO2決定部１９、脈拍数算出部２０、信頼度算出部２１および差分算出部２２を備え、そして、予め記憶された制御プログラムに従い、データ取得部１および出力部３０を当該機能に応じてそれぞれ制御する。
【００８５】
測定部１０は、データ取得部１によって取得された脈波の時系列データに対して、所定の前処理を行った後に、ノイズ成分の除去および生体情報の演算に必要な初期値を算出し、前記初期値を用いて、脈波の時系列データに含まれるノイズ成分を除去するノイズ成分除去処理を行い、ノイズ成分除去後の脈波の時系列データに基づき、生体情報等を測定し、その生体情報の測定値を出力部３０へ出力する装置である。測定部１０は、例えば、後述する演算処理を行う各演算処理プログラム等を記憶するＲＯＭ（Read Only Memory）、いわゆるワーキングメモリとして機能し一時的にデータを格納するＲＡＭ（Random Access Memory）および前記演算処理プログラム等を前記ＲＯＭから読み出して実行する中央処理装置（ＣＰＵ）およびその周辺回路等である。
【００８６】
AC/DC(R)変換部１１は、データ取得部１のセンサ（R）部２から入力される赤色光（R）の時系列データに対し、直流成分と交流成分との比であるデータに変換する回路である。なお、本実施形態では、AC/DC(R)変換部１１は、前記前処理として、暗電流による成分を除去するダーク処理を行う。
【００８７】
BPF(R)部１３は、AC/DC(R)変換部１１から赤色光（R）の時系列データが入力され、ノイズ成分を除去して赤色光（R）についての周波数成分を得るためのバンドパスフィルタ（ＢＰＦ）である。BPF(R)部１３は、フィルタリング後の時系列データであるR_signalをΣR_signal＊IR_signal算出部１５、Σ(R_signal)2算出部１６およびSvO2推定SpO2決定部１９へそれぞれ出力する。
【００８８】
AC/DC(IR)変換部１２は、データ取得部１のセンサ（IR）部３から入力される赤外光（IR）の時系列データに対し、直流成分と交流成分との比であるデータに変換する回路である。なお、本実施形態では、AC/DC(IR)変換部１２は、前記前処理として、暗電流による成分を除去するダーク処理を行う。
【００８９】
BPF(IR)部１４は、AC/DC(IR)変換部１２から赤外光（IR）の時系列データが入力され、ノイズ成分を除去して赤外光（IR）についての周波数成分を得るためのバンドパスフィルタである。BPF(IR)部１４は、フィルタリング後の時系列データであるIR_signalをΣR_signal＊IR_signal算出部１５、ΣIR_signal2算出部１７およびSvO2推定SpO2決定部１９へそれぞれ出力する。
【００９０】
ΣR_signal＊IR_signal算出部１５は、BPF(R)部１３から入力される赤色光(R)の時系列データと、BPF(IR)部１４から入力される赤外光(IR)の時系列データとを用いて、相互相関ΣR_signal＊IR_signalを算出し、算出した相互相関ΣR_signal＊IR_signalを初期値算出部１８および信頼度算出部２１に出力する回路である。
【００９１】
Σ(R_signal)2算出部１６は、BPF(R)部１３から入力される赤色光(R)の時系列データを用いて、自己相関Σ(R_signal)2を算出し、算出した自己相関Σ(R_signal)2を初期値算出部１８および信頼度算出部２１に出力する回路である。
【００９２】
Σ(IR_signal)2算出部１７は、BPF(IR)部１４から入力される赤外光(IR)の時系列データを用いて、自己相関Σ(IR_signal)2を算出し、算出した自己相関Σ(IR_signal)2を初期値算出部１８および信頼度算出部２１に出力する回路である。
【００９３】
差分算出部２２は、BPF(R)部１３から入力される赤色光(R)の時系列データと、BPF(IR)部１４から入力される赤外光(IR)の時系列データとを用いて、q(0,i)およびdkv_IR(i)（式（Ａ１）および式（Ａ７）参照）を算出し、算出したq(0,i)およびdkv_IR(i)をSvO2推定SpO2決定部１９に出力する回路である。
【００９４】
初期値算出部１８は、初期値および後の演算で使用する値を算出する回路であり、より具体的には、ΣR_signal＊IR_signal算出部１５から入力される相互相関ΣR_signal＊IR_signal、Σ(R_signal)2算出部１６から入力される自己相関Σ(R_signal)2、および、Σ(IR_signal)2算出部１７から入力される自己相関Σ(IR_signal)2を用いて、初期値および後の演算で使用する値（以下、まとめて「初期値」という場合がある。）を算出し、算出した初期値等をSvO2推定SpO2決定部１９へ出力する。
【００９５】
SvO2推定SpO2決定部１９は、周期性を有する信号成分を含む時系列データから、周期性を用いて信号成分に重畳したノイズ成分を除去し、ノイズ成分を除去した時系列データに基づいて、動脈の血中酸素飽和度を算出する回路である。SvO2推定SpO2決定部１９は、課題を解決するための手段で開示した信号生成部および推定部に相当する。SvO2推定SpO2決定部１９は、より具体的には、BPF(R)部１３から入力される時系列データR_signal、BPF(IR)部１４から入力される時系列データIR_signalおよび差分算出部２２から入力される時系列データq(0,i)およびdkv_IR(i)から、初期値算出部１８で算出された初期値を基に、動脈の血中酸素飽和度を算出する回路であり、必要に応じてノイズ成分に基づき静脈の血中酸素飽和度を算出し、算出した血中酸素飽和度をSpO2出力部３１へ出力する。
【００９６】
脈拍数算出部２０は、SvO2推定SpO2決定部１９でノイズ成分を除去した時系列データに基づいて、脈拍数を算出し、算出した脈拍数を脈拍数出力部３２へ出力する。
【００９７】
信頼度算出部２１は、測定した生体情報について、誤差の度合いを表す信頼度を算出する回路であり、算出した信頼度を信頼度出力部３３へ出力する。
【００９８】
なお、必要に応じて生体情報測定装置４０は、図略の外部記憶部をさらに備えてもよい。外部記憶部は、例えば、メモリカード、フレキシブルディスク、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Compact Disc Recordable）、ＤＶＤ−Ｒ（Digital Versatile Disc Recordable）およびブルーレイディスク（Blue-ray Disc）等の記憶媒体との間でデータを読み込みおよび／または書き込みを行う装置であり、例えば、メモリカードインタフェース、フレキシブルディスクドライブ、ＣＤ−ＲＯＭドライブ、ＣＤ−Ｒドライブ、ＤＶＤ−Ｒドライブおよびブルーレイディスクドライブ等である。
【００９９】
ここで、生体情報測定装置４０は、測定部１０に演算処理プログラム等が格納されていない場合には、これらプログラム等を記録した記録媒体から、前記外部記憶部を介して測定部１０にインストールされるように構成されてもよい。あるいは、生体情報測定装置４０は、検出された生体情報等のデータが前記外部記憶部を介して記録媒体に記録されるように構成されてもよい。
【０１００】
＜動作＞
次に、本実施形態の動作について、図２〜７を用いて説明する。
【０１０１】
図２の生体情報測定処理のフローチャート（ステップＳ１０〜ステップＳ２３）は、主に３つの処理を行っている。ステップＳ１０は初期処理を示し、ステップＳ１１〜ステップＳ２０は血中酸素飽和度検出処理を示し、ステップＳ２１は脈拍数検出処理を示す。ステップＳ２２は、算出した血中酸素飽和度、脈拍数等を表示する処理を示し、ステップＳ２３は、血中酸素飽和度等を測定周期で測定するよう制御する処理を示している。また、図３は、初期処理の詳細を示すフローチャートであり、図４〜図８は、血中酸素飽和度検出処理の中の一部の処理の詳細を示すフローチャートおよび説明のための図である。図９は、脈拍数検出処理の詳細を示すフローチャートである。
【０１０２】
生体情報測定装置４０は、例えば、その起動によって演算処理プログラムを実行する。この演算処理プログラムの実行によって、測定部１０の各部（１１〜２２）が機能的に構成される。そして、生体情報測定装置４０は、以下の動作によって、データ取得部１で取得された時系列データに基づいて、ノイズ成分を低減した例えば血中酸素飽和度等の生体情報を測定する。
【０１０３】
＜初期処理（ステップ１０）＞
生体情報測定装置４０は、まず初期処理を行う（ステップＳ１０）。初期処理では、主に、４つの処理を行う。１つは、所定数の時系列データを収集する処理である。２つ目は、k_vを求める処理で用いる値（後述するqの差分）等を算出して記憶する処理であり、３つ目は、演算で用いる固定値（SUMで始まる定数等）を算出して記憶する処理であり、４つ目は、初期値を算出する処理である。
【０１０４】
初期処理の詳細を、図３を用いて説明する。図３は、初期処理を示すフローチャートである。
【０１０５】
生体情報測定装置４０の測定が開始されると、測定部１０には、データ取得部１から、センサ（R）部２の暗電流による成分R_dark、および、センサ（IR）部３の暗電流による成分IR_darkがそれぞれ入力される。そして、測定部１０には、データ取得部１から、生体組織を透過または反射した波長IRおよび波長Rの光の強度変化に関するデータ時系列であるR_signal_and_dark(i)およびIR_signal_and_dark(i)がそれぞれ入力される。ここで、iは、１からN個まである時系列データのi番目であることを表し、iと時間との対応は、前記発明の原理で述べた通りである。Nは、酸素飽和度の算出に必要なデータ数が選ばれ、例えば２００等の値が用いられる。簡便のため、以降は、時系列データのiを省略して記載する場合がある。なお、波長Rに関するデータR_signal_and_darkには、センサ（R）部２の暗電流による成分R_darkが、波長IRに関するデータIR_signal_and_darkには、センサ（IR）部３の暗電流による成分IR_darkがそれぞれ含まれている。
【０１０６】
次に、AC/DC(R)部１１は、R_signal_and_darkからR_darkを差し引くダーク処理を行う。同様に、AC/DC(IR)部１２は、IR_signal_and_darkからIR_darkを差し引くダーク処理を行う（ステップＳ３０）。
【０１０７】
次に、ダーク処理された各波長の信号を直流成分と交流成分との比（交流成分／直流成分）にそれぞれ変換する（ステップＳ３１）。より具体的には、IR_signal_and_darkからIR_dark成分が除去され、生体由来の光の強度のみからなる直流成分と交流成分との比であるIR_signalが算出される。同様にR_signal_and_darkからR_dark成分が除去され、生体由来の光の強度のみからなる直流成分と交流成分との比であるR_signalが算出される。
【０１０８】
次に、BPF(R)部１３は、R_signalをフィルタリングし、不要な周波数成分を除去して所望の周波数成分のみを得るための処理を行う。同様にBPF(IR)部１４は、IR_signalをフィルタリングし、不要な周波数成分を除去して所望の周波数成分信号成分を得るための処理を行う（ステップＳ３２）。
【０１０９】
次に、差分算出部２２は、式（Ａ１）を用いて、q(0,i)を算出し、また、式（Ａ７）を用いて、IR_signal(i)からdkv_IR(i)を算出する（ステップＳ３３）。差分算出部２２は、算出したq(0,i)およびdkv_IR(i)を、SvO2推定SpO2決定部１９に出力する。
【０１１０】
また、ステップＳ３２で取り出されたデータR_signal、IR_signalに対し、ΣR_signal*IR_signal算出部１５、Σ(R_signal)2算出部１６、Σ(IR_signal)2算出部１７は、相互相関ΣR_signal＊IR_signal、波長Rについての自己相関ΣR_signal2および波長IRについての自己相関ΣIR_signal2をそれぞれを算出し、初期値算出部１８に出力する。なお、Σはiに関して１からNまで取られる総和であり、短い時間間隔で、k_vとk_aは、一定であるという仮定が成り立つ。
【０１１１】
次に、初期値算出部１８は、所定数（N個）のデータを収集したら（ステップＳ３４：Ｙｅｓ）、SUM_sq_R、SUM_sq_IR、および、SUM_prdct_R_IRを、式（Ａ４）、式（Ａ５）、および、式（Ａ６）を用いて算出する（ステップＳ３５）。
【０１１２】
次に、初期値算出部１８は、k_vの初期値であるk_v_0に、適当な値を代入する（ステップＳ３６）。なお、２度目の測定からは、前回求められたkv_minを記憶しておき、それを代入する。
【０１１３】
次に、ステップＳ３５およびステップＳ３６で求めたSUM_sq_R、SUM_sq_IR、SUM_prdct_R_IRおよびk_v_0を用いて、初期値算出部１８は、血中酸素飽和度検出フローで必要となる、k_aの仮値k_a_0およびk_nsの仮値であるk_ns_0をそれぞれ算出する（ステップＳ３７）。初期値算出部１８は、算出した各値を、SvO2推定SpO2決定部１９に出力する。
【０１１４】
SvO2推定SpO2決定部１９は、SUM_sq_R、SUM_sq_IR、SUM_prdct_R_IR、k_v_0、k_a_0、および、k_ns_0を前記メモリ等（図示せず）に記憶する。
【０１１５】
ここで、k_nsは、ノイズ成分を抽出するか否かを判定するための指標であり、ノイズ成分の強度と信号成分の強度との比等で表され、例えば、式（９）のように表される。
【０１１６】
【数９】
【０１１７】
ここで、ノイズ成分の強度に関する値であるn_squareと信号成分の強度に関する値であるs_squareは、式（８）にk_v_0を代入して等号が成り立つか否かによって異なる式で与えられる。
【０１１８】
式（８）にk_v_0を代入して、等号が成り立つ（すなわち左辺がゼロに等しくなる）ようなk_aの値が求められた場合には、そのk_aの値がk_a_0とされ、以下の式（１０）と式（１１）と（９）を用いてk_ns_0が算出される。
【０１１９】
【数１０】
【０１２０】
【数１１】
【０１２１】
一方、式（８）にk_v_0を代入し、等号が成り立つ生理学的に妥当なk_aが求まらない場合には、前回のk_aなど適当なk_aをk_a_0とし、以下の式（１２）と式（１３）と式（９）を用いてk_ns_0が算出される。また、動脈血と静脈血の酸素飽和度の差として例えば１０％と仮定して、それに相当する値をk_v_0から差し引いた値をk_a_0としても良い。
【０１２２】
【数１２】
【０１２３】
【数１３】
【０１２４】
＜血中酸素飽和度検出（ステップＳ１１〜ステップＳ２０）＞
上述の初期処理（図２のステップＳ１０）が終了した後、SvO2推定SpO2決定部１９によって血中酸素飽和度検出処理が実行される。
【０１２５】
SvO2推定SpO2決定部１９は、はじめに、ノイズ成分の指標であるk_ns_0と所定の値とを比較することにより、ノイズ成分状態を判定する（ステップＳ１１）。この判定結果により次の２通りに処理が分かれる。なお、所定の値は、除去すべきノイズ成分のレベルに応じて適宜に選定され、例えば、この実施例では0.05等である。
【０１２６】
ノイズ成分の指標が所定の値以下であるとSvO2推定SpO2決定部１９が判定した場合（ステップＳ１１：Ｎｏ）には、例えば、従来の技術として知られている、k_aを波長IRの吸光度の変化分の信号成分と波長Rの吸光度の変化分の信号成分との比の平均値として求める方法である式（１４）を用いてk_aが算出され（ステップＳ１２）、続いてステップＳ２０が実行される。
【０１２７】
【数１４】
【０１２８】
尚、式(８)から式(１４)はR_signal、IR_signalをそれぞれ時間差分に置き換えたものでも良い。
【０１２９】
一方、ノイズ成分の指標が所定の値より大きいとSvO2推定SpO2決定部１９が判定した場合（ステップＳ１１：Ｙｅｓ）には、kv_minが算出されたことがあるか、すなわち、生体情報測定装置４０が起動されて最初の測定処理か否かを判断する（ステップＳ１２）。
【０１３０】
この判断は、算出済フラグを参照して判断する。算出済フラグがＯＮであり「算出がされた」旨を示している場合には（ステップＳ１２：Ｙｅｓ）、保存しておいたkv_minを用いる（ステップＳ１５）。また、算出済フラグがＯＦＦであり「算出がされていない」旨を示している場合には（ステップＳ１２：Ｎｏ）、ステップＳ１６が実行され、kv_minが算出され、続いてステップＳ１８が実行される。kv_min算出処理は、後の＜kv_min算出処理＞の項で説明する。
【０１３１】
kv_minが求まると、メモリに記憶し、算出済フラグをＯＮにして「算出がされた」旨を設定する（ステップＳ１７）。
【０１３２】
次に、SvO2推定SpO2決定部１９は、ステップＳ１６で求めたkv_minを式（８）に代入してこの方程式を解くことでk_aを算出する（ステップＳ１８）。
【０１３３】
k_aを算出したSvO2推定SpO2決定部１９は、次に、式（９）を用いて、ノイズ成分の指標となるk_nsを算出する（ステップＳ１９）。k_nsの算出には、初期処理（ステップＳ１０）でメモリに記憶されたΣR_signal＊IR_signal、ΣR_signal2およびΣIR_signal2が用いられる。
【０１３４】
次に、SvO2推定SpO2決定部１９は、k_aの信頼度又はk_vの信頼度を算出し、k_aに基づき、血中酸素飽和度を算出する（ステップＳ２０）。ここで、k_aの信頼度とは、k_aに含まれる誤差の度合いを表す値である。同様に、k_vの信頼度とは、k_vに含まれる誤差の度合いを表す値である。これら信頼度の算出処理は、後の＜信頼度算出処理＞の項で説明する。
【０１３５】
SvO2推定SpO2決定部１９は、次回の測定に備えて、kv_minを記憶する。
【０１３６】
次に、脈拍数算出部２０は、脈拍数PRを算出する（ステップＳ２１）。脈拍数PR算出処理は、後の＜脈拍数PR算出処理＞の項で説明する。
【０１３７】
SpO2出力部３１、脈拍数出力部３２および信頼度出力部３３は、それぞれステップＳ２０およびステップＳ２１で検出した血中酸素飽和度と脈拍数PRとk_aの信頼度とを出力する（ステップＳ２２）。
【０１３８】
なお、信頼度出力部３３は、k_aの信頼度に代えてk_vの信頼度を出力してもよい。また、測定部１０は、ステップＳ１９で算出したk_nsと所定の値（ステップＳ１１の所定値）とを比較し、k_nsが所定の値（ステップＳ１１の所定値）を超えた場合には、k_nsが所定の値を超えた旨を警告するための警告をSpO2出力部３１に出力する。さらに、測定部１０は、k_aの信頼度と予め設定されているk_aの信頼度判定用の値とを比較し、k_aがその値を超えた場合には、k_aがその値を超えた旨を示す警告をさらに信頼度出力部３３に出力してもよい。また、測定部１０は、k_aの信頼度に代えてk_vの信頼度を用いて、k_vの信頼度と予め設定されているk_vの信頼度判定用の値とを比較し、k_vがその値を超えた場合には、k_vがその値を超えた旨を示す警告をさらに信頼度出力部３３に出力してもよい。
【０１３９】
測定部１０は、測定周期が経過したことを内部のタイマから検出するまで待ち（ステップＳ２３：Ｎｏ）、測定周期が経過したことを検出すると（ステップＳ２３：Ｙｅｓ）、k_v算出周期が経過したか否かを判断する（ステップＳ２４）。
【０１４０】
k_v算出周期が経過していた場合は（ステップＳ２４：Ｙｅｓ）、算出済みフラグをＯＦＦにし（ステップＳ２５）、経過していなかった場合は（ステップＳ２４：Ｎｏ）、算出済みフラグをＯＮのままにしておく。そして、ステップＳ１０から処理を開始する。
【０１４１】
ここで、k_v算出周期とは、kv_minを算出する周期であり、k_vが変化しないと想定される時間である。k_vは急激には変化しないことから、k_aおよびSpO2の算出周期よりもkv_minの算出周期を長くする。例えば、k_aおよびSpO2の算出が例えば１秒毎に実施されるとする（ステップ２３の測定周期に該当）。これに対して、kv_minを算出する周期を、k_aおよびSpO2の算出周期よりも長い周期(例えば２秒)とする（ステップＳ２４のk_v算出周期に該当）。このようにすることで、測定結果の正しさを変えることなく、kv_minの算出回数を減らして計算量を減らすことが可能となる。
【０１４２】
＜信頼度算出処理＞
一般に、算出された動脈血酸素飽和度の値の大小によって、それに含まれる誤差は異なることが知られており、例えば、動脈血酸素飽和度が約９５％の場合と７０％の場合とでは動脈血酸素飽和度の誤差の大きさが異なる。そこで、k_nsの所定値は、以下のように変化させる。k_nsの所定値は、動脈血酸素飽和度が所定の動脈血酸素飽和度より小さい場合（k_aとk_vとは大きく、静脈の血中酸素飽和度が小さい場合に相当)にはk_nsの所定値を小さく設定し、一方、逆に動脈血酸素飽和度が所定の動脈血酸素飽和度より大きい場合（k_aとk_vとは小さく、静脈の血中酸素飽和度が大きい場合に相当）にはk_nsの所定値を大きく設定するとよい。このように算出された動脈血酸素飽和度の値に応じて、k_nsの所定値を変化させることで、例えば、動脈血酸素飽和度が大きい場合は閾値を大きくし、動脈血酸素飽和度の誤差が動脈血酸素飽和度が小さい場合とほぼ同じ値で警告や出力禁止を行うことができる。ここで、所定の動脈血酸素飽和度および所定の動脈血酸素飽和度に応じたk_nsの所定値は、過去の測定データを元にして決定した値を測定部１０に予め記憶しておいてもよいし、生体や生体の状態に応じて適宜に変化させてもよい。求められたk_a又はk_vの信頼度の指標として前記k_ns以外に、例えば、式(１６)から式(２１)の何れかの式で与えられる値zを用いてもよい。zの絶対値が大きいと、k_aの信頼度は低い。また、zの絶対値が大きいと、k_vの信頼度は低い。
【０１４３】
【数１６】
【０１４４】
【数１７】
【０１４５】
【数１８】
【０１４６】
【数１９】
【０１４７】
【数２０】
【０１４８】
【数２１】
【０１４９】
さらに、k_aの信頼度の指標の代わりとして、以下の式(２２)で与えられる値wのバラツキ(標準偏差、最大値−最小値など)を用いてもよい。値wのバラツキが小さいことは、求めたk_aの信頼度は高いことを示している。
【０１５０】
【数２２】
【０１５１】
この他、R_signalをy、IR_signalをxとした場合における（x,y）の回帰直線に対するyのバラツキや最大値と最小値との差等からk_aの信頼性を評価してもよい。式(１６)から式(２２)はR_signal、IR_signalをそれぞれ時間差分に置き換えたものでも良い。
【０１５２】
＜kv_min算出処理＞
ステップＳ１６のkv_min算出処理について、図４を用いて説明する。
【０１５３】
上記＜本発明の前提となる従来技術＞で説明したように、式（４）の左辺に相当するq（以下の式（１５）参照）が最も強く周期性を持つように選択されたk_vが、最適なk_vである。
【０１５４】
【数１５】
【０１５５】
ここで、qは、あるk_vに対する、i番目のデータ系列である式（１５）の左辺についての値を表す。
【０１５６】
より具体的に、k_vの値に応じて、式（１５）は、どのように変化するかについて、図５を用いて以下に説明する。
【０１５７】
図５は、測定データによるR-kv*IRの波形を示す図である。図５（ａ）は、k_vを最適値より小さい値で与えた場合の、k_vに対する式（１５）で与えられるデータ系列による波形を示し、図５（ｂ）は、k_vを最適値に近似される値で与えた場合の、k_vに対する式（１５）で与えられるデータ系列による波形を示し、そして、図５（ｃ）は、k_vを最適値より大きい値で与えた場合の、k_vに対する式（１５）で与えられるデータ系列による波形を示す。なお、波形は、時系列データのデータ間で補間している。図５（ａ）から（ｃ）において、これら横軸は、所定の短い時間（データ系列の番号iでは１からNまでのデータ系列に相当）であり、縦軸は、式１５のqである。
【０１５８】
最適値より小さいk_vを式（１５）に与えた場合には、図５（ａ）に示す波形で、最適値に近似される値k_vを式（１５）に与えた場合には、図５（ｂ）に示す波形となっている。これら２つの図を比較すると分かるように、最適値に近似される値k_vが代入されている場合（図５（ｂ）では、信号成分にノイズ成分が重畳している場合でも波形に強い周期性が認められるが、最適値より小さいk_vを式（１５）に与えた場合には、信号成分にノイズ成分が重畳していると、波形は、乱れ、周期性が弱い。そして、最適値より大きいk_vを式（１５）に与えた場合には、図５（ｃ）に示す波形となっている。この図５（ｃ）を図５（ｂ）と比較すると分かるように、最適値より大きいk_vを式（１５）に与えた場合（図５（ｃ））には、信号成分にノイズ成分が重畳していると、波形は、乱れ、周期性が弱い。
【０１５９】
以上のk_vに対する式（１５）の変化を利用するため、式（１５）の周期性が検出される。
【０１６０】
図４は、kv_min算出処理のフローチャートである。
【０１６１】
実施形態のkv_min算出処理では、k_vの初期値をk_v_sとして、k_vを所定の範囲内（例えば、k_v_s-0.5＜k_v＜k_v_s+0.5）で、Δk_vずつ、プラス側とマイナス側にそれぞれM回ずつ変化させて、(2*M+1)個のq（式（１５）参照）を得る。そして、それぞれのqを２値化することによって、それぞれの波形の幅を算出する。
【０１６２】
以下、具体的に、kv_min算出処理を説明する。
【０１６３】
まず、SvO2推定SpO2決定部１９は、本処理における初期値を設定する（ステップＳ４０）。k_vを変化させた回数を示すカウンタmに０（ゼロ）を設定し、パルス幅のバラツキ指数の最小値を示す変数Var_PW_minに、バラツキ指数としてあり得ない値、例えば、１００を設定する。また、k_v_sに、初回は例えば１を、２回目からは前回のkv_minを設定する。
【０１６４】
次に、各qに適した閾値を算出する際に用いる値を算出し、メモリ等に記憶する（ステップＳ４１）。具体的には、式（Ａ２）を用いてQ(0)を、式（Ａ３）を用いてT(0)を算出し、式（Ａ８）を用いてUを算出する。これらの算出に必要な値は、初期処理で算出してメモリに記憶されているSUM_sq_R、SUM_sq_IRおよびSUM_prdct_R_IRを読み出し、式（Ａ９）を用いてdkv_SUM_prdct_R_IRを算出し、（Ａ１０）を用いてdkv_SUM_sq_IRを算出し、式（Ａ１１）を用いてsq_dkvSUM_sq_IRを算出する。
【０１６５】
次に、q(0,i)を、メモリから読み出す（ステップＳ４２）。なお、２度目の測定からは、後述のステップＳ６２においてk_v_minに対応するq_min(i)をストアしておき、q_min(i)を初期値q(0,i)とする。より、計算量の増加を防ぐことができる。
【０１６６】
次に、ステップＳ４３〜ステップＳ５０の処理、すなわち、k_vをΔk_vずつプラスしてバラツキを算出する処理をＭ回繰り返す。
【０１６７】
まず、Q(m)を、式（Ｂ２）で算出し、２値化処理で用いる閾値を算出する（ステップＳ４３）。本実施形態では、Q(m)から、以下の式で求められる閾値thresh_1と閾値thresh_2とを設定する。なお、f_1およびf_2は、例えば、それぞれ、０．２、−０．２である。
【０１６８】
Sqr_Q=Q(m)0.5
thresh_1=f_1*Sqr_Q
thresh_2=f_2*Sqr_Q
なお、ｍ＝０の場合は、ステップＳ４１で求めたQ(0)を用いる。
【０１６９】
次に、q(m,i)の２値化処理を行う（ステップＳ４４）。ｍ＝０の場合は、ステップＳ４２で求めたq(0,i)を用い、ｍ＞０の場合は、ステップＳ５０で求めたq(m,i)を用いる。
【０１７０】
次に、２値化処理によって得た情報から、バラツキ指数Var_PWを算出する（ステップＳ４４）。２値化処理については、後の＜２値化処理＞の項で説明し、バラツキ指数算出処理については、後の＜バラツキ指数算出処理＞の項で説明する。
【０１７１】
ステップＳ４５で求めたバラツキ指数Var_PWと、最小バラツキ指数Var_PW_minとを比較し、小さいほうの値を最小バラツキ指数Var_PW_minに設定する。また、最小バラツキ指数Var_PW_minが設定された場合には、設定されたVar_PWのk_vをkv_minに設定する（ステップＳ４６）。
【０１７２】
ステップＳ４３〜ステップＳ５０の処理をＭ回繰り返していない場合は（ステップＳ４７：Ｙｅｓ）、k_vにΔk_vをプラスして、次の処理で用いるk_vを求める（ステップＳ４８）。
【０１７３】
そして、式（Ｂ２）および式（Ｂ３）を用いて、Q(m)およびT(m)を算出する（ステップＳ４９）。また、式（Ｂ１）を用いて、q(m,i)を算出する（ステップＳ５０）。
【０１７４】
ステップＳ４３〜ステップＳ５０の処理をＭ回繰り返した場合には（ステップＳ４７：Ｎｏ）、Δk_vずつ、マイナス側に変化させて、バラツキを算出する処理に移る。
【０１７５】
まず、SvO2推定SpO2決定部１９は、k_vを変化させた回数を示すカウンタmに−１を設定し（ステップＳ５１）、式（Ｂ５）を用いてQ(-1)を、式（Ｂ６）を用いてT(-1)を算出し（ステップＳ５２）、式（Ｂ４）を用いてq(-1,i)を算出する（ステップＳ５３）。
【０１７６】
次に、ステップＳ５４〜ステップＳ６１の処理、すなわち、k_vをΔk_vずつマイナスしてバラツキを算出する処理をＭ回繰り返す。
【０１７７】
ステップＳ５４〜ステップＳ５７の処理は、ステップＳ４３〜ステップＳ４６の処理と同様である。
【０１７８】
ステップＳ５４〜ステップＳ６１の処理をＭ回繰り返していない場合は（ステップＳ５８：Ｙｅｓ）、k_vからΔk_vをマイナスして、次の処理で用いるk_vを求める（ステップＳ５９）。
【０１７９】
そして、式（Ｂ５）および式（Ｂ６）を用いて、Q(m)およびT(m)を算出する（ステップＳ６０）。また、式（Ｂ４）を用いて、q(m,i)を算出する（ステップＳ６１）。
【０１８０】
ステップＳ５４〜ステップＳ６１の処理をＭ回繰り返した場合には（ステップＳ５８：Ｎｏ）、最小バラツキ指数Var_PW_minの時のq(m,i)を、q_min(i)として記憶する（ステップＳ６２）。この時には最小バラツキ指数Var_PW_minには、最小のバラツキ指数Var_PWが設定されており、kv_minには、最小バラツキ指数Var_PW_minのk_vが設定されている。
【０１８１】
本実施形態では、固定値であるΔk_vを用いているが、まずは大きいΔk_vを用いて広いk_vの範囲で、２値化したq(m,i)のパルス幅のバラツキを最小にするkv_min_1を求め、次にkv_min_1の周辺のより狭いk_vの範囲で小さいΔk_vを用いて、２値化したq(m,i)のパルス幅のバラツキを最小にするkv_minを求めることとしてもよい。これにより、計算量を増やすことなく高い精度でkv_minを決定することが可能となる。
【０１８２】
＜２値化処理＞
ここで、ステップＳ４４の２値化処理について、図６及び図７を用いて説明する。
【０１８３】
図６は、２値化処理のフローチャートである。図７は、２値化処理を説明するためのｑの波形の拡大図である。
【０１８４】
本２値化処理では、ステップＳ４３で求めた閾値thresh_1と閾値thresh_2とを用いる。
【０１８５】
そして、q(m,i)が、閾値thresh_1を下回る値から超えた時のｉをテーブルt1_l_to_h(kr_1)に、閾値thresh_1を上回る値から下回った時のｉをテーブルt1_h_to_l(kf_1)に記憶する。同様に、q(m,i)が、閾値thresh_2を下回る値から超えた時のｉをテーブルt2_l_to_h(kr_2)に、閾値thresh_2を上回る値から下回った時のｉをテーブルt2_h_to_l(kf_2)に記憶する（図７参照）。kr_1、kf_1、kr_2及びkf_2は、０から始まるインデックスである。
【０１８６】
まず、２値化処理で用いる変数の初期化を行う（ステップＳ７０）。ここで、status_1は、q(m,i)が閾値thresh_1を超えている場合に１が設定され、以下の場合に０が設定されるフラグである。status_2は、q(m,i)が閾値thresh_2を超えている場合に１が設定され、以下の場合に０が設定されるフラグである。このフラグは、現在処理しているq(m,i)のひとつ前のq(m,i-1)の値の状態を示す。なお、閾値は、図４のステップＳ４３で算出されている。
【０１８７】
フラグstatus_1が０で且つq(m,i)が閾値thresh_1より大きい場合、すなわち、閾値thresh_1を下回る値から閾値thresh_1を超えた場合（ステップＳ７１：Ｙｅｓ）は、ｉをテーブルt1_l_to_h(kr_1)に記憶する（ステップＳ７２）。その他の場合は（ステップＳ７１：Ｎｏ）、何も記憶しない。
【０１８８】
フラグstatus_1が１で且つq(m,i)が閾値thresh_1以下である場合、すなわち、閾値thresh_1を上回る値から閾値thresh_1以下となった場合（ステップＳ７３：Ｙｅｓ）は、ｉをテーブルt1_h_to_l(kf_1)に記憶する（ステップＳ７４）。その他の場合は（ステップＳ７３：Ｎｏ）、何も記憶しない。
【０１８９】
そして、q(m,i)が閾値thresh_1を超えているか、以下かを判断し、フラグstatus_1に０又は１を設定する（ステップＳ７５）。
【０１９０】
次に同様に、フラグstatus_2が０で且つq(m,i)が閾値thresh_2より大きい場合、すなわち、閾値thresh_2を下回る値から閾値thresh_2を超えた場合（ステップＳ７６：Ｙｅｓ）は、ｉをテーブルt2_l_to_h(kr_2)に記憶する（ステップＳ７７）。その他の場合は（ステップＳ７７：Ｎｏ）、何も記憶しない。
【０１９１】
フラグstatus_2が１で且つq(m,i)が閾値thresh_2以下である場合、すなわち、閾値thresh_2を上回る値から閾値thresh_2以下となった場合（ステップＳ７８：Ｙｅｓ）は、ｉをテーブルt2_h_to_l(kf_2)に記憶する（ステップＳ７９）。その他の場合は（ステップＳ７８：Ｎｏ）、何も記憶しない。
【０１９２】
そして、q(m,i)が閾値thresh_2を超えているか、否かを判断し、フラグstatus_2に０又は１を設定する（ステップＳ８０）。
【０１９３】
ステップＳ７１〜ステップＳ８０の処理を、ｉが１〜Ｎまで繰り返し（ステップＳ８１、ステップＳ８２）、t1_l_to_h(kr_1)等のテーブルに記憶されたｉの数をkr_1_max等に記憶する（ステップＳ８３）。
【０１９４】
＜バラツキ指数算出処理＞
ここで、ステップＳ４５のバラツキ指数算出処理について、図８を用いて説明する。
【０１９５】
図８は、バラツキ指数算出処理のフローチャートである。
【０１９６】
本実施形態では、２値化処理で求めたテーブルt1_l_to_h(kr_1)、t1_h_to_l(kf_1)、t2_l_to_h(kr_2)およびt2_h_to_l(kf_2)から、バラツキ指数Var_PWを求める。
【０１９７】
具体的には、図７で示すように、すべての波形の山と谷のａ１、ａ２、ｂ１およびｂ３の長さを求め、ａ１、ａ２、ｂ１およびｂ３それぞれの最小値と最大値を求める。そして、その調和平均を求めてバラツキ指数とする。
【０１９８】
まず、テーブルt1_l_to_h(0)とテーブルt1_h_to_l(0)に登録されているｉの大きさを比べる。波形が山から始まっているのか、谷から始まっているのかを判断する為である。
【０１９９】
t1_l_to_h(0)＞t1_h_to_l(0)の場合は（ステップＳ９０：Ｙｅｓ）、谷から始まっていると判断し、テーブルt1_h_to_lのインデックスを１ずらして計算し、ａ１を求める。すべての山のａ１を算出し、その最大値t_1_h_maxと最小値t_1_h_minとを求める（ステップＳ９３）。
【０２００】
また、すべての山のｂ１を算出し、その最大値t_1_l_maxと最小値t_1_l_minとを求める（ステップＳ９４）。
【０２０１】
一方、t1_l_to_h(0)≦t1_h_to_l(0)の場合は（ステップＳ９０：Ｎｏ）、山から始まっていると判断し、すべての山のａ１を算出し、その最大値t_1_h_maxと最小値t_1_h_minとを求める（ステップＳ９１）。
【０２０２】
また、テーブルt1_h_to_lのインデックスを１ずらして計算し、すべての山のｂ１を算出し、その最大値t_1_l_maxと最小値t_1_l_minとを求める（ステップＳ９２）。
【０２０３】
次に、テーブルt2_l_to_h(0)とテーブルt2_h_to_l(0)に登録されているｉの大きさを比べる。波形が山から始まっているのか、谷から始まっているのかを判断する為である。
【０２０４】
t2_l_to_h(0)＞t2_h_to_l(0)の場合は（ステップＳ９５：Ｙｅｓ）、谷から始まっていると判断し、テーブルt1_h_to_lのインデックスを１ずらして計算し、ａ２を求める。すべての山のａ２を算出し、その最大値t_2_h_maxと最小値t_2_h_minとを求める（ステップＳ９８）。
【０２０５】
また、すべての山のｂ２を算出し、その最大値t_2_l_maxと最小値t_2_l_minとを求める（ステップＳ９９）。
【０２０６】
一方、t2_l_to_h(0)≦t2_h_to_l(0)の場合は（ステップＳ９５：Ｎｏ）、山から始まっていると判断し、すべての山のａ２を算出し、その最大値t_2_h_maxと最小値t_2_h_minとを求める（ステップＳ９６）。
【０２０７】
また、テーブルt2_h_to_lのインデックスを１ずらして計算し、すべての山のｂ２を算出し、その最大値t_2_l_maxと最小値t_2_l_minとを求める（ステップＳ９７）。
【０２０８】
次に、最大値t_1_l_max、最小値t_1_l_min、最大値t_1_h_max、最小値t_1_h_min、最大値t_2_l_max、最小値t_2_l_min、最大値t_2_h_max、および、最小値t_2_h_minからバラツキ指数Var_PWを求める（ステップＳ１００）。
【０２０９】
具体的には、まず、(t_1_l_max-t_1_l_min)と(t_2_l_max-t_2_l_min)との平均調和t_l_max_minを例えば、以下の式で求める。
【０２１０】
t_l_max_min = {(t_1_l_max-t_1_l_min)*(t_2_l_max-t_2_l_min)}
/{(t_1_l_max-t_1_l_min)+(t_2_l_max-t_2_l_min)}
また、同様に、(t_1_h_max-t_1_h_min)と(t_2_h_max-t_2_h_min)との平均調和t_h_max_minを求める。
【０２１１】
そして、t_l_max_minとt_h_max_minとの内、大きい方をバラツキ指数Var_PWとする。
【０２１２】
＜脈拍数PR算出処理（ステップＳ２１）＞
上述の血中酸素飽和度検出処理（図２のステップＳ１１〜ステップＳ２０）が終了した後、脈拍数算出部２０によって脈拍数PR算出処理が実行される。
【０２１３】
図９は、脈拍数PR算出処理のフローチャートである。
【０２１４】
実施形態では、図４のステップＳ６２で記憶したq_min(i)の状態変化の間隔の平均T_aveを算出し、脈拍数PRを求める。
【０２１５】
まず、q_min(i)に対し、２値化処理を行う（ステップＳ１１０〜ステップＳ１１８）。閾値を０（ゼロ）として、図６の２値化処理と同様に行う。この場合の閾値は、前記のk_vを決めるときに用いた閾値でも良い。
【０２１６】
詳細には、２値化処理で用いる変数の初期化を行う（ステップＳ１１０）。ここで、statusは、q_min(i)が閾値を超えている場合に１が設定され、閾値以下の場合に０が設定されるフラグである。このフラグは、現在処理しているq_min(i)のひとつ前のq_min(i-1)の値の状態を示す。
【０２１７】
フラグstatusが０で且つq_min(i)が閾値より大きい場合、すなわち、閾値を下回る値から閾値を超えた場合（ステップＳ１１１：Ｙｅｓ）は、ｉをテーブルt_l_to_h(kr)に記憶する（ステップＳ１１２）。その他の場合は（ステップＳ１１１：Ｎｏ）、何も記憶しない。
【０２１８】
フラグstatusが１で且つq_min(i)が閾値以下である場合、すなわち、閾値を上回る値から閾値以下となった場合（ステップＳ１１３：Ｙｅｓ）は、ｉをテーブルt_h_to_l(kf)に記憶する（ステップＳ１１４）。その他の場合は（ステップＳ１１３：Ｎｏ）、何も記憶しない。
【０２１９】
そして、q_min(i)が閾値を超えているか、閾値以下かを判断し、フラグstatusに０又は１を設定する（ステップＳ１１５）。
【０２２０】
ステップＳ１１１〜ステップＳ１１７の処理を、ｉが１〜Ｎまで繰り返し（ステップＳ１１６、ステップＳ１１７）、テーブルt_l_to_h(kr)およびt_h_to_l(kf)に記憶されたｉの数をkr_maxおよびkf_maxに記憶する（ステップＳ１１８）。
【０２２１】
次に、脈拍数算出部２０は、q_min(i)の所定時間内の周期T(j)の平均値T_aveを以下の式で算出する（ステップＳ１１９）。
T_ave＝((Σk=0kr_max-1|t_l_to_h(k)-t_l_to_h(k+1)|)
+(Σk=0kf_max-1|t_h_to_l(k)-t_h_to_l(k+1)|))/(kr_max+kf_max)
次に、脈拍数を算出する（ステップＳ１２０）。脈拍数は、前記周期の逆数として求められる。すなわち、脈拍数をPR(回／min)として、PR＝60／T_aveで表される。
＜変形例＞
実施形態では、kv_min算出処理（図４）において、q(m,i)の２値化処理を行う場合、測定データとして記憶しているすべてのデータ、すなわち、i=1〜i=Nまでのデータそれぞれについて２値化の処理を行っている。変形例では、このデータを間引いて、例えば、ｉが偶数のデータのみで２値化処理を行うことで、計算量を減らす。
【０２２２】
この場合、データを間引いたことによってパルス幅検出の精度を落とさないために、信号qが閾値をまたいだ場合には、またいだデータの間で間引いたデータについて２値化の処理を行い、閾値をまたいだ時刻を正確に取得するものとする。なお、データの間引きは、信号qのパルス幅のバラツキを検出できて、適切なkv_minを算出できる範囲で行うものとする。
【０２２３】
また、実施形態では、脈拍数PR算出処理（図９）において、q_min(i)から心拍数PRを求める場合、i=1〜i=Nまでのデータすべてを用いているが、変形例の２値化処理と同様に、このデータを間引いて、例えば、ｉが偶数のデータのみで２値化処理を行うことで、計算量を減らす。
【０２２４】
以下、変形例の２値化処理および脈拍数PR算出処理を説明する。
【０２２５】
＜２値化処理＞
図１０は、変形例の２値化処理のフローチャートである。ここでは、１つ置きにデータを間引くものとする。
【０２２６】
変形例の２値化処理は、図６を用いて説明した実施形態の２値化処理と、処理の流れはほぼ同様である。異なる点は、まず、q(m,i)のインデックスｉに２を加算する点（ステップＳ２１４）である。１つ置きのデータを処理するからである。また、実施形態では、q(m,i)が閾値をまたいだ時には、そのｉをテーブルに記憶している。しかし、変形例では、q(m,i)が閾値をまたいだ時には、その時飛ばしたデータq(m,i-1)が閾値をまたいでいるかを確認して、そのｉまたはi-1をテーブルに記憶する（ステップＳ２０２〜ステップＳ２０５、ステップＳ２０７〜ステップＳ２０９、ステップＳ２１３〜ステップＳ２１６、ステップＳ２１８〜ステップＳ２２１）。
【０２２７】
まず、２値化処理で用いる変数の初期化を行う（ステップＳ２００）。この変数は、図６と同じである。
【０２２８】
フラグstatus_1が０で且つq(m,i)が閾値thresh_1より大きい場合、すなわち、閾値thresh_1を下回る値から閾値thresh_1を超えた場合（ステップＳ２０１：Ｙｅｓ）は、q(m,i-1)の値を算出する（ステップＳ２０２）。すなわち、飛ばしたデータを算出する。
【０２２９】
そして、q(m,i-1)が閾値thresh_1より大きい場合、すなわち、飛ばしたデータが閾値thresh_1を下回る値から閾値thresh_1を超えた場合（ステップＳ２０３：Ｙｅｓ）は、ｉ-1をテーブルt1_l_to_h(kr_1)に記憶し（ステップＳ２０４）、超えない場合（ステップＳ２０３：Ｎｏ）は、ｉをテーブルt1_l_to_h(kr_1)に記憶する（ステップＳ２０５）。
【０２３０】
次に、フラグstatus_1が１で且つq(m,i)が閾値thresh_1以下である場合、すなわち、閾値thresh_1を上回る値から閾値thresh_1以下となった場合（ステップＳ２０６：Ｙｅｓ）も、q(m,i-1)の値を算出する（ステップＳ２０７）。
【０２３１】
そして、q(m,i-1)が閾値thresh_1以下の場合、すなわち、飛ばしたデータが閾値thresh_1を上回る値から閾値thresh_1以下となった場合（ステップＳ２０８：Ｙｅｓ）は、ｉ-1をテーブルt1_l_to_h(kr_1)に記憶し（ステップＳ２０９）、以下とならなかった場合（ステップＳ２０８：Ｎｏ）は、ｉをテーブルt1_h_to_l(kf_1)に記憶する（ステップＳ２１０）。
【０２３２】
そして、q(m,i)が閾値thresh_1を超えているか、以下かを判断し、フラグstatus_1に０又は１を設定する（ステップＳ２１１）。
【０２３３】
次に同様に、フラグstatus_2が０で且つq(m,i)が閾値thresh_2より大きい場合、すなわち、閾値thresh_2を下回る値から閾値thresh_2を超えた場合（ステップＳ２１２：Ｙｅｓ）は、q(m,i-1)の値を算出する（ステップＳ２１３）。すなわち、飛ばしたデータを算出する。
【０２３４】
そして、q(m,i-1)が閾値thresh_2より大きい場合、すなわち、飛ばしたデータが閾値thresh_2を下回る値から閾値thresh_2を超えた場合（ステップＳ２１４：Ｙｅｓ）は、ｉ-1をテーブルt2_l_to_h(kr_1)に記憶し（ステップＳ２１５）、超えない場合（ステップＳ２１４：Ｎｏ）は、ｉをテーブルt2_l_to_h(kr_1)に記憶する（ステップＳ２１６）。
【０２３５】
次に、フラグstatus_2が１で且つq(m,i)が閾値thresh_2以下である場合、すなわち、閾値thresh_2を上回る値から閾値thresh_2以下となった場合（ステップＳ２１７：Ｙｅｓ）も、q(m,i-1)の値を算出する（ステップＳ２１８）。
【０２３６】
そして、q(m,i-1)が閾値thresh_2より小さい場合、すなわち、飛ばしたデータが閾値thresh_2を上回る値から閾値thresh_2以下となった場合（ステップＳ２１９：Ｙｅｓ）は、ｉ-1をテーブルt2_h_to_l(kf_1)に記憶し（ステップＳ２２０）、以下とならなかった場合（ステップＳ２１９：Ｎｏ）は、ｉをテーブルt2_h_to_l(kf_1)に記憶する（ステップＳ２２１）。
【０２３７】
そして、q(m,i)が閾値thresh_2を超えているか、以下かを判断し、フラグstatus_2に０又は１を設定する（ステップＳ２２２）。
【０２３８】
ステップＳ２０１〜ステップＳ２２２の処理を、ｉが１〜Ｎまで繰り返し（ステップＳ２１３、ステップＳ２１４）、t1_l_to_h(kr_1)等のテーブルに記憶されたｉの数をkr_1_max等に記憶する（ステップＳ２１５）。
【０２３９】
＜脈拍数PR算出処理＞
図１１は、変形例の脈拍数PR算出処理のフローチャートである。
【０２４０】
図４のステップＳ６２で記憶したq_min(i)の状態変化の間隔の平均T_aveを算出し、脈拍数PRを求める。
【０２４１】
まず、q_min(i)に対し、２値化処理を行う（ステップＳ２２０〜ステップＳ２３４）。閾値を０（ゼロ）として、図１０の２値化処理と同様に行う。この場合の閾値は、前記のk_vを決めるときに用いた閾値でも良い。
【０２４２】
詳細には、２値化処理で用いる変数の初期化を行う（ステップＳ２２０）。この変数は、図９と同じである。
【０２４３】
フラグstatusが０で且つq_min(i)が閾値より大きい場合、すなわち、閾値を下回る値から閾値を超えた場合（ステップＳ２２１：Ｙｅｓ）は、q_min(i-1)の値を算出する（ステップＳ２２２）。すなわち、飛ばしたデータを算出する。
【０２４４】
そして、q_min(i-1)が閾値より大きい場合、すなわち、飛ばしたデータが閾値を下回る値から閾値を超えた場合（ステップＳ２２３：Ｙｅｓ）は、i-1をテーブルt_l_to_h(kr)に記憶し（ステップＳ２２４）、超えない場合（ステップＳ２２３：Ｎｏ）は、ｉをテーブルt_l_to_h(kr)に記憶する（ステップＳ２２５）。
【０２４５】
次に、フラグstatusが１で且つq_min(i)が閾値以下である場合、すなわち、閾値を上回る値から閾値以下となった場合（ステップＳ２２６：Ｙｅｓ）も、q_min(i-1)の値を算出する（ステップＳ２２７）。
【０２４６】
そして、q_min(i-1)が閾値以下の場合、すなわち、飛ばしたデータが閾値を上回る値から閾値以下となった場合（ステップＳ２２８：Ｙｅｓ）は、i-1をテーブルt_h_to_l(kf)に記憶し（ステップＳ２２９）、以下とならなかった場合（ステップＳ２２８：Ｎｏ）は、ｉをテーブルt_h_to_l(kf)に記憶する（ステップＳ２３０）。
【０２４７】
そして、q_min(i)が閾値を超えているか、以下かを判断し、フラグstatusに０又は１を設定する（ステップＳ２３１）。
【０２４８】
次に、q_min(i)の所定時間内の周期T(j)の平均値T_aveを、図９のステップＳ１１９と同様に算出する（ステップＳ２３５）。
【０２４９】
次に、脈拍数PRを、図９のステップＳ１２０と同様に算出する（ステップＳ２３６）。
【０２５０】
以上のような実施形態及び変形例によれば、漸化式を用いて周期性が最も良い信号成分を求めているので、従来技術より演算処理を簡単化することができ、演算処理量を低減することができる。変形例では、更に、より少ないデータに対して処理を行うため、より演算処理量を低減できる。そのため、演算処理に伴う消費電力を抑制することができる。この結果、生体情報測定装置４０に電池を使用することが可能となり、携帯に便利な小型軽量な生体情報測定装置を提供することができる。このため、例えば、登山者や在宅酸素療法患者が常時携帯しても体力的な負担を低下させることができる。
【０２５１】
また、本実施形態に係る生体情報測定装置４０は、信号成分の強度と前記ノイズ成分の強度との比に関する情報を有する指標に基づいて、ノイズ成分を抽出するか否かの判断を行うので、ノイズ成分が無視できるほど小さい場合においては、このノイズ成分の除去を行う演算を省略することができるので、生体情報を得るための演算時間が短縮され、消費電力を抑えて電池動作時間をより長くできるとともに、より早く生体情報を得ることができる。
【０２５２】
また、本実施形態に係る生体情報測定装置４０は、音声、文字および光等によって生体情報、信号信頼度およびノイズ信頼度を出力するので、これらを容易に確認することが可能となる。さらに、本実施形態に係る生体情報測定装置４０は、指標が所定の値を超えたことを示す警告、あるいは信号信頼度またはノイズ信頼度の低下を示す警告を出力するので、測定対象者や測定者等に生体情報の取り扱いに関して注意を喚起することもできる。そして、前記指標、前記ノイズ信頼度および前記信号信頼度に応じて、演算された生体情報の取り扱いを変更することが可能となり、例えば、指標を基に再計測等を行うことができる。
【０２５３】
また、本実施形態に係る生体情報測定装置４０は、脈拍数を算出する際に、２値化した波形を用いることで、例えばダブルピーク等による誤測定を低減することができる。
【０２５４】
また、本実施形態に係る生体情報測定方法は、同じ計算資源（ハードウェア資源）において、生体情報をより早い計算時間で算出することできる。
【０２５５】
上記実施形態では、一例として、周期性をもつ２波長の光の強度変化を検出して血中酸素飽和度を求める生体情報測定装置４０を例示した。この他、周期性をもつ１波長の光の強度変化に基づき、脈拍、不整脈（心房細動・期外収縮）の検出、除細動時のモニター、自律神経障害、血管年齢等の診断等を行う生体情報測定装置にも本発明を適用することができる。さらに、脈波波形以外に、心電等の生体情報を測定対象とする装置にも適用することができる。
【０２５６】
実施形態の２値化処理は、式（１５）が周期性を持つようなk_vを求めるための指標として、２値化処理後の波形の山と谷の周期を用いたが、これに限定されるものではない。２値化処理後の波形における、他の指標として、山が繰り返される周期のみ、谷が繰り返される周期のみを用いてもよい。また、山の周期の最大値‐山の周期の最小値、谷の周期の最大値‐谷の周期の最小値、山の周期の標準偏差、谷の周期の標準偏差、山の幅の最大値‐山の幅の最小値、谷の幅の最大値‐谷の幅の最小値、山の幅の標準偏差、谷の幅の標準偏差を用いても良い。さらに、指標として、R_signal - k_v * IR_signalの波形の極大値と極小値から求められる振幅の最大‐最小や標準偏差等を用いてもよい。
【０２５７】
実施形態では、閾値を２つ用いたが、３つ以上の閾値を設け、式（１５）に対応した信号がそれぞれの閾値を超えるか否かで２値化した３種以上の２値化パルスそれぞれについて前記の種々の周期性の指標を算出して、その平均値（単純平均、調和平均、相乗平均を含む）を最小にするk_vを最適値としてもよい。閾値の数を多くすると計算量は増えるが、周期性の評価の精度が向上する利点がある。
【０２５８】
いずれの指標においても、指標の標準偏差や指標の最大値と最小値との差等から、所定の時間内で、最も指標のばらつきが最小となるようなk_vを求めてもよい。複数の指標を用いた場合、それぞれの指標を最小にするk_vが異なる場合がある。その場合はそれぞれの指標を最小にするk_vの調和平均値を最適値とする。平均値としては単純平均値、相乗平均値でもよい。また、複数の指標を最小にする複数のk_vのうち前回のk_vに最も近いものを採用しても良い。また、複数の指標の平均値(単純平均値、調和平均値、相乗平均値など)を指標として評価しても良い。
【０２５９】
複数の指標に応じて複数のk_vを求めた場合には、各k_vに対して、式（８）の等号が成立するようなk_aをおのおの算出しても良い。複数のk_aから最終的なk_aを決める方法は、複数のk_aの平均値(単純平均値、調和平均値、相乗平均値など)を最終的なk_aとしてもよいし、複数のk_aのうち前回の最終的なk_aに最も近いものを今回の最終的なk_aとしてもよいし、複数のk_aのうち最大と最小を除いたものの平均値を最終的なk_aとしてもよい。さらに、前回のk_aとの差の絶対値に対応した重みを掛けて荷重平均しても良い。
【０２６０】
また、k_vは例えば２００組のR_signal、IR_signalについて前記の方法で決定し、k_aを求める際は２００組のデータを例えば４分割して４組の式(８)を解いて４個のk_aを算出してそれらの平均値(単純平均、調和平均、相乗平均値、前回k_aとの差の絶対値を重みとする荷重平均値など)を求めても良い。２００組のデータ全体で求めたk_aと前記の４個のk_aの平均値をさらに平均しても良い。
【０２６１】
この他、複数の指標からそれぞれに応じて複数のk_vが求められた場合には、連続するk_vの範囲の中央値、または直前の動脈血酸素飽和度演算で得られたk_vに最も近い値を最適値としてもよい。
【０２６２】
なお、今までの説明ではノイズ成分の指標が所定の値以下のときは式（１４）など従来の方法によりk_aを算出しているが、式（８）に前回のk_v(＝k_v_0)を代入してk_aを求めても良い。ノイズが小さい場合は式(８)の左辺は
Σ{IR_signal * k_v - R_signal} * {IR_signal * k_a - R_signal}
≒(k_v - k_a) * ΣIR_signal * {IR_signal * k_a - R_signal}
≒(k_v - k_a) * {ΣIR_signal2 * k_a - ΣIR_signal * R_signal}
なので、k_v としてk_aとは異なる値を代入すれば、方程式は
ΣIR_signal2 * k_a - ΣIR_signal * R_signal＝０
と等価である。この解は
k_a＝ΣIR_signal * R_signal ／ ΣIR_signal2
が得られる。ノイズが小さいときはR_signal≒真のk_a * IR_signalなので、上式の右辺は真のk_aになる。
【０２６３】
したがって、式(８)に前回のk_vなどk_aとは異なると思われる値を代入して解くことにより、真のk_aに十分近い値を得ることができる。
【０２６４】
また、ノイズが小さいときはR_signal − k_v * IR_signalは広い範囲のk_vに対して強い周期性を持つが、その場合でも前回のk_vとk_vの差の絶対値を前記の種々の指標に掛けた値を評価すれば、k_vの最適値はほぼ前回のk_vに近い値になるので、ノイズの指標が大きいときと同様の処理を行っても真のk_aに十分近い値を得ることができる。計算速度よりもプログラムの簡素化を重視する場合はノイズ指標の大きさによって処理を分けなくても良い。
【０２６５】
なお、周期性を利用する一例として、式（１５）の左辺が周期性を満たすようにk_vを求める例を挙げたが、式（１５）に限定されるものではなく、k_vが信号成分で表される数式であれば、式（１５）に代えて用いてもよい。
【０２６６】
本発明を表現するために、上述において図面を参照しながら実施形態を通して本発明を適切且つ十分に説明したが、当業者であれば上述の実施形態を変更および／または改良することは容易に為し得ることであると認識すべきである。したがって、当業者が実施する変更形態または改良形態が、請求の範囲に記載された請求項の権利範囲を離脱するレベルのものでない限り、当該変更形態または当該改良形態は、当該請求項の権利範囲に包括されると解釈される。
【符号の説明】
【０２６７】
１ データ取得部
２ センサ（R）部
３ センサ（IR）部
１０ 測定部
１１ AC/DC(R)変換部
１２ AC/DC(IR)変換部
１３ BPF(R)部
１４ BPF(IR)部
１５ ΣR_signal*IR_signal算出部
１６ Σ(R_signal)2算出部
１７ Σ(IR_signal)2算出部
１８ 初期値算出部
１９ SvO2推定SpO2決定部
２０ 脈拍数算出部
２１ 信頼度算出部
２２ 差分算出部
３０ 出力部
３１ SpO2出力部
３２ 脈拍数出力部
３３ 信頼度出力部
４０ 生体情報測定装置
【技術分野】
【０００１】
本発明は、時系列信号からノイズ成分を除去する信号処理装置、信号処理方法および生体情報測定装置に関する。
【背景技術】
【０００２】
従来、ノイズ成分が重畳した時系列データからノイズ成分を除去する信号処理に関する技術が、様々な信号処理装置に応用されてきた。特に時系列データが生体情報に関する情報を含んでいる場合には、上記の信号処理装置は、生体情報測定装置と呼ばれている。生体情報測定装置は、生体組織から生体情報を非侵襲で検出する装置であり、具体的には光電脈波計と呼ばれる生体の脈波波形および脈拍数を測定する測定装置や、パルスオキシメータと呼ばれる動脈血中酸素飽和濃度を測定する測定装置等である。これらの測定装置の原理は、生体組織を透過または反射した光を受光することによって得られる、生体組織の脈動による変動分に対応した信号成分に基づいて、血中における吸光物質の濃度等の生体情報を求めるものである。
【０００３】
一般に、生体組織を透過または反射した光を受光することによって得られる、生体情報の検出に必要なデータには様々なノイズ成分が重畳されている。ノイズ成分は、主に、生体情報測定装置を使用している際に、生体が体を動かす等の体動を行うことによるものである。ノイズ成分が信号成分に重畳すると生体情報の算出において誤差要因となるため、ノイズ成分を除去することが望まれている。
【０００４】
互いに波長の異なる複数の光を生体にそれぞれ照射した場合に、生体組織を透過または反射した光の強度における直流交流比に基づいて、生体情報を算出する技術が提案されてきた。特に、信号成分にノイズ成分が重畳している場合には、各波長についての直流交流比は、信号成分とノイズ成分とで表される。このように表されたノイズ成分を除去する技術として、測定されたデータに基づいて生成された周期性を有する信号成分を含むデータから、周期性を用いることによって信号成分を抽出し、データからノイズ成分を除去する技術が提案されている（特許文献１等参照）。
【０００５】
この技術によれば、信号成分とノイズ成分との相互相関を用いた従来の手法に比べ、ノイズ成分除去のための演算処理量を低減することができる。
【先行技術文献】
【特許文献】
【０００６】
【特許文献１】国際公開第２０１０／０７３９０８号明細書
【発明の概要】
【発明が解決しようとする課題】
【０００７】
しかしながら、特許文献１の技術によって演算処理量が低減したとはいえ、乗算回数が比較的多いため、消費電力が大きくなってしまうことは否めない。特に、携帯用の生体情報測定装置では、通常、電池で駆動されるため、消費電力は重大な問題である。
【０００８】
また、生体情報測定装置は、手術室、集中治療室等の病棟のみならず、呼吸不全患者、在宅酸素療法患者の日常生活中の呼吸状態のデータ収集や管理、睡眠時無呼吸症候群のスクリーニング、登山等のスポーツ分野等常時身に着ける用途にまで用いられつつある。このような状況に鑑みても、生体情報測定装置には、小型化、軽量化、省電力化および低価格化等も求められている。
【０００９】
例として取りあげた生体情報測定装置に限らず、一般に、時系列データからノイズ成分を除去する信号処理を行う信号処理装置においても、ノイズ成分を除去する信号処理は消費電力の点で重大な問題といえる。
【００１０】
本発明は、上述の事情に鑑みて為された発明であり、その目的は、信号処理量を低減し消費電力を抑えた信号処理装置、信号処理方法および生体情報測定装置を提供することである。
【課題を解決するための手段】
【００１１】
本発明者は、種々検討した結果、上記目的は、以下の本発明により達成されることを見出した。
【００１２】
本発明の一態様に係る信号処理装置は、周期性を有する第１信号成分と、第１ノイズ成分とを含む第１の時系列信号と、前記第１信号成分と所定の関係を有する第２信号成分および前記第１ノイズ成分と所定の関係を有する第２ノイズ成分を含む第２の時系列信号と、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係とに基づいて、前記第１の時系列信号と前記第２の時系列信号とから前記第１信号成分を含む信号であって、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を所定の差分ずつ異ならせた複数の信号を生成する信号生成部と、第１の所定時間範囲での前記信号生成部で生成した各信号の周期性を用いて、前記第１の所定時間範囲での前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を推定する推定部とを備え、前記信号生成部は、基準となる１つの信号を生成し、当該信号と、前記第１時系列信号と前記所定の差分とから算出したデータとを用いた漸化式によって他の信号を生成することを特徴とする。
【００１３】
また、上述の信号処理装置において、前記信号生成部が、生成した信号に２値化処理を行うことで２値化信号をさらに生成し、前記推定部が、前記２値化信号の周期性を用いて前記信号生成部で生成した信号の周期を推定することを特徴とする。
【００１４】
そして、本発明の他の一態様に係る信号処理方法は、信号処理方法であって、第１の所定時間範囲での、周期性を有する第１信号成分と第１ノイズ成分とを含む第１の時系列信号と、前記第１信号成分と所定の関係を有する第２信号成分および前記第１ノイズ成分と所定の関係を有する第２ノイズ成分を含む第２の時系列信号と、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係とに基づいて、前記第１の時系列信号と前記第２の時系列信号とから前記第１信号成分を含む信号であって、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を所定の差分ずつ異ならせた複数の信号を生成する信号生成ステップと、第１の所定時間範囲での前記生成ステップで生成した各信号の周期性を用いて、前記第１の所定時間範囲での前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を推定する推定ステップとを備え、前記信号生成ステップは、基準となる１つの信号を生成し、当該信号と、前記第１時系列信号と前記所定の差分とから算出したデータとを用いた漸化式によって他の信号を生成することを特徴とする。
【００１５】
このような構成の信号処理装置では、前記周期性を用いることによってノイズ成分を除去する際に必要な信号処理において、漸化式を用いることで、乗除算よりも加減算を多用するので、結果的に従来技術より信号処理量を低減することができる。このため、信号処理量の低減により、必要とされるＣＰＵ性能を低く抑えることが可能となり、また、ノイズ成分除去の演算に伴う消費電力も抑制することが可能となる。
【００１６】
また、上述の信号処理装置において、前記信号生成部は、基準となる信号の閾値を算出し、当該閾値と、前記第１の時系列信号又は前記第２の時系列信号と前記所定の差分とを用いて算出したデータとを用いて、漸化式によって他の信号の閾値を算出することを特徴とする。
【００１７】
このような構成の信号処理装置によれば、２値化処理を行う際の閾値を、生成信号毎に設定でき、且つ、閾値を求める際の信号処理量が従来よりも低減できるので、より正確で迅速な２値化処理が可能となる。
【００１８】
また、上述の信号処理装置において、前記信号生成部が生成した信号は、時系列信号であり、前記信号生成部は、生成した信号から所定数毎に間引かれた結果のデータを用いて２値化処理を行うことを特徴とする。
【００１９】
このような構成の信号処理装置によれば、測定データのうちの一部のデータを用いて信号処理を行うので、より信号処理量を低減することができる。
【００２０】
また、上述の信号処理装置において、前記閾値を挟む２つの前記間引かれた結果のデータの間の間引いたデータに、前記２値化処理を行うことを特徴とする。
【００２１】
このような構成の信号処理装置によれば、閾値近辺では、測定した間隔のデータを用いて２値化するので、迅速であり、且つ、より正確な２値化処理が可能となる。
【００２２】
また、本発明の一態様に係る生体情報測定装置は、互いに波長の異なる複数の光を生体へそれぞれ照射して前記生体を透過または反射した各光をそれぞれ受光することによって得られた少なくとも第１測定データまたは第２測定データに基づいて、前記生体の生体情報を測定する生体情報測定装置であって、周期性を有する第１信号成分と、第１ノイズ成分からなる第１測定データを測定する第１測定部と、前記第１信号成分と所定の関係を有する第２信号成分および第１ノイズ成分と所定の関係を有する第２ノイズ成分を含む第２測定データを測定する第２測定部と、前記第１測定データと、前記第２測定データと、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係とに基づいて、前記第１測定データと前記第２測定データとから前記第１信号成分を含む信号であって、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を所定の差分ずつ異ならせた複数の信号を生成する信号生成部と、前記信号生成部が生成した信号の周期性を用いて、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を推定する推定部とを備え、前記信号生成部は、基準となる１つの信号を生成し、当該信号と、前記第１測定データと前記所定の差分とから算出したデータとを用いた漸化式によって他の信号を生成することを特徴とする。
【００２３】
このような構成の生体情報測定装置によれば、演算処理量を低減し消費電力を抑制した生体の血中酸素飽和度を測定する生体情報測定装置を提供することができる。また、演算処理量が低減されていることから、高性能な処理機能を備える必要が少なくなり、低価格化、小型化、軽量化の実現が可能となる。
【００２４】
また、上述の生体情報測定装置において、前記推定部が推定した所定の関係に基づいて、所定の周期で血中酸素飽和度を算出する血中酸素飽和度算出部を備え、前記生成信号部は、前記所定の周期よりも長い周期で、前記複数の信号を生成して、２値化信号を生成し、前記推定部は、前記２値化信号の周期性に基づいて、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を推定することを特徴とする。
【００２５】
測定データから血中酸素飽和度を求める周期よりも長い周期で、前記第１ノイズ成分と前記第２ノイズ成分との比である所定の関係k_vを算出する。k_vは、急激には変化しないからである。従って、複数の信号を生成し２値化する等の処理、データを測定するごとk_vを求めるよりも信号処理量が減ることになり、より迅速な信号処理量が可能となる。
【００２６】
また、上述の生体情報測定装置において、前記推定部が推定した所定の関係を有する信号に２値化処理を行うことで２値化信号を生成し、当該２値化信号の周期性を用いて脈拍数を算出する脈拍数算出部を備えることを特徴とする。
【００２７】
このような構成の信号処理装置では、前記周期性を用いることによってノイズ成分を除去し、心拍数を推定する際に必要な信号処理において、漸化式を用いることで、乗除算よりも加減算を多用するので、結果的に従来技術より信号処理量を低減することができる。
【発明の効果】
【００２８】
本発明に係る信号処理装置、信号処理方法および生体情報測定装置は、ノイズ成分除去のために必要な信号処理量を低減することができ、これにより、ノイズ成分除去の信号処理に伴う消費電力を抑制することが可能となる。
【図面の簡単な説明】
【００２９】
【図１】実施形態における生体情報測定装置の構成を示すブロック図である。
【図２】実施形態の生体情報測定装置における生体情報測定処理を示すフローチャートである。
【図３】実施形態の生体情報測定装置における初期処理を示すフローチャートである。
【図４】実施形態の生体情報測定装置におけるkv_min算出処理を示すフローチャートである。
【図５】測定データによるR-kv*IRの波形を示す図である。
【図６】実施形態の生体情報測定装置における２値化処理を示すフローチャートである。
【図７】図６の２値化処理を説明するための図である。
【図８】実施形態の生体情報測定装置におけるバラツキ指数算出処理を示すフローチャートである。
【図９】実施形態の生体情報測定装置における脈拍数PR算出処理を示すフローチャートである。
【図１０】変形例の生体情報測定装置における２値化処理を示すフローチャートである。
【図１１】変形例の生体情報測定装置における脈拍数PR算出処理を示すフローチャートである。
【発明を実施するための形態】
【００３０】
以下に、本発明の前提となる従来技術および実施形態について述べる。なお、便宜上、信号処理装置のうち、特に生体情報測定装置を例として取りあげるが、生体情報測定装置のみならず、ノイズ成分を除去する信号処理装置にも本発明は適用可能である。
＜本発明の前提となる従来技術＞
まず、本発明の前提となる従来技術について説明する。この説明をするにあたり、一例として、互いに波長IR、Rの異なる複数の光を生体へそれぞれ照射して前記生体を透過または反射した各光をそれぞれ受光することによって得られた各測定データに基づいて前記生体の生体情報として血中酸素飽和度を測定する場合について説明する。
【００３１】
ランバート・ベールの法則によって、生体組織を透過または反射したある波長の光の強度における交流成分と直流成分との比は、その波長での生体組織の吸光度の変化分に等しいと近似される。
【００３２】
上記ランバート・ベールの法則による近似を用いることによって、赤外波長IRについての、透過光または反射光の強度の直流成分と交流成分との比である赤外直交比IR_signalは、波長IRについての生体組織の吸光度の変化分と等しいと見なすことができる。同様に、赤色波長Rについての、透過光または反射光の強度の直流成分と交流成分との比である赤色直交比R_signalも、波長Rについての生体組織の吸光度の変化分と等しいと見なすことができる。
【００３３】
前記赤外直交比IR_signalは、式（１）のように表される。
【００３４】
【数１】
【００３５】
ここで、sは、吸光度の変化分の信号成分で、nは、信号成分に重畳しているノイズ成分である。
【００３６】
そして、前記赤色直交比R_signalは、式（２）で表される。
【００３７】
【数２】
【００３８】
ここで、k_aは、波長IRにおける吸光度の変化分の信号成分sと波長Rの吸光度の変化分の信号成分との比であり、k_vは、波長IRの信号成分に重畳したノイズ成分nと波長Rの信号成分に重畳したノイズ成分との比である。
【００３９】
この式（２）のk_aと動脈血中酸素飽和度とは、一対一に対応することが知られており、k_aを求めることによって、動脈血中酸素飽和度を求めることができる。
【００４０】
また、式（１）にk_vを乗算することによって式（３）が得られ、式（３）から式（２）を減算することによって、以下の式（４）が得られる。
【００４１】
【数３】
【００４２】
【数４】
【００４３】
同様に、上の式（１）にk_aを乗算することによって式（５）が得られ、式（５）から（式２）を減算することによって、以下の式（６）が得られる。
【００４４】
【数５】
【００４５】
【数６】
【００４６】
ここで、信号成分sとノイズ成分nは、独立であるという関係、すなわち以下の関係式（７）を用い、かつ、短い時間内では、k_aおよびk_vが一定という条件の下で、上の（４）式と（６）式の相関をとることによって、式（８）が得られる。
【００４７】
【数７】
【００４８】
【数８】
【００４９】
ここで、Σは、k_aおよびk_vが一定であるような短い時間に関しての総和である。iは、光の強度の変化分の時系列データIR_signal、R_signalのデータ番号であり、データの測定時間間隔をΔt、測定開始時刻をt0として、t＝Δt＊i+ｔ0という関係で時間tと結ばれている。
式（８）にはk_vとk_aという未知数が２つ含まれているため、式（８）のみからk_aとk_vとを求めることができない。
【００５０】
ここで、k_vを求めるために式（４）の右辺は、ほぼ周期的であることに着目し、式（４）の左辺が周期性をもつようなk_vを求める。
【００５１】
このように求められたk_vを式（８）に代入して、式（８）を満たす場合のk_aを求める。このk_aに基づいて、ノイズ成分を低減した動脈血中酸素飽和度を求めることができる。そして、上述の演算方法では、例えば、演算処理量が比較的大きなフーリエ変換等を用いる必要がないので、その分演算処理量を低減することが可能となっている。
＜実施形態＞
本実施形態は、式（４）の左辺が周期性をもつようなk_vを求める際の演算処理量を上記式をそのまま演算してk_vを求める場合の演算処理量に比べて大幅に削減するものである。
【００５２】
ここで、k_vを求める際の演算処理の削減方法について説明する。
【００５３】
式（４）の左辺に相当するq（以下の式参照）が最も周期性を持つようなk_vを求める。
【００５４】
q(i)=IR_signal(i)*k_v-R_signal(i)
ここで、iは、１からN個まである時系列データのi番目であることを表し、iと時間との対応は、上記＜本発明の前提となる従来技術＞で述べた通りである。Nは、動脈酸素飽和度の算出に必要なデータ数が選ばれ、例えば２００等の値が用いられる。
【００５５】
動脈血中酸素飽和度を求めるには、まず、qを求め、qを２値化してそのパルス周期のバラツキを表わす指標が最小になるk_vを決定する。そして、決定したk_vをkv_minとして（以下、パルス周期のバラツキを表わす指標が最小になるk_vを「kv_min」という。）、
ΣR_signal(i)2-kv_min*ΣR_signal(i)*IR_signal(i)
-k_a*{ΣR_signal(i)*IR_signal(i)-kv_min*ΣIR_signal(i)2}=0
を解いて、
k_a={ΣR_signal(i)2-kv_min*ΣR_signal(i)*IR_signal(i)}
/{ΣR_signal(i)*IR_signal(i)-kv_min*ΣIR_signal(i)2}
よりk_aを決定する。そして、k_aに基づいて血中酸素飽和度を算出する。Σはiに関して１からNまで取られる総和である（以下の式でも同様である。）。また、ΣR_signal(i)2、ΣIR_signal(i)2、ΣR_signal(i)*IR_signal(i)を、それぞれR_signal2、ΣIR_signal2、ΣR_signal*IR_signalと記載することがある。
【００５６】
k_a={ΣR_signal(i)2-kv_min*ΣR_signal(i)*IR_signal(i)}
/{ΣR_signal(i)*IR_signal(i)-kv_min*ΣIR_signal(i)2}
において、R_signal、IR_signalをそれぞれ時間差分に置き換えたものでも良い。
【００５７】
kv_minを求める方法は、次の通りである。
【００５８】
まず、k_vをΔk_vずつ変化させて、それぞれのk_vについてqを求める。
【００５９】
例えば、k_vをある値からΔk_vずつM回プラスし、また、M回マイナスして、(2*M+1)通りのqを算出する。この場合、各k_vについての演算では、N回の乗算を行うことになるので、すべてのk_vについてのqを求めるには、N*(2*M+1)回の乗算を行うことになる。
【００６０】
次に、それぞれのqにおける周期性を判断するために、それぞれのqを２値化する。２値化信号を生成する場合に用いる閾値は、各k_v応じて閾値を設定することが好ましい。qの振幅の大きさは各k_vによって異なるからである。
【００６１】
閾値は例えば、以下の式
閾値=f*√(1/N)Σq(i)2
より求めることができる。fは、例えば０．２である。
【００６２】
(1/N)Σq(i)2をQと表わすと
Q=(1/N)*Σq2
=k_v*k_v*(1/N)*ΣIR_signal2
-2*k_v*(1/N)*ΣR_signal*IR_signal+(1/N)*ΣR_signal2
である。
【００６３】
従って、(1/N)*ΣIR_signal2、(1/N)*ΣR_signal*IR_signal、(1/N)*ΣR_signal2は一度求めておけばよい。しかし、各k_vのQの算出には４回の乗算が必要であるので、すべてのk_vについては、4*(2*M+1)回の乗算が必要であることになり、計算量が非常に多くなる。
【００６４】
さらに、より精度を求めて、k_vを細かく変化させてq(i)のパルス幅のバラツキを調べる場合には、計算量が大幅に増えてしまうことになる。
【００６５】
そこで、本実施形態では計算量を削減するために、以下に示す手順でqおよびQを求める。
【００６６】
m番目(mは-M〜Mの範囲の整数)のk_vに対応するqおよびQは、q(m,i)、Q(m)と表し、m=0のときのq(0,i)、Q(0)、及び、mが0でない場合のQ(m)を求めるために必要なT(0)は、以下の式（Ａ１）〜（Ａ３）で求める。
【００６７】
q(0,i)=R_signal(i)-k_v(0)*IR_signal(i) ・・・（Ａ１）
Q(0)=k_v(0)2*SUM_sq_IR-2*k_v(0)*SUM_prdct_R_IR+SUM_sq_R ・・・（Ａ２）
T(0)=2*k_v(0)*dkv_SUM_sq_IR-2*dkv_SUM_prdct_R_IR ・・・（Ａ３）
なお、以下の値は事前に計算しておく。
【００６８】
SUM_sq_R=(1/N)*ΣR_signal(i)2 ・・・（Ａ４）
SUM_sq_IR=(1/N)*ΣIR_signal(i)2 ・・・（Ａ５）
SUM_prdct_R_IR=(1/N)*ΣR_signal(i)*IR_signal(i) ・・・（Ａ６）
次に、1≦m≦Mのとき(m=0から増加させる)のq及びQは、以下の式（Ｂ１）〜（Ｂ３）で求める。
【００６９】
q(m,i)=q(m-1,i)+dkv_IR(i) ・・・（Ｂ１）
Q(m)=Q(m-1)+T(m-1) ・・・（Ｂ２）
T(m)=T(m-1)+U ・・・（Ｂ３）
また、-M≦m≦-1のとき(m=0から減少させる)のq及びQは、以下の式（Ｂ４）〜（Ｂ６）で求める。
【００７０】
q(m,i)=q(m+1,i)-dkv_IR(i) ・・・（Ｂ４）
Q(m)=Q(m+1)-T(m+1)+U ・・・（Ｂ５）
T(m)=T(m+1)-U ・・・（Ｂ６）
なお、以下の値は事前に計算しておく。
【００７１】
dkv_IR(i)=Δk_v*IR_signal(i) ・・・（Ａ７）
U=2*sq_dkvSUM_sq_IR ・・・（Ａ８）
dkv_SUM_prdct_R_IR=Δk_v*SUM_prdct_R_IR ・・・（Ａ９）
dkv_SUM_sq_IR=Δk_v*SUM_sq_IR ・・・（Ａ１０）
sq_dkvSUM_sq_IR=Δk_v*dkv_SUM_sq_IR ・・・（Ａ１１）
なお、dkv_IR(i)はIR_signal(i)を算出するごとに算出しておいてもよい。
【００７２】
このように、式（Ａ１）を用いてq(0,i)、式（Ａ７）を用いてdkv_IR(i)を算出しておけば、q(m,i)の算出は式（Ｂ１）および式（Ｂ４）で示すように加減算のみで実行できる。
【００７３】
また、同様に式（Ａ２）を用いてQ(0)、式（Ａ３）を用いてT(0)、式（Ａ８）を用いてUを算出しておけば、Q(m)の算出は式（Ｂ２）および式（Ｂ３）、並びに、式（Ｂ５）および式（Ｂ６）で示すように加減算のみで実行できる。従って、計算量を大幅に少なくすることができる。
【００７４】
ここで、式（Ｂ２）および式（Ｂ３）、式（Ｂ５）および式（Ｂ６）について、説明する。
mを増加させるときのQ(m)(式Ｂ(２))は以下のように求めることができる。
Q(m)={k_v'(m)}2*SUM_sq_IR
-2*k_v'(m)*SUM_prdct_R_IR+SUM_sq_IR
={k_v'(m-1)+Δk_v}2*SUM_sq_IR
-2*{k_v'(m-1)+Δk_v}*SUM_prdct_R_IR+SUM_sq_IR
={k_v'(m-1)}2*SUM_sq_IR+2*k_v'(m-1)*dkv_SUM_sq_IR+sq_dkv_SUM_sq_IR
-2*k_v'(m-1)*SUM_prdct_R_IR-2*dkv_SUM_prdct_R_IR+SUM_sq_IR
ここで
{k_v'(m-1)}2*SUM_sq_IR-2*k_v'(m-1)*SUM_prdct_R_IR+SUM_sq_IR=Q(m-1)だから、
Q(m)=Q(m-1)+2*k_v'(m-1)*dkv_SUM_sq_IR+sq_dkv_SUM_sq_IR-2*dkv_SUM_prdct_R_IR
2*k_v'(m)*dkv_SUM_sq_IR+sq_dkv_SUM_sq_IR-2*dkv_SUM_prdct_R_IR=T(m)と定義すると、
2*k_v'(m-1)*dkv_SUM_sq_IR+sq_dkv_SUM_sq_IR-2*dkv_SUM_prdct_R_IR=T(m-1)だから
Q(m)=Q(m-1)+T(m-1)
また、
T(m)=2*k_v'(m)*dkv_SUM_sq_IR+sq_dkv_SUM_sq_IR-2*dkv_SUM_prdct_R_IR
=2*k_v'(m-1)*dkv_SUM_sq_IR+sq_dkv_SUM_sq_IR-2*dkv_SUM_prdct_R_IR
+2*sq_dkv_SUM_sq_IR
=T(m-1)+2*sq_dkv_SUM_sq_IR
である。
【００７５】
2*sq_dkv_SUM_sq_IR=Uと表わすとT(m)は次の漸化式で求めることができる。
T(m)=T(m-1)+U
したがって、以下の２つの漸化式によりQ(m)を順次求めることができる。
Q(m)=Q(m-1)+T(m-1)
T(m)=T(m-1)+U
次に、mを減少させるときのQ(m)(式Ｂ(５))は以下のように求めることができる。
mを増加させるときのT(m)=T(m-1)+Uより
T(m-1)=T(m)-U
だからこれをQ(m)=Q(m-1)+T(m-1)に代入すると
Q(m)=Q(m-1)+T(m)-U
Q(m-1)=Q(m)-T(m)+U
mをm+1に置き換えると、以下の漸化式でmを-1から順次減少させてQ(m)を求めることができる。
Q(m)=Q(m+1)-T(m+1)+U
T(m)=T(m+1)-U
以下、本発明に係る実施の一形態を図面に基づいて説明する。なお、各図において同一の符号を付した構成は、同一の構成であることを示し、適宜、その説明を省略する。
【００７６】
＜構成＞
まず、本発明の実施形態の構成について説明する。図１は、実施形態における生体情報測定装置の構成を示すブロック図である。なお、図中の実線は、後述する脈派の時系列データに相当する電気信号成分の各ブロック間での流れを表す。
【００７７】
この生体情報測定装置４０は、例えば体動等によるノイズ成分の影響を除去した生体情報を測定するために、測定対象（被験体）の生体を透過または反射した光の強度の時系列データを用いて測定対象の体動等による生体情報測定値への影響を除去し、ノイズ成分が除去された時系列データに基づいて、例えば血中酸素飽和度等の生体情報を測定する測定装置である。
【００７８】
本発明の一態様に即してより具体的に言えば、生体情報測定装置４０の生体情報測定装置本体は、互いに波長の異なる複数の光を生体へそれぞれ照射して生体を透過または反射した各光をそれぞれ受光することによって得られた各測定データに基づいて生体の生体情報を測定する。該生体情報測定装置本体は、測定データに基づいて生成され周期性を有する信号成分を含むデータから、信号成分の周期性を用いることによってデータからノイズ成分を除去し、ノイズ成分を除去したデータに基づいて生体の生体情報を測定する測定部１０を備える。
【００７９】
すなわち、生体情報測定装置４０は、互いに波長の異なる複数の光を生体へそれぞれ照射して生体を透過または反射した各光をそれぞれ受光することによって各測定データを取得するデータ取得部１と、データ取得部１で得られた各測定データに基づいて生体の生体情報を測定する該生体情報測定装置本体とを備える。
【００８０】
以下、より詳細に、この生体情報測定装置４０の構成について述べる。生体情報測定装置４０は、例えば、図１に示すように、データ取得部１と、測定部１０と、出力部３０とを備えて構成される。
【００８１】
データ取得部１は、測定部１０での生体情報の測定に必要な、生体の脈動に関する時系列データを、所定の時間間隔（測定周期）で測定して取得するための装置である。ここで、脈動を検出する方法としては、各種の方法が採用可能であるが、例えば生体組織のヘモグロビンの吸光特性を利用する方法を好適に採用することができる。周知の通り、酸素は、ヘモグロビンによって生体の各細胞に運ばれるが、ヘモグロビンは、肺で酸素と結合して酸化ヘモグロビンとなり、生体の細胞で酸素が消費されるとヘモグロビンに戻る。酸素飽和度は、血中の酸化ヘモグロビンの割合をいう。これらヘモグロビンの吸光度および酸化ヘモグロビンの吸光度は、波長依存性を有しており、例えば、ヘモグロビンは、赤色領域の波長Rの赤色光に対し酸化ヘモグロビンよりも光を多く吸収するが、赤外線領域の波長IRの赤外光に対しては酸化ヘモグロビンよりも光の吸収が少ない。この方法は、このようなヘモグロビンと酸化ヘモグロビンとの赤色光と赤外光とに対する吸光特性の違いを利用して例えば血中酸素飽和度や脈拍数等の生体情報を求めるものである。データ取得部１は、例えば、図１に示すように、赤色光（以下、R）を所定の生体組織に照射する発光素子（R）、および、前記発光素子（R）で照射され測定対象の生体組織を透過または反射した光を受光する受光素子（R）を備えたセンサ（R）部２と、赤外光（以下、IR）を前記所定の生体組織に照射する発光素子（IR）、および、前記発光素子（IR）で照射され測定対象の生体組織を透過または反射した各光を受光する受光素子（IR）を備えたセンサ（IR）部３とを含む反射型若しくは透過型センサである。このような構成のデータ取得部１は、所定の生体組織にセットされ、前記受光素子（R、IR）によって各受光量をそれぞれモニタして、これら受光された各光を光強度に従って電気信号へそれぞれ光電変換することによって脈波に関する前記各時系列データをそれぞれ取得する。データ取得部１は、測定部１０に接続され、これら各時系列データを後述するAC/DC(R)部１１へ出力する。
【００８２】
出力部３０は、測定部１０に接続され、測定部１０で測定された生体に関する生体情報の測定値等を出力するための装置である。出力部３０は、例えば、液晶表示装置（ＬＣＤ；Liquid Crystal Display）、７セグメントＬＥＤ、有機フォトルミネセンス表示装置、ＣＲＴ（Cathode Ray Tube）表示装置およびプラズマ表示装置等の表示装置や、マイクやスピーカー等の音出力装置や、プリンタ等の印刷装置等である。例えば、脈波データのデータ解析結果等の各種測定情報は、光点灯（点消灯、点滅を含む）、文字、画像、音声あるいは印刷等の適宜に任意の形態で出力される。
【００８３】
この出力部３０は、例えば、図１に示すように、SpO2出力部３１と、脈拍数出力部３２と、信頼度出力部３３とを備えて構成される。SpO2出力部３１は、後述するSvO2推定SpO2決定部１９で算出された血中酸素飽和度を出力するための装置である。脈拍数出力部３２は、後述する脈拍数算出部２０で算出された脈拍数を出力するための装置である。そして、信頼度出力部３３は、後述する信頼度算出部２１で算出された信頼度を出力するための装置である。
【００８４】
この測定部１０は、例えば、図１に示すように、機能的に、AC/DC(R)変換部１１、AC/DC(IR)変換部１２、BPF(R)部１３、BPF(IR)部１４、ΣR_signal*IR_signal算出部１５、ΣR2算出部１６、ΣIR2算出部１７、初期値算出部１８、SvO2推定SpO2決定部１９、脈拍数算出部２０、信頼度算出部２１および差分算出部２２を備え、そして、予め記憶された制御プログラムに従い、データ取得部１および出力部３０を当該機能に応じてそれぞれ制御する。
【００８５】
測定部１０は、データ取得部１によって取得された脈波の時系列データに対して、所定の前処理を行った後に、ノイズ成分の除去および生体情報の演算に必要な初期値を算出し、前記初期値を用いて、脈波の時系列データに含まれるノイズ成分を除去するノイズ成分除去処理を行い、ノイズ成分除去後の脈波の時系列データに基づき、生体情報等を測定し、その生体情報の測定値を出力部３０へ出力する装置である。測定部１０は、例えば、後述する演算処理を行う各演算処理プログラム等を記憶するＲＯＭ（Read Only Memory）、いわゆるワーキングメモリとして機能し一時的にデータを格納するＲＡＭ（Random Access Memory）および前記演算処理プログラム等を前記ＲＯＭから読み出して実行する中央処理装置（ＣＰＵ）およびその周辺回路等である。
【００８６】
AC/DC(R)変換部１１は、データ取得部１のセンサ（R）部２から入力される赤色光（R）の時系列データに対し、直流成分と交流成分との比であるデータに変換する回路である。なお、本実施形態では、AC/DC(R)変換部１１は、前記前処理として、暗電流による成分を除去するダーク処理を行う。
【００８７】
BPF(R)部１３は、AC/DC(R)変換部１１から赤色光（R）の時系列データが入力され、ノイズ成分を除去して赤色光（R）についての周波数成分を得るためのバンドパスフィルタ（ＢＰＦ）である。BPF(R)部１３は、フィルタリング後の時系列データであるR_signalをΣR_signal＊IR_signal算出部１５、Σ(R_signal)2算出部１６およびSvO2推定SpO2決定部１９へそれぞれ出力する。
【００８８】
AC/DC(IR)変換部１２は、データ取得部１のセンサ（IR）部３から入力される赤外光（IR）の時系列データに対し、直流成分と交流成分との比であるデータに変換する回路である。なお、本実施形態では、AC/DC(IR)変換部１２は、前記前処理として、暗電流による成分を除去するダーク処理を行う。
【００８９】
BPF(IR)部１４は、AC/DC(IR)変換部１２から赤外光（IR）の時系列データが入力され、ノイズ成分を除去して赤外光（IR）についての周波数成分を得るためのバンドパスフィルタである。BPF(IR)部１４は、フィルタリング後の時系列データであるIR_signalをΣR_signal＊IR_signal算出部１５、ΣIR_signal2算出部１７およびSvO2推定SpO2決定部１９へそれぞれ出力する。
【００９０】
ΣR_signal＊IR_signal算出部１５は、BPF(R)部１３から入力される赤色光(R)の時系列データと、BPF(IR)部１４から入力される赤外光(IR)の時系列データとを用いて、相互相関ΣR_signal＊IR_signalを算出し、算出した相互相関ΣR_signal＊IR_signalを初期値算出部１８および信頼度算出部２１に出力する回路である。
【００９１】
Σ(R_signal)2算出部１６は、BPF(R)部１３から入力される赤色光(R)の時系列データを用いて、自己相関Σ(R_signal)2を算出し、算出した自己相関Σ(R_signal)2を初期値算出部１８および信頼度算出部２１に出力する回路である。
【００９２】
Σ(IR_signal)2算出部１７は、BPF(IR)部１４から入力される赤外光(IR)の時系列データを用いて、自己相関Σ(IR_signal)2を算出し、算出した自己相関Σ(IR_signal)2を初期値算出部１８および信頼度算出部２１に出力する回路である。
【００９３】
差分算出部２２は、BPF(R)部１３から入力される赤色光(R)の時系列データと、BPF(IR)部１４から入力される赤外光(IR)の時系列データとを用いて、q(0,i)およびdkv_IR(i)（式（Ａ１）および式（Ａ７）参照）を算出し、算出したq(0,i)およびdkv_IR(i)をSvO2推定SpO2決定部１９に出力する回路である。
【００９４】
初期値算出部１８は、初期値および後の演算で使用する値を算出する回路であり、より具体的には、ΣR_signal＊IR_signal算出部１５から入力される相互相関ΣR_signal＊IR_signal、Σ(R_signal)2算出部１６から入力される自己相関Σ(R_signal)2、および、Σ(IR_signal)2算出部１７から入力される自己相関Σ(IR_signal)2を用いて、初期値および後の演算で使用する値（以下、まとめて「初期値」という場合がある。）を算出し、算出した初期値等をSvO2推定SpO2決定部１９へ出力する。
【００９５】
SvO2推定SpO2決定部１９は、周期性を有する信号成分を含む時系列データから、周期性を用いて信号成分に重畳したノイズ成分を除去し、ノイズ成分を除去した時系列データに基づいて、動脈の血中酸素飽和度を算出する回路である。SvO2推定SpO2決定部１９は、課題を解決するための手段で開示した信号生成部および推定部に相当する。SvO2推定SpO2決定部１９は、より具体的には、BPF(R)部１３から入力される時系列データR_signal、BPF(IR)部１４から入力される時系列データIR_signalおよび差分算出部２２から入力される時系列データq(0,i)およびdkv_IR(i)から、初期値算出部１８で算出された初期値を基に、動脈の血中酸素飽和度を算出する回路であり、必要に応じてノイズ成分に基づき静脈の血中酸素飽和度を算出し、算出した血中酸素飽和度をSpO2出力部３１へ出力する。
【００９６】
脈拍数算出部２０は、SvO2推定SpO2決定部１９でノイズ成分を除去した時系列データに基づいて、脈拍数を算出し、算出した脈拍数を脈拍数出力部３２へ出力する。
【００９７】
信頼度算出部２１は、測定した生体情報について、誤差の度合いを表す信頼度を算出する回路であり、算出した信頼度を信頼度出力部３３へ出力する。
【００９８】
なお、必要に応じて生体情報測定装置４０は、図略の外部記憶部をさらに備えてもよい。外部記憶部は、例えば、メモリカード、フレキシブルディスク、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Compact Disc Recordable）、ＤＶＤ−Ｒ（Digital Versatile Disc Recordable）およびブルーレイディスク（Blue-ray Disc）等の記憶媒体との間でデータを読み込みおよび／または書き込みを行う装置であり、例えば、メモリカードインタフェース、フレキシブルディスクドライブ、ＣＤ−ＲＯＭドライブ、ＣＤ−Ｒドライブ、ＤＶＤ−Ｒドライブおよびブルーレイディスクドライブ等である。
【００９９】
ここで、生体情報測定装置４０は、測定部１０に演算処理プログラム等が格納されていない場合には、これらプログラム等を記録した記録媒体から、前記外部記憶部を介して測定部１０にインストールされるように構成されてもよい。あるいは、生体情報測定装置４０は、検出された生体情報等のデータが前記外部記憶部を介して記録媒体に記録されるように構成されてもよい。
【０１００】
＜動作＞
次に、本実施形態の動作について、図２〜７を用いて説明する。
【０１０１】
図２の生体情報測定処理のフローチャート（ステップＳ１０〜ステップＳ２３）は、主に３つの処理を行っている。ステップＳ１０は初期処理を示し、ステップＳ１１〜ステップＳ２０は血中酸素飽和度検出処理を示し、ステップＳ２１は脈拍数検出処理を示す。ステップＳ２２は、算出した血中酸素飽和度、脈拍数等を表示する処理を示し、ステップＳ２３は、血中酸素飽和度等を測定周期で測定するよう制御する処理を示している。また、図３は、初期処理の詳細を示すフローチャートであり、図４〜図８は、血中酸素飽和度検出処理の中の一部の処理の詳細を示すフローチャートおよび説明のための図である。図９は、脈拍数検出処理の詳細を示すフローチャートである。
【０１０２】
生体情報測定装置４０は、例えば、その起動によって演算処理プログラムを実行する。この演算処理プログラムの実行によって、測定部１０の各部（１１〜２２）が機能的に構成される。そして、生体情報測定装置４０は、以下の動作によって、データ取得部１で取得された時系列データに基づいて、ノイズ成分を低減した例えば血中酸素飽和度等の生体情報を測定する。
【０１０３】
＜初期処理（ステップ１０）＞
生体情報測定装置４０は、まず初期処理を行う（ステップＳ１０）。初期処理では、主に、４つの処理を行う。１つは、所定数の時系列データを収集する処理である。２つ目は、k_vを求める処理で用いる値（後述するqの差分）等を算出して記憶する処理であり、３つ目は、演算で用いる固定値（SUMで始まる定数等）を算出して記憶する処理であり、４つ目は、初期値を算出する処理である。
【０１０４】
初期処理の詳細を、図３を用いて説明する。図３は、初期処理を示すフローチャートである。
【０１０５】
生体情報測定装置４０の測定が開始されると、測定部１０には、データ取得部１から、センサ（R）部２の暗電流による成分R_dark、および、センサ（IR）部３の暗電流による成分IR_darkがそれぞれ入力される。そして、測定部１０には、データ取得部１から、生体組織を透過または反射した波長IRおよび波長Rの光の強度変化に関するデータ時系列であるR_signal_and_dark(i)およびIR_signal_and_dark(i)がそれぞれ入力される。ここで、iは、１からN個まである時系列データのi番目であることを表し、iと時間との対応は、前記発明の原理で述べた通りである。Nは、酸素飽和度の算出に必要なデータ数が選ばれ、例えば２００等の値が用いられる。簡便のため、以降は、時系列データのiを省略して記載する場合がある。なお、波長Rに関するデータR_signal_and_darkには、センサ（R）部２の暗電流による成分R_darkが、波長IRに関するデータIR_signal_and_darkには、センサ（IR）部３の暗電流による成分IR_darkがそれぞれ含まれている。
【０１０６】
次に、AC/DC(R)部１１は、R_signal_and_darkからR_darkを差し引くダーク処理を行う。同様に、AC/DC(IR)部１２は、IR_signal_and_darkからIR_darkを差し引くダーク処理を行う（ステップＳ３０）。
【０１０７】
次に、ダーク処理された各波長の信号を直流成分と交流成分との比（交流成分／直流成分）にそれぞれ変換する（ステップＳ３１）。より具体的には、IR_signal_and_darkからIR_dark成分が除去され、生体由来の光の強度のみからなる直流成分と交流成分との比であるIR_signalが算出される。同様にR_signal_and_darkからR_dark成分が除去され、生体由来の光の強度のみからなる直流成分と交流成分との比であるR_signalが算出される。
【０１０８】
次に、BPF(R)部１３は、R_signalをフィルタリングし、不要な周波数成分を除去して所望の周波数成分のみを得るための処理を行う。同様にBPF(IR)部１４は、IR_signalをフィルタリングし、不要な周波数成分を除去して所望の周波数成分信号成分を得るための処理を行う（ステップＳ３２）。
【０１０９】
次に、差分算出部２２は、式（Ａ１）を用いて、q(0,i)を算出し、また、式（Ａ７）を用いて、IR_signal(i)からdkv_IR(i)を算出する（ステップＳ３３）。差分算出部２２は、算出したq(0,i)およびdkv_IR(i)を、SvO2推定SpO2決定部１９に出力する。
【０１１０】
また、ステップＳ３２で取り出されたデータR_signal、IR_signalに対し、ΣR_signal*IR_signal算出部１５、Σ(R_signal)2算出部１６、Σ(IR_signal)2算出部１７は、相互相関ΣR_signal＊IR_signal、波長Rについての自己相関ΣR_signal2および波長IRについての自己相関ΣIR_signal2をそれぞれを算出し、初期値算出部１８に出力する。なお、Σはiに関して１からNまで取られる総和であり、短い時間間隔で、k_vとk_aは、一定であるという仮定が成り立つ。
【０１１１】
次に、初期値算出部１８は、所定数（N個）のデータを収集したら（ステップＳ３４：Ｙｅｓ）、SUM_sq_R、SUM_sq_IR、および、SUM_prdct_R_IRを、式（Ａ４）、式（Ａ５）、および、式（Ａ６）を用いて算出する（ステップＳ３５）。
【０１１２】
次に、初期値算出部１８は、k_vの初期値であるk_v_0に、適当な値を代入する（ステップＳ３６）。なお、２度目の測定からは、前回求められたkv_minを記憶しておき、それを代入する。
【０１１３】
次に、ステップＳ３５およびステップＳ３６で求めたSUM_sq_R、SUM_sq_IR、SUM_prdct_R_IRおよびk_v_0を用いて、初期値算出部１８は、血中酸素飽和度検出フローで必要となる、k_aの仮値k_a_0およびk_nsの仮値であるk_ns_0をそれぞれ算出する（ステップＳ３７）。初期値算出部１８は、算出した各値を、SvO2推定SpO2決定部１９に出力する。
【０１１４】
SvO2推定SpO2決定部１９は、SUM_sq_R、SUM_sq_IR、SUM_prdct_R_IR、k_v_0、k_a_0、および、k_ns_0を前記メモリ等（図示せず）に記憶する。
【０１１５】
ここで、k_nsは、ノイズ成分を抽出するか否かを判定するための指標であり、ノイズ成分の強度と信号成分の強度との比等で表され、例えば、式（９）のように表される。
【０１１６】
【数９】
【０１１７】
ここで、ノイズ成分の強度に関する値であるn_squareと信号成分の強度に関する値であるs_squareは、式（８）にk_v_0を代入して等号が成り立つか否かによって異なる式で与えられる。
【０１１８】
式（８）にk_v_0を代入して、等号が成り立つ（すなわち左辺がゼロに等しくなる）ようなk_aの値が求められた場合には、そのk_aの値がk_a_0とされ、以下の式（１０）と式（１１）と（９）を用いてk_ns_0が算出される。
【０１１９】
【数１０】
【０１２０】
【数１１】
【０１２１】
一方、式（８）にk_v_0を代入し、等号が成り立つ生理学的に妥当なk_aが求まらない場合には、前回のk_aなど適当なk_aをk_a_0とし、以下の式（１２）と式（１３）と式（９）を用いてk_ns_0が算出される。また、動脈血と静脈血の酸素飽和度の差として例えば１０％と仮定して、それに相当する値をk_v_0から差し引いた値をk_a_0としても良い。
【０１２２】
【数１２】
【０１２３】
【数１３】
【０１２４】
＜血中酸素飽和度検出（ステップＳ１１〜ステップＳ２０）＞
上述の初期処理（図２のステップＳ１０）が終了した後、SvO2推定SpO2決定部１９によって血中酸素飽和度検出処理が実行される。
【０１２５】
SvO2推定SpO2決定部１９は、はじめに、ノイズ成分の指標であるk_ns_0と所定の値とを比較することにより、ノイズ成分状態を判定する（ステップＳ１１）。この判定結果により次の２通りに処理が分かれる。なお、所定の値は、除去すべきノイズ成分のレベルに応じて適宜に選定され、例えば、この実施例では0.05等である。
【０１２６】
ノイズ成分の指標が所定の値以下であるとSvO2推定SpO2決定部１９が判定した場合（ステップＳ１１：Ｎｏ）には、例えば、従来の技術として知られている、k_aを波長IRの吸光度の変化分の信号成分と波長Rの吸光度の変化分の信号成分との比の平均値として求める方法である式（１４）を用いてk_aが算出され（ステップＳ１２）、続いてステップＳ２０が実行される。
【０１２７】
【数１４】
【０１２８】
尚、式(８)から式(１４)はR_signal、IR_signalをそれぞれ時間差分に置き換えたものでも良い。
【０１２９】
一方、ノイズ成分の指標が所定の値より大きいとSvO2推定SpO2決定部１９が判定した場合（ステップＳ１１：Ｙｅｓ）には、kv_minが算出されたことがあるか、すなわち、生体情報測定装置４０が起動されて最初の測定処理か否かを判断する（ステップＳ１２）。
【０１３０】
この判断は、算出済フラグを参照して判断する。算出済フラグがＯＮであり「算出がされた」旨を示している場合には（ステップＳ１２：Ｙｅｓ）、保存しておいたkv_minを用いる（ステップＳ１５）。また、算出済フラグがＯＦＦであり「算出がされていない」旨を示している場合には（ステップＳ１２：Ｎｏ）、ステップＳ１６が実行され、kv_minが算出され、続いてステップＳ１８が実行される。kv_min算出処理は、後の＜kv_min算出処理＞の項で説明する。
【０１３１】
kv_minが求まると、メモリに記憶し、算出済フラグをＯＮにして「算出がされた」旨を設定する（ステップＳ１７）。
【０１３２】
次に、SvO2推定SpO2決定部１９は、ステップＳ１６で求めたkv_minを式（８）に代入してこの方程式を解くことでk_aを算出する（ステップＳ１８）。
【０１３３】
k_aを算出したSvO2推定SpO2決定部１９は、次に、式（９）を用いて、ノイズ成分の指標となるk_nsを算出する（ステップＳ１９）。k_nsの算出には、初期処理（ステップＳ１０）でメモリに記憶されたΣR_signal＊IR_signal、ΣR_signal2およびΣIR_signal2が用いられる。
【０１３４】
次に、SvO2推定SpO2決定部１９は、k_aの信頼度又はk_vの信頼度を算出し、k_aに基づき、血中酸素飽和度を算出する（ステップＳ２０）。ここで、k_aの信頼度とは、k_aに含まれる誤差の度合いを表す値である。同様に、k_vの信頼度とは、k_vに含まれる誤差の度合いを表す値である。これら信頼度の算出処理は、後の＜信頼度算出処理＞の項で説明する。
【０１３５】
SvO2推定SpO2決定部１９は、次回の測定に備えて、kv_minを記憶する。
【０１３６】
次に、脈拍数算出部２０は、脈拍数PRを算出する（ステップＳ２１）。脈拍数PR算出処理は、後の＜脈拍数PR算出処理＞の項で説明する。
【０１３７】
SpO2出力部３１、脈拍数出力部３２および信頼度出力部３３は、それぞれステップＳ２０およびステップＳ２１で検出した血中酸素飽和度と脈拍数PRとk_aの信頼度とを出力する（ステップＳ２２）。
【０１３８】
なお、信頼度出力部３３は、k_aの信頼度に代えてk_vの信頼度を出力してもよい。また、測定部１０は、ステップＳ１９で算出したk_nsと所定の値（ステップＳ１１の所定値）とを比較し、k_nsが所定の値（ステップＳ１１の所定値）を超えた場合には、k_nsが所定の値を超えた旨を警告するための警告をSpO2出力部３１に出力する。さらに、測定部１０は、k_aの信頼度と予め設定されているk_aの信頼度判定用の値とを比較し、k_aがその値を超えた場合には、k_aがその値を超えた旨を示す警告をさらに信頼度出力部３３に出力してもよい。また、測定部１０は、k_aの信頼度に代えてk_vの信頼度を用いて、k_vの信頼度と予め設定されているk_vの信頼度判定用の値とを比較し、k_vがその値を超えた場合には、k_vがその値を超えた旨を示す警告をさらに信頼度出力部３３に出力してもよい。
【０１３９】
測定部１０は、測定周期が経過したことを内部のタイマから検出するまで待ち（ステップＳ２３：Ｎｏ）、測定周期が経過したことを検出すると（ステップＳ２３：Ｙｅｓ）、k_v算出周期が経過したか否かを判断する（ステップＳ２４）。
【０１４０】
k_v算出周期が経過していた場合は（ステップＳ２４：Ｙｅｓ）、算出済みフラグをＯＦＦにし（ステップＳ２５）、経過していなかった場合は（ステップＳ２４：Ｎｏ）、算出済みフラグをＯＮのままにしておく。そして、ステップＳ１０から処理を開始する。
【０１４１】
ここで、k_v算出周期とは、kv_minを算出する周期であり、k_vが変化しないと想定される時間である。k_vは急激には変化しないことから、k_aおよびSpO2の算出周期よりもkv_minの算出周期を長くする。例えば、k_aおよびSpO2の算出が例えば１秒毎に実施されるとする（ステップ２３の測定周期に該当）。これに対して、kv_minを算出する周期を、k_aおよびSpO2の算出周期よりも長い周期(例えば２秒)とする（ステップＳ２４のk_v算出周期に該当）。このようにすることで、測定結果の正しさを変えることなく、kv_minの算出回数を減らして計算量を減らすことが可能となる。
【０１４２】
＜信頼度算出処理＞
一般に、算出された動脈血酸素飽和度の値の大小によって、それに含まれる誤差は異なることが知られており、例えば、動脈血酸素飽和度が約９５％の場合と７０％の場合とでは動脈血酸素飽和度の誤差の大きさが異なる。そこで、k_nsの所定値は、以下のように変化させる。k_nsの所定値は、動脈血酸素飽和度が所定の動脈血酸素飽和度より小さい場合（k_aとk_vとは大きく、静脈の血中酸素飽和度が小さい場合に相当)にはk_nsの所定値を小さく設定し、一方、逆に動脈血酸素飽和度が所定の動脈血酸素飽和度より大きい場合（k_aとk_vとは小さく、静脈の血中酸素飽和度が大きい場合に相当）にはk_nsの所定値を大きく設定するとよい。このように算出された動脈血酸素飽和度の値に応じて、k_nsの所定値を変化させることで、例えば、動脈血酸素飽和度が大きい場合は閾値を大きくし、動脈血酸素飽和度の誤差が動脈血酸素飽和度が小さい場合とほぼ同じ値で警告や出力禁止を行うことができる。ここで、所定の動脈血酸素飽和度および所定の動脈血酸素飽和度に応じたk_nsの所定値は、過去の測定データを元にして決定した値を測定部１０に予め記憶しておいてもよいし、生体や生体の状態に応じて適宜に変化させてもよい。求められたk_a又はk_vの信頼度の指標として前記k_ns以外に、例えば、式(１６)から式(２１)の何れかの式で与えられる値zを用いてもよい。zの絶対値が大きいと、k_aの信頼度は低い。また、zの絶対値が大きいと、k_vの信頼度は低い。
【０１４３】
【数１６】
【０１４４】
【数１７】
【０１４５】
【数１８】
【０１４６】
【数１９】
【０１４７】
【数２０】
【０１４８】
【数２１】
【０１４９】
さらに、k_aの信頼度の指標の代わりとして、以下の式(２２)で与えられる値wのバラツキ(標準偏差、最大値−最小値など)を用いてもよい。値wのバラツキが小さいことは、求めたk_aの信頼度は高いことを示している。
【０１５０】
【数２２】
【０１５１】
この他、R_signalをy、IR_signalをxとした場合における（x,y）の回帰直線に対するyのバラツキや最大値と最小値との差等からk_aの信頼性を評価してもよい。式(１６)から式(２２)はR_signal、IR_signalをそれぞれ時間差分に置き換えたものでも良い。
【０１５２】
＜kv_min算出処理＞
ステップＳ１６のkv_min算出処理について、図４を用いて説明する。
【０１５３】
上記＜本発明の前提となる従来技術＞で説明したように、式（４）の左辺に相当するq（以下の式（１５）参照）が最も強く周期性を持つように選択されたk_vが、最適なk_vである。
【０１５４】
【数１５】
【０１５５】
ここで、qは、あるk_vに対する、i番目のデータ系列である式（１５）の左辺についての値を表す。
【０１５６】
より具体的に、k_vの値に応じて、式（１５）は、どのように変化するかについて、図５を用いて以下に説明する。
【０１５７】
図５は、測定データによるR-kv*IRの波形を示す図である。図５（ａ）は、k_vを最適値より小さい値で与えた場合の、k_vに対する式（１５）で与えられるデータ系列による波形を示し、図５（ｂ）は、k_vを最適値に近似される値で与えた場合の、k_vに対する式（１５）で与えられるデータ系列による波形を示し、そして、図５（ｃ）は、k_vを最適値より大きい値で与えた場合の、k_vに対する式（１５）で与えられるデータ系列による波形を示す。なお、波形は、時系列データのデータ間で補間している。図５（ａ）から（ｃ）において、これら横軸は、所定の短い時間（データ系列の番号iでは１からNまでのデータ系列に相当）であり、縦軸は、式１５のqである。
【０１５８】
最適値より小さいk_vを式（１５）に与えた場合には、図５（ａ）に示す波形で、最適値に近似される値k_vを式（１５）に与えた場合には、図５（ｂ）に示す波形となっている。これら２つの図を比較すると分かるように、最適値に近似される値k_vが代入されている場合（図５（ｂ）では、信号成分にノイズ成分が重畳している場合でも波形に強い周期性が認められるが、最適値より小さいk_vを式（１５）に与えた場合には、信号成分にノイズ成分が重畳していると、波形は、乱れ、周期性が弱い。そして、最適値より大きいk_vを式（１５）に与えた場合には、図５（ｃ）に示す波形となっている。この図５（ｃ）を図５（ｂ）と比較すると分かるように、最適値より大きいk_vを式（１５）に与えた場合（図５（ｃ））には、信号成分にノイズ成分が重畳していると、波形は、乱れ、周期性が弱い。
【０１５９】
以上のk_vに対する式（１５）の変化を利用するため、式（１５）の周期性が検出される。
【０１６０】
図４は、kv_min算出処理のフローチャートである。
【０１６１】
実施形態のkv_min算出処理では、k_vの初期値をk_v_sとして、k_vを所定の範囲内（例えば、k_v_s-0.5＜k_v＜k_v_s+0.5）で、Δk_vずつ、プラス側とマイナス側にそれぞれM回ずつ変化させて、(2*M+1)個のq（式（１５）参照）を得る。そして、それぞれのqを２値化することによって、それぞれの波形の幅を算出する。
【０１６２】
以下、具体的に、kv_min算出処理を説明する。
【０１６３】
まず、SvO2推定SpO2決定部１９は、本処理における初期値を設定する（ステップＳ４０）。k_vを変化させた回数を示すカウンタmに０（ゼロ）を設定し、パルス幅のバラツキ指数の最小値を示す変数Var_PW_minに、バラツキ指数としてあり得ない値、例えば、１００を設定する。また、k_v_sに、初回は例えば１を、２回目からは前回のkv_minを設定する。
【０１６４】
次に、各qに適した閾値を算出する際に用いる値を算出し、メモリ等に記憶する（ステップＳ４１）。具体的には、式（Ａ２）を用いてQ(0)を、式（Ａ３）を用いてT(0)を算出し、式（Ａ８）を用いてUを算出する。これらの算出に必要な値は、初期処理で算出してメモリに記憶されているSUM_sq_R、SUM_sq_IRおよびSUM_prdct_R_IRを読み出し、式（Ａ９）を用いてdkv_SUM_prdct_R_IRを算出し、（Ａ１０）を用いてdkv_SUM_sq_IRを算出し、式（Ａ１１）を用いてsq_dkvSUM_sq_IRを算出する。
【０１６５】
次に、q(0,i)を、メモリから読み出す（ステップＳ４２）。なお、２度目の測定からは、後述のステップＳ６２においてk_v_minに対応するq_min(i)をストアしておき、q_min(i)を初期値q(0,i)とする。より、計算量の増加を防ぐことができる。
【０１６６】
次に、ステップＳ４３〜ステップＳ５０の処理、すなわち、k_vをΔk_vずつプラスしてバラツキを算出する処理をＭ回繰り返す。
【０１６７】
まず、Q(m)を、式（Ｂ２）で算出し、２値化処理で用いる閾値を算出する（ステップＳ４３）。本実施形態では、Q(m)から、以下の式で求められる閾値thresh_1と閾値thresh_2とを設定する。なお、f_1およびf_2は、例えば、それぞれ、０．２、−０．２である。
【０１６８】
Sqr_Q=Q(m)0.5
thresh_1=f_1*Sqr_Q
thresh_2=f_2*Sqr_Q
なお、ｍ＝０の場合は、ステップＳ４１で求めたQ(0)を用いる。
【０１６９】
次に、q(m,i)の２値化処理を行う（ステップＳ４４）。ｍ＝０の場合は、ステップＳ４２で求めたq(0,i)を用い、ｍ＞０の場合は、ステップＳ５０で求めたq(m,i)を用いる。
【０１７０】
次に、２値化処理によって得た情報から、バラツキ指数Var_PWを算出する（ステップＳ４４）。２値化処理については、後の＜２値化処理＞の項で説明し、バラツキ指数算出処理については、後の＜バラツキ指数算出処理＞の項で説明する。
【０１７１】
ステップＳ４５で求めたバラツキ指数Var_PWと、最小バラツキ指数Var_PW_minとを比較し、小さいほうの値を最小バラツキ指数Var_PW_minに設定する。また、最小バラツキ指数Var_PW_minが設定された場合には、設定されたVar_PWのk_vをkv_minに設定する（ステップＳ４６）。
【０１７２】
ステップＳ４３〜ステップＳ５０の処理をＭ回繰り返していない場合は（ステップＳ４７：Ｙｅｓ）、k_vにΔk_vをプラスして、次の処理で用いるk_vを求める（ステップＳ４８）。
【０１７３】
そして、式（Ｂ２）および式（Ｂ３）を用いて、Q(m)およびT(m)を算出する（ステップＳ４９）。また、式（Ｂ１）を用いて、q(m,i)を算出する（ステップＳ５０）。
【０１７４】
ステップＳ４３〜ステップＳ５０の処理をＭ回繰り返した場合には（ステップＳ４７：Ｎｏ）、Δk_vずつ、マイナス側に変化させて、バラツキを算出する処理に移る。
【０１７５】
まず、SvO2推定SpO2決定部１９は、k_vを変化させた回数を示すカウンタmに−１を設定し（ステップＳ５１）、式（Ｂ５）を用いてQ(-1)を、式（Ｂ６）を用いてT(-1)を算出し（ステップＳ５２）、式（Ｂ４）を用いてq(-1,i)を算出する（ステップＳ５３）。
【０１７６】
次に、ステップＳ５４〜ステップＳ６１の処理、すなわち、k_vをΔk_vずつマイナスしてバラツキを算出する処理をＭ回繰り返す。
【０１７７】
ステップＳ５４〜ステップＳ５７の処理は、ステップＳ４３〜ステップＳ４６の処理と同様である。
【０１７８】
ステップＳ５４〜ステップＳ６１の処理をＭ回繰り返していない場合は（ステップＳ５８：Ｙｅｓ）、k_vからΔk_vをマイナスして、次の処理で用いるk_vを求める（ステップＳ５９）。
【０１７９】
そして、式（Ｂ５）および式（Ｂ６）を用いて、Q(m)およびT(m)を算出する（ステップＳ６０）。また、式（Ｂ４）を用いて、q(m,i)を算出する（ステップＳ６１）。
【０１８０】
ステップＳ５４〜ステップＳ６１の処理をＭ回繰り返した場合には（ステップＳ５８：Ｎｏ）、最小バラツキ指数Var_PW_minの時のq(m,i)を、q_min(i)として記憶する（ステップＳ６２）。この時には最小バラツキ指数Var_PW_minには、最小のバラツキ指数Var_PWが設定されており、kv_minには、最小バラツキ指数Var_PW_minのk_vが設定されている。
【０１８１】
本実施形態では、固定値であるΔk_vを用いているが、まずは大きいΔk_vを用いて広いk_vの範囲で、２値化したq(m,i)のパルス幅のバラツキを最小にするkv_min_1を求め、次にkv_min_1の周辺のより狭いk_vの範囲で小さいΔk_vを用いて、２値化したq(m,i)のパルス幅のバラツキを最小にするkv_minを求めることとしてもよい。これにより、計算量を増やすことなく高い精度でkv_minを決定することが可能となる。
【０１８２】
＜２値化処理＞
ここで、ステップＳ４４の２値化処理について、図６及び図７を用いて説明する。
【０１８３】
図６は、２値化処理のフローチャートである。図７は、２値化処理を説明するためのｑの波形の拡大図である。
【０１８４】
本２値化処理では、ステップＳ４３で求めた閾値thresh_1と閾値thresh_2とを用いる。
【０１８５】
そして、q(m,i)が、閾値thresh_1を下回る値から超えた時のｉをテーブルt1_l_to_h(kr_1)に、閾値thresh_1を上回る値から下回った時のｉをテーブルt1_h_to_l(kf_1)に記憶する。同様に、q(m,i)が、閾値thresh_2を下回る値から超えた時のｉをテーブルt2_l_to_h(kr_2)に、閾値thresh_2を上回る値から下回った時のｉをテーブルt2_h_to_l(kf_2)に記憶する（図７参照）。kr_1、kf_1、kr_2及びkf_2は、０から始まるインデックスである。
【０１８６】
まず、２値化処理で用いる変数の初期化を行う（ステップＳ７０）。ここで、status_1は、q(m,i)が閾値thresh_1を超えている場合に１が設定され、以下の場合に０が設定されるフラグである。status_2は、q(m,i)が閾値thresh_2を超えている場合に１が設定され、以下の場合に０が設定されるフラグである。このフラグは、現在処理しているq(m,i)のひとつ前のq(m,i-1)の値の状態を示す。なお、閾値は、図４のステップＳ４３で算出されている。
【０１８７】
フラグstatus_1が０で且つq(m,i)が閾値thresh_1より大きい場合、すなわち、閾値thresh_1を下回る値から閾値thresh_1を超えた場合（ステップＳ７１：Ｙｅｓ）は、ｉをテーブルt1_l_to_h(kr_1)に記憶する（ステップＳ７２）。その他の場合は（ステップＳ７１：Ｎｏ）、何も記憶しない。
【０１８８】
フラグstatus_1が１で且つq(m,i)が閾値thresh_1以下である場合、すなわち、閾値thresh_1を上回る値から閾値thresh_1以下となった場合（ステップＳ７３：Ｙｅｓ）は、ｉをテーブルt1_h_to_l(kf_1)に記憶する（ステップＳ７４）。その他の場合は（ステップＳ７３：Ｎｏ）、何も記憶しない。
【０１８９】
そして、q(m,i)が閾値thresh_1を超えているか、以下かを判断し、フラグstatus_1に０又は１を設定する（ステップＳ７５）。
【０１９０】
次に同様に、フラグstatus_2が０で且つq(m,i)が閾値thresh_2より大きい場合、すなわち、閾値thresh_2を下回る値から閾値thresh_2を超えた場合（ステップＳ７６：Ｙｅｓ）は、ｉをテーブルt2_l_to_h(kr_2)に記憶する（ステップＳ７７）。その他の場合は（ステップＳ７７：Ｎｏ）、何も記憶しない。
【０１９１】
フラグstatus_2が１で且つq(m,i)が閾値thresh_2以下である場合、すなわち、閾値thresh_2を上回る値から閾値thresh_2以下となった場合（ステップＳ７８：Ｙｅｓ）は、ｉをテーブルt2_h_to_l(kf_2)に記憶する（ステップＳ７９）。その他の場合は（ステップＳ７８：Ｎｏ）、何も記憶しない。
【０１９２】
そして、q(m,i)が閾値thresh_2を超えているか、否かを判断し、フラグstatus_2に０又は１を設定する（ステップＳ８０）。
【０１９３】
ステップＳ７１〜ステップＳ８０の処理を、ｉが１〜Ｎまで繰り返し（ステップＳ８１、ステップＳ８２）、t1_l_to_h(kr_1)等のテーブルに記憶されたｉの数をkr_1_max等に記憶する（ステップＳ８３）。
【０１９４】
＜バラツキ指数算出処理＞
ここで、ステップＳ４５のバラツキ指数算出処理について、図８を用いて説明する。
【０１９５】
図８は、バラツキ指数算出処理のフローチャートである。
【０１９６】
本実施形態では、２値化処理で求めたテーブルt1_l_to_h(kr_1)、t1_h_to_l(kf_1)、t2_l_to_h(kr_2)およびt2_h_to_l(kf_2)から、バラツキ指数Var_PWを求める。
【０１９７】
具体的には、図７で示すように、すべての波形の山と谷のａ１、ａ２、ｂ１およびｂ３の長さを求め、ａ１、ａ２、ｂ１およびｂ３それぞれの最小値と最大値を求める。そして、その調和平均を求めてバラツキ指数とする。
【０１９８】
まず、テーブルt1_l_to_h(0)とテーブルt1_h_to_l(0)に登録されているｉの大きさを比べる。波形が山から始まっているのか、谷から始まっているのかを判断する為である。
【０１９９】
t1_l_to_h(0)＞t1_h_to_l(0)の場合は（ステップＳ９０：Ｙｅｓ）、谷から始まっていると判断し、テーブルt1_h_to_lのインデックスを１ずらして計算し、ａ１を求める。すべての山のａ１を算出し、その最大値t_1_h_maxと最小値t_1_h_minとを求める（ステップＳ９３）。
【０２００】
また、すべての山のｂ１を算出し、その最大値t_1_l_maxと最小値t_1_l_minとを求める（ステップＳ９４）。
【０２０１】
一方、t1_l_to_h(0)≦t1_h_to_l(0)の場合は（ステップＳ９０：Ｎｏ）、山から始まっていると判断し、すべての山のａ１を算出し、その最大値t_1_h_maxと最小値t_1_h_minとを求める（ステップＳ９１）。
【０２０２】
また、テーブルt1_h_to_lのインデックスを１ずらして計算し、すべての山のｂ１を算出し、その最大値t_1_l_maxと最小値t_1_l_minとを求める（ステップＳ９２）。
【０２０３】
次に、テーブルt2_l_to_h(0)とテーブルt2_h_to_l(0)に登録されているｉの大きさを比べる。波形が山から始まっているのか、谷から始まっているのかを判断する為である。
【０２０４】
t2_l_to_h(0)＞t2_h_to_l(0)の場合は（ステップＳ９５：Ｙｅｓ）、谷から始まっていると判断し、テーブルt1_h_to_lのインデックスを１ずらして計算し、ａ２を求める。すべての山のａ２を算出し、その最大値t_2_h_maxと最小値t_2_h_minとを求める（ステップＳ９８）。
【０２０５】
また、すべての山のｂ２を算出し、その最大値t_2_l_maxと最小値t_2_l_minとを求める（ステップＳ９９）。
【０２０６】
一方、t2_l_to_h(0)≦t2_h_to_l(0)の場合は（ステップＳ９５：Ｎｏ）、山から始まっていると判断し、すべての山のａ２を算出し、その最大値t_2_h_maxと最小値t_2_h_minとを求める（ステップＳ９６）。
【０２０７】
また、テーブルt2_h_to_lのインデックスを１ずらして計算し、すべての山のｂ２を算出し、その最大値t_2_l_maxと最小値t_2_l_minとを求める（ステップＳ９７）。
【０２０８】
次に、最大値t_1_l_max、最小値t_1_l_min、最大値t_1_h_max、最小値t_1_h_min、最大値t_2_l_max、最小値t_2_l_min、最大値t_2_h_max、および、最小値t_2_h_minからバラツキ指数Var_PWを求める（ステップＳ１００）。
【０２０９】
具体的には、まず、(t_1_l_max-t_1_l_min)と(t_2_l_max-t_2_l_min)との平均調和t_l_max_minを例えば、以下の式で求める。
【０２１０】
t_l_max_min = {(t_1_l_max-t_1_l_min)*(t_2_l_max-t_2_l_min)}
/{(t_1_l_max-t_1_l_min)+(t_2_l_max-t_2_l_min)}
また、同様に、(t_1_h_max-t_1_h_min)と(t_2_h_max-t_2_h_min)との平均調和t_h_max_minを求める。
【０２１１】
そして、t_l_max_minとt_h_max_minとの内、大きい方をバラツキ指数Var_PWとする。
【０２１２】
＜脈拍数PR算出処理（ステップＳ２１）＞
上述の血中酸素飽和度検出処理（図２のステップＳ１１〜ステップＳ２０）が終了した後、脈拍数算出部２０によって脈拍数PR算出処理が実行される。
【０２１３】
図９は、脈拍数PR算出処理のフローチャートである。
【０２１４】
実施形態では、図４のステップＳ６２で記憶したq_min(i)の状態変化の間隔の平均T_aveを算出し、脈拍数PRを求める。
【０２１５】
まず、q_min(i)に対し、２値化処理を行う（ステップＳ１１０〜ステップＳ１１８）。閾値を０（ゼロ）として、図６の２値化処理と同様に行う。この場合の閾値は、前記のk_vを決めるときに用いた閾値でも良い。
【０２１６】
詳細には、２値化処理で用いる変数の初期化を行う（ステップＳ１１０）。ここで、statusは、q_min(i)が閾値を超えている場合に１が設定され、閾値以下の場合に０が設定されるフラグである。このフラグは、現在処理しているq_min(i)のひとつ前のq_min(i-1)の値の状態を示す。
【０２１７】
フラグstatusが０で且つq_min(i)が閾値より大きい場合、すなわち、閾値を下回る値から閾値を超えた場合（ステップＳ１１１：Ｙｅｓ）は、ｉをテーブルt_l_to_h(kr)に記憶する（ステップＳ１１２）。その他の場合は（ステップＳ１１１：Ｎｏ）、何も記憶しない。
【０２１８】
フラグstatusが１で且つq_min(i)が閾値以下である場合、すなわち、閾値を上回る値から閾値以下となった場合（ステップＳ１１３：Ｙｅｓ）は、ｉをテーブルt_h_to_l(kf)に記憶する（ステップＳ１１４）。その他の場合は（ステップＳ１１３：Ｎｏ）、何も記憶しない。
【０２１９】
そして、q_min(i)が閾値を超えているか、閾値以下かを判断し、フラグstatusに０又は１を設定する（ステップＳ１１５）。
【０２２０】
ステップＳ１１１〜ステップＳ１１７の処理を、ｉが１〜Ｎまで繰り返し（ステップＳ１１６、ステップＳ１１７）、テーブルt_l_to_h(kr)およびt_h_to_l(kf)に記憶されたｉの数をkr_maxおよびkf_maxに記憶する（ステップＳ１１８）。
【０２２１】
次に、脈拍数算出部２０は、q_min(i)の所定時間内の周期T(j)の平均値T_aveを以下の式で算出する（ステップＳ１１９）。
T_ave＝((Σk=0kr_max-1|t_l_to_h(k)-t_l_to_h(k+1)|)
+(Σk=0kf_max-1|t_h_to_l(k)-t_h_to_l(k+1)|))/(kr_max+kf_max)
次に、脈拍数を算出する（ステップＳ１２０）。脈拍数は、前記周期の逆数として求められる。すなわち、脈拍数をPR(回／min)として、PR＝60／T_aveで表される。
＜変形例＞
実施形態では、kv_min算出処理（図４）において、q(m,i)の２値化処理を行う場合、測定データとして記憶しているすべてのデータ、すなわち、i=1〜i=Nまでのデータそれぞれについて２値化の処理を行っている。変形例では、このデータを間引いて、例えば、ｉが偶数のデータのみで２値化処理を行うことで、計算量を減らす。
【０２２２】
この場合、データを間引いたことによってパルス幅検出の精度を落とさないために、信号qが閾値をまたいだ場合には、またいだデータの間で間引いたデータについて２値化の処理を行い、閾値をまたいだ時刻を正確に取得するものとする。なお、データの間引きは、信号qのパルス幅のバラツキを検出できて、適切なkv_minを算出できる範囲で行うものとする。
【０２２３】
また、実施形態では、脈拍数PR算出処理（図９）において、q_min(i)から心拍数PRを求める場合、i=1〜i=Nまでのデータすべてを用いているが、変形例の２値化処理と同様に、このデータを間引いて、例えば、ｉが偶数のデータのみで２値化処理を行うことで、計算量を減らす。
【０２２４】
以下、変形例の２値化処理および脈拍数PR算出処理を説明する。
【０２２５】
＜２値化処理＞
図１０は、変形例の２値化処理のフローチャートである。ここでは、１つ置きにデータを間引くものとする。
【０２２６】
変形例の２値化処理は、図６を用いて説明した実施形態の２値化処理と、処理の流れはほぼ同様である。異なる点は、まず、q(m,i)のインデックスｉに２を加算する点（ステップＳ２１４）である。１つ置きのデータを処理するからである。また、実施形態では、q(m,i)が閾値をまたいだ時には、そのｉをテーブルに記憶している。しかし、変形例では、q(m,i)が閾値をまたいだ時には、その時飛ばしたデータq(m,i-1)が閾値をまたいでいるかを確認して、そのｉまたはi-1をテーブルに記憶する（ステップＳ２０２〜ステップＳ２０５、ステップＳ２０７〜ステップＳ２０９、ステップＳ２１３〜ステップＳ２１６、ステップＳ２１８〜ステップＳ２２１）。
【０２２７】
まず、２値化処理で用いる変数の初期化を行う（ステップＳ２００）。この変数は、図６と同じである。
【０２２８】
フラグstatus_1が０で且つq(m,i)が閾値thresh_1より大きい場合、すなわち、閾値thresh_1を下回る値から閾値thresh_1を超えた場合（ステップＳ２０１：Ｙｅｓ）は、q(m,i-1)の値を算出する（ステップＳ２０２）。すなわち、飛ばしたデータを算出する。
【０２２９】
そして、q(m,i-1)が閾値thresh_1より大きい場合、すなわち、飛ばしたデータが閾値thresh_1を下回る値から閾値thresh_1を超えた場合（ステップＳ２０３：Ｙｅｓ）は、ｉ-1をテーブルt1_l_to_h(kr_1)に記憶し（ステップＳ２０４）、超えない場合（ステップＳ２０３：Ｎｏ）は、ｉをテーブルt1_l_to_h(kr_1)に記憶する（ステップＳ２０５）。
【０２３０】
次に、フラグstatus_1が１で且つq(m,i)が閾値thresh_1以下である場合、すなわち、閾値thresh_1を上回る値から閾値thresh_1以下となった場合（ステップＳ２０６：Ｙｅｓ）も、q(m,i-1)の値を算出する（ステップＳ２０７）。
【０２３１】
そして、q(m,i-1)が閾値thresh_1以下の場合、すなわち、飛ばしたデータが閾値thresh_1を上回る値から閾値thresh_1以下となった場合（ステップＳ２０８：Ｙｅｓ）は、ｉ-1をテーブルt1_l_to_h(kr_1)に記憶し（ステップＳ２０９）、以下とならなかった場合（ステップＳ２０８：Ｎｏ）は、ｉをテーブルt1_h_to_l(kf_1)に記憶する（ステップＳ２１０）。
【０２３２】
そして、q(m,i)が閾値thresh_1を超えているか、以下かを判断し、フラグstatus_1に０又は１を設定する（ステップＳ２１１）。
【０２３３】
次に同様に、フラグstatus_2が０で且つq(m,i)が閾値thresh_2より大きい場合、すなわち、閾値thresh_2を下回る値から閾値thresh_2を超えた場合（ステップＳ２１２：Ｙｅｓ）は、q(m,i-1)の値を算出する（ステップＳ２１３）。すなわち、飛ばしたデータを算出する。
【０２３４】
そして、q(m,i-1)が閾値thresh_2より大きい場合、すなわち、飛ばしたデータが閾値thresh_2を下回る値から閾値thresh_2を超えた場合（ステップＳ２１４：Ｙｅｓ）は、ｉ-1をテーブルt2_l_to_h(kr_1)に記憶し（ステップＳ２１５）、超えない場合（ステップＳ２１４：Ｎｏ）は、ｉをテーブルt2_l_to_h(kr_1)に記憶する（ステップＳ２１６）。
【０２３５】
次に、フラグstatus_2が１で且つq(m,i)が閾値thresh_2以下である場合、すなわち、閾値thresh_2を上回る値から閾値thresh_2以下となった場合（ステップＳ２１７：Ｙｅｓ）も、q(m,i-1)の値を算出する（ステップＳ２１８）。
【０２３６】
そして、q(m,i-1)が閾値thresh_2より小さい場合、すなわち、飛ばしたデータが閾値thresh_2を上回る値から閾値thresh_2以下となった場合（ステップＳ２１９：Ｙｅｓ）は、ｉ-1をテーブルt2_h_to_l(kf_1)に記憶し（ステップＳ２２０）、以下とならなかった場合（ステップＳ２１９：Ｎｏ）は、ｉをテーブルt2_h_to_l(kf_1)に記憶する（ステップＳ２２１）。
【０２３７】
そして、q(m,i)が閾値thresh_2を超えているか、以下かを判断し、フラグstatus_2に０又は１を設定する（ステップＳ２２２）。
【０２３８】
ステップＳ２０１〜ステップＳ２２２の処理を、ｉが１〜Ｎまで繰り返し（ステップＳ２１３、ステップＳ２１４）、t1_l_to_h(kr_1)等のテーブルに記憶されたｉの数をkr_1_max等に記憶する（ステップＳ２１５）。
【０２３９】
＜脈拍数PR算出処理＞
図１１は、変形例の脈拍数PR算出処理のフローチャートである。
【０２４０】
図４のステップＳ６２で記憶したq_min(i)の状態変化の間隔の平均T_aveを算出し、脈拍数PRを求める。
【０２４１】
まず、q_min(i)に対し、２値化処理を行う（ステップＳ２２０〜ステップＳ２３４）。閾値を０（ゼロ）として、図１０の２値化処理と同様に行う。この場合の閾値は、前記のk_vを決めるときに用いた閾値でも良い。
【０２４２】
詳細には、２値化処理で用いる変数の初期化を行う（ステップＳ２２０）。この変数は、図９と同じである。
【０２４３】
フラグstatusが０で且つq_min(i)が閾値より大きい場合、すなわち、閾値を下回る値から閾値を超えた場合（ステップＳ２２１：Ｙｅｓ）は、q_min(i-1)の値を算出する（ステップＳ２２２）。すなわち、飛ばしたデータを算出する。
【０２４４】
そして、q_min(i-1)が閾値より大きい場合、すなわち、飛ばしたデータが閾値を下回る値から閾値を超えた場合（ステップＳ２２３：Ｙｅｓ）は、i-1をテーブルt_l_to_h(kr)に記憶し（ステップＳ２２４）、超えない場合（ステップＳ２２３：Ｎｏ）は、ｉをテーブルt_l_to_h(kr)に記憶する（ステップＳ２２５）。
【０２４５】
次に、フラグstatusが１で且つq_min(i)が閾値以下である場合、すなわち、閾値を上回る値から閾値以下となった場合（ステップＳ２２６：Ｙｅｓ）も、q_min(i-1)の値を算出する（ステップＳ２２７）。
【０２４６】
そして、q_min(i-1)が閾値以下の場合、すなわち、飛ばしたデータが閾値を上回る値から閾値以下となった場合（ステップＳ２２８：Ｙｅｓ）は、i-1をテーブルt_h_to_l(kf)に記憶し（ステップＳ２２９）、以下とならなかった場合（ステップＳ２２８：Ｎｏ）は、ｉをテーブルt_h_to_l(kf)に記憶する（ステップＳ２３０）。
【０２４７】
そして、q_min(i)が閾値を超えているか、以下かを判断し、フラグstatusに０又は１を設定する（ステップＳ２３１）。
【０２４８】
次に、q_min(i)の所定時間内の周期T(j)の平均値T_aveを、図９のステップＳ１１９と同様に算出する（ステップＳ２３５）。
【０２４９】
次に、脈拍数PRを、図９のステップＳ１２０と同様に算出する（ステップＳ２３６）。
【０２５０】
以上のような実施形態及び変形例によれば、漸化式を用いて周期性が最も良い信号成分を求めているので、従来技術より演算処理を簡単化することができ、演算処理量を低減することができる。変形例では、更に、より少ないデータに対して処理を行うため、より演算処理量を低減できる。そのため、演算処理に伴う消費電力を抑制することができる。この結果、生体情報測定装置４０に電池を使用することが可能となり、携帯に便利な小型軽量な生体情報測定装置を提供することができる。このため、例えば、登山者や在宅酸素療法患者が常時携帯しても体力的な負担を低下させることができる。
【０２５１】
また、本実施形態に係る生体情報測定装置４０は、信号成分の強度と前記ノイズ成分の強度との比に関する情報を有する指標に基づいて、ノイズ成分を抽出するか否かの判断を行うので、ノイズ成分が無視できるほど小さい場合においては、このノイズ成分の除去を行う演算を省略することができるので、生体情報を得るための演算時間が短縮され、消費電力を抑えて電池動作時間をより長くできるとともに、より早く生体情報を得ることができる。
【０２５２】
また、本実施形態に係る生体情報測定装置４０は、音声、文字および光等によって生体情報、信号信頼度およびノイズ信頼度を出力するので、これらを容易に確認することが可能となる。さらに、本実施形態に係る生体情報測定装置４０は、指標が所定の値を超えたことを示す警告、あるいは信号信頼度またはノイズ信頼度の低下を示す警告を出力するので、測定対象者や測定者等に生体情報の取り扱いに関して注意を喚起することもできる。そして、前記指標、前記ノイズ信頼度および前記信号信頼度に応じて、演算された生体情報の取り扱いを変更することが可能となり、例えば、指標を基に再計測等を行うことができる。
【０２５３】
また、本実施形態に係る生体情報測定装置４０は、脈拍数を算出する際に、２値化した波形を用いることで、例えばダブルピーク等による誤測定を低減することができる。
【０２５４】
また、本実施形態に係る生体情報測定方法は、同じ計算資源（ハードウェア資源）において、生体情報をより早い計算時間で算出することできる。
【０２５５】
上記実施形態では、一例として、周期性をもつ２波長の光の強度変化を検出して血中酸素飽和度を求める生体情報測定装置４０を例示した。この他、周期性をもつ１波長の光の強度変化に基づき、脈拍、不整脈（心房細動・期外収縮）の検出、除細動時のモニター、自律神経障害、血管年齢等の診断等を行う生体情報測定装置にも本発明を適用することができる。さらに、脈波波形以外に、心電等の生体情報を測定対象とする装置にも適用することができる。
【０２５６】
実施形態の２値化処理は、式（１５）が周期性を持つようなk_vを求めるための指標として、２値化処理後の波形の山と谷の周期を用いたが、これに限定されるものではない。２値化処理後の波形における、他の指標として、山が繰り返される周期のみ、谷が繰り返される周期のみを用いてもよい。また、山の周期の最大値‐山の周期の最小値、谷の周期の最大値‐谷の周期の最小値、山の周期の標準偏差、谷の周期の標準偏差、山の幅の最大値‐山の幅の最小値、谷の幅の最大値‐谷の幅の最小値、山の幅の標準偏差、谷の幅の標準偏差を用いても良い。さらに、指標として、R_signal - k_v * IR_signalの波形の極大値と極小値から求められる振幅の最大‐最小や標準偏差等を用いてもよい。
【０２５７】
実施形態では、閾値を２つ用いたが、３つ以上の閾値を設け、式（１５）に対応した信号がそれぞれの閾値を超えるか否かで２値化した３種以上の２値化パルスそれぞれについて前記の種々の周期性の指標を算出して、その平均値（単純平均、調和平均、相乗平均を含む）を最小にするk_vを最適値としてもよい。閾値の数を多くすると計算量は増えるが、周期性の評価の精度が向上する利点がある。
【０２５８】
いずれの指標においても、指標の標準偏差や指標の最大値と最小値との差等から、所定の時間内で、最も指標のばらつきが最小となるようなk_vを求めてもよい。複数の指標を用いた場合、それぞれの指標を最小にするk_vが異なる場合がある。その場合はそれぞれの指標を最小にするk_vの調和平均値を最適値とする。平均値としては単純平均値、相乗平均値でもよい。また、複数の指標を最小にする複数のk_vのうち前回のk_vに最も近いものを採用しても良い。また、複数の指標の平均値(単純平均値、調和平均値、相乗平均値など)を指標として評価しても良い。
【０２５９】
複数の指標に応じて複数のk_vを求めた場合には、各k_vに対して、式（８）の等号が成立するようなk_aをおのおの算出しても良い。複数のk_aから最終的なk_aを決める方法は、複数のk_aの平均値(単純平均値、調和平均値、相乗平均値など)を最終的なk_aとしてもよいし、複数のk_aのうち前回の最終的なk_aに最も近いものを今回の最終的なk_aとしてもよいし、複数のk_aのうち最大と最小を除いたものの平均値を最終的なk_aとしてもよい。さらに、前回のk_aとの差の絶対値に対応した重みを掛けて荷重平均しても良い。
【０２６０】
また、k_vは例えば２００組のR_signal、IR_signalについて前記の方法で決定し、k_aを求める際は２００組のデータを例えば４分割して４組の式(８)を解いて４個のk_aを算出してそれらの平均値(単純平均、調和平均、相乗平均値、前回k_aとの差の絶対値を重みとする荷重平均値など)を求めても良い。２００組のデータ全体で求めたk_aと前記の４個のk_aの平均値をさらに平均しても良い。
【０２６１】
この他、複数の指標からそれぞれに応じて複数のk_vが求められた場合には、連続するk_vの範囲の中央値、または直前の動脈血酸素飽和度演算で得られたk_vに最も近い値を最適値としてもよい。
【０２６２】
なお、今までの説明ではノイズ成分の指標が所定の値以下のときは式（１４）など従来の方法によりk_aを算出しているが、式（８）に前回のk_v(＝k_v_0)を代入してk_aを求めても良い。ノイズが小さい場合は式(８)の左辺は
Σ{IR_signal * k_v - R_signal} * {IR_signal * k_a - R_signal}
≒(k_v - k_a) * ΣIR_signal * {IR_signal * k_a - R_signal}
≒(k_v - k_a) * {ΣIR_signal2 * k_a - ΣIR_signal * R_signal}
なので、k_v としてk_aとは異なる値を代入すれば、方程式は
ΣIR_signal2 * k_a - ΣIR_signal * R_signal＝０
と等価である。この解は
k_a＝ΣIR_signal * R_signal ／ ΣIR_signal2
が得られる。ノイズが小さいときはR_signal≒真のk_a * IR_signalなので、上式の右辺は真のk_aになる。
【０２６３】
したがって、式(８)に前回のk_vなどk_aとは異なると思われる値を代入して解くことにより、真のk_aに十分近い値を得ることができる。
【０２６４】
また、ノイズが小さいときはR_signal − k_v * IR_signalは広い範囲のk_vに対して強い周期性を持つが、その場合でも前回のk_vとk_vの差の絶対値を前記の種々の指標に掛けた値を評価すれば、k_vの最適値はほぼ前回のk_vに近い値になるので、ノイズの指標が大きいときと同様の処理を行っても真のk_aに十分近い値を得ることができる。計算速度よりもプログラムの簡素化を重視する場合はノイズ指標の大きさによって処理を分けなくても良い。
【０２６５】
なお、周期性を利用する一例として、式（１５）の左辺が周期性を満たすようにk_vを求める例を挙げたが、式（１５）に限定されるものではなく、k_vが信号成分で表される数式であれば、式（１５）に代えて用いてもよい。
【０２６６】
本発明を表現するために、上述において図面を参照しながら実施形態を通して本発明を適切且つ十分に説明したが、当業者であれば上述の実施形態を変更および／または改良することは容易に為し得ることであると認識すべきである。したがって、当業者が実施する変更形態または改良形態が、請求の範囲に記載された請求項の権利範囲を離脱するレベルのものでない限り、当該変更形態または当該改良形態は、当該請求項の権利範囲に包括されると解釈される。
【符号の説明】
【０２６７】
１ データ取得部
２ センサ（R）部
３ センサ（IR）部
１０ 測定部
１１ AC/DC(R)変換部
１２ AC/DC(IR)変換部
１３ BPF(R)部
１４ BPF(IR)部
１５ ΣR_signal*IR_signal算出部
１６ Σ(R_signal)2算出部
１７ Σ(IR_signal)2算出部
１８ 初期値算出部
１９ SvO2推定SpO2決定部
２０ 脈拍数算出部
２１ 信頼度算出部
２２ 差分算出部
３０ 出力部
３１ SpO2出力部
３２ 脈拍数出力部
３３ 信頼度出力部
４０ 生体情報測定装置
【特許請求の範囲】
【請求項１】
周期性を有する第１信号成分と、第１ノイズ成分とを含む第１の時系列信号と、前記第１信号成分と所定の関係を有する第２信号成分および前記第１ノイズ成分と所定の関係を有する第２ノイズ成分を含む第２の時系列信号と、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係とに基づいて、前記第１の時系列信号と前記第２の時系列信号とから前記第１信号成分を含む信号であって、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を所定の差分ずつ異ならせた複数の信号を生成する信号生成部と、
第１の所定時間範囲での前記信号生成部で生成した各信号の周期性を用いて、前記第１の所定時間範囲での前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を推定する推定部とを備え、
前記信号生成部は、基準となる１つの信号を生成し、当該信号と、前記第１時系列信号と前記所定の差分とから算出したデータとを用いた漸化式によって他の信号を生成する
ことを特徴とする信号処理装置。
【請求項２】
請求項１において、
前記信号生成部が、生成した信号に２値化処理を行うことで２値化信号をさらに生成し、
前記推定部が、前記２値化信号の周期性を用いて前記信号生成部で生成した信号の周期を推定すること
を特徴とする信号処理装置。
【請求項３】
請求項２において、
前記信号生成部は、基準となる信号の閾値を算出し、当該閾値と、前記第１の時系列信号又は前記第２の時系列信号と前記所定の差分とを用いて算出したデータとを用いて、漸化式によって他の信号の閾値を算出すること
を特徴とする信号処理装置。
【請求項４】
請求項２または請求項３において、
前記信号生成部が生成した信号は、時系列信号であり、
前記信号生成部は、生成した信号から所定数毎に間引かれた結果のデータを用いて２値化処理を行うこと
を特徴とする信号処理装置。
【請求項５】
請求項４において、
前記閾値を挟む２つの前記間引かれた結果のデータの間の間引いたデータに、前記２値化処理を行うこと
を特徴とする信号処理装置。
【請求項６】
信号処理方法であって、
第１の所定時間範囲での、周期性を有する第１信号成分と第１ノイズ成分とを含む第１の時系列信号と、前記第１信号成分と所定の関係を有する第２信号成分および前記第１ノイズ成分と所定の関係を有する第２ノイズ成分を含む第２の時系列信号と、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係とに基づいて、前記第１の時系列信号と前記第２の時系列信号とから前記第１信号成分を含む信号であって、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を所定の差分ずつ異ならせた複数の信号を生成する信号生成ステップと、
第１の所定時間範囲での前記生成ステップで生成した各信号の周期性を用いて、前記第１の所定時間範囲での前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を推定する推定ステップとを備え、
前記信号生成ステップは、基準となる１つの信号を生成し、当該信号と、前記第１時系列信号と前記所定の差分とから算出したデータとを用いた漸化式によって他の信号を生成する
ことを特徴とする信号処理方法。
【請求項７】
互いに波長の異なる複数の光を生体へそれぞれ照射して前記生体を透過または反射した各光をそれぞれ受光することによって得られた少なくとも第１測定データまたは第２測定データに基づいて、前記生体の生体情報を測定する生体情報測定装置であって、
周期性を有する第１信号成分と、第１ノイズ成分からなる第１測定データを測定する第１測定部と、
前記第１信号成分と所定の関係を有する第２信号成分および第１ノイズ成分と所定の関係を有する第２ノイズ成分を含む第２測定データを測定する第２測定部と、
前記第１測定データと、前記第２測定データと、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係とに基づいて、前記第１測定データと前記第２測定データとから前記第１信号成分を含む信号であって、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を所定の差分ずつ異ならせた複数の信号を生成する信号生成部と、
前記信号生成部が生成した信号の周期性を用いて、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を推定する推定部とを備え、
前記信号生成部は、基準となる１つの信号を生成し、当該信号と、前記第１測定データと前記所定の差分とから算出したデータとを用いた漸化式によって他の信号を生成すること
を特徴とする生体情報測定装置。
【請求項８】
請求項７において、
前記推定部が推定した所定の関係に基づいて、所定の周期で血中酸素飽和度を算出する血中酸素飽和度算出部を備え、
前記生成信号部は、前記所定の周期よりも長い周期で、前記複数の信号を生成して、２値化信号を生成し、
前記推定部は、前記２値化信号の周期性に基づいて、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を推定すること
を特徴とする生体情報測定装置。
【請求項９】
請求項７または請求項８において、
前記推定部が推定した所定の関係を有する信号に２値化処理を行うことで２値化信号を生成し、当該２値化信号の周期性を用いて脈拍数を算出する脈拍数算出部を備えること
を特徴とする生体情報測定装置。
【請求項１】
周期性を有する第１信号成分と、第１ノイズ成分とを含む第１の時系列信号と、前記第１信号成分と所定の関係を有する第２信号成分および前記第１ノイズ成分と所定の関係を有する第２ノイズ成分を含む第２の時系列信号と、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係とに基づいて、前記第１の時系列信号と前記第２の時系列信号とから前記第１信号成分を含む信号であって、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を所定の差分ずつ異ならせた複数の信号を生成する信号生成部と、
第１の所定時間範囲での前記信号生成部で生成した各信号の周期性を用いて、前記第１の所定時間範囲での前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を推定する推定部とを備え、
前記信号生成部は、基準となる１つの信号を生成し、当該信号と、前記第１時系列信号と前記所定の差分とから算出したデータとを用いた漸化式によって他の信号を生成する
ことを特徴とする信号処理装置。
【請求項２】
請求項１において、
前記信号生成部が、生成した信号に２値化処理を行うことで２値化信号をさらに生成し、
前記推定部が、前記２値化信号の周期性を用いて前記信号生成部で生成した信号の周期を推定すること
を特徴とする信号処理装置。
【請求項３】
請求項２において、
前記信号生成部は、基準となる信号の閾値を算出し、当該閾値と、前記第１の時系列信号又は前記第２の時系列信号と前記所定の差分とを用いて算出したデータとを用いて、漸化式によって他の信号の閾値を算出すること
を特徴とする信号処理装置。
【請求項４】
請求項２または請求項３において、
前記信号生成部が生成した信号は、時系列信号であり、
前記信号生成部は、生成した信号から所定数毎に間引かれた結果のデータを用いて２値化処理を行うこと
を特徴とする信号処理装置。
【請求項５】
請求項４において、
前記閾値を挟む２つの前記間引かれた結果のデータの間の間引いたデータに、前記２値化処理を行うこと
を特徴とする信号処理装置。
【請求項６】
信号処理方法であって、
第１の所定時間範囲での、周期性を有する第１信号成分と第１ノイズ成分とを含む第１の時系列信号と、前記第１信号成分と所定の関係を有する第２信号成分および前記第１ノイズ成分と所定の関係を有する第２ノイズ成分を含む第２の時系列信号と、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係とに基づいて、前記第１の時系列信号と前記第２の時系列信号とから前記第１信号成分を含む信号であって、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を所定の差分ずつ異ならせた複数の信号を生成する信号生成ステップと、
第１の所定時間範囲での前記生成ステップで生成した各信号の周期性を用いて、前記第１の所定時間範囲での前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を推定する推定ステップとを備え、
前記信号生成ステップは、基準となる１つの信号を生成し、当該信号と、前記第１時系列信号と前記所定の差分とから算出したデータとを用いた漸化式によって他の信号を生成する
ことを特徴とする信号処理方法。
【請求項７】
互いに波長の異なる複数の光を生体へそれぞれ照射して前記生体を透過または反射した各光をそれぞれ受光することによって得られた少なくとも第１測定データまたは第２測定データに基づいて、前記生体の生体情報を測定する生体情報測定装置であって、
周期性を有する第１信号成分と、第１ノイズ成分からなる第１測定データを測定する第１測定部と、
前記第１信号成分と所定の関係を有する第２信号成分および第１ノイズ成分と所定の関係を有する第２ノイズ成分を含む第２測定データを測定する第２測定部と、
前記第１測定データと、前記第２測定データと、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係とに基づいて、前記第１測定データと前記第２測定データとから前記第１信号成分を含む信号であって、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を所定の差分ずつ異ならせた複数の信号を生成する信号生成部と、
前記信号生成部が生成した信号の周期性を用いて、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を推定する推定部とを備え、
前記信号生成部は、基準となる１つの信号を生成し、当該信号と、前記第１測定データと前記所定の差分とから算出したデータとを用いた漸化式によって他の信号を生成すること
を特徴とする生体情報測定装置。
【請求項８】
請求項７において、
前記推定部が推定した所定の関係に基づいて、所定の周期で血中酸素飽和度を算出する血中酸素飽和度算出部を備え、
前記生成信号部は、前記所定の周期よりも長い周期で、前記複数の信号を生成して、２値化信号を生成し、
前記推定部は、前記２値化信号の周期性に基づいて、前記第１ノイズ成分と前記第２ノイズ成分との前記所定の関係を推定すること
を特徴とする生体情報測定装置。
【請求項９】
請求項７または請求項８において、
前記推定部が推定した所定の関係を有する信号に２値化処理を行うことで２値化信号を生成し、当該２値化信号の周期性を用いて脈拍数を算出する脈拍数算出部を備えること
を特徴とする生体情報測定装置。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図７】
【図８】
【図９】
【図１０】
【図１１】
【公開番号】特開２０１２−２３５８９０（Ｐ２０１２−２３５８９０Ａ）
【公開日】平成２４年１２月６日（２０１２．１２．６）
【国際特許分類】
【出願番号】特願２０１１−１０６５４６（Ｐ２０１１−１０６５４６）
【出願日】平成２３年５月１１日（２０１１．５．１１）
【出願人】（３０３０５０１６０）コニカミノルタオプティクス株式会社 (175)
【Ｆターム（参考）】
【公開日】平成２４年１２月６日（２０１２．１２．６）
【国際特許分類】
【出願日】平成２３年５月１１日（２０１１．５．１１）
【出願人】（３０３０５０１６０）コニカミノルタオプティクス株式会社 (175)
【Ｆターム（参考）】
[ Back to top ]