JP3646702B2 - Frequency analyzer - Google Patents

Description

である。ただし、ω＝２πｐ／Ｎτである。これより、
【数２】

である。（３）式は入力信号ｓ（ｎ）について、周波数ｐ／Ｎτにおける周波数成分を表すものである。
【００１０】
ｇ（ｎ）は窓関数と呼ばれるもので、例えば、次のようなハニング窓関数などが一般的に用いられる。
ｇ（ｎ）＝１／２｛１＋ｃｏｓ２π／Ｎ（ｎ−Ｎ／２）｝
ただし、ｎ＝０〜Ｎ−１、それ以外のｎではｇ（ｎ）＝０
【００１１】
図１９は、以上の原理を実現する短時間フーリエ変換を用いた周波数分析装置の構成である。例えば、無線信号が受信されると、適当な中間周波数信号に変換され帯域フィルタ１に入力される。帯域フィルタ１では入力信号について、分析対象とする周波数帯の信号のみを取り出し、Ａ／Ｄ変換器２ではこの取り出した信号について適当な時間間隔τでサンプリングして、デジタル信号に変換する。このデジタル変換した信号列が上述のｓ（ｎ）であり、ｓ（ｎ）は順次メモリ３に記憶される。なお、以下、Ａ／Ｄ変換については、ＩおよびＱチャネルを使用するなどにより、信号が複素数として取得されることを想定するが、これは、信号が実数としてサンプリングされる場合に比べると、サンプリング間隔が２倍で良い他は本質的な違いはない。
【００１２】
重み付け手段４では、メモリ２に記録された信号を順次読み出し、信号列ｓ（ｎ）に窓関数ｇ（ｎ）を掛け合わせ（１）式に表されたような演算を行い、高速フーリエ変換手段５では、重み付けの行われた信号ｓ_ｋ（ｎ）に対し、順次（３）式の処理を行い、各時点での信号に含まれる周波数成分を分析し、その結果は表示装置６に表示する。
【００１３】
以上のように、従来の周波数分析装置は、短時間フーリエ変換を使用して受信される電波信号の周波数分布を求めている。ここで、短時間フーリエ変換の性質を考えると、（３）式の性質から以下の特徴が導き出せる。
【００１４】
短時間フーリエ変換を行った結果、周波数ｆ＝ｐ／Ｎτ（ｐは整数）となり、離散的な値を持つ。その内、最も高い周波数はｐ＝Ｎ−１のときの（Ｎ−１）／Ｎτである。即ち、サンプリング定理により、エイリアジングの発生しない周波数範囲は１／τである。即ち、分析しようとする最大周波数をｆ_ｗとすると、ｆ_ｗ＜１／τならば、周波数の特定が可能であることを意味する。周波数１／τはナイキスト周波数ともいう。
【００１５】
また、（３）式の右辺は、
【数３】

のような関係が成り立つ。ここでｒは整数である。これは、真の周波数ｐ／Ｎτからｒ／τだけずれた周波数に同一の値を持った成分が現れることを意味する。即ち、図２１に示したように、本来の信号の周波数の他に、そこから１／τごとにエイリアジングと呼ばれる周波数成分が発生する。
【００１６】
帯域フィルタ１によって分析する周波数成分の範囲を１／τ以下に制限すればこのようなエイリアジングの発生を防ぐことができる。しかし、上述のように周波数分析装置には分析周波数の広帯域化が望まれるが、これにより、帯域フィルタの通過させる周波数範囲を１／τ（τはサンプリング間隔）より広帯域化すれば、それだけ多くのエイリアジングが生じ、信号の真の周波数と分離することができない。以上が、エイリアジング発生による問題点である。
【００１７】
次に、上述したように、有限の時間範囲で窓関数によって畳込み演算を行うことで原理的にスペクトルが拡散する誤差について説明する。即ち、ｇ（ｎ）のフーリエ変換をＧ（２πｐ／Ｎτ）とすると、

また、（１）式および（２）式から、ｓ（ｎ）のフーリエ変換をＳ（ω）とすると、
【数４】

ここで、便宜上ｍ＝ｎ−ｋＮのように置換すると、
【数５】

従って、Ｇ（ω）が複数の値を持つことから、短時間フーリエ変換で計算される周波数分布は、真の周波数の周囲に、窓関数の性質にしたがって拡散する。この拡散の様子を示したものを図２２に示す。図２２（ａ）は信号の元の周波数分布であり、（ｂ）は窓関数によって拡散した周波数分布を表す。拡散した各成分の強度の和は、（８）式の畳込み演算の性質から元の信号の強度に一致する。また、その拡散の程度は窓関数ｇ（ｎ）の性質によるが、性質の良いハニング窓関数であっても、真の周波数の高低に３つの拡散成分が発生することは免れない。
【００１８】
この拡散成分を抑え真のスペクトルを算出する周波数分析装置の一例として、例えば、特開平６−１６０４４５号公報に示された周波数分析装置がある。図２３に特開平６−１６０４４５号公報の周波数分析装置を示す。図中、１はＡＤ変換器、２はメモリ、３は窓関数重み付け手段、４は高速フーリエ変換手段、５は高速フーリエ変換器４の出力側に設けられたスペクトル補間判定手段、６は表示器、７は補間判定手段５の判定結果によりスペクトル補間を実行するスペクトル補間手段、８はスペクトル補間手段７によって算出したスペクトルの周波数、振幅、位相データを時間軸データに逆変換する逆フーリエ変換器、９は逆フーリエ変換器８で逆フーリエ変換した時間軸データを累積する累積器、１１は累積器９によって累積された累積結果をメモリ２から読み出される時間データから減算し、この減算結果を再び窓関数重み付け手段３に入力する減算器、１２はスペクトル補間手段７によって算出した周波数領域データを累積する累積器、１３はこの累積器１２に累積した周波数領域データを減算器１１で減算した残差データを再度高速フーリエ変換して得られた周波数領域データＳ（ｔ）に加算する加算器である。
【００１９】
この周波数分析装置では、フーリエ変換して得られた離散周波数スペクトルの中で３本の周波数スペクトルが互いに隣接して存在し、しかもその上位と下側のスペクトルの位相がほぼ同じであったら、この３本のスペクトルは真のスペクトルが分散されているものと判定し、補間処理によりその真のスペクトルの周波数と振幅、位相を求める。
【００２０】
この真のスペクトルデータを累積器１２に累積するとともに、逆フーリエ変換して時間軸データに逆変換し、この時間軸データを累積器９で累積し、その累積データを入力データ列から減算するから、補間処理により真のスペクトル線を算出したあとでは残差エネルギーが低下し、線スペクトルは広がりを持つことはなく、線スペクトルとして表示される。この結果、スペクトル分解能が向上する。また、補間処理により振幅、位相も真値近くに近似して算出するから振幅値および位相も正確に測定することができる。
【００２１】
しかし、この場合、フィードバックループにより累積したデータを繰り返し減算して真の周波数成分を取り出していく構成であるため、値が収束するまでに何回もサンプリングを行う必要があり、時間がかかり、また、フィードバックループを無くした場合は、微少信号を抽出できない。
【００２２】
また、上記の例では真の周波数の周辺の周波数を正しく３つ選択できることが必要だが、真の周波数が、離散周波数値のほぼ中間に存在した場合、例えば、ｎ／Ｎτと（ｎ＋１）／Ｎτの中間の場合、真の周波数に最も近い周波数を（ｎ＋１）／Ｎτと誤判定しやすく、真の周波数の周辺の３つの周波数を間違えてｎ／Ｎτ、（ｎ＋１）／Ｎτ、（ｎ＋２）／Ｎτのように選んでしまうことが起こりやすいため、実際の周波数算出は不安定になりやすい。
【００２３】
さらに、この場合、真の周波数の算出方法は特開平６−１６０４４５号公報の（１４）式などのように引き算が主体であり、ｓｉｎやｃｏｓが０の付近では計算結果が不安定であり、周波数算出に誤差が出やすい。
【００２４】
【発明が解決しようとする課題】
本発明は、以上のような問題点を解決するためになされたもので、短時間フーリエ変換の際に、受信信号を広帯域化するに伴う計算量の増大を抑え、周波数分布を効率よく正確に算出することを目的とする。
【００２５】
【課題を解決するための手段】
本発明に係わる周波数分析装置は、周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、上記帯域制限されたアナログ入力信号を第一の時間間隔で順次デジタルデータ列に変換するＡＤ変換器と、上記ＡＤ変換器によって変換された入力データ列を順次記憶するメモリと、上記メモリに記憶された入力データ列から上記第一の時間間隔よりも長い第二の時間間隔を持った第一の分配データ列、およびこの第一の分配データ列から上記第一の時間間隔ずつ遅れた第二の分配データ列とを抽出するデータ分配手段と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、上記第一または第二の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、信号が検出された場合、上記第一および第二の周波数領域データの位相を比較して真の周波数を算出するエイリアジング判定手段とを備えるものである。
【００２６】
また、周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、上記帯域制限されたアナログ入力信号を２つに分配する信号分配器と、上記分配された一方のアナログ入力信号を一定の時間間隔でデジタル信号に変換し、第一の分配データ列を生成する第一のＡ／Ｄ変換器と、上記分配された他方のアナログ入力信号を、上記第一のＡ／Ｄ変換器のデジタル変換より上記一定の時間間隔より短い時間間隔だけ遅れて、一定の時間間隔でデジタル信号に変換し、第二の分配データ列を生成する第二のＡ／Ｄ変換器と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、上記第一または第二の周波数領域データを各周波数ごとに信号が存在するかを検出する信号検出手段と、信号が検出された場合、上記第一および第二の周波数領域データを比較することで、エイリアジングの影響を排除し、真の周波数を算出するエイリアジング判定手段とを備えるものである。
【００２７】
また、周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、上記帯域制限されたアナログ入力信号を第一の時間間隔で順次デジタルデータ列に変換するＡＤ変換器と、上記ＡＤ変換器によって変換された入力データ列を順次記憶するメモリと、上記メモリに記憶された入力データ列から上記第一の時間間隔よりも長い第二の時間間隔を持った第一の分配データ列と、この第一の分配データ列から上記第一の時間間隔ずつ遅れた第二の分配データ列と、この第一の分配データ列から上記第一および第二の時間間隔とは異なる第三の時間間隔ずつ遅れた第三の分配データ列とを抽出するデータ分配手段と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、上記第三の分配データ列に所定の窓関数を乗算し、重み付けを行う第三の重み付け手段と、この重み付けされた第三の分配データ列をフーリエ変換し、第三の周波数領域データに変換する第三の高速フーリエ変換手段と、上記第一または第二または第三の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、信号が検出された場合、上記第一、第二、第三の周波数領域データのいずれか２つの単位時間当たりの位相差を比較して、信号が検出された周波数における信号成分が単一か否かを判定し、信号成分が単一の場合のみ真の周波数を算出するエイリアジング判定手段とを備えるものである。
【００２８】
また、周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、上記帯域制限されたアナログ入力信号を第一の時間間隔で順次デジタルデータ列に変換するＡＤ変換器と、上記ＡＤ変換器によって変換された入力データ列を順次記憶するメモリと、上記メモリに記憶された入力データ列から上記第一の時間間隔よりも長い第二の時間間隔を持った第一の分配データ列と、この第一の分配データ列から上記第一の時間間隔ずつ遅れた第二の分配データ列と、この第一の分配データ列から上記第一および第二の時間間隔とは異なる第三の時間間隔ずつ遅れた第三の分配データ列とを抽出するデータ分配手段と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、上記第三の分配データ列に所定の窓関数を乗算し、重み付けを行う第三の重み付け手段と、この重み付けされた第三の分配データ列をフーリエ変換し、第三の周波数領域データに変換する第三の高速フーリエ変換手段と、上記第一または第二または第三の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、信号が検出された場合、上記第一、第二、第三の各周波数における周波数領域データを要素とする行列に、当該信号が検出された周波数から上記第二の時間間隔の逆数の整数倍ずつずれた複数の周波数を要素とする行列の逆行列を乗算することで真の周波数を算出するエイリアジング判定手段とを備えるものである。
【００２９】
また、周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、上記帯域制限されたアナログ入力信号を第一の時間間隔であらかじめ決められた一定回数だけ順次デジタルデータに変換するＡＤ変換器と、上記ＡＤ変換器によって変換された入力データを順次記憶するメモリと、上記メモリに記憶された入力データから上記第一の時間間隔を持ったデータを抽出してなる第一の分配データ列、および上記第一の時間間隔に上記デジタル変換した一定の回数を掛けた数の二分の一より短い第二の時間間隔ずつ上記第一の分配データ列より遅れた第二の分配データ列とを抽出するデータ分配手段と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、上記第一または第二の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、信号が検出された場合、上記第一および第二の周波数領域データの位相を比較して、この位相が同じの場合、真の周波数を算出し、算出した周波数において当該信号強度を足し合せていく処理を行う畳込み効果判定手段とを備えるものである。
【００３０】
また、周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、上記帯域制限されたアナログ入力信号を第一の時間間隔であらかじめ定められた一定回数だけ順次デジタルデータ列に変換するＡＤ変換器と、上記ＡＤ変換器によって変換された入力データ列を順次記憶するメモリと、上記メモリに記憶された入力データから上記第一の時間間隔を持ったデータを抽出してなる第一の分配データ列、および上記第一の時間間隔に上記デジタル変換した一定の回数を掛けた数の二分の一より短い第二の時間間隔ずつ上記第一の分配データ列より遅れた第二の分配データ列と、上記第一の時間間隔に上記デジタル変換した一定の回数を掛けた数の二分の一より短かく、かつ上記第一および第二の時間間隔とは異なる第三の時間間隔ずつ上記第一の分配データ列より遅れた第三の分配データ列とを抽出するデータ分配手段と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、上記第三の分配データ列に所定の窓関数を乗算し、重み付けを行う第三の重み付け手段と、この重み付けされた第三の分配データ列をフーリエ変換し、第三の周波数領域データに変換する第三の高速フーリエ変換手段と、上記第一または第二または第三の周波数領域データいずれかの各周波数ごとに信号が存在するかを検出する信号検出手段と、信号が検出された場合、上記第一、第二、第三の周波数領域データのいずれか２つの単位時間当たりの位相差を比較して、信号が検出された周波数における信号成分が単一か否かを判定し、信号成分が単一の場合のみ上記第一および第二の周波数領域データの位相を比較して、この位相が同じの場合、真の周波数を算出し、算出した周波数において当該信号強度を足し合せていく処理を行う畳込み効果判定手段とを備えるものである。
【００３１】
【発明の実施の形態】
実施の形態１．
本発明の実施の形態１について説明する。本実施の形態１に係わる周波数分析装置は、Ａ／Ｄ変換された信号から微少タイミングずらしてサンプルして得られた２系列の信号列をそれぞれ短時間フーリエ変換処理して局所的な周波数分布を求めるものである。また、その際生じるエイリアジングを微少タイミングずらしたことにより生じる信号列間の位相差を検出することで真の周波数を特定するものである。
【００３２】
以下、詳細について説明する。図１は本実施の形態１に係わる周波数分析装置の構成を表すブロック図である。図１において、図１９と同じ構成要素には同じ符号を付す。７はメモリ３内に記録した信号のデータからサンプリングタイミングのわずかに異なる２つのデータ列を抽出、分配するデータ分配手段、８はフーリエ変換結果から各周波数に関して信号の存在を判定する信号検出手段、９は信号検出手段８で判定された信号の存在する周波数に関して、２つの短時間フーリエ変換の結果を比較してエイリアジングの影響を排除し、正しい周波数を算出するエイリアジング判定手段である。
【００３３】
図２は本実施の形態１に係わる周波数分析装置のサンプリング方法を表す図である。図２に示した通り、時間間隔Δτでサンプリングしたデータの中から、適当な時間間隔τの信号列ｓ（ｎ）と、当該ｓ（ｎ）の各信号からΔτだけ遅れた信号列ｓ_Δ（ｎ）を考える。信号の周波数をｆ（＝ｐ／Ｎτ，ｐは整数）とすると、信号列ｓ（ｎ）とｓ_Δ（ｎ）とは、位相遅れ２πｆΔτが発生するため、
【数６】

