JP6536038B2 - 周期推定装置、周期推定方法及びプログラム。 - Google Patents

周期推定装置、周期推定方法及びプログラム。 Download PDF

Info

Publication number
JP6536038B2
JP6536038B2 JP2015008008A JP2015008008A JP6536038B2 JP 6536038 B2 JP6536038 B2 JP 6536038B2 JP 2015008008 A JP2015008008 A JP 2015008008A JP 2015008008 A JP2015008008 A JP 2015008008A JP 6536038 B2 JP6536038 B2 JP 6536038B2
Authority
JP
Japan
Prior art keywords
frequency
signal
period
dimensional
unit
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
JP2015008008A
Other languages
English (en)
Other versions
JP2016131713A (ja
Inventor
山本　康平
康平 山本
素子 加賀谷
素子 加賀谷
前野　蔵人
蔵人 前野
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Oki Electric Industry Co Ltd
Original Assignee
Oki Electric Industry Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Oki Electric Industry Co Ltd filed Critical Oki Electric Industry Co Ltd
Priority to JP2015008008A priority Critical patent/JP6536038B2/ja
Priority to US14/955,994 priority patent/US20160206230A1/en
Publication of JP2016131713A publication Critical patent/JP2016131713A/ja
Application granted granted Critical
Publication of JP6536038B2 publication Critical patent/JP6536038B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/103—Detecting, measuring or recording devices for testing the shape, pattern, colour, size or movement of the body or parts thereof, for diagnostic purposes
    • A61B5/11—Measuring movement of the entire body or parts thereof, e.g. head or hand tremor, mobility of a limb
    • A61B5/1126—Measuring movement of the entire body or parts thereof, e.g. head or hand tremor, mobility of a limb using a particular sensing technique
    • A61B5/1128—Measuring movement of the entire body or parts thereof, e.g. head or hand tremor, mobility of a limb using a particular sensing technique using image analysis
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7235—Details of waveform analysis
    • A61B5/7246—Details of waveform analysis using correlation, e.g. template matching or determination of similarity
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/08—Detecting, measuring or recording devices for evaluating the respiratory organs
    • A61B5/0816—Measuring devices for examining respiratory frequency
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/103—Detecting, measuring or recording devices for testing the shape, pattern, colour, size or movement of the body or parts thereof, for diagnostic purposes
    • A61B5/11—Measuring movement of the entire body or parts thereof, e.g. head or hand tremor, mobility of a limb
    • A61B5/113—Measuring movement of the entire body or parts thereof, e.g. head or hand tremor, mobility of a limb occurring during breathing
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/103—Detecting, measuring or recording devices for testing the shape, pattern, colour, size or movement of the body or parts thereof, for diagnostic purposes
    • A61B5/11—Measuring movement of the entire body or parts thereof, e.g. head or hand tremor, mobility of a limb
    • A61B5/113—Measuring movement of the entire body or parts thereof, e.g. head or hand tremor, mobility of a limb occurring during breathing
    • A61B5/1135—Measuring movement of the entire body or parts thereof, e.g. head or hand tremor, mobility of a limb occurring during breathing by monitoring thoracic expansion
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7225—Details of analog processing, e.g. isolation amplifier, gain or sensitivity adjustment, filtering, baseline or drift compensation
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7271—Specific aspects of physiological measurement analysis
    • A61B5/7278—Artificial waveform generation or derivation, e.g. synthesising signals from measured signals
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/48—Diagnostic techniques
    • A61B8/488—Diagnostic techniques involving Doppler signals
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B8/00—Diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/52—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves
    • A61B8/5207—Devices using data or image processing specially adapted for diagnosis using ultrasonic, sonic or infrasonic waves involving processing of raw data to produce diagnostic data, e.g. for generating an image
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/05—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves 
    • A61B5/0507—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves  using microwaves or terahertz waves

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Molecular Biology (AREA)
  • Animal Behavior & Ethology (AREA)
  • Veterinary Medicine (AREA)
  • Public Health (AREA)
  • Physics & Mathematics (AREA)
  • General Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Pathology (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Surgery (AREA)
  • Physiology (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Signal Processing (AREA)
  • Radiology & Medical Imaging (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Psychiatry (AREA)
  • Artificial Intelligence (AREA)
  • Dentistry (AREA)
  • Oral & Maxillofacial Surgery (AREA)
  • Pulmonology (AREA)
  • Power Engineering (AREA)
  • Measurement Of The Respiration, Hearing Ability, Form, And Blood Characteristics Of Living Organisms (AREA)
  • Radar Systems Or Details Thereof (AREA)

Description

本発明は、周期推定装置、周期推定方法及びプログラムに関する。
近年、被験者の健康状態や安否情報等を確認するために、当該被験者の呼吸などの運動を検出し、当該運動を評価するシステムに関する技術開発が進められている。
例えば、下記特許文献１で開示された技術は、ドップラーセンサを用いて被験者から得られた出力信号に基づき、呼吸数を計測している。より詳細に説明すると、下記特許文献１では、ドップラーセンサによって被験者から得られたドップラー信号をフーリエ解析し、最もピークの高い周波数を呼吸周波数として抽出し、抽出された呼吸周波数に基づき呼吸数を算出する技術が開示されている。
また、下記特許文献２で開示された技術は、睡眠を評価するシステムを実現するために、呼吸波形に基づく波形情報を利用して、呼吸の安定度を算出している。より詳細に説明すると、下記特許文献２では、睡眠中の呼吸波形を抽出し、得られた波形をフーリエ解析して振幅スペクトル上で振幅がピークとなる成分の周波数を検出し、当該周波数値の時間変化の統計量である平均値や標準偏差に基づき、呼吸の安定度を評価する技術が開示されている。
特開２０１３−２１１７７９号公報 国際公開第２０１１／０１９０９１号
しかしながら、例えば上記特許文献においては、呼吸波形の周期は、複数の周期を含む波形をフーリエ解析することにより得られた周波数スペクトルに基づいて算出されるが、呼吸のような１周期毎に微小に変化する信号の周期を、１周期毎に正確に算出することは困難であるという問題があった。より具体的に説明すると、例えば、フーリエ解析として離散フーリエ変換（ＤＦＴ：Ｄｉｓｃｒｅｔｅ Ｆｏｕｒｉｅｒ Ｔｒａｎｓｆｏｒｍ）を用いた場合、０．１Ｈｚの分解能を得るためには１０秒間の区間が必要であり、得られる呼吸周期はその平均的な値となってしまう。そのため、呼吸の１周期毎の微小な揺らぎを把握することが難しい。
そこで、本発明は、上記問題に鑑みてなされたものであり、本発明の目的とするところは、信号の周期をより高い精度で１周期毎に算出することが可能な、新規かつ改良された周期推定装置、周期推定方法及びプログラムを提供することにある。
上記課題を解決するために、本発明のある観点によれば、物体の運動の周期の情報を含む一次元信号と参照信号との位相差および前記位相差の時間変化に基づいて、前記一次元信号についての第１の周波数を推定する第１の周波数推定部と、前記第１の周波数に対応する第１の周期に相当する区間長を有する比較信号を前記一次元信号から切り出し、前記比較信号と前記区間長を有する前記一次元信号との相関係数を算出することにより、前記一次元信号についての第２の周波数を推定する第２の周波数推定部と、を備える、周期推定装置が提供される。
前記第１の周波数推定部は、推定された前記第１の周波数に応じた周波数を有する前記参照信号を生成してもよい
前記第１の周波数推定部は、推定された前記第１の周波数が所定の周波数帯域に存在しない場合、前記第１の周波数推定部で用いられた前記参照信号の周波数を、前記一次元信号についての前記第１の周波数として推定してもよい。
前記第２の周波数推定部は、前記第１の周期により定められる範囲内において前記相関係数が最大となる際のタイムラグに基づき、前記一次元信号についての第２の周波数を推定してもよい。
前記周期推定装置は、前記第１の周波数、前記第２の周波数および前記第２の周波数が推定された際の前記相関係数に基づき、前記第１の周波数または前記第２の周波数のいずれかを前記一次元信号の推定周波数として判定する判定部をさらに備えてもよい。
前記周期推定装置は、前記物体に対して放射波を放射し、前記放射波の周波数と、前記放射波が前記物体により反射した反射波の周波数との差分の周波数を有するビート信号を出力するドップラーセンサと、前記ドップラーセンサから出力された前記ビート信号を取得するビート信号取得部と、取得された前記ビート信号を、前記一次元信号に変換する信号変換部と、をさらに備えてもよい。
前記周期推定装置は、前記ビート信号から低周波成分を減少させ、前記低周波成分を減少させて得られた補正信号を出力するフィルタ部をさらに備えてもよい。
前記信号変換部は、二次元ベクトルとして表現される前記補正信号について主成分分析を行い、前記補正信号を、前記主成分分析から得られた固有ベクトルにより変換された前記補正信号の主成分方向の値を出力値とする前記一次元信号に変換してもよい。
前記信号変換部は、前記補正信号を、二次元平面上における、前記補正信号の強度と、前記補正信号の偏角の時間微分との積を変数とする連続関数の解を出力値とする前記一次元信号に変換してもよい。
前記物体の運動は、生体の呼吸運動であってもよい。
前記周期推定装置は、前記推定周波数に対応する周期に相当する区間長を有する基準信号を、前記一次元信号から切り出し、前記基準信号と、前記推定周波数および所定の初期位相を有する余弦波との位相差に基づいて、前記一次元信号のピーク位置を推定する、ピーク位置推定部をさらに備えてもよい。
前記周期推定装置は、前記ピーク位置推定部により前記ピーク位置を推定した際に用いた前記基準信号の区間において過去に推定された各ピーク位置を抽出し、抽出された前記各ピーク位置の分布に基づき、前記ピーク位置を確定する、ピーク位置確定部をさらに備えてもよい。
前記周期推定装置は、連続する前記確定された２つのピーク位置の差を、ピーク間隔として算出する、ピーク間隔算出部をさらに備えてもよい。
また、上記課題を解決するために、本発明の別の観点によれば、物体の運動の周期の情報を含む一次元信号と参照信号との位相差および前記位相差の時間変化に基づいて、前記一次元信号についての第１の周波数を推定するステップと、前記第１の周波数に対応する第１の周期に相当する区間長を有する比較信号を前記一次元信号から切り出し、前記比較信号と前記区間長を有する前記一次元信号との相関係数を算出することにより、前記一次元信号についての第２の周波数を推定するステップと、を含む、周期推定方法が提供される。
また、上記課題を解決するために、本発明の別の観点によれば、コンピュータを、物体の運動の周期の情報を含む一次元信号と参照信号との位相差および前記位相差の時間変化に基づいて、前記一次元信号についての第１の周波数を推定する第１の周波数推定部と、前記第１の周波数に対応する第１の周期に相当する区間長を有する比較信号を前記一次元信号から切り出し、前記比較信号と前記区間長を有する前記一次元信号との相関係数を算出することにより、前記一次元信号についての第２の周波数を推定する第２の周波数推定部と、として機能させるための、プログラムが提供される。
以上説明したように本発明によれば、信号の周期をより高い精度で１周期毎に算出することが可能である。
本発明の実施形態に係る周期推定装置の概要を示した図である。 同実施形態に係る周期推定装置の構成を示したブロック図である。 同実施形態に係る周波数推定部の構成を示したブロック図である。 同実施形態に係る周期推定装置の動作例を示した流れ図である。 同実施形態に係る周期推定装置が受信するビート信号を二次元平面上に示した図である。 同実施形態に係るフィルタ処理後に得られた補正信号の軌跡の一例を示した図である。 同実施形態に係るフィルタ処理後に得られた補正信号の軌跡の一例を示した図である。 同実施形態に係る第２の周波数推定部による周波数の推定方法の一例を示した図である。 同実施形態に係る判定部の動作例を示すフローチャートである。 同実施形態に係るピーク位置確定部によるピーク位置の確定方法の一例を示した図である。 同実施形態に係る周期推定装置のハードウェア構成を示したブロック図である。
以下に添付図面を参照しながら、本発明の好適な実施の形態について詳細に説明する。なお、本明細書及び図面において、実質的に同一の機能構成を有する構成要素については、同一の符号を付することにより重複説明を省略する。
＜１．概要＞
まず初めに、本発明の一実施形態に係る周期推定装置の概要について、図１を参照しながら説明する。図１は、本発明の一実施形態に係る周期推定装置の概要について説明する説明図である。
図１に示されているように、本実施形態に係る周期推定装置２０は、被験者１０に対して光、電磁波、または音波などの放射波を放射し、放射波を受けた被験者１０から反射した反射波を取得する。周期推定装置２０はドップラーセンサを有し、ドップラーセンサは、放射波の周波数と反射波の周波数との差分の周波数を有するビート信号（ドップラー信号）を出力する。周期推定装置２０は、この出力されたビート信号の周期を推定する。
周期推定装置２０は、例えば、光、電磁波、または音波などの放射波を用いるため、被験者１０に直接接触することなく、呼吸などの体動の周期を推定することができる。そのため、例えば周期推定装置２０は、被験者１０の呼吸などの生体情報の変化を、放射波の届く範囲において検出することが可能である。
なお、ドップラーセンサは、対象物（ここでは被験者１０）の移動速度に比例して反射波の周波数がシフトする現象（ドップラー効果）を利用して、対象物の運動を検出するセンサである。ドップラーセンサから出力されるビート信号は、ドップラー効果によりシフトされた周波数を有するため、ビート信号を解析することにより、対象物の運動情報を検出することが可能である。しかし、ここで出力されるビート信号は、例えば検出対象である呼吸だけではなく、身体の揺れ、心拍および周囲の環境雑音等を含む。そのため、例えば呼吸を検出対象とするとき、呼吸の波形情報のみをビート信号から検出することが必要となる。
例えば、上記特許文献１では、出力されたビート信号を高速フーリエ変換処理によりフーリエ解析を行い、呼吸に相当する周波数帯において最もピークの高い周波数を呼吸の周波数とし、その周波数に基づき呼吸数を検出する、という発明が記載されている。しかしフーリエ解析を用いて呼吸を検出する場合、フーリエ変換に必要なウインドウサイズが、少なくとも呼吸数回分の長さであるため、呼吸一回分のみの情報を検出することができなかった。また、フーリエ解析では、呼吸一回ごとの波形のピーク位置を特定することができないため、呼吸一回ごとに変化する周期の変化を検出することができなかった。
そこで、上記事情を一着眼点にして、周期推定装置２０を創作するに至った。本発明の実施形態に係る周期推定装置２０は、例えば呼吸などの信号波形の周期をより高い精度で１周期毎に算出することが可能である。以下、このような本実施形態に係る周期推定装置２０の構成についてより詳細に説明する。
＜２．構成＞
図２は、本実施形態に係る周期推定装置２０の構成を示したブロック図である。図２に示したように、周期推定装置２０は、ドップラーセンサ２１、信号加工部３０、および周期推定部４０を主に備える。信号加工部３０は、ビート信号取得部３１、フィルタ部３２、および信号変換部３３を含む。また、周期推定部４０は、周波数推定部４１、ピーク位置推定部４２、ピーク位置確定部４３、およびピーク間隔算出部４４を含む。
ドップラーセンサ２１は、被験者１０に対して、所定の周波数を有する放射波を放射し、放射波を受けた被験者１０から反射した反射波を検出する装置である。反射波は、ドップラー効果により、被験者１０の運動に応じて、放射波の周波数からシフトが生じた周波数を有する。ドップラーセンサ２１は、内部ミキサによって反射波と放射波とを積算し、積算した波を低域通過処理することにより、放射波の周波数と反射波の周波数の差分の周波数を有するビート信号を出力する。ドップラーセンサ２１は直交検波方式を搭載してもよく、その場合、ドップラーセンサ２１は、余弦波成分（Ｉ成分）および正弦波成分（Ｑ成分）の２種類のビート信号を出力する。なお、ここで用いられた添字のＩとＱとは、Ｉｎ−ｐｈａｓｅ：同相、Ｑｕａｄｒａｔｕｒｅ：直角位相を示す。
ビート信号取得部３１は、ドップラーセンサ２１から出力されたビート信号を取得する機能を有する。取得されたビート信号は、フィルタ部３２へと出力される。
フィルタ部３２は、ビート信号取得部３１より出力されたビート信号に含まれる直流成分などの低周波成分を減少させるフィルタ機能を有する。そして、フィルタ部３２は、当該フィルタ機能により低周波成分が減少したビート信号を、補正信号として信号変換部３３に出力する。なお、フィルタ部３２は、例えばハイパスフィルタやバンドパスフィルタ、ＩＩＲフィルタなどの多様なフィルタを用いることができ、またこれらのフィルタの組み合わせを用いてもよい。
信号変換部３３は、フィルタ部３２から出力された補正信号を、一次元信号に変換する機能を有する。補正信号を変換して得られた一次元信号は、周波数推定部４１へ出力される。信号変換部３３の具体的な機能については、動作例において詳しく説明する。なお、ここで用いられる一次元信号とは、例えば、時系列信号である。より詳細に説明すると、一次元信号とは、時間変化に伴い変化する単一の信号である。そして、単一の信号は、実数・虚数・複素数などを含み、以下の説明においては、複数の一次元信号が、まとめて一次元信号と表現されることもある。
周波数推定部４１は、一次元信号の有する周波数を推定する機能を有する。図３は周波数推定部４１の構成を示すブロック図を説明する図である。周波数推定部４１は、第１の周波数推定部４１１、第２の周波数推定部４１２、および判定部４１３を含む。まず、一次元信号は、第１の周波数推定部４１１、および第２の周波数推定部４１２に入力される。そして、判定部４１３は、第１の周波数推定部４１１で推定された第１の周波数、または第２の周波数推定部４１２で推定された第２の周波数のいずれかを、一次元信号の推定周波数として判定し、ピーク位置推定部４２へ出力する。また、第１の周波数推定部４１１で推定された第１の周波数は、第１の周波数推定部４１１および第２の周波数推定部４１２へとフィードバックされる。具体的な機能については、動作例において詳しく説明する。
ピーク位置推定部４２は、周波数推定部４１から出力された上記推定周波数を用いて一次元信号のピーク位置を推定する機能を有する。具体的な機能については、動作例において詳しく説明する。
ピーク位置確定部４３は、ピーク位置推定部４２により推定されたピーク位置について、過去の推定ピーク位置の分布を利用して、一次元信号のピーク位置を確定する機能を有する。具体的な機能については、動作例において詳しく説明する。
ピーク間隔算出部４４は、ピーク位置確定部４３において得られた連続する２つの確定ピーク位置の差を算出する機能を有する。ピーク間隔算出部４４で得られたピーク位置の差は、物体の運動の周期として出力される。
なお、周期推定装置２０は、上述した各機能部において実施された処理により得られた信号や周波数などの情報を記憶する記憶部（不図示）を有してもよい。具体的には、記憶部は、各機能部において得られた信号や周波数等の情報を記憶し、また、各機能部は、その機能を発揮する際に必要な情報を記憶部から呼び出す。例えば、信号変換部３３およびピーク位置確定部４３は、後述する動作例において、記憶部に記憶された過去の信号やピーク位置の情報を用いることにより、ビート信号を補正信号に変換し、および一次元信号のピーク位置を確定することができる。
以上、本実施形態に係る周期推定装置２０の構成を説明した。本実施形態に係る周期推定装置２０は、ドップラーセンサにより得られたビート信号を一次元信号へと変換し、一次元信号の周波数を２種類の手段で推定し、その後一次元信号のピーク位置を確定し、一周期ごとの周期を算出する。これにより、周期推定装置２０は、呼吸などの運動の周期を一周期ごとに精度よく検出することができる。
＜３．動作例＞
次に、図４〜１０を参照しながら、本実施形態に係る周期推定装置２０の動作例について説明する。図４は、本実施形態に係る周期推定装置２０の動作例を示す流れ図である。以下、「ビート信号の取得と変換」、「周波数の推定」、「ピーク位置の推定」の３段階に分けて、本実施形態に係る周期推定装置２０の動作例について説明する。
［ビート信号の取得と変換］
まず、ドップラーセンサ２１は、放射波と、放射波を受けた被験者１０から反射した反射波に基づきビート信号Ｄ（ｔ）を生成し、当該ビート信号Ｄ（ｔ）を出力する（Ｓ１０１）。ドップラーセンサ２１は、光、電磁波、または音波等の放射波を放射し、放射波を受けた被験者１０から反射された反射波を受信する。ドップラーセンサ２１は、受信した反射波と放射波をミキシングすることにより、放射波の周波数と反射波の周波数との差分の周波数を有するビート信号Ｄ（ｔ）を生成する。ビート信号Ｄ（ｔ）はＩ成分とＱ成分の２波の成分を有し、振幅をＡ（ｔ）、波長をλ、時刻ｔにおけるドップラーセンサ２１と被験者１０との距離をｄ（ｔ）、初期位相をφ ０ 、直流成分をＯ、ノイズ成分をｗとすると、ビート信号Ｄ（ｔ）は、下記の数式１のように表現される。
Figure 0006536038
続いて、ビート信号取得部３１は、ドップラーセンサ２１より出力されたビート信号Ｄ（ｔ）を取得する（Ｓ１０２）。取得されたビート信号Ｄ（ｔ）は、フィルタ部３２へと出力される。なお、ビート信号取得部３１は、取得されたビート信号Ｄ（ｔ）を時系列に記憶部に記憶してもよい。
図５は、ビート信号Ｄ（ｔ）をＩ成分とＱ成分からなる二次元平面上に示した図である。ビート信号Ｄ（ｔ）の位相変化は、直流成分Ｏを中心とした円１００の回転で表される。ここで、ベクトル１０１は、２次元平面における直流成分Ｏを始点としたビート信号Ｄ（ｔ）のベクトルである。被験者１０の周囲に存在する壁などの静止物で反射した反射波は、放射波と同一の周波数を有するため、ドップラー効果が生じない。そのため、放射波と反射波の定常の位相差のみを有する直流成分Ｏが、ビート信号Ｄ（ｔ）のバックグラウンドとして存在し、Ｉ平面とＱ平面からなる二次元平面上において、ベクトル１０１の始点として示されている。被験者１０の呼吸運動により、ベクトル１０１は、円弧１０２上に沿って振動する。円弧１０２の長さは、被験者１０とドップラーセンサ２１との距離の変化の大きさに比例する。また、ベクトル１０１の変化の速さは、呼吸運動の速度に比例する。
次に、フィルタ部３２は、ビート信号取得部３１において取得されたビート信号Ｄ（ｔ）に含まれる直流成分Ｏなどの低周波成分を減少させ、低周波成分を減少させた補正信号を出力する（Ｓ１０３）。フィルタ部３２は、ビート信号Ｄ（ｔ）から低周波成分等を減少させることにより、環境雑音、身体の体動など、呼吸運動とは無関係の成分を除去する。フィルタ部３２で用いられるフィルタは、例えばハイパスフィルタやバンドパスフィルタ、ＩＩＲフィルタなどの多様なフィルタであり、またこれらのフィルタを組み合わせたものを使用してもよい。なお、補正信号は、詳しくは後述するが、例えば、図６や図７で示されているように、Ｉ成分とＱ成分からなる二次元平面上において、様々な軌跡を描く。
次に、信号変換部３３は、フィルタ部３２から出力された補正信号を、一次元信号ｒ（ｔ）に変換する（Ｓ１０４）。より具体的には、信号変換部３３は、出力された補正信号に応じて、次に示す第１の変換手段、または第２の変換手段を用いて、もしくは第１の変換手段と第２の変換手段を組み合わせて、一次元信号ｒ（ｔ）を生成する。二次元の信号である補正信号を一次元化することにより、後処理である周期推定をより容易に実施できるという効果がある。
図６および図７は、いずれもフィルタ部３２から出力されたビート信号Ｄ（ｔ）の補正信号の軌跡を、Ｉ成分とＱ成分からなる二次元平面上で表現した図である。図６は、呼吸運動が小さい場合に、つまり、図５で示されている円弧１０２が円１００の円周上を一周しない場合に出力される補正信号の軌跡１０３を示している。軌跡１０３は、原点を通過する８の字形を形成している。補正信号が図６で示されている軌跡１０３を描く場合、信号変換部３３は、第１の変換手段として、補正信号の軌跡１０３を、原点を基準とした二次元ベクトルとして捉え、前記補正信号について主成分分析を行う。より詳細に説明すると、信号変換部３３は、軌跡１０３上に存在する過去の呼吸１周期分程度の補正信号を記憶部から呼び出し、各補正信号の二次元ベクトルから分散共分散行列を構成し、当該行列の固有値・固有ベクトルを算出する。そして、信号変換部３３は、補正信号と上記固有ベクトルとを積算して得られた信号の出力値を、一次元信号ｒ（ｔ）として出力する。図６においては、軌跡１０３の主成分方向が、破線の軸方向であることが示されている。信号変換部３３は、この軌跡１０３の主成分方向の値を連続的に出力し、一次元信号ｒ（ｔ）を生成する。なお、主成分分析において、主成分方向は、第一主成分と第二主成分の２種類存在するが、信号変換部３３は、各主成分方向に対応する固有値の大きい方を、一次元信号ｒ（ｔ）の生成に用いる主成分方向として採用する。
なお、主成分分析の方法は、補正信号のベクトルの分散共分散行列における固有値・固有ベクトルを算出する方法以外の分析方法でもよく、例えば、ベクトルのＩ成分およびＱ成分において分散が最大となるような主成分を算出する方法でもよい。
一方、図７は、呼吸運動が大きい場合に出力される補正信号の軌跡１０４を示している。深い呼吸など、呼吸運動が大きい場合、図５におけるベクトル１０１が描く円弧１０２は、円１００に沿って一周以上回転する場合が存在する。その場合、例えば、図７に示されているような補正信号の軌跡１０４が得られる。この軌跡１０４に対して上述した主成分分析を実施しても、主成分を正確に検出することができない。なぜなら、補正信号が軌跡１０４のような円状の分布を示す場合、主成分分析で得られる固有値が低くなり、主成分方向の出力値の重みが小さくなるためである。そこで、信号変換部３３は、主成分分析により一次元信号ｒ（ｔ）の正確な波形を得ることが困難である場合、第２の変換手段として、補正信号を、補正信号の二次元平面上における信号強度ｐ（ｔ）および偏角θ（ｔ）の時間微分θ’（ｔ）の積を変数とする関数を用いて、一次元信号ｒ（ｔ）に変換する。具体的には、信号変換部３３は、補正信号を、以下の数式２における連続関数（本実施形態においては、シグモイド関数）を利用した関数を用いて一次元信号ｒ（ｔ）を出力値として変換する。なお、数式２におけるａ、Ｋ、およびｃは定数であるが、必要に応じて変更が可能である。
Figure 0006536038
信号強度ｐ（ｔ）と偏角θ（ｔ）の時間微分θ’（ｔ）の積は、図７で示されている領域１０５の面積速度に相当する。面積速度は呼吸運動の大きさに相当する。時間微分θ’（ｔ）は、呼吸運動の方向により正負の値を取るため、呼吸の振動方向を表現できる。そのため、ｐ（ｔ）とθ’（ｔ）の積により、呼吸運動の状態を的確に表現することが可能となる。
信号変換部３３は、上記２つの変換手段を用いて、補正信号を一次元信号ｒ（ｔ）に変換する。ここで、信号変換部３３は、２つの変換手段のいずれか一方だけを用いてもよいし、２つの変換手段を組み合わせて用いてもよい。具体的には、本実施形態においては、信号変換部３３がいずれの変換手段を用いるかを、呼吸運動の大きさ、つまり、図５に示されている円弧１０２の長さによって判断してもよい。例えば、呼吸運動が小さい場合、つまり、円弧１０２の軌跡の長さが、円１００の円周上を一周しない程度の長さである場合は、補正信号の信号強度ｐ（ｔ）および偏角θ（ｔ）の時間微分θ’（ｔ）がいずれも小さくなるため、第２の変換手段では一次元信号ｒ（ｔ）を的確に表現できない。そのため、呼吸運動が小さい場合は、第１の変換手段を用いる。一方、呼吸運動が大きい場合、つまり、ベクトル１０１で示されている円弧１０２の軌跡が、円１００の円周上を一周以上回転して振動する場合は、第１の変換手段である主成分分析を用いて一次元信号ｒ（ｔ）を正確に得ることが困難であるため、第２の変換手段を用いることにより一次元信号ｒ（ｔ）が得られる。信号変換部３３における、第１の変換手段と第２の変換手段のどちらを用いるかの判断は、例えば、ビート信号Ｄ（ｔ）の軌跡を示す円弧１０２の長さに基づいて行われてもよいし、第１の変換手段で得られる固有値に基づいて行われてもよい。
［周波数の推定］
次に、本実施形態に係る一次元信号ｒ（ｔ）の周波数の推定に関する動作例について説明する。本動作例では、周波数推定部４１が、一次元信号ｒ（ｔ）の周波数を推定する（Ｓ１０５）。まず、信号変換部３３から出力された一次元信号ｒ（ｔ）は、周波数推定部４１に入力される。具体的には、一次元信号ｒ（ｔ）は、図３に示した第１の周波数推定部４１１および第２の周波数推定部４１２へと入力される。
まず、第１の周波数推定部４１１について説明する。ここで、周期を推定する具体的なタイミング時刻はｔｋとする。第１の周波数推定部４１１は、時刻ｔｋにおける入力された一次元信号ｒ（ｔｋ）と、参照信号とを積算し、積算された信号の低周波成分を、ローパスフィルタ等を用いて抽出することにより、一次元信号ｒ（ｔk）と参照信号との位相差ψ（ｔk）を算出する。なお、ここで用いる参照信号とは、１回の呼吸による信号変化を示すモデルとなる信号であり、例えば、理想的な呼吸信号としてのモデルや、近似的な呼吸信号としてのモデルがあり得る。本実施形態においては、位相変化による周期推定を容易とするため、位相とその時間変化をパラメータにもつ正弦波、あるいは余弦波が参照信号として用いられる。すなわち、本実施形態の参照信号は、第１の周波数推定部４１１によって時刻ｔk−１に推定された第１の推定周波数ｆＡ（ｔk−１）を有する正弦波、あるいは余弦波の信号である。
ここで、第１の周波数推定部４１１は、位相差ψ（ｔk）と、時刻ｔk−１に算出された位相差ψ（ｔk−１）との変化分を算出する。位相差の時間変化分ψ（ｔk）−ψ（ｔk−１）は、一次元信号ｒ（ｔ）の周波数の変化によって生じたものと仮定して、第１の周波数推定部４１１は、参照信号の周波数ｆＡ（ｔk−１）に、位相差の時間変化分ψ（ｔk）−ψ（ｔk−１）に相当する周波数を加算して、第１の推定周波数ｆＡ（ｔk）を出力する。
なお、上記位相差の時間変化分が所定の閾値以上であった場合、周波数のフィードバック値が大きく見積もられやすいため、第１の推定周波数ｆＡ（ｔk）が振動を伴って収束しないといった問題がある。その問題を参酌して、参照信号の周波数ｆＡ（ｔk−１）に加算する周波数は、位相差の時間変化分ψ（ｔk）−ψ（ｔk−１）に相当する周波数に、１未満の係数を積算したものでもよい。この処理によって、比較的短い時間でｆＡ（ｔｋ）が収束することが可能である。
また、算出された第１の推定周波数ｆＡ（ｔk）が、所定の周波数帯域に存在しない場合は、参照信号の周波数ｆＡ（ｔk−１）を、第１の推定周波数ｆＡ（ｔk）として出力してもよい。ここで所定の周波数帯域とは、例えば、本実施形態においては、呼吸運動の周期に相当する周波数帯域を意味する。具体的には、１分間の呼吸数が１２から６０であると仮定すると、第１の推定周波数ｆＡ（ｔk）は、０．２Ｈｚ＜ｆＡ（ｔk）＜１Ｈｚを満たす必要がある。
なお、この第１の周波数推定部４１１は、いわゆる位相同期回路（ＰＬＬ：ｐｈａｓｅ ｌｏｃｋｅｄ ｌｏｏｐ）により実現されてもよい。
続いて、第２の周波数推定部４１２について説明する。第２の周波数推定部４１２は、入力された一次元信号ｒ（ｔ）の相関分析を行うことにより、第１の周波数推定部４１１で得られた第１の周波数ｆＡ（ｔk）より精度の高い第２の周波数ｆＢ（ｔk）を推定する。
図８を参照しながら、第２の周波数推定部４１２における周波数推定手段について説明する。図８における曲線２００は、一次元信号ｒ（ｔ）である。まず、第２の周波数推定部４１２は、一次元信号ｒ（ｔ）から、時刻ｔk−ＴＡから時刻ｔkまでの区間を切り出し、これを比較信号として扱う。ＴＡは、先に推定された第１の周波数ｆＡ（ｔk−１）に対応する周期（第１の周期）である。次に、第２の周波数推定部４１２は、切り出された比較信号を用いて、相関係数を算出する。図８に示されているように、相関係数算出範囲２０１は、ｔk−ａ２ＴＡから、ｔk−ａ１ＴＡの範囲に相当する。ａ１およびａ２の値は任意に設定が可能であるが、ｔk−２ＴＡからｔk−ＴＡの範囲を含むことが望ましい。第２の周波数推定部４１２は、上記の相関係数算出範囲２０１において、区間長ＴＡを有する一次元信号ｒ（ｔ）と、比較信号２０２との相関係数Ｒ（τ）が最大となるタイムラグτｍａｘを算出し、また、τ＝τｍａｘのときの相関係数Ｒ（τｍａｘ）を記憶部に記憶する。なお、相関係数Ｒ（τ）は、以下の数式３を用いて算出される。なお、ＴｃおよびＴｄは比較信号の始点および終点の時刻であり、Ｔｄ −Ｔｃ＝ＴＡである。
Figure 0006536038
第２の周波数推定部４１２は、相関係数算出範囲２０１において算出される相関係数が最大となるときのタイムラグτｍａｘから、第２の推定周波数ｆＢ（ｔｋ）を算出する。ｆＢ（ｔk）は、以下の数式４を用いて算出される。
Figure 0006536038
ここで得られる第２の推定周波数ｆＢ（ｔk）は、先に得られた第１の推定周波数ｆＡ（ｔk）に相当する周期に基づく範囲内において、一次元信号ｒ（ｔ）の波形の相関分析により推定されるため、第２の周波数推定部４１２は、単純な正弦波、あるいは余弦波との位相差に基づき得られた第１の推定周波数ｆＡ（ｔk）より、高い精度で一次元信号ｒ（ｔ）の周波数を推定することが可能である。しかし、例えば、呼吸の波形が呼吸一回ごとに著しく変動する場合、第２の周波数推定部４１２において得られた第２の推定周波数ｆＢ（ｔk）は、第１の推定周波数ｆＡ（ｔk）と乖離している可能性があり、また、得られた第２の推定周波数ｆＢ（ｔk）が、第１の推定周波数ｆＡ（ｔk）と近い値であっても、第２の周波数推定部４１２において得られるＲ（τｍａｘ）が低い値を示す場合、第２の周波数推定部４１２における相関分析が、誤っている可能性も存在する。
そこで、判定部４１３は、第１の推定周波数ｆＡ（ｔk）、第２の推定周波数ｆＢ（ｔk）、および第２の推定周波数ｆＢ（ｔk）が推定された際の相関係数Ｒ（τｍａｘ）から、ｆＡ（ｔk）もしくはｆＢ（ｔk）のいずれかを推定周波数ｆＥ（ｔk）として採用するかを判定する。
図９は、判定部４１３の動作例を示したフローチャートである。まず、第１の推定周波数ｆＡ（ｔk）および第２の推定周波数ｆＢ（ｔk）が各周波数推定部にて推定される（Ｓ２０１）。推定された各周波数は、判定部４１３へ入力される。
次に、判定部４１３は、ｆＡ（ｔk）とｆＢ（ｔk）の差が、第１の閾値未満であるかどうかを判定する（Ｓ２０２）。ｆＡ（ｔk）とｆＢ（ｔk）の差が、第１の閾値未満でない場合（Ｓ２０２／ＮＯ）、判定部４１３は、時刻ｔｋにおけるＲ（τｍａｘ）を記憶部から呼び出し、Ｒ（τｍａｘ）が第２の閾値より大きいかどうかを判定する（Ｓ２０３）。ステップＳ２０２およびステップＳ２０３のいずれかにおいて判定結果がＹＥＳである場合は、判定部４１３は、ｆＢ（ｔk）が所定の周波数帯域に存在しているかどうか判定する（Ｓ２０４）。なお、この所定の周波数帯域は、上述したように、例えば、本実施形態においては、呼吸運動の周期に相当する周波数帯域を意味する。また、第１の閾値および第２の閾値は定数であるが、これらの定数は、必要に応じて変更が可能である。




ステップＳ２０２またはステップＳ２０３でＹＥＳと判定され、かつ、ステップＳ２０４でＹＥＳと判定された場合、判定部から出力される推定周波数ｆＥ（ｔk）には第２の推定周波数ｆＢ（ｔk）が代入される（Ｓ２０５）。逆に、ステップＳ２０２およびステップＳ２０３でＮＯと判定され、もしくは、ステップＳ２０４でＮＯと判定された場合、判定部から出力される推定周波数ｆＥ（ｔk）には第１の推定周波数ｆＡ（ｔk）が代入される（Ｓ２０６）。
なお、ステップＳ２０２およびステップＳ２０３において判定結果がいずれもＮＯであることは、第２の推定周波数ｆＢ（ｔk）が、第１の推定周波数ｆＡ（ｔk）と大きく乖離しており、かつ相関分析が適切に行われていないため、第２の推定周波数ｆＢ（ｔk）が実際の周波数と大きく異なっている可能性が高いことを意味する。その場合、判定部４１３は、第１の推定周波数ｆＡ（ｔk）を採用する。つまり、周波数推定部４１は、呼吸運動のような揺らぎを多く生じる波形の周期を推定するために、精度は低いがおおよその周波数を推定できる第１の周波数推定部４１１と、イレギュラーな変動には追従できないが、精度の高い周波数を推定できる第２の周波数推定部４１２とを相補的に用いることにより、呼吸波形を１周期ごとに高い精度で推定することを可能としている。
また、ステップＳ１０５で推定された第１の推定周波数ｆＡ（ｔk）は、記憶部に記憶される。そして、図３のブロック図に示されているように、記憶されたｆＡ（ｔk）は、時刻ｔｋ＋１において入力された一次元信号ｒ（ｔｋ＋１）の第１の推定周波数ｆＡ（ｔk＋１）を推定する際に、参照信号の周波数として、第１の周波数推定部４１１へフィードバックされる。さらに、記憶されたｆＡ（ｔk）は、第２の推定周波数ｆＢ（ｔk＋１）を推定する際に、相関分析を行うための情報として、第２の周波数推定部４１２へフィードバックされる。
なお、上記フィードバックに係る参照信号の周波数として、第１の推定周波数ｆＡ（ｔk）に１未満の係数を積算した値を用いてもよい。また、このような構成以外に、例えば、第２の推定周波数ｆＢ（ｔk）に１未満の係数を積算した値が、参照信号の周波数として用いられてもよい。また、上記の積算値のみでなく、例えば、第１の推定周波数ｆＡ（ｔk）や第２の推定周波数ｆＢ（ｔk）の前段や後段の処理によって得られる周波数が、参照信号の周波数として代用されてもよい。さらに、例えば、当該フィードバックされる周波数は、フーリエ変換など別の手段によってビート信号から直接的に求められた周波数でもよく、さらに加算処理および積算処理などによる上記周波数の組み合わせにより得られる周波数でもよい。
［ピーク位置の推定］
以上、本実施形態に係る周波数の推定に関する動作例について説明した。続いて、周波数推定部４１で得られた推定周波数ｆＥ（ｔk）に基づき、本実施形態に係る一次元信号ｒ（ｔ）のピーク位置を推定に関する動作例について説明する。
周波数推定部４１から出力された推定周波数ｆＥ（ｔk）は、続いてピーク位置推定部４２に入力される。ピーク位置推定部４２は、推定周波数ｆＥ（ｔk）に基づきピーク位置を推定する（Ｓ１０６）。具体的には、まず、ピーク位置推定部４２は、推定時刻ｔｋを終点とする、推定周波数ｆＥ（ｔk）に対応する周期ＴＥを区間長とする基準信号を、一次元信号ｒ（ｔ）から切り出す。ピーク位置推定部４２は、切り出された基準信号のピーク位置を推定する。より詳細に説明すると、ピーク位置推定部４２は、基準信号と同一の周期ＴＥを有する所定の初期位相（本実施形態においては、初期位相を０とする）の余弦波を生成し、当該余弦波と基準信号との位相差を算出する。ここで算出された位相差と時刻ｔｋから、一次元信号ｒ（ｔ）の基準信号区間内におけるピーク位置を推定することが可能である。
なお、ピーク位置推定部４２は、推定したピーク位置を記憶部に記憶することが可能である。記憶されたピーク位置は、次のピーク位置確定部４３で用いられる。
ここで、ピーク位置推定部４２において推定されるピーク位置は、必ずしも正確ではなく、誤差を含む場合が多い。そのため、推定されたピーク位置を補正し、ピーク位置を確定するための手段が必要となる。
ピーク位置確定部４３は、上記基準信号の区間内において推定されたピーク位置の分布から、統計手段によりピーク位置を確定する（Ｓ１０７）。より詳細に説明すると、ピーク位置確定部４３は、ピーク位置の推定を実施した時刻ｔｋを終点とする基準信号の区間ＴＥの範囲内において過去に推定された各ピーク位置を記憶部から抽出し、各ピーク位置の分布を得る。抽出された各ピーク位置の分布は、単位区間で仕切られるヒストグラムにおいて表現される。
図１０は、ピーク位置確定部４３によるピーク位置の確定方法の一例を説明するための説明図である。曲線３００は、一次元信号ｒ（ｔ）である。まず、ピーク位置推定部４２は、時刻ｔｋの時点において、位相差φ（ｔｋ）３０１を推定した、と仮定する。このときに推定された推定ピーク位置３０２は、他時刻に推定されたピーク位置３０３とは異なる位置が推定されたと仮定する。推定ピーク位置３０２に補正をかけなければ、推定ピーク位置３０２がこのまま出力される。
そこで、ピーク位置確定部４３は、ｔｋを終点とする基準信号の区間ＴＥ内において推定されたピーク位置の情報を、記憶部から抽出する。具体的には、ピーク位置確定部４３は、ｔｋ−ＴＥからｔｋまでの区間において、ピーク位置が推定された時刻ｔｋ、ｔｋ−１、ｔｋ−２、ｔｋ−３、・・・における各推定ピーク位置の情報を記憶部から抽出し、抽出された各推定ピーク位置に基づくヒストグラム３０４を作成する。ここで、上記区間におけるヒストグラム３０４において、ピーク位置３０３が最頻値となっていると仮定する。
その後、ピーク位置確定部４３は、各推定ピーク位置の平均値を算出し、平均値を中心とし、任意の分散値を有するガウス関数と、ヒストグラム３０４を乗算することにより、ヒストグラム３０４を重み付けする。ここで重み付けられたヒストグラム３０４の最頻値が、確定されたピーク位置となる。ピーク位置確定部４３により、推定されたピーク位置に誤りがあったとしても、最も確からしいピーク位置が出力されることが可能である。なお、得られた確定ピーク位置は、記憶部に記憶される。
ピーク位置が確定した後、ピーク間隔算出部４４は、得られた確定ピーク位置に基づき、ピーク間隔を算出する（Ｓ１０８）。具体的には、ピーク間隔算出部４４は、ピーク位置確定部４３において得られた連続する２つの確定ピーク位置の差を算出する。確定ピーク位置は、記憶部から呼び出される。ここで得られたピーク間隔が、呼吸運動の１周期の大きさとして出力される。
ピーク間隔の出力後、引き続き周期推定を行う場合はステップＳ１０１へ戻り（Ｓ１０９／ＹＥＳ）、周期推定を行わない場合は装置の使用を終了する（Ｓ１０９／ＮＯ）。
以上、各図を参照しながら、本発明の実施形態に係る周期推定装置２０の動作例について説明した。
なお、上述した本発明の実施形態の動作例は、周期推定装置２０が取得された信号の周期をリアルタイムで推定することを前提に説明されたが、例えば、周期推定装置２０は、リアルタイム処理ではなく、記憶部に記憶された信号や周波数等の情報に基づき、バッチ処理として当該記憶された過去の信号の周期を推定することも可能である。
例えば、周期推定装置２０は、ステップＳ１０２において取得されたビート信号Ｄ（ｔ）を記憶部に時系列で記憶しておき、バッチ処理としてビート信号Ｄ（ｔ）の周期を推定することも可能である。また、他の具体例としては、周期推定装置２０は、バッチ処理として、ステップＳ１０５において、時刻ｔｋにおける第２の推定周波数ｆＢ（ｔｋ）を推定する際に、相関係数算出範囲２０１を、ｔｋよりも未来の時刻にまで広げることが可能である。
また、上述した本発明の実施形態の動作例は、検出対象を呼吸運動としたが、本発明はかかる例に限定されない。例えば、周期推定装置２０は、心拍や体動などのその他の生体運動の情報を検出してもよいし、様々な物体に生じる微小振動を検出することも可能である。
＜５．ハードウェア構成＞
次に、図１１を参照しながら、本発明の実施形態に係る周期推定装置２０のハードウェア構成について説明する。本発明の実施形態に係る周期推定装置２０による情報処理は、ソフトウェアと、周期推定装置２０との協働により実現される。
周期推定装置２０は、主に、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）５０１と、ＲＯＭ（Ｒｅａｄ Ｏｎｌｙ Ｍｅｍｏｒｙ）５０２と、ＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）５０３と、ホストバス５０４と、を備える。また、周期推定装置２０は、ブリッジ５０５と、外部バス５０６と、インタフェース５０７と、入力装置５０８と、出力装置５０９と、ストレージ装置５１０と、ドライブ５１１と、ネットワークインタフェース５１２と、を備える。
ＣＰＵ５０１は、演算処理装置および制御装置として機能し、各種プログラムに従って周期推定装置２０内の動作全般を制御する。また、ＣＰＵ５０１は、マイクロプロセッサであってもよい。ＲＯＭ５０２は、ＣＰＵ５０１が使用するプログラムや演算パラメータ等を記憶する。ＲＡＭ５０３は、ＣＰＵ５０１の実行において使用するプログラムや、その実行において適宜変化するパラメータ等を一時記憶する。これらはＣＰＵバスなどから構成されるホストバス５０４により相互に接続されている。
ホストバス５０４は、ブリッジ５０５を介して、ＰＣＩ（Ｐｅｒｉｐｈｅｒａｌ Ｃｏｍｐｏｎｅｎｔ Ｉｎｔｅｒｃｏｎｎｅｃｔ／Ｉｎｔｅｒｆａｃｅ）バスなどの外部バス５０６に接続されている。なお、必ずしもホストバス５０４、ブリッジ５０５および外部バス５０６を分離構成する必要はなく、１つのバスにこれらの機能を実装してもよい。
入力装置５０８は、マウス、キーボード、タッチパネル、ボタン、マイクロフォン、スイッチおよびレバーなどユーザが情報を入力するための入力手段と、ユーザによる入力に基づいて入力信号を生成し、ＣＰＵ５０１に出力する入力制御回路などから構成されている。周期推定装置２０のユーザは、当該入力装置５０８を操作することにより、周期推定装置２０に対して各種のデータを入力したり処理動作を指示したりすることができる。
出力装置５０９は、例えば、ＣＲＴディスプレイ装置、液晶ディスプレイ（ＬＣＤ）装置、ＯＬＥＤ装置およびランプなどの表示装置を含む。さらに、出力装置５０９は、スピーカ及びヘッドホンなどの音声出力装置を含む。出力装置５０９は、例えば、再生されたコンテンツを出力する。具体的には、表示装置は再生された映像データ等の各種情報をテキストまたはイメージで表示する。一方、音声出力装置は、再生された音声データや表示装置に表示されたテキストデータ等を音声に変換して出力する。
ストレージ装置５１０は、本実施形態に係る周期推定装置２０におけるデータ格納用の装置である。ストレージ装置５１０は、記憶媒体、記憶媒体にデータを記録する記録装置、記憶媒体からデータを読み出す読み出し装置および記憶媒体に記憶されたデータを削除する削除装置などを含んでも良い。ストレージ装置は、例えば、ＨＤＤ（Ｈａｒｄ Ｄｉｓｃ Ｄｒｉｖｅ）やＳＳＤ（Ｓｏｌｉｄ Ｓｔａｔｅ Ｄｒｉｖｅ）で構成される。このストレージ装置５１０は、ＣＰＵ５０１が実行するプログラムや各種データを格納する。
ドライブ５１１は、記憶媒体用リーダライタであり、周期推定装置２０に内蔵、あるいは外付けされる。ドライブ５１１は、装着されている磁気ディスク、光ディスク、光磁気ディスク、または半導体メモリ等のリムーバブル記憶媒体５６に記録されている情報を読みだして、ＲＡＭ５０３に出力する。また、ドライブ５１１は、リムーバブル記憶媒体５６に情報を書き込むこともできる。
ネットワークインタフェース５１２は、例えば、他の装置に接続するための通信デバイス等で構成された通信インタフェースである。また、ネットワークインタフェース５１２は、例えば、無線ＬＡＮ（Ｌｏｃａｌ Ａｒｅａ Ｎｅｔｗｏｒｋ）対応通信装置、ＬＴＥ（Ｌｏｎｇ Ｔｅｒｍ Ｅｖｏｌｕｔｉｏｎ）対応通信装置、近距離無線通信技術対応通信装置、および有線による通信を行うワイヤー通信装置であってもよい。
なお、本発明の実施形態においては、周期推定装置２０は、例えば、ＣＰＵ５０１、ＲＯＭ５０２、ＲＡＭ５０３、ストレージ装置５１０等により実現される。
以上、本発明の実施形態に係る周期推定装置２０の機能を実現可能なハードウェア構成の一例を示した。上記の各構成要素は、汎用的な部材を用いて構成されていてもよいし、各構成要素の機能に特化したハードウェアにより構成されていてもよい。従って、本実施形態を実施する時々の技術レベルに応じて、適宜、利用するハードウェア構成を変更することが可能である。
＜５．まとめ＞
参照信号との比較による周波数の推定では、呼吸信号を近似的に表現しているため、正確な周期は特定できないものの、波形の概形が維持できていることから、大きくは外れない周波数と位相の推定が期待できる。また相関による周期の推定では、周期が変化する環境下では位相の概念を持たせることが困難で、多数倍の周期への誤検知や、微細構造へのフィッティングなどによる周波数誤検知が発生しうるが、周期間の波形形状が十分に安定している場合において正しい周期で相関係数がピークを形成することが期待できる。そのため、この両者を組み合わせた本発明の実施形態によれば、一次元信号と参照信号との位相差および当該位相差の時間変化に基づいて推定される第１の推定周波数と、第１の推定周波数に対応する周期に相当する区間長を有する比較信号と一次元信号との相関分析によって算出される第２の推定周波数とを相補的に採用することにより、誤検知が低減された高精度な周期推定が可能となり、呼吸運動のような揺らぎが多く生じる波形の周期を、１周期ごとに高精度に推定することが可能となる。
以上、添付図面を参照しながら本発明の好適な実施形態について詳細に説明したが、本発明はかかる例に限定されない。本発明の属する技術の分野における通常の知識を有する者であれば、特許請求の範囲に記載された技術的思想の範疇内において、各種の変更例または修正例に想到し得ることは明らかであり、これらについても、当然に本発明の技術的範囲に属するものと了解される。
また、周期推定装置２０に内蔵されるＣＰＵ５０１、ＲＯＭ５０２およびＲＡＭ５０３などのハードウェアを、上述した周期推定装置２０の各構成と同等の機能を発揮させるためのコンピュータプログラムも作成可能である。また、該コンピュータプログラムを記憶させた記憶媒体も提供される。
なお、本明細書において、流れ図やフローチャートに記載されたステップは、記載された順序に沿って時系列的に行われる処理はもちろん、必ずしも時系列的に処理されなくとも、並列的にまたは個別的に実行される処理をも含む。また時系列的に処理されるステップでも、場合によっては適宜順序を変更することが可能であることは言うまでもない。
２０ 周期推定装置
２１ ドップラーセンサ
３０ 信号加工部
３１ ビート信号取得部
３２ フィルタ部
３３ 信号変換部
４０ 周期推定部
４１ 周波数推定部
４２ ピーク位置推定部
４３ ピーク位置確定部
４４ ピーク間隔算出部
４１１ 第１の周波数推定部
４１２ 第２の周波数推定部
４１３ 判定部

Claims (15)

  1. 物体の運動の周期の情報を含む一次元信号と参照信号との位相差および前記位相差の時間変化に基づいて、前記一次元信号についての第１の周波数を推定する第１の周波数推定部と、
    前記第１の周波数に対応する第１の周期に相当する区間長を有する比較信号を前記一次元信号から切り出し、前記比較信号と前記区間長を有する前記一次元信号との相関係数を算出することにより、前記一次元信号についての第２の周波数を推定する第２の周波数推定部と、
    を備える、周期推定装置。
  2. 前記第１の周波数推定部は、推定された前記第１の周波数に応じた周波数を有する前記参照信号を生成する、請求項１に記載の周期推定装置。
  3. 前記第１の周波数推定部は、推定された前記第１の周波数が所定の周波数帯域に存在しない場合、前記第１の周波数推定部で用いられた前記参照信号の周波数を、前記一次元信号についての前記第１の周波数として推定する、請求項１または２に記載の周期推定装置。
  4. 前記第２の周波数推定部は、前記第１の周期により定められる範囲内において前記相関係数が最大となる際のタイムラグに基づき、前記一次元信号についての第２の周波数を推定する、請求項１〜３のいずれか１項に記載の周期推定装置。
  5. 前記周期推定装置は、前記第１の周波数、前記第２の周波数および前記第２の周波数が推定された際の前記相関係数に基づき、前記第１の周波数または前記第２の周波数のいずれかを前記一次元信号の推定周波数として判定する判定部をさらに備える、請求項１〜４のいずれか１項に記載の周期推定装置。
  6. 前記周期推定装置は、前記物体に対して放射波を放射し、前記放射波の周波数と、前記放射波が前記物体により反射した反射波の周波数との差分の周波数を有するビート信号を出力するドップラーセンサと、
    前記ドップラーセンサから出力された前記ビート信号を取得するビート信号取得部と、
    取得された前記ビート信号を、前記一次元信号に変換する信号変換部と、
    をさらに備える、請求項１〜５のいずれか１項に記載の周期推定装置。
  7. 前記周期推定装置は、前記ビート信号から低周波成分を減少させ、前記低周波成分を減少させて得られた補正信号を出力するフィルタ部をさらに備える、請求項６に記載の周期推定装置。
  8. 前記信号変換部は、二次元ベクトルとして表現される前記補正信号について主成分分析を行い、前記補正信号を、前記主成分分析から得られた固有ベクトルにより変換された前記補正信号の主成分方向の値を出力値とする前記一次元信号に変換する、請求項７に記載の周期推定装置。
  9. 前記信号変換部は、前記補正信号を、二次元平面上における、前記補正信号の強度と、前記補正信号の偏角の時間微分との積を変数とする連続関数の解を出力値とする前記一次元信号に変換する、請求項７に記載の周期推定装置。
  10. 前記物体の運動は、生体の呼吸運動である、請求項１〜９のいずれか１項に記載の周期推定装置。
  11. 前記周期推定装置は、前記推定周波数に対応する周期に相当する区間長を有する基準信号を、前記一次元信号から切り出し、前記基準信号と、前記推定周波数および所定の初期位相を有する余弦波との位相差に基づいて、前記一次元信号のピーク位置を推定する、ピーク位置推定部をさらに備える、請求項５に記載の周期推定装置。
  12. 前記周期推定装置は、前記ピーク位置推定部により前記ピーク位置を推定した際に用いた前記基準信号の区間において過去に推定された各ピーク位置を抽出し、抽出された前記各ピーク位置の分布に基づき、前記ピーク位置を確定する、ピーク位置確定部をさらに備える、請求項１１に記載の周期推定装置。
  13. 前記周期推定装置は、連続する前記確定された２つのピーク位置の差を、ピーク間隔として算出する、ピーク間隔算出部をさらに備える、請求項１２に記載の周期推定装置。
  14. 物体の運動の周期の情報を含む一次元信号と参照信号との位相差および前記位相差の時間変化に基づいて、前記一次元信号についての第１の周波数を推定するステップと、
    前記第１の周波数に対応する第１の周期に相当する区間長を有する比較信号を前記一次元信号から切り出し、前記比較信号と前記区間長を有する前記一次元信号との相関係数を算出することにより、前記一次元信号についての第２の周波数を推定するステップと、
    を含む、周期推定方法。
  15. コンピュータを、
    物体の運動の周期の情報を含む一次元信号と参照信号との位相差および前記位相差の時間変化に基づいて、前記一次元信号についての第１の周波数を推定する第１の周波数推定部と、
    前記第１の周波数に対応する第１の周期に相当する区間長を有する比較信号を前記一次元信号から切り出し、前記比較信号と前記区間長を有する前記一次元信号との相関係数を算出することにより、前記一次元信号についての第２の周波数を推定する第２の周波数推定部と、
    として機能させるための、プログラム。
JP2015008008A 2015-01-19 2015-01-19 周期推定装置、周期推定方法及びプログラム。 Active JP6536038B2 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
JP2015008008A JP6536038B2 (ja) 2015-01-19 2015-01-19 周期推定装置、周期推定方法及びプログラム。
US14/955,994 US20160206230A1 (en) 2015-01-19 2015-12-01 Period estimation apparatus, period estimation method and storage medium

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2015008008A JP6536038B2 (ja) 2015-01-19 2015-01-19 周期推定装置、周期推定方法及びプログラム。

Publications (2)

Publication Number Publication Date
JP2016131713A JP2016131713A (ja) 2016-07-25
JP6536038B2 true JP6536038B2 (ja) 2019-07-03

Family

ID=56406890

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2015008008A Active JP6536038B2 (ja) 2015-01-19 2015-01-19 周期推定装置、周期推定方法及びプログラム。

Country Status (2)

Country Link
US (1) US20160206230A1 (ja)
JP (1) JP6536038B2 (ja)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6558343B2 (ja) * 2016-10-17 2019-08-14 オムロン株式会社 生体センサ
JP6558350B2 (ja) * 2016-11-29 2019-08-14 オムロン株式会社 体動検出装置
JP7267002B2 (ja) * 2018-12-11 2023-05-01 日本無線株式会社 呼吸周期測定装置及び呼吸周期測定プログラム

Family Cites Families (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6937696B1 (en) * 1998-10-23 2005-08-30 Varian Medical Systems Technologies, Inc. Method and system for predictive physiological gating
CA2464029A1 (en) * 2004-04-08 2005-10-08 Valery Telfort Non-invasive ventilation monitor
JP4586465B2 (ja) * 2004-09-07 2010-11-24 横浜ゴム株式会社 作業性評価装置、作業性評価方法および作業性評価プログラム
JP4389821B2 (ja) * 2005-03-22 2009-12-24 ソニー株式会社 体動検出装置、コンテンツ再生装置、体動検出方法およびコンテンツ再生方法
JP2011519288A (ja) * 2008-04-03 2011-07-07 カイ メディカル、 インコーポレイテッド 非接触の生理的運動センサおよびその使用方法
EP2285276B1 (en) * 2008-05-09 2013-12-25 Philips Intellectual Property & Standards GmbH Contactless respiration monitoring of a patient
CA2789521A1 (en) * 2009-02-25 2010-09-02 Xanthia Global Limited Wireless physiology monitor
JP5140891B2 (ja) * 2009-06-09 2013-02-13 国立大学法人九州大学 信号ピーク測定システム
US20110021928A1 (en) * 2009-07-23 2011-01-27 The Boards Of Trustees Of The Leland Stanford Junior University Methods and system of determining cardio-respiratory parameters
JP5467395B2 (ja) * 2009-09-02 2014-04-09 株式会社産学連携機構九州 生体情報測定システム
JP5718126B2 (ja) * 2011-03-31 2015-05-13 沖電気工業株式会社 微細振動特徴量算出装置、微細振動特徴量算出方法及びプログラム
JP2013072865A (ja) * 2011-09-29 2013-04-22 Oki Electric Ind Co Ltd 検出装置、検出方法、およびプログラム
JP5974512B2 (ja) * 2012-02-01 2016-08-23 沖電気工業株式会社 情報処理装置、情報処理方法、及びプログラム
JP6035813B2 (ja) * 2012-03-28 2016-11-30 富士通株式会社 生体モニタリング装置、装置の制御方法、及び装置の制御プログラム
JP5477424B2 (ja) * 2012-07-02 2014-04-23 沖電気工業株式会社 物体検知装置、物体検知方法及びプログラム
JP6244178B2 (ja) * 2013-11-12 2017-12-06 沖電気工業株式会社 情報処理装置、情報処理方法及びプログラム

Also Published As

Publication number Publication date
JP2016131713A (ja) 2016-07-25
US20160206230A1 (en) 2016-07-21

Similar Documents

Publication Publication Date Title
JP6477199B2 (ja) 振動状態推定装置、振動状態推定方法、およびプログラム
JP6515670B2 (ja) 睡眠深度推定装置、睡眠深度推定方法、およびプログラム
US20170007137A1 (en) Method of estimating blood pressure based on image
JP6056389B2 (ja) 心拍推定装置、心拍推定方法及びプログラム
JP2013072865A (ja) 検出装置、検出方法、およびプログラム
US20100138191A1 (en) Method and system for acquiring and transforming ultrasound data
JP5718126B2 (ja) 微細振動特徴量算出装置、微細振動特徴量算出方法及びプログラム
JP6536038B2 (ja) 周期推定装置、周期推定方法及びプログラム。
JP2014171589A (ja) 心房細動解析装置およびプログラム
JP2017136164A (ja) センサ情報処理装置、センサユニット、及び、センサ情報処理プログラム
JP4326486B2 (ja) 周期算出装置及びプログラム
JP2008073077A (ja) データ処理装置，データ処理方法及びデータ処理プログラム
JP2017127398A (ja) 情報処理装置、情報処理システム、情報処理方法及びプログラム
JP6060563B2 (ja) 心房細動判定装置、心房細動判定方法およびプログラム
JP5974512B2 (ja) 情報処理装置、情報処理方法、及びプログラム
JP2017205145A (ja) 脈拍推定装置、脈拍推定システム、脈拍推定方法および脈拍推定プログラム
JP2017006540A (ja) 心拍間隔特定プログラム、心拍間隔特定装置、及び心拍間隔特定方法
JP5998516B2 (ja) 拍動検出装置、電子機器及びプログラム
US20130211273A1 (en) Method and apparatus for heart rate measurement
JP6316063B2 (ja) 情報処理装置、情報処理システム、情報処理方法、およびプログラム
JP6299172B2 (ja) 情報処理装置、情報処理方法及びプログラム
JP5035370B2 (ja) 運動検出装置、運動検出方法、及びプログラム
JP6642055B2 (ja) センサ情報処理装置、センサユニット、及び、センサ情報処理プログラム
JP4945745B2 (ja) 自己組織化マップを用いた波形解析システム、波形解析のための自己組織化マップの作成方法及び作成プログラム並びに脈波解析システム
JP5800776B2 (ja) 生体動情報検出装置

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20171120

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20180914

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20181016

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20181114

RD03 Notification of appointment of power of attorney

Free format text: JAPANESE INTERMEDIATE CODE: A7423

Effective date: 20190327

RD04 Notification of resignation of power of attorney

Free format text: JAPANESE INTERMEDIATE CODE: A7424

Effective date: 20190328

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20190507

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20190520

R150 Certificate of patent or registration of utility model

Ref document number: 6536038

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150