(57)【特許請求の範囲】
【請求項１】
入力信号が一定時間間隔でサンプリングされてＡＤ変換されたサンプルデータを用いて前記入力信号の周波数を求める信号周波数算出方法であって、
前記サンプルデータに所定の窓関数で演算を行う窓関数演算ステップと、
前記所定の窓関数で演算した結果に対して高速フーリエ変換演算を行い、サンプル点数分の離散スペクトル（ｋｎ：ｎは１からサンプル点数までの自然数）を算出するＦＦＴステップと、
前記離散スペクトルのうち、振幅が最大値となるピークスペクトルの周波数（ｋｍａｘ）と振幅値（Ｓｉｍａｘ）を抽出する抽出ステップと、
前記ピークスペクトルに隣接する両側の離散スペクトルの振幅値（Ｓｉｍａｘ−１，Ｓｉｍａｘ＋１）の差と前記ピークスペクトルの振幅値との比から周波数の補正値を算出するステップと、
前記ピークスペクトルと、前記補正値から前記入力信号の周波数を算出するステップと、を有し、前記入力信号の周波数を求めることを特徴とする信号周波数算出方法。
【請求項２】
前記所定の窓関数がハニング窓であることを特徴とする請求項１に記載の信号周波数算出方法。
【請求項３】
前記入力信号の周波数（ｆ）は、次の式（１）から求めることを特徴とする請求項２に記載の信号周波数算出方法。
ｆ＝ｋｍａｘ±｛−３＋√（９＋４ｄ２）｝／｜ｄ｜・・・式（１）
但し、ｄ＞０のとき＋、ｄ＜０のとき−。
また、ｄ＝ΔＰ／｜Ｓｉｍａｘ｜
ΔＰ＝｜Ｓｉｍａｘ−１｜−｜Ｓｉｍａｘ＋１｜
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、信号解析技術に関し、特に、入力された信号を周波数解析してピーク周波数を推定する周波数算出方法に関する。
【背景技術】
【０００２】
従来から、信号解析処理技術として、入力された信号の周波数を決定するために観測信号を離散化して周波数変換を行う、ピーク離散フーリエ変換や高速フーリエ変換（ＦＦＴ）を用いる方法が一般的に知られている。このような信号解析方法を用いた場合、ピーク周波数の分解能は、観測時間の逆数に制限されてしまうので、ピーク周波数の精度の向上を図る周波数算出方法に関する技術が提案されている。【０００３】
サンプリング点数で決定される周波数刻み以上の精度で周波数を決定する方法として補間法があるが、この補間法として非特許文献１（従来例１）に記載の方式がよく知られている。【０００４】
従来例１の周波数算出方法では、入力信号系列に対してハニング窓をかけてＦＦＴを行って離散スペクトル列Ｇ（ｋ）を求める。求めた離散スペクトルの中で振幅値が最大のスペクトルＧ（ｋｍａｘ）と、その低周波側に隣接するスペクトルＧ（ｋｍａｘ−１）または高周波側に隣接するスペクトルＧ（ｋｍａｘ＋１）を用いて、下記の式（４）または式（６）の演算を行うことにより入力信号の周波数ｆを求めるものである。【０００５】
（数１）
ｆ＝ｋｍａｘ＋（１−２ｒ）／（１＋ｒ）              ・・・式（４）
但し、ｒ＝｜Ｇ（ｋｍａｘ−１）｜／｜Ｇ（ｋｍａｘ）｜・・・式（５）
ｆ＝ｋｍａｘ−（１−２ｓ）／（１＋ｓ）              ・・・式（６）
但し、ｓ＝｜Ｇ（ｋｍａｘ＋１）｜／｜Ｇ（ｋｍａｘ）｜・・・式（７）【０００６】
また、特許文献１（従来例２）では、図５に示すように、窓関数をかけたビート信号の離散周波数スペクトルを求め、そのピーク周波数を挟む前後２つの離散周波数の信号強度Ｐ−１，Ｐ１  を検出し、その信号強度の比Δｐから、ビート信号の離散ピーク周波数と窓関数のピーク周波数との周波数ずれΔｆを求める。これによりビート信号の真のピーク周波数を求める周波数算出方法に関する技術が開示されている。
【先行技術文献】
【特許文献】
【０００７】
【特許文献１】特開２００３‐２４０８４２号公報
【非特許文献】
【０００８】
【非特許文献１】電子情報通信学会論文誌Ａ  １９８７年５月  Ｖｏｌ．Ｊ７０−Ａ  Ｎｏ．５  ｐ．７９８−８０５
【発明の概要】
【発明が解決しようとする課題】
【０００９】
しかしながら、上述した従来技術の周波数算出方法においては、信号強度の比を対数比で求めているため多くの演算時間が掛ってしまうという課題があった。【００１０】
本発明は、上述した課題を解決するもので、演算時間を短縮することができる周波数算出方法を提供することを目的とする。
【課題を解決するための手段】
【００１１】
この課題を解決するために、本発明の信号周波数算出方法は、入力信号が一定時間間隔でサンプリングされてＡＤ変換されたサンプルデータを用いて前記入力信号の周波数を求める信号周波数算出方法であって、前記サンプルデータに所定の窓関数で演算を行う窓関数演算ステップと、前記所定の窓関数で演算した結果に対して高速フーリエ変換演算を行い、サンプル点数分の離散スペクトル（ｋｎ：ｎは１からサンプル点数までの自然数）を算出するＦＦＴステップと、前記離散スペクトルのうち、振幅が最大値となるピークスペクトルの周波数（ｋｍａｘ）と振幅値（Ｓｉｍａｘ）を抽出する抽出ステップと、前記ピークスペクトルに隣接する両側の離散スペクトルの振幅値（Ｓｉｍａｘ−１，Ｓｉｍａｘ＋１）の差と前記ピークスペクトルの振幅値との比から周波数の補正値を算出するステップと、前記ピークスペクトルと、前記補正値から前記入力信号の周波数を算出するステップと、を有し、前記入力信号の周波数を求めることを特徴とする。【００１２】
これによれば、ピークスペクトルに隣接する両側の離散スペクトル（Ｓｉｍａｘ−１，Ｓｉｍａｘ＋１）の振幅の差とピークスペクトルの振幅値（Ｓｉｍａｘ）との比から周波数を算出することができるので、対数比の演算行うことがない。このため短い演算時間で、入力信号の周波数を算出する方法を提供することができる。【００１３】
また、本発明の信号周波数算出方法は、前記所定の窓関数がハニング窓であることを特徴とする。【００１４】
これによれば、窓関数にハニング窓を使用することで、離散スペクトルのサイドローブが低く抑えられるので、より精度良く入力信号の周波数を算出することができる。【００１５】
また、本発明の信号周波数算出方法は、前記入力信号の周波数は、次の式（１）から求めることを特徴とする。【００１６】
数（２）
ｆ＝ｋｍａｘ±｛−３＋√（９＋４ｄ２）｝／｜ｄ｜      ・・・式（１）
但し、ｄ＞０のとき＋、ｄ＜０のとき−。
また、ｄ＝ΔＰ／｜Ｓｉｍａｘ｜                      ・・・式（２）
ΔＰ＝｜Ｓｉｍａｘ−１｜−｜Ｓｉｍａｘ＋１｜        ・・・式（３）【００１７】
これによれば、式（１）に示される簡易な計算によって信号周波数を算出することができるので、演算量が少なく高速に入力信号の周波数を求めることができる。
【発明の効果】
【００１８】
本発明の信号周波数算出方法によれば、演算時間を短縮することができる周波数算出方法を提供することができる。
【図面の簡単な説明】
【００１９】
【図１】本発明の信号周波数算出方法を適用することができるセンサ装置の概略構成を示す図である。
【図２】本発明の実施形態に係る周波数算出方法を示すフローチャートである。
【図３】本発明の実施形態に係る周波数算出方法の手順について示す図である。
【図４】高速フーリエ変換で得られる離散周波数とその振幅値を示す図である。
【図５】特許文献１に記載の周波数算出方法を説明する図である。
【発明を実施するための形態】
【００２０】
［第１実施形態］
以下に第１実施形態における周波数算出方法について説明する。
【００２１】
図１は、本発明の信号周波数算出方法を適用することができるセンサ装置１００の概略構成を示す図で、図１（ａ）はセンサ装置１００の動作概要を示す図であり、図１（ｂ）はセンサ装置１００の構成を示すブロック図である。
【００２２】
センサ装置１００は図１（ａ）に示すように、送信波を放射し、検出対象８０からの反射波を検出し検出対象の検知を行う。
【００２３】
センサ装置１００は図１（ｂ）に示すように、アンテナ１０と、送信部２０と、検波部３０と、信号処理部４０と、制御部５０と、を備えている。
【００２４】
アンテナ１０は送信部２０と検波部３０と、に接続されており、送信部２０から出力される送信信号を送信波として空間に放射し、送信波が検出対象８０で反射された反射波を受信する。
【００２５】
検波部３０は、送信部２０からアンテナ１０を介して送信波を送信している間、検出対象８０で反射された反射波を送信信号の一部を参照信号として検波し、信号処理部４０へ出力する。
【００２６】
信号処理部４０は、検波部３０からの入力信号を増幅する機能、増幅した信号を一定時間間隔でサンプリングしてデジタル信号であるサンプルデータへ変換するＡＤ（Ａｎａｌｏｇ  ｔｏ  Ｄｉｇｉｔａｌ）変換機能、を有している。
【００２７】
信号処理部４０は、検波部３０から入力された信号を増幅し、ＡＤ変換によってサンプルデータに変換して制御部５０へ出力する。
【００２８】
制御部５０は、信号処理部４０から入力されたサンプルデータを用いて入力信号の周波数を求める。また制御部５０は、求めた周波数及びその変動から検出対象８０の有無や、その動きなどを検出し、検出した結果を表示したり、図示しない外部機器に制御信号を出力したりすることで、検出対象８０の検出結果を出力する。
【００２９】
次に制御部５０がサンプルデータを用いて入力信号の周波数を求める手順について図１から図４を用いて説明する。図２は、本発明の実施形態に係る周波数算出方法を示すフローチャートである。尚、図２に示す周波数算出方法の手順は、図１に示すセンサ装置１００の制御部５０が実行するものとして説明を進める。図３は、本発明の実施形態に係る周波数算出方法の手順について示す図である。図３（ａ）は信号処理部４０から入力したサンプルデータを時間波形として示している。図３（ｂ）は、窓関数（ハニング窓）のレスポンスを示す図である。図３（ｃ）は、図３（ａ）に示すサンプルデータに対して、図３（ｂ）に示す窓関数で演算を行うことにより、求められたデータ列を示す図である。
図４は、高速フーリエ変換（Ｆａｓｔ  Ｆｏｕｒｉｅｒ  Ｔｒａｎｓｆｏｒｍ以下ＦＦＴと記す）で得られる離散周波数とその振幅値の関係を示す図である。
【００３０】
図２に示す周波数算出の手順を示すフローチャートの手順Ｓ１は、サンプルデータ取得ステップである。制御部５０には信号処理部４０から、図３（ａ）に示すサンプルデータが入力されるので、制御部５０は図２に示す手順Ｓ１でこのサンプルデータを取得する。
【００３１】
手順Ｓ２は、窓関数演算ステップである。手順Ｓ２で制御部５０は、手順Ｓ１で取り込んだサンプルデータに、図３（ｂ）に示すあらかじめ記憶している所定の窓関数で演算を行うことで、図３（ｃ）に示すデータ列を求める。本実施形態では、所定の窓関数はハニング窓と称される窓関数を用いるものとして説明を行う。
【００３２】
手順Ｓ３はＦＦＴステップである。手順Ｓ３で制御部５０は、手順Ｓ２で演算した結果であるデータ列に対してＦＦＴ演算を行い、サンプル点数分の離散スペクトル（ｋｎ：ｎは１からサンプル点数までの自然数）を算出する。窓関数で演算したサンプルデータをＦＦＴ演算することにより、図４に示すような離散スペクトルを求めることができる。
【００３３】
手順Ｓ４は抽出ステップである。手順Ｓ３で求めた離散スペクトルのうち、振幅が最大値となるピークスペクトルの周波数（ｋｍａｘ）と振幅値（Ｓｉｍａｘ）と、ピークスペクトルに隣接する両側の離散スペクトルの振幅値（Ｓｉｍａｘ−１，Ｓｉｍａｘ＋１）を抽出する。尚図４では、ピークスペクトルの周波数ｋｍａｘが”０”となるよう作図しており、その振幅値Ｓｉｍａｘ及びｋｍａｘに隣接する両側の離散スペクトルの振幅値Ｓｉｍａｘ−１，Ｓｉｍａｘ＋１の関係が示されている。
【００３４】
手順Ｓ５は補正値を算出するステップである。手順Ｓ４で抽出したピークスペクトルに隣接する両側の離散スペクトルの振幅値（Ｓｉｍａｘ−１，Ｓｉｍａｘ＋１）の差ΔＰと、ピークスペクトルの振幅値（Ｓｉｍａｘ）との比から周波数の補正値を算出する。
【００３５】
手順Ｓ６は入力信号の周波数を算出するステップである。手順Ｓ４で抽出したピークスペクトルの周波数（ｋｍａｘ）と手順Ｓ５で算出した補正値とから入力信号の周波数を求める。
【００３６】
補正値は、前述の式（４）と式（６）からｆ−ｋｍａｘについて整理すると、次の式（８）を得るので、入力信号の周波数ｆは式（１０）で求めることができる。
【００３７】
（数３）
ｆ−ｋｍａｘ＝｛３±√（９＋４ｄ２）｝／ｄ      ・・・式（８）
但しｄ＝ｒ−ｓ      ・・・式（９）
ｆ＝ｋｍａｘ±｛−３＋√（９＋４ｄ２）｝／｜ｄ｜・・・式（１０）
但し、ｄ＞０のとき＋、ｄ＜０のとき−。
【００３８】
また、式（２）及び式（４）からΔＰを次の式（１１）に置き換えるとｄは式（１２）で表すことができる。
【００３９】
（数４）
ΔＰ＝｜Ｓｉｍａｘ−１｜−｜Ｓｉｍａｘ＋１｜    ・・・式（１１）
ｄ＝ΔＰ／｜Ｓｉｍａｘ｜                        ・・・式（１２）
【００４０】
次に、本実施形態としたことによる効果について説明する。
【００４１】
本実施形態の周波数算出方法では、サンプルデータに所定の窓関数で演算を行う窓関数演算ステップと、所定の窓関数で演算した結果に対して高速フーリエ変換演算を行い、サンプル点数分の離散スペクトル（ｋｎ：ｎは１からサンプル点数までの自然数）を算出するＦＦＴステップと、離散スペクトルのうち、振幅が最大値となるピークスペクトルの周波数（ｋｍａｘ）と振幅値（Ｓｉｍａｘ）を抽出する抽出ステップと、ピークスペクトルに隣接する両側の離散スペクトルの振幅値（Ｓｉｍａｘ−１，Ｓｉｍａｘ＋１）の差とピークスペクトルの振幅値との比から周波数の補正値を算出するステップと、ピークスペクトルと、補正値から前記入力信号の周波数を算出するステップと、を有し、前記入力信号の周波数を求めるようにした。
【００４２】
これにより、ピークスペクトルに隣接する両側の離散スペクトル（Ｓｉｍａｘ−１，Ｓｉｍａｘ＋１）の振幅の差ΔＰとピークスペクトルの振幅値（Ｓｉｍａｘ）との比から周波数を算出することができるので、対数での演算行うことがない。このため短い演算時間で、入力信号の周波数を算出する方法を提供することができる。
【００４３】
また、本実施形態の周波数算出方法では、窓関数にハニング窓を用いて周波数を求めるようにした。
【００４４】
これにより、離散スペクトルのサイドローブが低く抑えられるので、より精度良く入力信号の周波数を算出することができる。
【００４５】
また、本実施形態の周波数算出方法では、入力信号の周波数は、式（１）から求めるとした。
【００４６】
これにより、式（１）に示される簡易な計算によって信号周波数を算出することができるので、演算量が少なく高速に入力信号の周波数を求めることができる。
【００４７】
以上のように、本発明の実施形態に係る周波数算出方法を具体的に説明したが、本発明は上記の実施形態に限定されるものではなく、要旨を逸脱しない範囲で種々変更して実施することが可能である。例えば次のように変形して実施することができ、これらの実施形態も本発明の技術的範囲に属する。
【００４８】
＜変形例＞
（１）本実施形態において、センサ装置１００に備えられた制御部５０がサンプルデータを用いて入力信号の周波数を求める例を示して説明を行ったが、センサ装置とは別の外部機器がサンプルデータを取得して周波数を算出するようにしても良い。また、周波数を算出する手順をプログラムとして装置にプログラムの一部として実装しても良く、プログラムの供給形態に関わらず実施することが可能である。
【００４９】
（２）窓関数はハニング窓を用いて説明を行ったが、窓関数は必要に応じて変更しても良い。この場合、補正値は変更した窓関数に合わせて適宜変更して求めることが望ましい。
【００５０】
（３）式（１）の平方根内の式を展開した多項式にて演算することにより、平方根の演算が不要となるためさらに演算時間を短縮することが可能である。
【符号の説明】
【００５１】
１０    アンテナ
２０    送信部
３０    検波部
４０    信号処理部
５０    制御部
８０    検出対象
１００  センサ装置
【図１】
0006004510.tif 000002
【図２】
0006004510.tif 000003
【図３】
0006004510.tif 000004
【図４】
0006004510.tif 000005
【図５】
0006004510.tif 000006