の関係が成立する。ここで、係数ｅｘｐ（２πｊｆΔτ）は定数であるため、重み付け、フーリエ変換とあった短時間フーリエ変換の各プロセスによって変化することはなく、信号列ｓ_Δ（ｎ）は、ｓ（ｎ）に対して一様な周波数ｆによる位相差を持つ。このため、ｓ_Δ（ｎ）のフーリエ変換をＳＴＦＴ_Δ（２πｐ／Ｎτ）とすると、
【数７】

（１０）式より、ＳＴＦＴ（２πｐ／Ｎτ）≠０となる周波数ｆに対して、元のデータ列ｓ（ｎ）と、微少時間Δτだけずらしてサンプリングを行ったデータ列ｓ_Δ（ｎ）のそれぞれの短時間フーリエ変換を比較すると、係数ｅｘｐ（２πｊｆΔτ）の違いとなる。従って、（１０）式を変形すると、
【数８】

であり、（１１）式右辺の各ＦＦＴ成分については、それぞれ間隔τのデータを使用しているため、ｐに対する周波数ｆは、
ｆ＝ｐ／Ｎτ，（ｐ＋Ｎ）／Ｎτ，（ｐ＋２Ｎ）／Ｎτ…
のいずれかである。従って、ｒを整数としてｆ＝（ｐ＋ｒＮ）／Ｎτのように表現でき、（１１）式は、
【数９】

である。その際、−２πｊｒＮΔτ／Ｎτの項によって位相が一回転しない場合は、ｒを一意に求めることができる。即ち、
【数１０】

を満たす場合である。（１２）式を変形すると、
１／Δτ＞ｒ_Ｍ／τ （１４）の関係が得られる。１／Δτは本実施の形態１の周波数分析装置の特定可能な周波数の範囲である。ｒ_Ｍが余りに大きい場合には（１４）式を満たすことができなくなるため、帯域フィルタ１で適当な周波数帯域を選択する。
【００３４】
図３に（１３）式の関係を表すグラフを示す。図３は縦軸を振幅、横軸を周波数とした周波数分析結果の図である。周波数の値は離散値でその最小単位は１／Ｎτである。また、１／Δτがこの周波数分析装置で分析できる最大の周波数であり、１／Δτより上の周波数が入力されないように帯域フィルタ１で帯域制限する。一方、サンプリング間隔τの信号列ｓ（ｎ）を短時間フーリエ変換するのであるから、エイリアジングの発生しない周波数範囲は１／τである関係上、真の信号周波数から１／τずつずれたところにエイリアジングが発生する。
【００３５】
図では例えば、周波数（ｐ＋２Ｎ）／Ｎτに真の信号成分がある場合である。信号列ｓ（ｎ）およびｓ_Δ（ｎ）から見るとエイリアジングの発生しない周波数範囲は１／τであるのに対し、帯域フィルタ１から入力される受信信号はその周波数範囲を超えているので、（１１）式の右辺を実測によって求め、そこから左辺のｐおよびｒを求める。
【００３６】
以上が、本実施の形態１に係わる周波数分析装置の原理である。図１の周波数分析装置は以上の原理に基いて短時間フーリエ変換を行うものである。受信した電波を適当な中間周波数に変換した受信信号が、帯域フィルタ１に入力され、帯域制限される。Ａ／Ｄ変換器では、選択された信号をナイキスト周波数の逆数以下である適当な短い時間間隔Δτでサンプリングし、メモリ３に順次蓄積する。データ分配手段７では、メモリに蓄積した信号からそれぞれ適当な間隔τの２系列の信号列ｓ（ｎ）およびｓ_Δ（ｎ）を選択し、それぞれ重み付け手段４ａおよび４ｂに分配する。重み付け手段４ａおよび４ｂにおいて窓関数による重み付けがなされた後、フーリエ変換手段５ａおよび５ｂに入力される。フーリエ変換手段５ａおよび５ｂでは、重み付けされたｓ（ｎ）およびｓ_Δ（ｎ）のフーリエ変換を行い、（２）式のように、ＳＴＦＴ（ω）およびＳＴＦＴ_Δ（ω）を求める。
【００３７】
信号検出手段８ではフーリエ変換手段５ａの出力ＳＴＦＴ（ω）の結果を調べ、信号の存在する周波数を判定し、それらの周波数をエイリアジング判定手段９に通知する。エイリアジング判定手段９では、ＳＴＦＴ（ω）およびＳＴＦＴ_Δ（ω）を比較し、（１０）式の左辺ｅｘｐ（２πｊｆΔτ）を算出する。
【００３８】
信号検出手段８およびエイリアジング判定手段９で行う処理を図４に示したフローチャートにしたがって説明する。まず、信号検出手段８で離散値の各周波数ｆ＝（ｐ＋ｒＮ）／Ｎτをｐ＝０（ステップ１００）から順番に信号が存在するかの判定が行われる（ステップ１０１）。そして、信号が未検出ならば、そのときのｐ（今の場合、ｐ＝０）に対して、全てのｒについて受信信号の周波数成分が０である（ステップ１０２）。即ち、ｐ＝０で信号が未検出であれば、そのエイリアジングの可能性がある周波数ｆ＝Ｎ／Ｎτ，２Ｎ／Ｎτ，３Ｎ／Ｎτ，…における振幅＝０ということである。
【００３９】
一方、ｐ＝０で信号が検出された場合は、ＳＴＦＴ（ω）およびＳＴＦＴ_Δ（ω）から（１１）式を満足するｒを算出する（ステップ１０３）。そして、算出された値以外のｒについての受信信号の周波数成分は０である（ステップ１０４）。例えば、ｒ＝３と算出された場合は、それ以外の周波数ｆ＝（３＋Ｎ）／Ｎτ，（３＋２Ｎ）／Ｎτ，（３＋４Ｎ）／Ｎτ，…における振幅＝０ということである。
【００４０】
そして、ｐ＝０のとき、信号が検出、または未検出いずれかの場合の処理が終了すると、次の周波数、即ちｐ＝１（ステップ１０５）の場合について、同様の処理を行う。その際、ｐ≧Ｎであれば、演算する部分が重なってしまうため、ｐ≧Ｎであれば分析を終了し、ｐ＜Ｎの場合のみステップ１０１に戻る（ステップ１０６）。
【００４１】
このように、ｐ＝０からＮまで、即ち離散値の周波数ｆを０から１／τまで信号の有無を検出するだけで、実際は図３の０からｒ_Ｍ／τまでの信号帯域を分析することになるので、全ての周波数を分析する必要がなく、計算量が大幅に軽減することができる。
【００４２】
例えば、サンプル数Ｎ＝１２８の場合に、周波数帯域を１６倍に広げることを考えると、
（１）従来の周波数分析装置で周波数帯域を１６倍に広げる場合
フーリエ変換を行うデータ数が１６倍であるため、計算量は（１２８×１６）ｌｏｇ_２（１２８×１６）／１６ｌｏｇ_２１２８倍、即ち２５倍になる。
（２）本実施の形態１の周波数分析装置で周波数帯域を１６倍に広げる場合
フーリエ変換を行うデータ数自体は変化がない。そして、２つの信号列ｓ（ｎ）とｓ_Δ（ｎ）をフーリエ変換するため計算を２回行うためほぼ２倍である。ただし、（１１）式の計算も行う必要があるが、これは高々Ｎ回（今の場合１２８回）であり、フーリエ変換の計算量と比べれば１／１０程度であるため、計算量の増加にはほとんど影響しない。
このように、本実施の形態１の周波数分析装置を用いれば、約１２倍の計算負荷が軽減できる。
【００４３】
以上のように、本実施の形態１では、時間間隔τでサンプリングした信号列ｓ（ｎ）と、この信号列ｓ（ｎ）から微少時間Δτだけずらして時間間隔τでサンプリングした信号列ｓ_Δ（ｎ）とをそれぞれ短時間フーリエ変換し、比較することでΔτだけずらしたことにより生じる位相差ｅｘｐ（２πｊｆΔτ）を算出する。この位相差によって、受信信号を広帯域化に伴って生じるエイリアジングを利用して（１１）式によりｒを求めることで、帯域フィルタ１から出力される周波数帯域全てを走査しなくても受信信号の存在する周波数を算出することができるため、広帯域化しても計算量の増大を抑え、エイリアジングの効果を排除し、正確な周波数分析結果を得ることができる。
【００４４】
実施の形態２．
上記実施の形態１では、２つの信号列ｓ（ｎ）、ｓ_Δ（ｎ）を得るために、一律にΔτでサンプリングしたデータ群の中から適当な時間間隔τとなるデータを抽出してデータを分配する例を示したが、例えば、図５のように、タイミングのΔτだけ異なるクロックでサンプリングを行う２つのＡ／Ｄ変換器を使用して実現してもよい。
【００４５】
図６に本実施の形態２に係わる周波数分析装置の装置構成を表すブロック図を示す。図１と同じ構成要素には同じ符号を付す。１０は信号分配器である。本実施の形態２の周波数分析装置は図のように２つのＡ／Ｄ変換器２ａ、２ｂを使用するが、無駄な信号をサンプリングすることがなくなるため、メモリ３ａ、３ｂの容量を小さくて済む。
【００４６】
以下、動作について説明する。帯域フィルタ１は入力された受信信号を帯域制限して分析対象の周波数信号のみ選択する。信号分配器１０では、信号を２つのＡ／Ｄ変換器２ａ、２ｂに分配する。各Ａ／Ｄ変換器２ａ、２ｂではそれぞれ時間間隔τでサンプリングを行うが、例えば、Ａ／Ｄ変換器２ｂの方がＡ／Ｄ変換器２ａよりもΔτだけ遅れたタイミングでサンプリングを行う。これによって、実施の形態１で示したような２つの信号列ｓ（ｎ）およびｓ_Δ（ｎ）が生成され、順次メモリ３ａ、３ｂに記憶される。以下の動作は上記実施の形態１の周波数分析装置の動作と同様である。それぞれの信号列の短時間フーリエ変換ＳＴＦＴ（ω）、ＳＴＦＴ_Δ（ω）を比較し、（１１）式からｅｘｐ（２πｊｆΔτ）を算出する。
【００４７】
以上のように、本実施の形態２では、帯域制限した信号を予め分配してＡ／Ｄ変換器に入力するため、無駄な信号をサンプリングすることがなく、メモリの容量が小さくて済む。
【００４８】
実施の形態３．
上記（１１）式はそれぞれのｐについてＳＴＦＴ（ω）に含まれる信号が単一である場合のみ成立する。従って、上記実施の形態１および２では、受信信号の周波数分布が単一もしくは比較的疎な場合のみ分析が可能であった。本実施の形態３では、上記実施の形態１または２の周波数分析装置を用いた場合に互いにエイリアジング成分となるような複数の周波数成分が受信信号に含れていた場合、信号が複数のエイリアジング成分を含むか含まないかの判定を行い、含まないと判定された場合にのみ分析を行う周波数分析装置について説明する。
【００４９】
図７は本実施の形態３の周波数分析装置が対象とする信号の周波数スペクトル図である。点線は実際に測定されるスペクトル、実線は信号の真の周波数成分を表す。図７のように上記実施の形態１および２の周波数分析装置を用いた場合に、互いにエイリアジング成分となるような独立した信号成分ｗ_１、ｗ_３が受信信号に含まれていた場合、その周波数をｆ_１、ｆ_３とすると、
【数１１】

また、
【数１２】

これらの式をそれぞれフーリエ変換すると、
【数１３】

また、
【数１４】

である。
【００５０】
実際の信号から求まるのはＳＴＦＴ_Δ（２πｐ／Ｎτ）およびＳＴＦＴ（２πｐ／Ｎτ）の値のみである。従って、上記実施の形態１または２のようにＳＴＦＴ_Δ（２πｐ／Ｎτ）／ＳＴＦＴ（２πｐ／Ｎτ）を求めても、ｗ_１、ｗ_３が混ざっているため正しいｆ_１、ｆ_３の値を求めることができない。
【００５１】
そのために、本実施の形態３では、図８に示したように周波数分析装置を、時間間隔τを持った一つの信号系列ｓ（ｎ）と、そのｓ（ｎ）に対してそれぞれ２つの異なる微少な時間差Δτ_１、Δτ_２を持った信号系列ｓ_Δ１（ｎ）、ｓ_Δ２（ｎ）を生成する構成とし、それぞれ短時間フーリエ変換を行う。図において９ａは信号が単一か否かを判定しながら真の周波数を算出するエイリアジング判定手段である。
【００５２】
ここで、エイリアジングの周波数をｆ_１＝ｐ／Ｎτ、ｆ_２＝（ｐ＋Ｎ）／Ｎτ、ｆ_３＝（ｐ＋２Ｎ）／Ｎτ、…ｆ_ｒ＋１＝（ｐ＋ｒＮ）／Ｎτ、また各エイリアジング周波数における各信号強度をｗ_１、ｗ_２、ｗ_３、…ｗ_ｒとする。すると、信号が複数含まれている場合を一般化して表現すると、
【数１５】

【数１６】

