JP6004510B2

JP6004510B2 - 信号周波数算出方法

Info

Publication number: JP6004510B2
Application number: JP2015522711A
Authority: JP
Inventors: 大滝　幸夫; 幸夫大滝; 遠藤　誠; 誠遠藤
Original assignee: Alps Electric Co Ltd
Current assignee: Alps Alpine Co Ltd
Priority date: 2013-06-17
Filing date: 2014-05-30
Publication date: 2016-10-12
Anticipated expiration: 2034-05-30
Also published as: WO2014203708A1; JPWO2014203708A1

Description

本発明は、信号解析技術に関し、特に、入力された信号を周波数解析してピーク周波数を推定する周波数算出方法に関する。

従来から、信号解析処理技術として、入力された信号の周波数を決定するために観測信号を離散化して周波数変換を行う、ピーク離散フーリエ変換や高速フーリエ変換（ＦＦＴ）を用いる方法が一般的に知られている。このような信号解析方法を用いた場合、ピーク周波数の分解能は、観測時間の逆数に制限されてしまうので、ピーク周波数の精度の向上を図る周波数算出方法に関する技術が提案されている。

サンプリング点数で決定される周波数刻み以上の精度で周波数を決定する方法として補間法があるが、この補間法として非特許文献１（従来例１）に記載の方式がよく知られている。

従来例１の周波数算出方法では、入力信号系列に対してハニング窓をかけてＦＦＴを行って離散スペクトル列Ｇ（ｋ）を求める。求めた離散スペクトルの中で振幅値が最大のスペクトルＧ（ｋｍａｘ）と、その低周波側に隣接するスペクトルＧ（ｋｍａｘ−１）または高周波側に隣接するスペクトルＧ（ｋｍａｘ＋１）を用いて、下記の式（４）または式（６）の演算を行うことにより入力信号の周波数ｆを求めるものである。

（数１）
ｆ＝ｋｍａｘ＋（１−２ｒ）／（１＋ｒ）・・・式（４）
但し、ｒ＝｜Ｇ（ｋｍａｘ−１）｜／｜Ｇ（ｋｍａｘ）｜・・・式（５）
ｆ＝ｋｍａｘ−（１−２ｓ）／（１＋ｓ）・・・式（６）
但し、ｓ＝｜Ｇ（ｋｍａｘ＋１）｜／｜Ｇ（ｋｍａｘ）｜・・・式（７）

また、特許文献１（従来例２）では、図５に示すように、窓関数をかけたビート信号の離散周波数スペクトルを求め、そのピーク周波数を挟む前後２つの離散周波数の信号強度Ｐ−１，Ｐ１を検出し、その信号強度の比Δｐから、ビート信号の離散ピーク周波数と窓関数のピーク周波数との周波数ずれΔｆを求める。これによりビート信号の真のピーク周波数を求める周波数算出方法に関する技術が開示されている。

特開２００３‐２４０８４２号公報

電子情報通信学会論文誌Ａ１９８７年５月Ｖｏｌ．Ｊ７０−ＡＮｏ．５ｐ．７９８−８０５

しかしながら、上述した従来技術の周波数算出方法においては、信号強度の比を対数比で求めているため多くの演算時間が掛ってしまうという課題があった。

本発明は、上述した課題を解決するもので、演算時間を短縮することができる周波数算出方法を提供することを目的とする。

この課題を解決するために、本発明の信号周波数算出方法は、入力信号が一定時間間隔でサンプリングされてＡＤ変換されたサンプルデータを用いて前記入力信号の周波数を求める信号周波数算出方法であって、前記サンプルデータに所定の窓関数で演算を行う窓関数演算ステップと、前記所定の窓関数で演算した結果に対して高速フーリエ変換演算を行い、サンプル点数分の離散スペクトル（ｋｎ：ｎは１からサンプル点数までの自然数）を算出するＦＦＴステップと、前記離散スペクトルのうち、振幅が最大値となるピークスペクトルの周波数（ｋｍａｘ）と振幅値（Ｓｉｍａｘ）を抽出する抽出ステップと、前記ピークスペクトルに隣接する両側の離散スペクトルの振幅値（Ｓｉｍａｘ−１，Ｓｉｍａｘ＋１）の差と前記ピークスペクトルの振幅値との比から周波数の補正値を算出するステップと、前記ピークスペクトルと、前記補正値から前記入力信号の周波数を算出するステップと、を有し、前記入力信号の周波数を求めることを特徴とする。

これによれば、ピークスペクトルに隣接する両側の離散スペクトル（Ｓｉｍａｘ−１，Ｓｉｍａｘ＋１）の振幅の差とピークスペクトルの振幅値（Ｓｉｍａｘ）との比から周波数を算出することができるので、対数比の演算行うことがない。このため短い演算時間で、入力信号の周波数を算出する方法を提供することができる。

また、本発明の信号周波数算出方法は、前記所定の窓関数がハニング窓であることを特徴とする。

これによれば、窓関数にハニング窓を使用することで、離散スペクトルのサイドローブが低く抑えられるので、より精度良く入力信号の周波数を算出することができる。

また、本発明の信号周波数算出方法は、前記入力信号の周波数は、次の式（１）から求めることを特徴とする。

数（２）
ｆ＝ｋｍａｘ±｛−３＋√（９＋４ｄ^２）｝／｜ｄ｜・・・式（１）
但し、ｄ＞０のとき＋、ｄ＜０のとき−。
また、ｄ＝ΔＰ／｜Ｓｉｍａｘ｜・・・式（２）
ΔＰ＝｜Ｓｉｍａｘ−１｜−｜Ｓｉｍａｘ＋１｜・・・式（３）

これによれば、式（１）に示される簡易な計算によって信号周波数を算出することができるので、演算量が少なく高速に入力信号の周波数を求めることができる。

本発明の信号周波数算出方法によれば、演算時間を短縮することができる周波数算出方法を提供することができる。

本発明の信号周波数算出方法を適用することができるセンサ装置の概略構成を示す図である。本発明の実施形態に係る周波数算出方法を示すフローチャートである。本発明の実施形態に係る周波数算出方法の手順について示す図である。高速フーリエ変換で得られる離散周波数とその振幅値を示す図である。特許文献１に記載の周波数算出方法を説明する図である。

［第１実施形態］
以下に第１実施形態における周波数算出方法について説明する。

図１は、本発明の信号周波数算出方法を適用することができるセンサ装置１００の概略構成を示す図で、図１（ａ）はセンサ装置１００の動作概要を示す図であり、図１（ｂ）はセンサ装置１００の構成を示すブロック図である。

センサ装置１００は図１（ａ）に示すように、送信波を放射し、検出対象８０からの反射波を検出し検出対象の検知を行う。

センサ装置１００は図１（ｂ）に示すように、アンテナ１０と、送信部２０と、検波部３０と、信号処理部４０と、制御部５０と、を備えている。

アンテナ１０は送信部２０と検波部３０と、に接続されており、送信部２０から出力される送信信号を送信波として空間に放射し、送信波が検出対象８０で反射された反射波を受信する。

検波部３０は、送信部２０からアンテナ１０を介して送信波を送信している間、検出対象８０で反射された反射波を送信信号の一部を参照信号として検波し、信号処理部４０へ出力する。

信号処理部４０は、検波部３０からの入力信号を増幅する機能、増幅した信号を一定時間間隔でサンプリングしてデジタル信号であるサンプルデータへ変換するＡＤ（ＡｎａｌｏｇｔｏＤｉｇｉｔａｌ）変換機能、を有している。

信号処理部４０は、検波部３０から入力された信号を増幅し、ＡＤ変換によってサンプルデータに変換して制御部５０へ出力する。

制御部５０は、信号処理部４０から入力されたサンプルデータを用いて入力信号の周波数を求める。また制御部５０は、求めた周波数及びその変動から検出対象８０の有無や、その動きなどを検出し、検出した結果を表示したり、図示しない外部機器に制御信号を出力したりすることで、検出対象８０の検出結果を出力する。

次に制御部５０がサンプルデータを用いて入力信号の周波数を求める手順について図１から図４を用いて説明する。図２は、本発明の実施形態に係る周波数算出方法を示すフローチャートである。尚、図２に示す周波数算出方法の手順は、図１に示すセンサ装置１００の制御部５０が実行するものとして説明を進める。図３は、本発明の実施形態に係る周波数算出方法の手順について示す図である。図３（ａ）は信号処理部４０から入力したサンプルデータを時間波形として示している。図３（ｂ）は、窓関数（ハニング窓）のレスポンスを示す図である。図３（ｃ）は、図３（ａ）に示すサンプルデータに対して、図３（ｂ）に示す窓関数で演算を行うことにより、求められたデータ列を示す図である。
図４は、高速フーリエ変換（ＦａｓｔＦｏｕｒｉｅｒＴｒａｎｓｆｏｒｍ以下ＦＦＴと記す）で得られる離散周波数とその振幅値の関係を示す図である。