である。ｓ_Δ２（ｎ）も（１６）式のΔτ_１を、それぞれΔτ_２と置き換えたものと同様である。
【００５３】
ここで、例えば、周波数ｆ_３に真の信号成分が１つのみあるとすると、ｗ_３のみが０以外の値を持ち、かつｗ_３以外のｗ_１、ｗ_２、ｗ_４…ｗ_ｒは全て０である。
【００５４】
従って、その場合にｓ（ｎ）、ｓ_Δ１（ｎ）、ｓ_Δ２（ｎ）をフーリエ変換すると、それぞれ、以下のようになる。
【数１７】

【００５５】
逆に、信号成分が単一でなく、例えば、周波数ｆ_３の他に周波数ｆ_１に信号成分ｗ_１が存在していた場合は、そのフーリエ変換は（１７）、（１８）、（１９）式の代わりに、
【数１８】

【００５６】
ここで、各信号列間の単位時間当たりの位相差として以下の量を考える。
【数１９】

ただし、（２０）式において∠は位相を表す。
例えば、信号成分が周波数ｆ_３に一つの場合、（１７）〜（１９）式の右辺を（２１）式に代入すると、全て２πｆ_３が得られるが、信号成分が周波数ｆ_３およびｆ_１にあるような場合、（１７ａ）〜（１９ａ）式の右辺を（２０）式に代入しても一定の値は得られない。そのため、この（２０）式を用いて信号成分が複数あるか否かを判定する。
【００５７】
図８は本実施の形態３に係わる周波数分析装置の構成を表すブロック図である。図１と同じ構成要素には同じ符号を付す。１つの信号系列ｓ（ｎ）およびそのｓ（ｎ）から時間Δτ_１およびΔτ_２だけずれた２つの信号系列ｓ_Δ１（ｎ）およびｓ_Δ２（ｎ）の計３つの信号系列をサンプリングし、それぞれ短時間フーリエ変換する。そこから、信号検出手段８ｂおよびエイリアジング判定手段９ｂによって、上記（２０）式のような単位時間当たりの位相差を求め、信号が一つか否かを判定しながら真の周波数を求める。
【００５８】
以下、図９に示したフローチャートに従って、信号検出手段８ｂおよびエイリアジング判定手段９ｂの動作を説明する。まず、信号検出手段８ｂで離散値の各周波数ｆ＝（ｐ＋ｒＮ）／Ｎτをｐ＝０（ステップ２００）から順番に信号が存在するかの判定が行われる（ステップ２０１）。そして、信号が未検出ならば、そのときのｐ（今の場合、ｐ＝０）に対して、全てのｒについて受信信号の周波数成分が０である（ステップ２０２）。即ち、上記実施の形態１の場合と同様に、ｐ＝０で信号が未検出であれば、そのエイリアジングの可能性がある周波数ｆ＝Ｎ／Ｎτ，２Ｎ／Ｎτ，３Ｎ／Ｎτ，…についての周波数成分も０ということである。
【００５９】
一方、ｐ＝０で信号が検出された場合は、ＳＴＦＴ（ω）、ＳＴＦＴ_Δ１（ω）、ＳＴＦＴ_Δ２（ω）、ＳＴＦＴ_Δ３（ω）から（２０）式を計算し（ステップ２０４）、全て同じ値を算出するか判定する（ステップ２０５）。もし、（２１）式の計算結果が一致しないならば、上述したように信号成分が複数存在することになるので、当該周波数の成分は不定である（ステップ２０３）。一方、（２０）式の結果が同じ値ならば信号は１つであり、周波数の分析が可能であるので実施の形態１と同様に（１１）式を満足するｒ_Ｔを算出する（ステップ２０６）。そして、算出したｒ_Ｔ以外のｒに対し、周波数（ｐ＋ｒＮ）／Ｎτの成分が０である（ステップ２０７）。
【００６０】
以上、いずれかの処理が終わったら周波数（ｐ＋ｒＮ）／Ｎτのｐをｐ＝ｐ＋１として、次の周波数成分について同様の処理を行う（ステップ２０８）。その際、ｐ≧Ｎであれば、演算する部分が重なってしまうため、ｐ≧Ｎであれば分析を終了し、ｐ＜Ｎの場合のみステップ２０１に戻る（ステップ２０９）。
【００６１】
このように、ｐ＝０からＮ−１まで、即ち離散値の周波数ｆを０から１／τまで信号の有無を検出するだけで、全ての周波数を分析する必要がなく、計算量が大幅に軽減することができる。
【００６２】
また、信号がいくつあるか判定しながら周波数分析を行うため、信号成分が複数ある場合でも、誤った周波数を算出しない。
【００６３】
実施の形態４．
本実施の形態４では、それぞれ異なる時間だけずれた２つ以上の信号系列を短時間フーリエ変換して、その離散値の各周波数に含まれる周波数成分が１つか否かを判定し、複数の信号が含まれていてもエイリアジングの影響を除くことができる周波数分析装置について説明する。
【００６４】
ここでは上記実施の形態３の図７のように互いのエイリアジング成分となるような周波数に独立した複数の信号がある場合に、３つ信号系列を短時間フーリエ変換する例を示す。
【００６５】
図１０は本実施の形態４の周波数分析装置の構成を表すブロック図である。図において、図８と同じ構成要素には同じ符号を付す。９ｂは本実施の形態３に係わるエイリアジング判定手段である。前記実施の形態３と同様に、信号系列ｓ（ｎ）に対し複数の時間差Δτ_１、Δτ_２、Δτ_３を持った信号系列ｓ_Δ１（ｎ）、ｓ_Δ２（ｎ）、ｓ_Δ３（ｎ）をサンプリングして生成し、短時間フーリエ変換を行う。ただし、ここで、Δτ_１＜Δτ_２＜Δτ_３であって、図７のｒ／τ＜１／Δτ_１であるとする。これら４つの信号系列の短時間フーリエ変換後のスペクトルＳＴＦＴ（２πｐ／Ｎτ）、ＳＴＦＴ_Δ１（２πｐ／Ｎτ）、ＳＴＦＴ_Δ２（２πｐ／Ｎτ）、ＳＴＦＴ_Δ３（２πｐ／Ｎτ）それぞれについて（１１）式のような位相差を考える。
【００６６】
ここで、エイリアジングの性質として、真の周波数から１／τごとにエイリアジング成分が生じることから、エイリアジングの生じる可能性のある周波数は、
【数２０】

のように表現でき、この各周波数における振幅をｗ_１、ｗ_２、…ｗ_ｒ＋１とすれば、４つの信号系列は以下のような行列にて表現できる。
【数２１】