図２に示す周波数算出の手順を示すフローチャートの手順Ｓ１は、サンプルデータ取得ステップである。制御部５０には信号処理部４０から、図３（ａ）に示すサンプルデータが入力されるので、制御部５０は図２に示す手順Ｓ１でこのサンプルデータを取得する。

手順Ｓ２は、窓関数演算ステップである。手順Ｓ２で制御部５０は、手順Ｓ１で取り込んだサンプルデータに、図３（ｂ）に示すあらかじめ記憶している所定の窓関数で演算を行うことで、図３（ｃ）に示すデータ列を求める。本実施形態では、所定の窓関数はハニング窓と称される窓関数を用いるものとして説明を行う。

手順Ｓ３はＦＦＴステップである。手順Ｓ３で制御部５０は、手順Ｓ２で演算した結果であるデータ列に対してＦＦＴ演算を行い、サンプル点数分の離散スペクトル（ｋｎ：ｎは１からサンプル点数までの自然数）を算出する。窓関数で演算したサンプルデータをＦＦＴ演算することにより、図４に示すような離散スペクトルを求めることができる。

手順Ｓ４は抽出ステップである。手順Ｓ３で求めた離散スペクトルのうち、振幅が最大値となるピークスペクトルの周波数（ｋｍａｘ）と振幅値（Ｓｉｍａｘ）と、ピークスペクトルに隣接する両側の離散スペクトルの振幅値（Ｓｉｍａｘ−１，Ｓｉｍａｘ＋１）を抽出する。尚図４では、ピークスペクトルの周波数ｋｍａｘが”０”となるよう作図しており、その振幅値Ｓｉｍａｘ及びｋｍａｘに隣接する両側の離散スペクトルの振幅値Ｓｉｍａｘ−１，Ｓｉｍａｘ＋１の関係が示されている。

手順Ｓ５は補正値を算出するステップである。手順Ｓ４で抽出したピークスペクトルに隣接する両側の離散スペクトルの振幅値（Ｓｉｍａｘ−１，Ｓｉｍａｘ＋１）の差ΔＰと、ピークスペクトルの振幅値（Ｓｉｍａｘ）との比から周波数の補正値を算出する。

手順Ｓ６は入力信号の周波数を算出するステップである。手順Ｓ４で抽出したピークスペクトルの周波数（ｋｍａｘ）と手順Ｓ５で算出した補正値とから入力信号の周波数を求める。

補正値は、前述の式（４）と式（６）からｆ−ｋｍａｘについて整理すると、次の式（８）を得るので、入力信号の周波数ｆは式（１０）で求めることができる。

（数３）
ｆ−ｋｍａｘ＝｛３±√（９＋４ｄ^２）｝／ｄ・・・式（８）
但しｄ＝ｒ−ｓ・・・式（９）
ｆ＝ｋｍａｘ±｛−３＋√（９＋４ｄ^２）｝／｜ｄ｜・・・式（１０）
但し、ｄ＞０のとき＋、ｄ＜０のとき−。

また、式（２）及び式（４）からΔＰを次の式（１１）に置き換えるとｄは式（１２）で表すことができる。

（数４）
ΔＰ＝｜Ｓｉｍａｘ−１｜−｜Ｓｉｍａｘ＋１｜・・・式（１１）
ｄ＝ΔＰ／｜Ｓｉｍａｘ｜・・・式（１２）

次に、本実施形態としたことによる効果について説明する。

本実施形態の周波数算出方法では、サンプルデータに所定の窓関数で演算を行う窓関数演算ステップと、所定の窓関数で演算した結果に対して高速フーリエ変換演算を行い、サンプル点数分の離散スペクトル（ｋｎ：ｎは１からサンプル点数までの自然数）を算出するＦＦＴステップと、離散スペクトルのうち、振幅が最大値となるピークスペクトルの周波数（ｋｍａｘ）と振幅値（Ｓｉｍａｘ）を抽出する抽出ステップと、ピークスペクトルに隣接する両側の離散スペクトルの振幅値（Ｓｉｍａｘ−１，Ｓｉｍａｘ＋１）の差とピークスペクトルの振幅値との比から周波数の補正値を算出するステップと、ピークスペクトルと、補正値から前記入力信号の周波数を算出するステップと、を有し、前記入力信号の周波数を求めるようにした。

これにより、ピークスペクトルに隣接する両側の離散スペクトル（Ｓｉｍａｘ−１，Ｓｉｍａｘ＋１）の振幅の差ΔＰとピークスペクトルの振幅値（Ｓｉｍａｘ）との比から周波数を算出することができるので、対数での演算行うことがない。このため短い演算時間で、入力信号の周波数を算出する方法を提供することができる。

また、本実施形態の周波数算出方法では、窓関数にハニング窓を用いて周波数を求めるようにした。

これにより、離散スペクトルのサイドローブが低く抑えられるので、より精度良く入力信号の周波数を算出することができる。

また、本実施形態の周波数算出方法では、入力信号の周波数は、式（１）から求めるとした。

これにより、式（１）に示される簡易な計算によって信号周波数を算出することができるので、演算量が少なく高速に入力信号の周波数を求めることができる。

以上のように、本発明の実施形態に係る周波数算出方法を具体的に説明したが、本発明は上記の実施形態に限定されるものではなく、要旨を逸脱しない範囲で種々変更して実施することが可能である。例えば次のように変形して実施することができ、これらの実施形態も本発明の技術的範囲に属する。

＜変形例＞
（１）本実施形態において、センサ装置１００に備えられた制御部５０がサンプルデータを用いて入力信号の周波数を求める例を示して説明を行ったが、センサ装置とは別の外部機器がサンプルデータを取得して周波数を算出するようにしても良い。また、周波数を算出する手順をプログラムとして装置にプログラムの一部として実装しても良く、プログラムの供給形態に関わらず実施することが可能である。

（２）窓関数はハニング窓を用いて説明を行ったが、窓関数は必要に応じて変更しても良い。この場合、補正値は変更した窓関数に合わせて適宜変更して求めることが望ましい。

（３）式（１）の平方根内の式を展開した多項式にて演算することにより、平方根の演算が不要となるためさらに演算時間を短縮することが可能である。

１０アンテナ
２０送信部
３０検波部
４０信号処理部
５０制御部
８０検出対象
１００センサ装置

Claims

入力信号が一定時間間隔でサンプリングされてＡＤ変換されたサンプルデータを用いて前記入力信号の周波数を求める信号周波数算出方法であって、
前記サンプルデータに所定の窓関数で演算を行う窓関数演算ステップと、
前記所定の窓関数で演算した結果に対して高速フーリエ変換演算を行い、サンプル点数分の離散スペクトル（ｋｎ：ｎは１からサンプル点数までの自然数）を算出するＦＦＴステップと、
前記離散スペクトルのうち、振幅が最大値となるピークスペクトルの周波数（ｋｍａｘ）と振幅値（Ｓｉｍａｘ）を抽出する抽出ステップと、
前記ピークスペクトルに隣接する両側の離散スペクトルの振幅値（Ｓｉｍａｘ−１，Ｓｉｍａｘ＋１）の差と前記ピークスペクトルの振幅値との比から周波数の補正値を算出するステップと、
前記ピークスペクトルと、前記補正値から前記入力信号の周波数を算出するステップと、を有し、前記入力信号の周波数を求めることを特徴とする信号周波数算出方法。
前記所定の窓関数がハニング窓であることを特徴とする請求項１に記載の信号周波数算出方法。
前記入力信号の周波数（ｆ）は、次の式（１）から求めることを特徴とする請求項２に記載の信号周波数算出方法。
ｆ＝ｋｍａｘ±｛−３＋√（９＋４ｄ^２）｝／｜ｄ｜・・・式（１）
但し、ｄ＞０のとき＋、ｄ＜０のとき−。
また、ｄ＝ΔＰ／｜Ｓｉｍａｘ｜
ΔＰ＝｜Ｓｉｍａｘ−１｜−｜Ｓｉｍａｘ＋１｜