として、
Ｘ＝Ｓ・Ａ （２２）
である。
【００６７】
ここで、行列Ｘの各要素については４つの信号系列を実際に短時間フーリエ変換することで、その値が得られ、また、行列Ｓの各要素についても、上述のようにエイリアジング成分の性質から各周波数が定められるため、（２１）各式中のｐを固定すればその値が定まる。従って、Ｓの一般逆行列を用いて、
Ａ＝Ｓ^−１Ｘ＝（Ｓ^＊Ｓ）^−１・Ｓ^＊・Ｘ （２２ａ）
を計算することで、行列Ａの各要素を求めることができる。
【００６８】
エイリアジング判定手段９ｂでは、以上のような原理に基づいて周波数０から１／τまでのｐについて信号の有無を判定し、信号がある場合にはそのｐについて上記（２３ａ）式によって行列Ａの各要素を求める。
【００６９】
この過程をフローチャートにしたものが図１１である。以下、図に従って説明する。まず、（２１）式のｐ＝０から分析を開始する（ステップ３００）。そして、信号の有無を確認するが（ステップ３０１）、例えば、信号がｐ＝０のとき発見されなければ、（２１）式の各周波数において振幅＝０、即ち、信号はない（ステップ３０２）として、ｐを１だけ増やして（ステップ３０４）再びステップ３０１に戻って、信号の有無を確認する。一方、信号が発見された場合には、上述したように、（２３ａ）式を用いて３つの信号系列から行列Ａの各要素、即ち、各周波数における振幅を計算することで、信号成分が複数あっても正しく算出することができる。そして、算出できたらｐを１だけ増やして（ステップ３０４）再び信号の分析を行う。そして、ｐ＝Ｎ−１まで分析したら、周波数０からｒ_Ｍ／τまで分析が終わったことになるので、分析を終了する。
【００７０】
このように、周波数を０から１／τまで分析するだけで、周波数０からｒ_Ｍ／τまで分析を行うことになるため、信号を広帯域化しても効率のよい分析を行うことができる。
【００７１】
また、エイリアジング成分と複数の真の周波数成分とが重なりあっている場合でも、互いに異なる微少時間だけずれた信号系列を３つ生成し、位相差を比較することで、複数の周波数成分を分離して正しい周波数を算出することができるため、正確な分析を行うことができる。
【００７２】
以上の説明では、全ての信号成分を（２２）式を用いて算出する例について示したが、（２２）式は（１１）式に比べて計算量が多くなるため、実施の形態３に示したように、信号を単一か否かを判定するステップを挿入し、信号が単一の場合には（１１）式を用いて周波数を算出するようにしてもよい。
【００７３】
図１２に信号成分が単一か否かを判定しながら、信号成分を算出する場合のフローチャートを示す。図１１と異なる点は、ステップ４０４にて実施の形態３で行ったように（２０）式を用いて信号が単一か否かを判定する点である。信号成分が単一ならば実施の形態１と同様に算出し、複数ある場合は本実施の形態４で上述したように（２２）式を用いて全ての信号成分を算出する。
【００７４】
なお、本実施の形態４では、生成する信号系列を３つとしたため、周波数成分が２つ以内の場合に対応できるが、重み付け手段と、高速フーリエ変換手段を更に増やせば、さらに多くの周波数成分がエイリアジング成分と重なっても正確に分離することができる。
【００７５】
以上のように、本実施の形態４では、生成した３つの信号系列を短時間フーリエ変換し、エイリアジング成分の性質を利用して行列の要素を決定し、位相差を比較して信号成分を算出するため、複数の信号が同じ周波数に重なっている場合でも、互いに分離して算出することができる。
【００７６】
実施の形態５．
以上、実施の形態１〜４では、広帯域化に伴って発生するエイリアジング成分を排除しつつ、計算量の増大を抑え効率のよい周波数分析装置について説明したが、本実施の形態５では、窓関数の効果による周波数成分の拡散による誤差を排除しつつ、短時間のサンプリングで正確な周波数分析が可能な周波数分析装置について説明する。
【００７７】
図１３は本実施の形態５に係わる周波数分析装置の装置構成を表すブロック図である。図において図１と同じ構成要素には同じ符号を付す。１１は畳込み効果判定手段である。
【００７８】
以下、動作について説明する。帯域フィルタ１で信号の周波数帯を適当な帯域に制限し、Ａ／Ｄ変換器２で適当な時間間隔τでサンプリングする。サンプリングしたデジタル信号はメモリ３に記憶し、データ分配手段７によって、データを重み付け手段４ａ、４ｂに分配し、一方の信号系列ｓ（ｎ）に対して、他方ｓ（ｎ）に対してΔτだけずれたｓ_Δ（ｎ）が生成するようにする。それぞれの重み付け手段４ａ、４ｂでは窓関数による重み付けを行い、高速フーリエ変換手段５ａ、５ｂにおいてフーリエ変換を行う。信号検出手段８ではフーリエ変換後の各周波数をチェックして信号の有無を確認し、そして、畳込み効果判定手段１１では、このフーリエ変換後の周波数スペクトル同士の位相を比較することで、真の周波数を特定する。
【００７９】
このように構成することで、短時間のサンプリングで畳込みによる周波数成分の拡散を抑えることができ、連続波のみならず、断続波の場合でも、正確な周波数スペクトルを得る。
【００８０】
ここで、窓関数による畳込みの効果には以下に示すような特徴がある。即ち、（Ａ）畳込みにより拡散する各周波数の信号強度の和は元の信号強度にほぼ一致する。（Ｂ）１つの信号成分（周波数ｆ）から畳込みにより拡散した各周波数成分について、（１１）式は同じ値を持つ。（Ｃ）ハニング窓関数を用いることで、１つの信号の畳込みによる拡散の影響は、真の周波数周辺の３つの周波数値に限定することができる。
【００８１】
（Ａ）は、図２３および、（８）式の説明で示したように、畳込み演算の性質によるものである。この場合、周波数は離散値を取るため、近似的に一致する。また、（Ｂ）については、ＳＴＦＴ、ＳＴＦＴ_Δとも同じ元の周波数成分を窓関数ｇ（ｎ）により拡散したものであるから相互の位相の関係は同じになり、即ち、位相を表す（１０）式が異なる値を取ることはないからである。
【００８２】
従って、これらから畳込みによる周波数分布の拡散を排除して正しい周波数分布を求めるには、検出された周波数の信号成分が畳込みによる拡散の結果であることを（１０）式により判定し、畳込みによる拡散によると判定された場合には、拡散した各信号成分の信号強度を足し合せ、（１０）式のより求められた周波数に足し合せた信号強度の信号成分があるとすればよい。
【００８３】
なお、（Ｃ）で示したように、周波数の拡散はハニング窓関数を用いることで、真の周波数周辺の３つの周波数値に限定することができる。そのため、Δτについて（１０）式のｆが一意の値を持つためには、３つの周波数の幅の内で、（１０）式左辺の位相が１回転しないような値に限定する必要がある。即ち、ｅｘｐ（２πｊｆΔτ）が、（ｐ−１）／Ｎτ≦ｆ＜（ｐ＋１）／Ｎτで一意となるには、
Δτ＜Ｎτ／２ （２３）
であることが必要である。従って、本実施の形態５では、（２３）式のようなΔτとなるよう、データ分配手段７を動作させる。
【００８４】
図１４に本実施の形態５の信号のサンプリング方法を表す。また、フーリエ変換後の周波数スペクトルを図１５に示す。図１５のように、本実施の形態５はエイリアジングの影響を排除するための前記実施の形態１から４までと異なり、Δτはτより長い時間間隔となっている。（２３）式の条件からすれば、Δτはτより短くても問題はないが、その場合、あまり短く取りすぎると、雑音の影響が出てくるため、Δτの間隔は適用した各装置の特性に応じて定められる。
【００８５】
図１６に本実施の形態５の信号検出手段８と畳込み判定手段１１の動作のフローチャートを示す。まず、ｐ＝０から分析を開始し（ステップ５００）、信号の有無を検出する（ステップ５０１）。信号が検出された場合は、（１０）式から左辺のｆを算出する（ステップ５０２）。このｆにおける、信号強度に｜ＳＴＦＴ（２πｐ／Ｎτ）｜^２を加算する（ステップ５０３）。これが周波数ｐ／Ｎτにおける１サイクルであり、以上が終了したらｐの値を１増やして（ステップ５０４）、再び信号の有無を確認してから（１０）式を用いてｆを計算する。
【００８６】
上述したように、畳込みによって１の周波数が３つの周波数に拡散するが、（１０）式のｆは同じ値が出てくるため、同じｆならば畳込みによる拡散で生じたスペクトルであるとして、ステップ５０２および５０３の手順によって畳込みの影響を除去する。このように従来例のように真の周波数周辺の周波数を３つ選択するのではなく、（１０）式によってあらかじめ真の周波数を求めて、この周波数における信号強度を足し合せていくため、真の周波数が離散周波数値のほぼ中間にあっても、精度よく真の周波数を算出できる。これをｐ＝０〜Ｎまで行って（ステップ５０５）、分析可能範囲をすべて分析する。
【００８７】
以上のように、本実施の形態５では、サンプリングした信号から同時に２つの信号列を生成し、それらの位相を比較することで、畳込みによってスペクトルが拡散する効果を抑え、また、加算を主体にして信号強度を算出するため、正確な周波数分析を行う周波数分析装置を得る。
【００８８】
実施の形態６．
上記実施の形態５では、２つの信号列を得るために、時間間隔τでサンプリングを行った結果を記録し、その中から、適当なタイミングのデータをそれぞれ抽出し、２つの信号列の時間差がΔτとなるようにしたものを示したが、本実施の形態６では、図１７のように、タイミングのΔτだけ異なるクロックでサンプリングを行う２つのＡ／Ｄ変換器を使用して実現してもよい。
【００８９】
図１７に本実施の形態６に係わる周波数分析装置の装置構成を表すブロック図で示す。図６と同じ構成要素には同じ符号を付す。１１は上記実施の形態５と同様の畳込み効果判定手段である。本実施の形態６の周波数分析装置は図のように２つのＡ／Ｄ変換器２ａ、２ｂを使用するが、無駄な信号をサンプリングすることがなくなるため、メモリ３ａ、３ｂの容量を小さくて済む。
【００９０】
以下、動作について説明する。帯域フィルタ１は入力された受信信号を帯域制限して分析対象の周波数信号のみ選択する。信号分配器１０では、信号を２つのＡ／Ｄ変換器２ａ、２ｂに分配する。各Ａ／Ｄ変換器２ａ、２ｂではそれぞれ時間間隔τでサンプリングを行うが、例えば、Ａ／Ｄ変換器２ｂの方がＡ／Ｄ変換器２ａよりもΔτだけ遅れたタイミングでサンプリングを行う。これによって、実施の形態５で示したような２つの信号列ｓ（ｎ）およびｓ_Δ（ｎ）が生成され、順次メモリ３ａ、３ｂに記憶される。以下の動作は上記実施の形態５の周波数分析装置の動作と同様である。
【００９１】
以上のように、本実施の形態６では、帯域制限した信号を予め分配してＡ／Ｄ変換器に入力するため、無駄な信号をサンプリングすることがなく、メモリの容量が小さくて済む。
【００９２】
実施の形態７．
上記実施の形態５および６では、受信信号の周波数分布が単一もしくは比較的疎な場合のみ分析が可能であったが、本実施の形態７では同じ周波数に畳込みによって拡散した異なる信号成分が重なりあった場合でも、信号成分が単一か否かを判定し、信号成分が単一の場合のみ分析を行う周波数分析装置について説明する。
【００９３】
図１８は本実施の形態７の周波数分析装置の装置構成を表すブロック図である。本実施の形態７では、信号系列を３つ用いてるため、実施の形態３の図８と同じ構成要素には同じ符号を付す。１１ａは本実施の形態７の特徴である、各周波数における信号成分が単一か否かを判定しながら畳込みの効果を抑えて正しい周波数を算出する畳込み効果判定手段である。
【００９４】
次に、動作について説明する。１つの信号系列ｓ（ｎ）およびそのｓ（ｎ）から時間Δτ_１およびΔτ_２だけずれた２つの信号系列ｓ_Δ１（ｎ）およびｓ_Δ２（ｎ）の計３つの信号系列をサンプリングし、それぞれ短時間フーリエ変換する。そこから、信号検出手段８および畳込み効果判定手段１１ａによって、上記実施の形態３と同様に（２１）式のような単位時間当たりの位相差を求め、信号が一つか否かを判定しながら真の周波数を求める。
【００９５】
以下、信号検出手段８および畳込み効果判定手段１１ａの動作については、図１９に示したフローチャートにしたがって説明する。まず、信号検出手段８ｂで離散値の各周波数ｆ＝ｐ／Ｎτをｐ＝０（ステップ６００）から順番に信号が存在するかの判定が行われる（ステップ６０１）。そして、信号が未検出ならば、ｐ＝ｐ＋１（ステップ６０６）とし、次のｐについて再び信号が存在するか判定を行う。
【００９６】
一方、ｐ＝０で信号が検出された場合は、ＳＴＦＴ（ω）、ＳＴＦＴ_Δ１（ω）、ＳＴＦＴ_Δ２（ω）、ＳＴＦＴ_Δ３（ω）から（２０）式を計算し（ステップ６０２）、全て同じ値を算出するか判定する（ステップ６０３）。もし、（２１）式の計算結果が一致しないならば、上述したように信号成分が複数存在することになるので、当該周波数の成分は不定であるとして、次のｐの値について信号の存在の有無を判定する。一方、（２０）式がすべて同じ値ならば信号は１つであり、周波数の分析が可能であるので実施の形態５と同様に（１１）式を満足する周波数ｆを算出する（ステップ６０４）。そして、算出した周波数ｆにおける信号強度に補正を加える（ステップ６０５）。
【００９７】
以上、いずれかの処理が終わったら周波数ｐ／Ｎτのｐをｐ＝ｐ＋１として、次の周波数成分について同様の処理を行う（ステップ６０６）。その際、ｐ≧Ｎであれば、演算する部分が重なってしまうため、ｐ≧Ｎであれば分析を終了し、ｐ＜Ｎの場合のみステップ５０１に戻る（ステップ６０７）。
【００９８】
以上のように、本実施の形態７では、信号がいくつあるか判定しながら周波数分析を行うため、信号の周波数成分が複数ある場合でも、誤った周波数を算出せず、短時間のサンプリングで周波数分析が終了し、また、加算を主体にして信号強度を算出するため、正確な周波数分析を行うことが可能な周波数分析装置を得る。
【００９９】
【発明の効果】
以上のように、本発明の周波数分析装置は、サンプリングしたデジタルデータから同時に２つの信号系列を生成し、この２つの信号系列を短時間フーリエ変換する。そして、それぞれの位相差を算出し、エイリアジングを利用して、周波数帯域全てを走査しなくても受信信号の存在する周波数を算出することができるため、分析すべき信号の周波数帯が広帯域化しても計算量の増大を抑え、エイリアジングの効果を排除し、正確な周波数分析結果を得ることができる。
【０１００】
また、帯域制限した信号を予め分配してＡ／Ｄ変換器に入力するため、無駄な信号をサンプリングすることがなく、メモリの容量が小さくて済む。
【０１０１】
また、生成した３つの信号系列を短時間フーリエ変換し、各信号系列の単位時間当たりの位相差を算出して各周波数において信号成分が単一か否かを判定しながら分析を行うため、信号成分が複数ある場合でも、誤った周波数を算出しない。
【０１０２】
また、生成した３つの信号系列を短時間フーリエ変換し、エイリアジング成分の性質を利用して行列の要素を決定し、位相差を比較して信号成分を算出するため、複数の信号が同じ周波数に重なっている場合でも、互いに分離して算出することができる。
【０１０３】
また、加算によって信号強度を算出する畳込み効果判定手段を設けたため、窓関数の効果による周波数成分の拡散による誤差を排除し、正確な周波数分析結果を得る。
【図面の簡単な説明】
【図１】 本発明の実施の形態１に係わる周波数分析装置の構成を表すブロック図である。
【図２】 本発明の実施の形態１に係わる周波数分析装置のサンプリング方法を表す図である。
【図３】 本発明の実施の形態１に係わる周波数分析装置の周波数分析方法を表す図である。
【図４】 本発明の実施の形態１に係わる周波数分析装置の周波数算出方法のフローチャートである。
【図５】 本発明の実施の形態２に係わる周波数分析装置のサンプリング方法を表す図である。
【図６】 本発明の実施の形態２に係わる周波数分析装置の装置構成を表すブロック図である。
【図７】 本発明の実施の形態３が対象とする周波数分布を表す図である。
【図８】 本発明の実施の形態３に係わる周波数分析装置の装置構成を表すブロック図である。
【図９】 本発明の実施の形態３に係わる周波数分析装置の周波数算出方法のフローチャートである。
【図１０】 本発明の実施の形態４に係わる周波数分析装置の装置構成を表すブロック図である。
【図１１】 本発明の実施の形態４に係わる周波数分析装置の周波数算出方法のフローチャートである。
【図１２】 本発明の実施の形態４に係わる周波数分析装置の別の周波数算出方法のフローチャートである。
【図１３】 本発明の実施の形態５に係わる周波数分析装置の装置構成を表すブロック図である。
【図１４】 本発明の実施の形態５に係わる周波数分析装置のサンプリング方法を表す図である。
【図１５】 本発明の実施の形態５に係わる周波数分析結果の一例である。
【図１６】 本発明の実施の形態５に係わる周波数分析装置の周波数算出方法のフローチャートである。
【図１７】 本発明の実施の形態６に係わる周波数分析装置の装置構成を表すブロック図である。
【図１８】 本発明の実施の形態７に係わる周波数分析装置の装置構成を表すブロック図である。
【図１９】 本発明の実施の形態７に係わる周波数分析装置の周波数算出方法のフローチャートである。
【図２０】 従来の周波数分析装置の装置構成を表すブロック図である。
【図２１】 短時間フーリエ変換の原理を表す図である。（ａ）サンプリングにより生成された信号列である。（ｂ）窓関数の形状を表す図である。（ｃ）上記生成された信号列と窓関数とを畳込み演算によって重み付けした後の図である。（ｄ）それぞれの時間における周波数分布を表す図である。（ｅ）各時間における周波数分布を並べて時間変化を見る図である。
【図２２】 エイリアジングの効果を表す図である。
【図２３】 窓関数との畳込みの効果を表す図である。
【図２４】 従来の周波数分析装置を表すブロック図である。
【符号の説明】
１ 帯域フィルタ、 ２、２ａ、２ｂ Ａ／Ｄ変換器、
３、３ａ、３ｂ メモリ、 ４ａ、４ｂ 重み付け手段、
５ａ、５ｂ 高速フーリエ変換手段、 ６ 表示器、
７ データ分配手段、 ８ 信号検出手段、
９、９ａ エイリアジング判定手段、 １０ 信号分配器、
１１、１１ａ 畳込み効果判定手段。[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a frequency analyzer that analyzes a frequency distribution of a received signal by short-time Fourier transform.
[0002]
[Prior art]
In a frequency analyzer that receives a radio signal from a radar or the like and analyzes the frequency component of the signal by short-time Fourier transform, it is desired that the frequency band to be analyzed be wide in order to cope with frequency hopping of the radar transmission wave. .
[0003]
Incidentally, it is generally known that a signal AD-converted at a time interval τ has a signal component of f + r / τ (r is an integer) in addition to the frequency f of the original signal. This component is called aliasing. In order to widen the frequency band while avoiding the influence of this aliasing, it is necessary to handle a signal obtained by converting data at fine time intervals in proportion to the reciprocal of the band. This means that the amount of data to be handled needs to be processed for a large amount of data in proportion to the frequency band. For example, when the frequency band is increased 16 times, the time interval of the AD-converted signal needs to be reduced by 1/16, and the number of data handled is increased 16 times. In the FFT algorithm that performs short-time Fourier transform efficiently, if the number of data is Z, the amount of calculation is Zlog for the number of data Z. ₂ It is known that it increases in proportion to Z. When the number of data Z is multiplied by Y, the amount of calculation is (ZYlog ₂ ZY) / (ZlogZ) = Y × (1+ (log ₂ Y / log ₂ Z)). As an example, when Z = 128 and Y = 16, about 25 times the amount of calculation is required, and this increase in processing load makes it difficult to analyze a wideband signal.
[0004]
In addition, the short-time Fourier transform has a problem that a true frequency component is diffused by a convolution operation using a window function and a frequency cannot be calculated with high accuracy.
[0005]
Hereinafter, a conventional frequency analyzer using short-time Fourier transform will be described, and first, problems due to aliasing will be clarified. FIG. 19 shows the configuration of a conventional frequency analyzer. In the figure, 1 is a bandpass filter that selects only a frequency component to be monitored from a received radio wave signal, 2 is an analog-to-digital converter that samples an analog signal at a constant time interval, and converts it into a digital signal. A memory for sequentially recording the input signals converted into signals, 4 is a weighting means for reading the signals recorded in the memory 3 by a certain length and multiplying them by a window function, and 5 is a weighted constant length. A Fourier transform means for sequentially performing Fourier transform on the signals of No. 1 and performing frequency analysis, and 6 is a display device for outputting the analysis results.
[0006]
Next, the operation will be described. The frequency analyzer shown in FIG. 19 calculates the frequency distribution of a received signal using the principle of short-time Fourier transform. For the short-time Fourier transform, see Hitoshi Kiya, “Multi-rate signal processing” (Shokodo) pp. Although it is described in documents such as 149 to 153, it is as follows.
[0007]
As shown in FIG. 20, short-time Fourier transform is a processing technique in which a signal is cut out by multiplying a window function having a finite length, and is subjected to Fourier transform. FIG. 20A shows an actual signal sequence obtained by sampling, FIG. 20B shows the shape of the window function, and FIG. 20C shows the weighted signal obtained by the convolution operation of the signal sequence and the window function. (D) is the figure which carried out the Fourier transform of each, and was converted into the frequency component, (e) is a figure showing the time change which arranged the frequency component in each time. The window function is shifted according to time, and the time and frequency information is examined by repeating the process.
[0008]
The short-time Fourier transform can be expressed mathematically as follows.
Assume that the signal sampled at the clock period τ is s (n) and the window function is g (n) (g (n) is n = 0, τ, 2τ,..., (N−1) τ and is not 0). It is assumed that s (n) has a waveform as shown in FIG. 20 (a) and g (n) has a waveform as shown in FIG. 20 (b). Here, a function g (n−kN) obtained by shifting g (n) by Nτ is considered, and by multiplying g (n−kN) and s (n), a signal of the k-th block is obtained as shown in Expression (1). s _k (N) is extracted.
s _k (N) = g (n−kN) s (n) (1)
It is. In FIG. _k (N) is shown.
[0009]
S in equation (1) _k The process of Fourier transforming (n) is short-time Fourier transform. Therefore, s _N The short-time Fourier transform of (n) is
[Expression 1]

It is. However, ω = 2πp / Nτ. Than this,
[Expression 2]

It is. Equation (3) represents the frequency component at the frequency p / Nτ for the input signal s (n).
[0010]
g (n) is called a window function. For example, the following Hanning window function is generally used.
g (n) = 1/2 {1 + cos2π / N (n−N / 2)}
However, n = 0 to N-1, and g (n) = 0 for other n
[0011]
FIG. 19 shows the configuration of a frequency analyzer using short-time Fourier transform that realizes the above principle. For example, when a radio signal is received, it is converted into an appropriate intermediate frequency signal and input to the band filter 1. The band filter 1 extracts only the signal in the frequency band to be analyzed from the input signal, and the A / D converter 2 samples the extracted signal at an appropriate time interval τ and converts it to a digital signal. This digitally converted signal sequence is the above-described s (n), and s (n) is sequentially stored in the memory 3. In the following, for A / D conversion, it is assumed that the signal is acquired as a complex number by using the I and Q channels, but this is different from the case where the signal is sampled as a real number. There is no essential difference except that the interval may be doubled.
[0012]
The weighting means 4 sequentially reads out the signals recorded in the memory 2, multiplies the signal sequence s (n) by the window function g (n), performs an operation represented by the equation (1), and performs fast Fourier transform means. 5, the weighted signal s _k The processing of equation (3) is sequentially performed on (n), the frequency component included in the signal at each time point is analyzed, and the result is displayed on the display device 6.
[0013]
As described above, the conventional frequency analyzer obtains the frequency distribution of the radio signal received using the short-time Fourier transform. Here, considering the properties of the short-time Fourier transform, the following features can be derived from the properties of equation (3).
[0014]
As a result of performing the short-time Fourier transform, the frequency f = p / Nτ (p is an integer) and has a discrete value. Among them, the highest frequency is (N−1) / Nτ when p = N−1. That is, according to the sampling theorem, the frequency range in which aliasing does not occur is 1 / τ. That is, the maximum frequency to be analyzed is f _w Then f _w If <1 / τ, it means that the frequency can be specified. The frequency 1 / τ is also called the Nyquist frequency.
[0015]
Also, the right side of equation (3) is
[Equation 3]

The following relationship holds. Here, r is an integer. This means that a component having the same value appears at a frequency shifted by r / τ from the true frequency p / Nτ. That is, as shown in FIG. 21, in addition to the frequency of the original signal, a frequency component called aliasing is generated every 1 / τ.
[0016]
If the range of the frequency component analyzed by the band filter 1 is limited to 1 / τ or less, such aliasing can be prevented. However, as described above, it is desired that the frequency analysis device has a wider analysis frequency, and as a result, if the frequency range to be passed through the bandpass filter is wider than 1 / τ (τ is the sampling interval), the more the frequency is. Aliasing occurs and cannot be separated from the true frequency of the signal. The above is a problem caused by aliasing.
[0017]
Next, as described above, an error in which a spectrum is diffused in principle by performing a convolution operation with a window function in a finite time range will be described. That is, if the Fourier transform of g (n) is G (2πp / Nτ),

Further, from the equations (1) and (2), when the Fourier transform of s (n) is S (ω),
[Expression 4]

Here, if it replaces like m = n-kN for convenience,
[Equation 5]

Therefore, since G (ω) has a plurality of values, the frequency distribution calculated by the short-time Fourier transform spreads around the true frequency according to the nature of the window function. FIG. 22 shows the state of this diffusion. FIG. 22A shows the original frequency distribution of the signal, and FIG. 22B shows the frequency distribution spread by the window function. The sum of the intensities of the diffused components matches the intensity of the original signal due to the nature of the convolution operation of equation (8). The degree of diffusion depends on the nature of the window function g (n). However, even with a good Hanning window function, it is inevitable that three diffusion components are generated at high and low true frequencies.
[0018]
An example of a frequency analyzer that suppresses this diffusion component and calculates a true spectrum is, for example, a frequency analyzer disclosed in Japanese Patent Laid-Open No. 6-160445. FIG. 23 shows a frequency analyzer disclosed in Japanese Patent Laid-Open No. 6-160445. In the figure, 1 is an AD converter, 2 is a memory, 3 is a window function weighting means, 4 is a fast Fourier transform means, 5 is a spectrum interpolation determination means provided on the output side of the fast Fourier transformer 4, and 6 is a display. , 7 is a spectrum interpolation unit that performs spectrum interpolation based on the determination result of the interpolation determination unit 5, and 8 is an inverse Fourier transformer that inversely converts the frequency, amplitude, and phase data of the spectrum calculated by the spectrum interpolation unit 7 into time axis data, 9 is an accumulator for accumulating the time axis data inverse Fourier transformed by the inverse Fourier transformer 8; 11 is a subtracter of the accumulated result accumulated by the accumulator 9 from the time data read out from the memory 2; A subtractor that is input to the function weighting means 3, 12 is an accumulator that accumulates the frequency domain data calculated by the spectrum interpolation means 7, and 13 is this accumulator. Vessel 12 an adder for adding the accumulated frequency-domain data obtained by subtracting the re fast Fourier transform residual data by the subtractor 11 has the frequency domain data S (t) to.
[0019]
In this frequency analyzer, if three frequency spectra exist adjacent to each other in the discrete frequency spectrum obtained by Fourier transform and the phases of the upper and lower spectra are almost the same, It is determined that the true spectrum is dispersed among the three spectra, and the frequency, amplitude, and phase of the true spectrum are obtained by interpolation processing.
[0020]
This true spectrum data is accumulated in the accumulator 12, and inverse Fourier transformed to inversely transform into time axis data. This time axis data is accumulated in the accumulator 9, and the accumulated data is subtracted from the input data string. After the true spectral line is calculated by the interpolation process, the residual energy is reduced and the line spectrum is not spread and is displayed as a line spectrum. As a result, the spectral resolution is improved. In addition, since the amplitude and the phase are approximated and calculated close to the true value by the interpolation processing, the amplitude value and the phase can be accurately measured.
[0021]
However, in this case, the data accumulated by the feedback loop is repeatedly subtracted to extract the true frequency component, so it is necessary to sample many times before the value converges, and it takes time. When the feedback loop is eliminated, a minute signal cannot be extracted.
[0022]
Further, in the above example, it is necessary to be able to correctly select three frequencies around the true frequency, but when the true frequency exists in the middle of the discrete frequency value, for example, n / Nτ and (n + 1) / Nτ. , It is easy to misjudge the frequency closest to the true frequency as (n + 1) / Nτ, and mistaken three frequencies around the true frequency are n / Nτ, (n + 1) / Nτ, (n + 2) / Since it is likely to be selected like Nτ, actual frequency calculation tends to be unstable.
[0023]
Further, in this case, the true frequency is calculated mainly by subtraction as in equation (14) of Japanese Patent Laid-Open No. 6-160445, and the calculation result is unstable when sin and cos are around 0, Error is likely to occur in frequency calculation.
[0024]
[Problems to be solved by the invention]
The present invention has been made to solve the above-described problems, and suppresses an increase in the amount of calculation associated with widening a received signal in a short-time Fourier transform, thereby efficiently and accurately frequency distribution. The purpose is to calculate.
[0025]
[Means for Solving the Problems]
A frequency analysis apparatus according to the present invention converts a band-pass filter that limits an analog input signal to be frequency-analyzed to a certain frequency band, and sequentially converts the band-limited analog input signal into a digital data string at a first time interval. An AD converter, a memory for sequentially storing the input data sequence converted by the AD converter, and a second time interval longer than the first time interval from the input data sequence stored in the memory A data distribution means for extracting a first distribution data string and a second distribution data string delayed by the first time interval from the first distribution data string; and a predetermined distribution in the first distribution data string First weighting means for multiplying and weighting the window function, and a first high-speed that performs Fourier transform on the weighted first distribution data string and converts it to first frequency domain data A Fourier transform means, a second weighting means for multiplying the second distribution data string by a predetermined window function and performing weighting, and a Fourier transform of the weighted second distribution data string to obtain a second frequency Second fast Fourier transform means for converting to region data, signal detection means for detecting whether a signal exists for each frequency of the first or second frequency region data, and when a signal is detected, And an aliasing determining means for calculating the true frequency by comparing the phases of the first and second frequency domain data.
[0026]
A band-pass filter for limiting an analog input signal to be analyzed to a certain frequency band; a signal distributor for distributing the band-limited analog input signal into two; and the one of the distributed analog input signals. A first A / D converter that converts a digital signal at a constant time interval to generate a first distribution data string, and the other analog input signal that has been distributed is converted to the first A / D converter. A second A / D converter for generating a second distribution data string by converting the digital signal at a constant time interval, delayed by a time interval shorter than the constant time interval, A first weighting unit that multiplies the distribution data string by a predetermined window function and performs weighting, and a Fourier transform of the weighted first distribution data string to convert the weighted first distribution data string into first frequency domain data. High speed A second transforming unit, a second weighting unit that multiplies the second distribution data string by a predetermined window function and performs weighting, and a Fourier transform of the weighted second distribution data string to obtain a second frequency Second fast Fourier transform means for transforming into area data, signal detecting means for detecting whether a signal exists for each frequency of the first or second frequency domain data, and when a signal is detected, By comparing the first and second frequency domain data, an aliasing determining unit that eliminates the influence of aliasing and calculates a true frequency is provided.
[0027]
A band-pass filter that limits an analog input signal to be frequency-analyzed to a certain frequency band; an AD converter that sequentially converts the band-limited analog input signal into a digital data sequence at a first time interval; and the AD A memory for sequentially storing the input data strings converted by the converter, and a first distribution data string having a second time interval longer than the first time interval from the input data string stored in the memory; A second distribution data string delayed from the first distribution data string by the first time interval, and a third time different from the first and second time intervals from the first distribution data string. Data distribution means for extracting a third distribution data string delayed by an interval, first weighting means for multiplying the first distribution data string by a predetermined window function, and performing weighting, and weighting A first fast Fourier transform means for performing Fourier transform on the first distribution data string and converting it to first frequency domain data; and a second for performing weighting by multiplying the second distribution data string by a predetermined window function. Weighting means, a second fast Fourier transform means for Fourier transforming the weighted second distribution data string to convert it into second frequency domain data, and a predetermined window function for the third distribution data string A third weighting means for performing weighting, a third fast Fourier transforming means for Fourier transforming the weighted third distribution data string and converting it to third frequency domain data, and the first Or a signal detection means for detecting whether a signal exists for each frequency of the second or third frequency domain data, and, if a signal is detected, the first, second, or third frequency domain data. Z The phase difference per unit time is compared to determine whether or not the signal component at the frequency at which the signal is detected is single, and aliasing determination that calculates the true frequency only when the signal component is single Means.
[0028]
A band-pass filter that limits an analog input signal to be frequency-analyzed to a certain frequency band; an AD converter that sequentially converts the band-limited analog input signal into a digital data sequence at a first time interval; and the AD A memory for sequentially storing the input data strings converted by the converter, and a first distribution data string having a second time interval longer than the first time interval from the input data string stored in the memory; A second distribution data string delayed from the first distribution data string by the first time interval, and a third time different from the first and second time intervals from the first distribution data string. Data distribution means for extracting a third distribution data string delayed by an interval, first weighting means for multiplying the first distribution data string by a predetermined window function, and performing weighting, and weighting A first fast Fourier transform means for performing Fourier transform on the first distribution data string and converting it to first frequency domain data; and a second for performing weighting by multiplying the second distribution data string by a predetermined window function. Weighting means, a second fast Fourier transform means for Fourier transforming the weighted second distribution data string to convert it into second frequency domain data, and a predetermined window function for the third distribution data string A third weighting means for performing weighting, a third fast Fourier transforming means for Fourier transforming the weighted third distribution data string and converting it to third frequency domain data, and the first Or signal detecting means for detecting whether a signal exists for each frequency of the second or third frequency domain data, and, if a signal is detected, the frequency at each of the first, second, and third frequencies. The true frequency is obtained by multiplying the matrix whose elements are domain data by the inverse matrix of a matrix whose elements are multiple frequencies shifted by an integer multiple of the inverse of the second time interval from the frequency at which the signal is detected. And an aliasing determining means for calculating
[0029]
Also, a band filter for limiting the analog input signal to be frequency-analyzed to a certain frequency band, and an AD for sequentially converting the band-limited analog input signal into digital data a predetermined number of times at a first time interval. A converter, a memory for sequentially storing the input data converted by the AD converter, and a first distribution data obtained by extracting the data having the first time interval from the input data stored in the memory A second distribution data sequence that is delayed from the first distribution data sequence by a second time interval that is shorter than one-half the number obtained by multiplying the first time interval by the digitally converted fixed number of times. A data distribution means for extracting the first distribution data string, a first weighting means for multiplying the first distribution data string by a predetermined window function, and weighting, and a weighted first distribution. A first fast Fourier transform means for Fourier transforming the data string to convert the data string into first frequency domain data; and a second weighting means for multiplying the second distribution data string by a predetermined window function and performing weighting. , A second fast Fourier transform means for Fourier transforming the weighted second distribution data string to convert it into second frequency domain data, and a signal for each frequency of the first or second frequency domain data. When the signal is detected and the signal is detected, the phase of the first and second frequency domain data is compared, and if this phase is the same, the true frequency is calculated, And a convolution effect determining means for performing processing for adding the signal intensities at the calculated frequency.
[0030]
In addition, a bandpass filter that limits an analog input signal to be frequency-analyzed to a certain frequency band, and the band-limited analog input signal is sequentially converted into a digital data string a predetermined number of times at a first time interval. An AD converter, a memory for sequentially storing an input data string converted by the AD converter, and a first data obtained by extracting data having the first time interval from the input data stored in the memory Distribution data string and second distribution data delayed from the first distribution data string by a second time interval shorter than one half of the first time interval multiplied by the digitally converted fixed number of times And a third time interval that is shorter than one half of the first time interval multiplied by the digitally converted fixed number of times and different from the first and second time intervals. A data distribution means for extracting a third distribution data string that is delayed from the first distribution data string by an interval, and a first window that performs weighting by multiplying the first distribution data string by a predetermined window function Weighting means, first fast Fourier transform means for Fourier transforming the weighted first distribution data string and converting it to first frequency domain data, and a predetermined window function for the second distribution data string A second weighting means for multiplying and weighting; a second fast Fourier transform means for Fourier transforming the weighted second distribution data string to convert it to second frequency domain data; A third weighting unit that multiplies the distribution data sequence by a predetermined window function and performs weighting, and a third high-speed that performs Fourier transform on the weighted third distribution data sequence and converts it to third frequency domain data. -Lier conversion means, signal detection means for detecting whether there is a signal for each frequency of the first, second or third frequency domain data, and if a signal is detected, the first, first When the phase difference per unit time of any two of the second and third frequency domain data is compared to determine whether or not the signal component at the frequency where the signal is detected is single, and the signal component is single Only when the phases of the first and second frequency domain data are compared and the phases are the same, the true frequency is calculated and the convolution effect is performed to add the signal strengths at the calculated frequency. Determination means.
[0031]
DETAILED DESCRIPTION OF THE INVENTION
Embodiment 1 FIG.
Embodiment 1 of the present invention will be described. The frequency analysis apparatus according to the first embodiment performs a short-time Fourier transform process on each of two series of signal sequences obtained by sampling with a slight timing shift from the A / D converted signal, thereby obtaining a local frequency distribution. It is what you want. In addition, the true frequency is specified by detecting the phase difference between the signal sequences generated by shifting the aliasing generated at that time by a slight timing.
[0032]
Details will be described below. FIG. 1 is a block diagram showing the configuration of the frequency analyzer according to the first embodiment. 1, the same components as those in FIG. 19 are denoted by the same reference numerals. 7 is a data distribution means for extracting and distributing two data sequences having slightly different sampling timings from the data of the signal recorded in the memory 3, and 8 is a signal detection means for determining the presence of the signal for each frequency from the Fourier transform result. Reference numeral 9 denotes aliasing determination means for comparing the results of two short-time Fourier transforms with respect to the frequency at which the signal determined by the signal detection means 8 is present, eliminating the influence of aliasing, and calculating the correct frequency.
[0033]
FIG. 2 is a diagram illustrating a sampling method of the frequency analyzer according to the first embodiment. As shown in FIG. 2, from the data sampled at the time interval Δτ, the signal sequence s (n) at an appropriate time interval τ and the signal sequence s delayed by Δτ from each signal of the s (n). _Δ Consider (n). If the frequency of the signal is f (= p / Nτ, p is an integer), the signal sequences s (n) and s _Δ (N) means that a phase delay of 2πfΔτ occurs.
[Formula 6]

The relationship is established. Here, since the coefficient exp (2πjfΔτ) is a constant, the coefficient exp (2πjfΔτ) is not changed by each process of the short-time Fourier transform including weighting and Fourier transform, and the signal sequence s _Δ (N) has a phase difference with a uniform frequency f with respect to s (n). For this reason, s _Δ (N) Fourier transform STFT _Δ If (2πp / Nτ),
[Expression 7]

From the equation (10), with respect to the frequency f at which STFT (2πp / Nτ) ≠ 0, the original data string s (n) and the data string s sampled by shifting by a minute time Δτ. _Δ When the short-time Fourier transforms of (n) are compared, the difference is in the coefficient exp (2πjfΔτ). Therefore, if equation (10) is transformed,
[Equation 8]

For each FFT component on the right side of equation (11), the data of the interval τ is used, so the frequency f with respect to p is
f = p / Nτ, (p + N) / Nτ, (p + 2N) / Nτ...
One of them. Therefore, r can be expressed as f = (p + rN) / Nτ, where n is an integer.
[Equation 9]

It is. At that time, if the phase does not rotate once due to the term of −2πjrNΔτ / Nτ, r can be uniquely obtained. That is,
[Expression 10]

This is the case. When the equation (12) is transformed,
1 / Δτ> r _M The relationship / τ (14) is obtained. 1 / Δτ is a frequency range that can be specified by the frequency analysis apparatus according to the first embodiment. r _M Is too large, the expression (14) cannot be satisfied, and the bandpass filter 1 selects an appropriate frequency band.
[0034]
FIG. 3 shows a graph representing the relationship of equation (13). FIG. 3 is a diagram of frequency analysis results with the vertical axis representing amplitude and the horizontal axis representing frequency. The frequency value is a discrete value and its minimum unit is 1 / Nτ. Further, 1 / Δτ is the maximum frequency that can be analyzed by this frequency analyzer, and the band filter 1 limits the band so that a frequency above 1 / Δτ is not input. On the other hand, since the signal sequence s (n) at the sampling interval τ is Fourier-transformed for a short time, the frequency range in which aliasing does not occur is 1 / τ, and therefore is shifted by 1 / τ from the true signal frequency. Aliasing occurs.
[0035]
In the figure, for example, there is a case where there is a true signal component at the frequency (p + 2N) / Nτ. Signal sequences s (n) and s _Δ As seen from (n), the frequency range in which aliasing does not occur is 1 / τ, whereas the received signal input from the bandpass filter 1 exceeds the frequency range, so the right side of equation (11) is measured. And p and r on the left side are obtained therefrom.
[0036]
The above is the principle of the frequency analyzer according to the first embodiment. The frequency analyzer of FIG. 1 performs short-time Fourier transform based on the above principle. A reception signal obtained by converting the received radio wave into an appropriate intermediate frequency is input to the band filter 1 and band-limited. In the A / D converter, the selected signal is sampled at an appropriate short time interval Δτ that is equal to or less than the reciprocal of the Nyquist frequency, and sequentially stored in the memory 3. In the data distribution means 7, two series of signal sequences s (n) and s with appropriate intervals τ are respectively determined from the signals stored in the memory. _Δ (N) is selected and distributed to the weighting means 4a and 4b, respectively. The weighting means 4a and 4b are weighted by the window function and then input to the Fourier transform means 5a and 5b. In the Fourier transform means 5a and 5b, weighted s (n) and s _Δ (N) Fourier transform is performed, and STFT (ω) and STFT _Δ Find (ω).
[0037]
The signal detection unit 8 examines the result of the output STFT (ω) of the Fourier transform unit 5a, determines the frequency at which the signal exists, and notifies the aliasing determination unit 9 of these frequencies. In the aliasing determination means 9, STFT (ω) and STFT _Δ (Ω) is compared, and the left side exp (2πjfΔτ) of equation (10) is calculated.
[0038]
Processing performed by the signal detection means 8 and the aliasing determination means 9 will be described with reference to the flowchart shown in FIG. First, the signal detection means 8 determines whether or not there is a signal in order from the p = 0 (step 100) for each discrete frequency f = (p + rN) / Nτ (step 101). If no signal is detected, the frequency component of the received signal is 0 for all r with respect to p at that time (in this case, p = 0) (step 102). That is, if no signal is detected at p = 0, the amplitude at frequencies f = N / Nτ, 2N / Nτ, 3N / Nτ,...
[0039]
On the other hand, if a signal is detected at p = 0, STFT (ω) and STFT _Δ From (ω), r that satisfies the expression (11) is calculated (step 103). The frequency component of the received signal for r other than the calculated value is 0 (step 104). For example, when r = 3 is calculated, the amplitude at other frequencies f = (3 + N) / Nτ, (3 + 2N) / Nτ, (3 + 4N) / Nτ,...
[0040]
When p = 0, when the signal is detected or not detected, the same processing is performed for the next frequency, that is, p = 1 (step 105). At this time, if p ≧ N, the parts to be calculated overlap, so if p ≧ N, the analysis is terminated, and only when p <N, the process returns to step 101 (step 106).
[0041]
In this way, from p = 0 to N, that is, by simply detecting the presence / absence of a signal from a discrete value frequency f from 0 to 1 / τ, actually from 0 to r in FIG. _M Since the signal band up to / τ is analyzed, it is not necessary to analyze all the frequencies, and the calculation amount can be greatly reduced.
[0042]
For example, when the number of samples is N = 128, considering that the frequency band is expanded 16 times,
(1) When expanding the frequency band 16 times with a conventional frequency analyzer
Since the number of data to be subjected to Fourier transform is 16 times, the amount of calculation is (128 × 16) log ₂ (128 × 16) / 16log ₂ It becomes 128 times, that is, 25 times.
(2) When the frequency band is expanded 16 times with the frequency analyzer of the first embodiment
There is no change in the number of data to be subjected to Fourier transform. And two signal sequences s (n) and s _Δ Since the calculation is performed twice in order to Fourier transform (n), it is almost double. However, although it is also necessary to perform the calculation of equation (11), this is at most N times (128 times in this case), which is about 1/10 of the calculation amount of the Fourier transform, so the increase in the calculation amount Has little effect.
Thus, if the frequency analyzer of the first embodiment is used, the calculation load of about 12 times can be reduced.
[0043]
As described above, in the first embodiment, the signal sequence s (n) sampled at the time interval τ and the signal sequence s sampled at the time interval τ shifted from the signal sequence s (n) by the minute time Δτ. _Δ The phase difference exp (2πjfΔτ) generated by shifting by (Δτ) is calculated by performing a short-time Fourier transform on (n) and comparing them. Due to this phase difference, r is obtained by the equation (11) using aliasing that occurs as the received signal becomes wider, so that it is possible to scan the received signal without scanning all the frequency bands output from the band filter 1. Since the existing frequency can be calculated, an increase in the amount of calculation can be suppressed even when the bandwidth is widened, the effect of aliasing can be eliminated, and an accurate frequency analysis result can be obtained.
[0044]
Embodiment 2. FIG.
In the first embodiment, two signal sequences s (n) and s _Δ In order to obtain (n), an example in which data having an appropriate time interval τ is extracted from a group of data uniformly sampled at Δτ and distributed is shown. For example, as shown in FIG. May be realized by using two A / D converters that sample with different clocks by Δτ.
[0045]
FIG. 6 is a block diagram showing the device configuration of the frequency analyzer according to the second embodiment. The same components as those in FIG. 1 are denoted by the same reference numerals. Reference numeral 10 denotes a signal distributor. The frequency analysis apparatus according to the second embodiment uses two A / D converters 2a and 2b as shown in the figure, but since unnecessary signals are not sampled, the capacity of the memories 3a and 3b can be reduced. .
[0046]
The operation will be described below. The band filter 1 limits the band of the input received signal and selects only the frequency signal to be analyzed. The signal distributor 10 distributes the signal to the two A / D converters 2a and 2b. Each A / D converter 2a, 2b samples at a time interval τ. For example, the A / D converter 2b samples at a timing delayed by Δτ from the A / D converter 2a. As a result, two signal sequences s (n) and s as shown in the first embodiment are used. _Δ (N) is generated and sequentially stored in the memories 3a and 3b. The following operation is the same as the operation of the frequency analyzer of the first embodiment. Short-time Fourier transform STFT (ω), STFT of each signal sequence _Δ (Ω) is compared, and exp (2πjfΔτ) is calculated from the equation (11).
[0047]
As described above, in the second embodiment, the band-limited signal is distributed in advance and inputted to the A / D converter, so that useless signals are not sampled and the memory capacity is small.
[0048]
Embodiment 3 FIG.
The above equation (11) holds only when there is a single signal included in STFT (ω) for each p. Therefore, in the first and second embodiments, analysis is possible only when the frequency distribution of the received signal is single or relatively sparse. In the third embodiment, when a plurality of frequency components that are mutually aliasing components are included in the received signal when the frequency analyzer of the first or second embodiment is used, the signal includes a plurality of aliases. A description will be given of a frequency analyzer that determines whether or not a ging component is included and performs an analysis only when it is determined that the sing component is not included.
[0049]
FIG. 7 is a frequency spectrum diagram of a signal targeted by the frequency analyzer of the third embodiment. The dotted line represents the spectrum actually measured, and the solid line represents the true frequency component of the signal. When the frequency analyzers of the first and second embodiments are used as shown in FIG. ₁ , W ₃ Is included in the received signal, its frequency is changed to f ₁ , F ₃ Then,
[Expression 11]

Also,
[Expression 12]

When these equations are Fourier transformed respectively,
[Formula 13]

Also,
[Expression 14]

It is.
[0050]
What is found from the actual signal is STFT _Δ Only the values of (2πp / Nτ) and STFT (2πp / Nτ). Therefore, as in the first or second embodiment, the STFT _Δ Even if (2πp / Nτ) / STFT (2πp / Nτ) is obtained, w ₁ , W ₃ F is correct because it is mixed ₁ , F ₃ Cannot be obtained.
[0051]
For this reason, in the third embodiment, as shown in FIG. 8, the frequency analyzer is divided into one signal sequence s (n) having a time interval τ and two different s (n). Minute time difference Δτ ₁ , Δτ ₂ Signal series with _Δ1 (N), s _Δ2 (N) is generated and short-time Fourier transform is performed for each. In the figure, reference numeral 9a denotes aliasing determination means for calculating the true frequency while determining whether or not the signal is single.
[0052]
Where the aliasing frequency is f ₁ = P / Nτ, f ₂ = (P + N) / Nτ, f ₃ = (P + 2N) / Nτ, ... f _{r + 1} = (P + rN) / Nτ, and each signal intensity at each aliasing frequency is w ₁ , W ₂ , W ₃ ... w _r And Then, generalizing the case where multiple signals are included,
[Expression 15]

[Expression 16]

It is. s _Δ2 (N) is also expressed as Δτ in equation (16). ₁ Respectively, Δτ ₂ It is the same as that replaced.
[0053]
Here, for example, the frequency f ₃ If there is only one true signal component in w ₃ Only has a non-zero value and w ₃ W other than ₁ , W ₂ , W ₄ ... w _r Are all zero.
[0054]
Therefore, in that case, s (n), s _Δ1 (N), s _Δ2 When (n) is Fourier-transformed, the results are as follows.
[Expression 17]

[0055]
Conversely, the signal component is not single, for example, the frequency f ₃ Frequency f ₁ Signal component w ₁ Is present in place of the equations (17), (18) and (19),
[Expression 18]

[0056]
Here, the following amount is considered as a phase difference per unit time between each signal sequence.
[Equation 19]

In Equation (20), ∠ represents a phase.
For example, the signal component has a frequency f ₃ If the right side of equations (17) to (19) is substituted into equation (21), all 2πf ₃ Is obtained, but the signal component has a frequency f. ₃ And f ₁ If the right side of equations (17a) to (19a) is substituted into equation (20), a constant value cannot be obtained. Therefore, it is determined whether there are a plurality of signal components by using this equation (20).
[0057]
FIG. 8 is a block diagram showing the configuration of the frequency analyzer according to the third embodiment. The same components as those in FIG. 1 are denoted by the same reference numerals. One signal sequence s (n) and its time Δτ from s (n) ₁ And Δτ ₂ Two signal sequences s shifted by _Δ1 (N) and s _Δ2 A total of three signal sequences (n) are sampled and short-time Fourier transformed. From there, the signal detection means 8b and aliasing determination means 9b determine the phase difference per unit time as in the above equation (20), and determine the true frequency while determining whether there is one signal.
[0058]
Hereinafter, the operations of the signal detection means 8b and the aliasing determination means 9b will be described with reference to the flowchart shown in FIG. First, the signal detection means 8b determines whether or not there is a signal in order from p = 0 (step 200) for each discrete frequency f = (p + rN) / Nτ (step 201). If no signal is detected, the frequency components of the received signal are 0 for all r with respect to p at that time (in this case, p = 0) (step 202). That is, as in the case of the first embodiment, if p = 0 and no signal is detected, the frequencies f = N / Nτ, 2N / Nτ, 3N / Nτ,. The frequency component of is also zero.
[0059]
On the other hand, if a signal is detected at p = 0, STFT (ω), STFT _Δ1 (Ω), STFT _Δ2 (Ω), STFT _Δ3 The equation (20) is calculated from (ω) (step 204), and it is determined whether all the same values are calculated (step 205). If the calculation result of the equation (21) does not match, there are a plurality of signal components as described above, and the frequency component is indefinite (step 203). On the other hand, if the result of the equation (20) is the same value, there is one signal and the frequency can be analyzed, so that r satisfies the equation (11) as in the first embodiment. _T Is calculated (step 206). And the calculated r _T The component of frequency (p + rN) / Nτ is 0 for r other than (step 207).
[0060]
As described above, when any of the processes is completed, the frequency (p + rN) / Nτ is set to p = p + 1, and the same process is performed for the next frequency component (step 208). At this time, if p ≧ N, the parts to be calculated overlap, so if p ≧ N, the analysis is terminated, and the process returns to step 201 only when p <N (step 209).
[0061]
In this way, it is only necessary to detect the presence or absence of a signal from p = 0 to N−1, that is, a discrete frequency f from 0 to 1 / τ, and it is not necessary to analyze all frequencies, and the calculation amount is greatly Can be reduced.
[0062]
Further, since frequency analysis is performed while determining how many signals are present, an erroneous frequency is not calculated even when there are a plurality of signal components.
[0063]
Embodiment 4 FIG.
In the fourth embodiment, two or more signal sequences shifted by different times are subjected to a short-time Fourier transform to determine whether or not there is one frequency component included in each frequency of the discrete value, and a plurality of signals A frequency analysis apparatus that can eliminate the influence of aliasing even if it is included will be described.
[0064]
Here, an example of performing short-time Fourier transform on three signal sequences when there are a plurality of independent signals having frequencies that become mutual aliasing components as shown in FIG. 7 of the third embodiment.
[0065]
FIG. 10 is a block diagram showing the configuration of the frequency analyzer according to the fourth embodiment. In the figure, the same components as those in FIG. Reference numeral 9b denotes aliasing determination means according to the third embodiment. Similar to the third embodiment, a plurality of time differences Δτ with respect to the signal sequence s (n). ₁ , Δτ ₂ , Δτ ₃ Signal series with _Δ1 (N), s _Δ2 (N), s _Δ3 (N) is generated by sampling and short-time Fourier transform is performed. Where Δτ ₁ <Δτ ₂ <Δτ ₃ Where r / τ <1 / Δτ in FIG. ₁ Suppose that Spectral STFT (2πp / Nτ) after short-time Fourier transform of these four signal sequences, STFT _Δ1 (2πp / Nτ), STFT _Δ2 (2πp / Nτ), STFT _Δ3 For each (2πp / Nτ), a phase difference such as the expression (11) is considered.
[0066]
Here, as an aliasing property, since an aliasing component is generated every 1 / τ from the true frequency, the frequency at which aliasing may occur is
[Expression 20]

The amplitude at each frequency is expressed as w ₁ , W ₂ ... w _{r + 1} Then, the four signal sequences can be expressed by the following matrix.
[Expression 21]

As
X = SA (22)
It is.
[0067]
Here, for each element of the matrix X, four signal sequences are actually subjected to a short-time Fourier transform to obtain their values, and for each element of the matrix S, the properties of the aliasing component as described above are obtained. (21) If p in each equation is fixed, the value is determined. Therefore, using the general inverse of S,
A = S ^-1 X = (S ^* S) ^-1 ・ S ^* ・ X (22a)
By calculating, each element of the matrix A can be obtained.
[0068]
The aliasing determining means 9b determines the presence / absence of a signal for p having frequencies from 0 to 1 / τ based on the above principle, and if there is a signal, the matrix A is expressed by the above equation (23a). Find each element.
[0069]
FIG. 11 is a flowchart showing this process. Hereinafter, it demonstrates according to a figure. First, the analysis is started from p = 0 in the equation (21) (step 300). Then, the presence / absence of the signal is confirmed (step 301). For example, if the signal is not found when p = 0, the amplitude = 0 in each frequency of equation (21), that is, there is no signal (step 302). , P is incremented by 1 (step 304), and the process returns to step 301 to check the presence or absence of a signal. On the other hand, when a signal is found, as described above, each element of the matrix A, that is, the amplitude at each frequency is calculated from the three signal sequences using the equation (23a), so that a plurality of signal components are obtained. Even if it is, it can be calculated correctly. When the calculation is completed, p is increased by 1 (step 304), and the signal is analyzed again. Then, when analyzing up to p = N−1, frequency 0 to r _M Since the analysis is finished up to / τ, the analysis is terminated.
[0070]
In this way, simply by analyzing the frequency from 0 to 1 / τ, the frequency 0 to r _M Since the analysis is performed up to / τ, an efficient analysis can be performed even if the signal is widened.
[0071]
Even if aliasing components and multiple true frequency components are overlapped, three signal sequences that are shifted by a minute time that differ from each other are generated, and phase differences are compared to separate multiple frequency components. Since the correct frequency can be calculated, accurate analysis can be performed.
[0072]
In the above description, an example is shown in which all signal components are calculated using Equation (22). However, since Equation (22) requires more computation than Equation (11), it is shown in Embodiment 3. As described above, a step of determining whether or not the signal is single may be inserted, and when the signal is single, the frequency may be calculated using the equation (11).
[0073]
FIG. 12 shows a flowchart for calculating a signal component while determining whether or not the signal component is single. 11 differs from FIG. 11 in that it is determined in step 404 whether or not the signal is single using the equation (20) as in the third embodiment. If there is a single signal component, the calculation is performed in the same manner as in the first embodiment. If there are a plurality of signal components, all the signal components are calculated using equation (22) as described above in the fourth embodiment.
[0074]
In the fourth embodiment, since three signal sequences are generated, it is possible to cope with the case where the number of frequency components is two or less. However, if the weighting means and the fast Fourier transform means are further increased, more frequency components are generated. Even if it overlaps with the aliasing component, it can be accurately separated.
[0075]
As described above, in the fourth embodiment, the generated three signal sequences are Fourier-transformed for a short time, matrix elements are determined using the properties of aliasing components, and signal components are compared by comparing phase differences. In order to calculate, even when a plurality of signals overlap the same frequency, they can be calculated separately from each other.
[0076]
Embodiment 5 FIG.
As described above, the first to fourth embodiments have described the frequency analysis apparatus that is efficient while suppressing the increase in the amount of calculation while eliminating the aliasing components that are generated with the increase in the bandwidth. A frequency analyzer capable of performing accurate frequency analysis by sampling in a short time while eliminating errors due to diffusion of frequency components due to function effects will be described.
[0077]
FIG. 13 is a block diagram showing a device configuration of the frequency analyzer according to the fifth embodiment. In the figure, the same components as those in FIG. Reference numeral 11 denotes a convolution effect determination means.
[0078]
The operation will be described below. The bandpass filter 1 limits the frequency band of the signal to an appropriate band, and the A / D converter 2 samples at an appropriate time interval τ. The sampled digital signal is stored in the memory 3, and the data distribution means 7 distributes the data to the weighting means 4a and 4b. For one signal series s (n), only Δτ for the other s (n). Shifted _Δ (N) is generated. Each of the weighting means 4a and 4b performs weighting by a window function, and the fast Fourier transform means 5a and 5b performs Fourier transform. The signal detection means 8 checks each frequency after the Fourier transform to confirm the presence or absence of the signal, and the convolution effect judgment means 11 compares the phases of the frequency spectra after the Fourier transform, Specify the frequency.
[0079]
With this configuration, it is possible to suppress diffusion of frequency components due to convolution with short-time sampling, and an accurate frequency spectrum is obtained not only for continuous waves but also for intermittent waves.
[0080]
Here, the convolution effect by the window function has the following characteristics. That is, (A) the sum of the signal intensities of each frequency spread by convolution substantially matches the original signal intensity. (B) For each frequency component diffused by convolution from one signal component (frequency f), equation (11) has the same value. (C) By using the Hanning window function, the influence of spreading due to convolution of one signal can be limited to three frequency values around the true frequency.
[0081]
(A) is due to the nature of the convolution operation, as shown in FIG. 23 and the description of equation (8). In this case, since the frequency takes a discrete value, it approximately matches. For (B), STFT, STFT _Δ In both cases, since the same original frequency component is diffused by the window function g (n), the mutual phase relationship is the same, that is, the expression (10) representing the phase does not take different values. .
[0082]
Therefore, in order to eliminate the spread of the frequency distribution due to convolution from these and obtain the correct frequency distribution, it is determined by Equation (10) that the signal component of the detected frequency is the result of diffusion due to the convolution, If it is determined that the signal is spread due to error, the signal intensities of the diffused signal components are added together, and there may be a signal component having the signal intensity added to the frequency obtained from equation (10).
[0083]
As shown in (C), frequency spreading can be limited to three frequency values around the true frequency by using a Hanning window function. Therefore, in order for f in Equation (10) to have a unique value for Δτ, it is necessary to limit the value of the left side of Equation (10) to a value that does not rotate once within the three frequency widths. That is, for exp (2πjfΔτ) to be unique when (p−1) / Nτ ≦ f <(p + 1) / Nτ,
Δτ <Nτ / 2 (23)
It is necessary to be. Therefore, in the fifth embodiment, the data distribution unit 7 is operated so as to satisfy Δτ as expressed by the equation (23).
[0084]
FIG. 14 shows a signal sampling method according to the fifth embodiment. Further, FIG. 15 shows the frequency spectrum after Fourier transform. As shown in FIG. 15, the fifth embodiment is different from the first to fourth embodiments for eliminating the influence of aliasing, and Δτ is a time interval longer than τ. According to the condition of equation (23), there is no problem even if Δτ is shorter than τ, but in that case, if it is made too short, the influence of noise appears, so the interval of Δτ is the characteristic of each device applied. It is decided according to.
[0085]
FIG. 16 shows a flowchart of operations of the signal detection means 8 and the convolution determination means 11 of the fifth embodiment. First, analysis is started from p = 0 (step 500), and the presence or absence of a signal is detected (step 501). If a signal is detected, f on the left side is calculated from equation (10) (step 502). The signal intensity at this f is | STFT (2πp / Nτ) | ² Are added (step 503). This is one cycle at the frequency p / Nτ. When the above is completed, the value of p is incremented by 1 (step 504), the presence / absence of a signal is confirmed again, and f is calculated using equation (10).
[0086]
As described above, the frequency of 1 is spread to three frequencies by convolution, but f in equation (10) has the same value, so if it is the same f, it is assumed that it is a spectrum generated by spreading by convolution The influence of convolution is removed by the procedure of steps 502 and 503. In this way, instead of selecting three frequencies around the true frequency as in the conventional example, the true frequency is obtained in advance by the equation (10), and the signal intensity at this frequency is added together. Even if the frequency is approximately in the middle of the discrete frequency value, the true frequency can be calculated with high accuracy. This is performed until p = 0 to N (step 505), and the entire analysis possible range is analyzed.
[0087]
As described above, in the fifth embodiment, two signal sequences are generated simultaneously from the sampled signals, and the phases are compared, thereby suppressing the effect of spectrum spreading due to convolution and adding mainly. Thus, in order to calculate the signal intensity, a frequency analyzer that performs accurate frequency analysis is obtained.
[0088]
Embodiment 6 FIG.
In the fifth embodiment, in order to obtain two signal sequences, the result of sampling at a time interval τ is recorded, and data at an appropriate timing is extracted from the result, and the time difference between the two signal sequences is calculated. In the sixth embodiment, as shown in FIG. 17, it can be realized by using two A / D converters that perform sampling with different clocks by Δτ as shown in FIG. Good.
[0089]
FIG. 17 is a block diagram showing the device configuration of the frequency analyzer according to the sixth embodiment. The same components as those in FIG. 6 are denoted by the same reference numerals. Reference numeral 11 denotes convolution effect determination means similar to that of the fifth embodiment. The frequency analysis apparatus according to the sixth embodiment uses two A / D converters 2a and 2b as shown in the figure, but since unnecessary signals are not sampled, the capacity of the memories 3a and 3b can be reduced. .
[0090]
The operation will be described below. The band filter 1 limits the band of the input received signal and selects only the frequency signal to be analyzed. The signal distributor 10 distributes the signal to the two A / D converters 2a and 2b. Each A / D converter 2a, 2b samples at a time interval τ. For example, the A / D converter 2b samples at a timing delayed by Δτ from the A / D converter 2a. As a result, two signal sequences s (n) and s as shown in the fifth embodiment are used. _Δ (N) is generated and sequentially stored in the memories 3a and 3b. The following operation is the same as that of the frequency analyzer of the fifth embodiment.
[0091]
As described above, in the sixth embodiment, since the band-limited signal is distributed in advance and inputted to the A / D converter, unnecessary signals are not sampled and the memory capacity can be reduced.
[0092]
Embodiment 7 FIG.
In the fifth and sixth embodiments, analysis is possible only when the frequency distribution of the received signal is single or relatively sparse, but in the seventh embodiment, different signal components spread by convolution to the same frequency are present. A description will be given of a frequency analyzer that determines whether or not there is a single signal component even when they overlap, and performs analysis only when the signal component is single.
[0093]
FIG. 18 is a block diagram showing a device configuration of the frequency analyzer of the seventh embodiment. In the seventh embodiment, since three signal sequences are used, the same components as those in FIG. 8 of the third embodiment are denoted by the same reference numerals. Reference numeral 11a denotes a convolution effect determining unit that calculates the correct frequency while suppressing the effect of convolution while determining whether there is a single signal component at each frequency, which is a feature of the seventh embodiment.
[0094]
Next, the operation will be described. One signal sequence s (n) and its time Δτ from s (n) ₁ And Δτ ₂ Two signal sequences s shifted by _Δ1 (N) and s _Δ2 A total of three signal sequences (n) are sampled and short-time Fourier transformed. From there, the signal detection means 8 and the convolution effect determination means 11a obtain the phase difference per unit time as in the equation (21) as in the third embodiment, and determine whether or not there is only one signal. Find the true frequency.
[0095]
Hereinafter, the operations of the signal detection means 8 and the convolution effect determination means 11a will be described with reference to the flowchart shown in FIG. First, the signal detection means 8b determines whether or not there is a signal in order from each of the discrete frequencies f = p / Nτ from p = 0 (step 600) (step 601). If no signal is detected, p = p + 1 (step 606), and it is determined again whether the signal exists for the next p.
[0096]
On the other hand, if a signal is detected at p = 0, STFT (ω), STFT _Δ1 (Ω), STFT _Δ2 (Ω), STFT _Δ3 The equation (20) is calculated from (ω) (step 602), and it is determined whether all the same values are calculated (step 603). If the calculation result of the equation (21) does not match, there are a plurality of signal components as described above. Therefore, it is assumed that the frequency component is indefinite, and the presence of the signal exists for the next p value. Determine presence or absence. On the other hand, if all of the equations (20) have the same value, the number of signals is one and the frequency can be analyzed. Therefore, the frequency f satisfying the equation (11) is calculated as in the fifth embodiment (step 604). . Then, the signal intensity at the calculated frequency f is corrected (step 605).
[0097]
As described above, when any of the processes is completed, the frequency p / Nτ is set to p = p + 1, and the same process is performed for the next frequency component (step 606). At this time, if p ≧ N, the parts to be calculated overlap, so if p ≧ N, the analysis is terminated, and only when p <N, the process returns to step 501 (step 607).
[0098]
As described above, in the seventh embodiment, frequency analysis is performed while determining how many signals are present. Therefore, even when there are a plurality of frequency components of the signal, an erroneous frequency is not calculated, and the frequency is obtained by short-time sampling. Since the analysis is completed and the signal intensity is calculated mainly by addition, a frequency analyzer capable of performing accurate frequency analysis is obtained.
[0099]
【The invention's effect】
As described above, the frequency analyzer of the present invention generates two signal sequences simultaneously from sampled digital data, and performs a short-time Fourier transform on the two signal sequences. Since each phase difference is calculated and aliasing can be used to calculate the frequency at which the received signal exists without scanning the entire frequency band, the frequency band of the signal to be analyzed becomes wider. However, an increase in the amount of calculation can be suppressed, the effect of aliasing can be eliminated, and an accurate frequency analysis result can be obtained.
[0100]
Further, since the band-limited signal is distributed in advance and inputted to the A / D converter, useless signals are not sampled, and the memory capacity can be reduced.
[0101]
In addition, short-time Fourier transform is performed on the three generated signal sequences, the phase difference per unit time of each signal sequence is calculated, and analysis is performed while determining whether or not there is a single signal component at each frequency. Even if there are multiple components, the wrong frequency is not calculated.
[0102]
In addition, the three generated signal sequences are Fourier-transformed for a short time, matrix elements are determined using the properties of aliasing components, and signal components are calculated by comparing phase differences. Even if they overlap, they can be calculated separately from each other.
[0103]
In addition, since the convolution effect determining means for calculating the signal intensity by addition is provided, an error due to the diffusion of the frequency component due to the effect of the window function is eliminated, and an accurate frequency analysis result is obtained.
[Brief description of the drawings]
FIG. 1 is a block diagram showing a configuration of a frequency analyzer according to a first embodiment of the present invention.
FIG. 2 is a diagram illustrating a sampling method of the frequency analyzer according to the first embodiment of the present invention.
FIG. 3 is a diagram illustrating a frequency analysis method of the frequency analyzer according to the first embodiment of the present invention.
FIG. 4 is a flowchart of a frequency calculation method of the frequency analyzer according to the first embodiment of the present invention.
FIG. 5 is a diagram illustrating a sampling method of the frequency analyzer according to the second embodiment of the present invention.
FIG. 6 is a block diagram showing a device configuration of a frequency analyzer according to a second embodiment of the present invention.
FIG. 7 is a diagram illustrating a frequency distribution targeted by a third embodiment of the present invention.
FIG. 8 is a block diagram showing a device configuration of a frequency analyzer according to a third embodiment of the present invention.
FIG. 9 is a flowchart of a frequency calculation method of the frequency analyzer according to the third embodiment of the present invention.
FIG. 10 is a block diagram showing a device configuration of a frequency analyzer according to a fourth embodiment of the present invention.
FIG. 11 is a flowchart of a frequency calculation method of the frequency analyzer according to the fourth embodiment of the present invention.
FIG. 12 is a flowchart of another frequency calculation method of the frequency analyzer according to the fourth embodiment of the present invention.
FIG. 13 is a block diagram showing a device configuration of a frequency analyzer according to a fifth embodiment of the present invention.
FIG. 14 is a diagram illustrating a sampling method of the frequency analyzer according to the fifth embodiment of the present invention.
FIG. 15 is an example of a frequency analysis result according to the fifth embodiment of the present invention.
FIG. 16 is a flowchart of a frequency calculation method of the frequency analyzer according to the fifth embodiment of the present invention.
FIG. 17 is a block diagram showing a device configuration of a frequency analyzer according to a sixth embodiment of the present invention.
FIG. 18 is a block diagram showing a device configuration of a frequency analyzer according to a seventh embodiment of the present invention.
FIG. 19 is a flowchart of a frequency calculation method of the frequency analyzer according to the seventh embodiment of the present invention.
FIG. 20 is a block diagram showing a device configuration of a conventional frequency analyzer.
FIG. 21 is a diagram illustrating the principle of short-time Fourier transform. (A) A signal sequence generated by sampling. (B) It is a figure showing the shape of a window function. (C) It is a figure after weighting the said produced | generated signal sequence and window function by the convolution calculation. (D) It is a figure showing the frequency distribution in each time. (E) It is a figure which arranges frequency distribution in each time, and looks at a time change.
FIG. 22 is a diagram illustrating the effect of aliasing.
FIG. 23 is a diagram illustrating the effect of convolution with a window function.
FIG. 24 is a block diagram showing a conventional frequency analyzer.
[Explanation of symbols]
1 band filter, 2, 2a, 2b A / D converter,
3, 3a, 3b memory, 4a, 4b weighting means,
5a, 5b Fast Fourier transform means, 6 display,
7 data distribution means, 8 signal detection means,
9, 9a aliasing judging means, 10 signal distributor,
11, 11a Convolution effect determination means.

Claims (6)

A bandpass filter that limits the analog input signal to be frequency analyzed to a certain frequency band;
An AD converter that sequentially converts the band-limited analog input signal into a digital data sequence at a first time interval;
A memory for sequentially storing input data strings converted by the AD converter;
A first distribution data string having a second time interval longer than the first time interval from the input data string stored in the memory, and the first time interval from the first distribution data string. Data distribution means for extracting the delayed second distribution data string;
A first weighting means for multiplying the first distribution data string by a predetermined window function and performing weighting;
A first fast Fourier transform means for Fourier transforming the weighted first distribution data string and transforming it into first frequency domain data;
A second weighting means for multiplying the second distribution data string by a predetermined window function and performing weighting;
A second fast Fourier transform means for Fourier transforming the weighted second distribution data string and transforming it into second frequency domain data;
Signal detection means for detecting whether a signal exists for each frequency of the first or second frequency domain data;
A frequency analyzer comprising: an aliasing determination unit that calculates a true frequency by comparing phases of the first and second frequency domain data when a signal is detected. A bandpass filter that limits the analog input signal to be frequency analyzed to a certain frequency band;
A signal distributor for distributing the band-limited analog input signal into two;
A first A / D converter that converts one of the distributed analog input signals into a digital signal at regular time intervals to generate a first distributed data string;
The other distributed analog input signal is converted into a digital signal at a fixed time interval, delayed by a time interval shorter than the fixed time interval from the digital conversion of the first A / D converter. A second A / D converter that generates a distribution data sequence of the first weighting means, a first weighting unit that performs weighting by multiplying the first distribution data sequence by a predetermined window function,
A first fast Fourier transform means for Fourier transforming the weighted first distribution data string and transforming it into first frequency domain data;
A second weighting means for multiplying the second distribution data string by a predetermined window function and performing weighting;
A second fast Fourier transform means for Fourier transforming the weighted second distribution data string and transforming it into second frequency domain data;
Signal detecting means for detecting whether a signal exists for each frequency in the first or second frequency domain data;
A frequency comprising: an aliasing determination means for calculating a true frequency by removing the influence of aliasing by comparing the first and second frequency domain data when a signal is detected Analysis equipment. A bandpass filter that limits the analog input signal to be frequency analyzed to a certain frequency band;
An AD converter that sequentially converts the band-limited analog input signal into a digital data sequence at a first time interval;
A memory for sequentially storing input data strings converted by the AD converter;
A first distribution data string having a second time interval longer than the first time interval from the input data string stored in the memory, and the first time interval from the first distribution data string. Data distribution for extracting a delayed second distribution data string and a third distribution data string delayed from the first distribution data string by a third time interval different from the first and second time intervals. Means,
A first weighting means for multiplying the first distribution data string by a predetermined window function and performing weighting;
A first fast Fourier transform means for Fourier transforming the weighted first distribution data string and transforming it into first frequency domain data;
A second weighting means for multiplying the second distribution data string by a predetermined window function and performing weighting;
A second fast Fourier transform means for Fourier transforming the weighted second distribution data string and transforming it into second frequency domain data;
A third weighting means for multiplying the third distribution data string by a predetermined window function and performing weighting;
A third fast Fourier transform means for Fourier transforming the weighted third distribution data string and transforming it into third frequency domain data;
Signal detection means for detecting whether a signal exists for each frequency of the first or second or third frequency domain data;
When a signal is detected, the phase difference per unit time of any one of the first, second, and third frequency domain data is compared to determine whether the signal component at the frequency at which the signal is detected is single. And an aliasing determination means for calculating a true frequency only when the signal component is single. A bandpass filter that limits the analog input signal to be frequency analyzed to a certain frequency band;
An AD converter that sequentially converts the band-limited analog input signal into a digital data sequence at a first time interval;
A memory for sequentially storing input data strings converted by the AD converter;
A first distribution data string having a second time interval longer than the first time interval from the input data string stored in the memory, and the first time interval from the first distribution data string. Data distribution for extracting a delayed second distribution data string and a third distribution data string delayed from the first distribution data string by a third time interval different from the first and second time intervals. Means,
A first weighting means for multiplying the first distribution data string by a predetermined window function and performing weighting;
A first fast Fourier transform means for Fourier transforming the weighted first distribution data string and transforming it into first frequency domain data;
A second weighting means for multiplying the second distribution data string by a predetermined window function and performing weighting;
A second fast Fourier transform means for Fourier transforming the weighted second distribution data string and transforming it into second frequency domain data;
A third weighting means for multiplying the third distribution data string by a predetermined window function and performing weighting;
A third fast Fourier transform means for Fourier transforming the weighted third distribution data string and transforming it into third frequency domain data;
Signal detection means for detecting whether a signal exists for each frequency of the first or second or third frequency domain data;
When a signal is detected, a matrix having frequency domain data at the first, second, and third frequencies as elements, and an integer multiple of the reciprocal of the second time interval from the frequency at which the signal is detected. A frequency analyzer comprising: an aliasing determination unit that calculates a true frequency by multiplying an inverse matrix of a matrix having a plurality of shifted frequencies as elements. A bandpass filter that limits the analog input signal to be frequency analyzed to a certain frequency band;
An AD converter that sequentially converts the band-limited analog input signal into digital data a predetermined number of times in a first time interval;
A memory for sequentially storing the input data converted by the AD converter;
A first distribution data string obtained by extracting the data having the first time interval from the input data stored in the memory, and a number obtained by multiplying the first time interval by the predetermined number of digital conversions. A data distribution means for extracting a second distribution data string delayed from the first distribution data string by a second time interval shorter than one half of
A first weighting means for multiplying the first distribution data string by a predetermined window function and performing weighting;
A first fast Fourier transform means for Fourier transforming the weighted first distribution data string and transforming it into first frequency domain data;
A second weighting means for multiplying the second distribution data string by a predetermined window function and performing weighting;
A second fast Fourier transform means for Fourier transforming the weighted second distribution data string and transforming it into second frequency domain data;
Signal detection means for detecting whether a signal exists for each frequency of the first or second frequency domain data;
When a signal is detected, the phase of the first and second frequency domain data is compared, the true frequency is calculated, and the process of adding the signal intensity at the calculated frequency is performed. A frequency analyzer. A bandpass filter that limits the analog input signal to be frequency analyzed to a certain frequency band;
An AD converter that sequentially converts the band-limited analog input signal into a digital data sequence a predetermined number of times in a first time interval;
A memory for sequentially storing input data strings converted by the AD converter;
A first distribution data string obtained by extracting the data having the first time interval from the input data stored in the memory, and a number obtained by multiplying the first time interval by the predetermined number of digital conversions. A second distribution data sequence that is delayed from the first distribution data sequence by a second time interval shorter than one half of the first time interval, and a half of the number obtained by multiplying the first time interval by the predetermined number of times of digital conversion. Data distribution means for extracting a third distribution data sequence that is shorter than the first and second time intervals different from the first and second time intervals by a third distribution data sequence that is delayed from the first distribution data sequence;
A first weighting means for multiplying the first distribution data string by a predetermined window function and performing weighting;
A first fast Fourier transform means for Fourier transforming the weighted first distribution data string and transforming it into first frequency domain data;
A second weighting means for multiplying the second distribution data string by a predetermined window function and performing weighting;
A second fast Fourier transform means for Fourier transforming the weighted second distribution data string and transforming it into second frequency domain data;
A third weighting means for multiplying the third distribution data string by a predetermined window function and performing weighting;
A third fast Fourier transform means for Fourier transforming the weighted third distribution data string and transforming it into third frequency domain data;
Signal detection means for detecting whether a signal exists for each frequency of the first or second or third frequency domain data;
When a signal is detected, the phase difference per unit time of any one of the first, second, and third frequency domain data is compared to determine whether the signal component at the frequency at which the signal is detected is single. Only when there is a single signal component, the phase of the first and second frequency domain data is compared, the true frequency is calculated, and the signal strength is added at the calculated frequency. A frequency analysis apparatus comprising: a convolution effect determination means for performing

Images