(57)【特許請求の範囲】
【請求項１】
周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、
上記帯域制限されたアナログ入力信号を第一の時間間隔で順次デジタルデータ列に変換するＡＤ変換器と、
上記ＡＤ変換器によって変換された入力データ列を順次記憶するメモリと、
上記メモリに記憶された入力データ列から上記第一の時間間隔よりも長い第二の時間間隔を持った第一の分配データ列、およびこの第一の分配データ列から上記第一の時間間隔ずつ遅れた第二の分配データ列とを抽出するデータ分配手段と、
上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、
この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、
上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、
この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、
上記第一または第二の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、
信号が検出された場合、上記第一および第二の周波数領域データの位相を比較して真の周波数を算出するエイリアジング判定手段とを備えることを特徴とする周波数分析装置。
【請求項２】
周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、
上記帯域制限されたアナログ入力信号を２つに分配する信号分配器と、
上記分配された一方のアナログ入力信号を一定の時間間隔でデジタル信号に変換し、第一の分配データ列を生成する第一のＡ／Ｄ変換器と、
上記分配された他方のアナログ入力信号を、上記第一のＡ／Ｄ変換器のデジタル変換より上記一定の時間間隔より短い時間間隔だけ遅れて、一定の時間間隔でデジタル信号に変換し、第二の分配データ列を生成する第二のＡ／Ｄ変換器と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、
この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、
上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、
この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、
上記第一または第二の周波数領域データを各周波数ごとに信号が存在するかを検出する信号検出手段と、
信号が検出された場合、上記第一および第二の周波数領域データを比較することで、エイリアジングの影響を排除し、真の周波数を算出するエイリアジング判定手段とを備えることを特徴とする周波数分析装置。
【請求項３】
周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、
上記帯域制限されたアナログ入力信号を第一の時間間隔で順次デジタルデータ列に変換するＡＤ変換器と、
上記ＡＤ変換器によって変換された入力データ列を順次記憶するメモリと、
上記メモリに記憶された入力データ列から上記第一の時間間隔よりも長い第二の時間間隔を持った第一の分配データ列と、この第一の分配データ列から上記第一の時間間隔ずつ遅れた第二の分配データ列と、この第一の分配データ列から上記第一および第二の時間間隔とは異なる第三の時間間隔ずつ遅れた第三の分配データ列とを抽出するデータ分配手段と、
上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、
この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、
上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、
この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、
上記第三の分配データ列に所定の窓関数を乗算し、重み付けを行う第三の重み付け手段と、
この重み付けされた第三の分配データ列をフーリエ変換し、第三の周波数領域データに変換する第三の高速フーリエ変換手段と、
上記第一または第二または第三の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、
信号が検出された場合、上記第一、第二、第三の周波数領域データのいずれか２つの単位時間当たりの位相差を比較して、信号が検出された周波数における信号成分が単一か否かを判定し、信号成分が単一の場合のみ真の周波数を算出するエイリアジング判定手段とを備えることを特徴とする周波数分析装置。
【請求項４】
周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、
上記帯域制限されたアナログ入力信号を第一の時間間隔で順次デジタルデータ列に変換するＡＤ変換器と、
上記ＡＤ変換器によって変換された入力データ列を順次記憶するメモリと、
上記メモリに記憶された入力データ列から上記第一の時間間隔よりも長い第二の時間間隔を持った第一の分配データ列と、この第一の分配データ列から上記第一の時間間隔ずつ遅れた第二の分配データ列と、この第一の分配データ列から上記第一および第二の時間間隔とは異なる第三の時間間隔ずつ遅れた第三の分配データ列とを抽出するデータ分配手段と、
上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、
この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、
上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、
この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、
上記第三の分配データ列に所定の窓関数を乗算し、重み付けを行う第三の重み付け手段と、
この重み付けされた第三の分配データ列をフーリエ変換し、第三の周波数領域データに変換する第三の高速フーリエ変換手段と、
上記第一または第二または第三の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、
信号が検出された場合、上記第一、第二、第三の各周波数における周波数領域データを要素とする行列に、当該信号が検出された周波数から上記第二の時間間隔の逆数の整数倍ずつずれた複数の周波数を要素とする行列の逆行列を乗算することで真の周波数を算出するエイリアジング判定手段とを備えることを特徴とする周波数分析装置。
【請求項５】
周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、
上記帯域制限されたアナログ入力信号を第一の時間間隔であらかじめ決められた一定回数だけ順次デジタルデータに変換するＡＤ変換器と、
上記ＡＤ変換器によって変換された入力データを順次記憶するメモリと、
上記メモリに記憶された入力データから上記第一の時間間隔を持ったデータを抽出してなる第一の分配データ列、および上記第一の時間間隔に上記デジタル変換した一定の回数を掛けた数の二分の一より短い第二の時間間隔ずつ上記第一の分配データ列より遅れた第二の分配データ列とを抽出するデータ分配手段と、
上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、
この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、
上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、
この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、
上記第一または第二の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、
信号が検出された場合、上記第一および第二の周波数領域データの位相を比較して、真の周波数を算出し、算出した周波数において当該信号強度を足し合せていく処理を行う畳込み効果判定手段とを備えることを特徴とする周波数分析装置。
【請求項６】
周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、
上記帯域制限されたアナログ入力信号を第一の時間間隔であらかじめ定められた一定回数だけ順次デジタルデータ列に変換するＡＤ変換器と、
上記ＡＤ変換器によって変換された入力データ列を順次記憶するメモリと、
上記メモリに記憶された入力データから上記第一の時間間隔を持ったデータを抽出してなる第一の分配データ列、および上記第一の時間間隔に上記デジタル変換した一定の回数を掛けた数の二分の一より短い第二の時間間隔ずつ上記第一の分配データ列より遅れた第二の分配データ列と、上記第一の時間間隔に上記デジタル変換した一定の回数を掛けた数の二分の一より短かく、かつ上記第一および第二の時間間隔とは異なる第三の時間間隔ずつ上記第一の分配データ列より遅れた第三の分配データ列とを抽出するデータ分配手段と、
上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、
この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、
上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、
この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、
上記第三の分配データ列に所定の窓関数を乗算し、重み付けを行う第三の重み付け手段と、
この重み付けされた第三の分配データ列をフーリエ変換し、第三の周波数領域データに変換する第三の高速フーリエ変換手段と、
上記第一または第二または第三の周波数領域データいずれかの各周波数ごとに信号が存在するかを検出する信号検出手段と、
信号が検出された場合、上記第一、第二、第三の周波数領域データのいずれか２つの単位時間当たりの位相差を比較して、信号が検出された周波数における信号成分が単一か否かを判定し、信号成分が単一の場合のみ上記第一および第二の周波数領域データの位相を比較して、真の周波数を算出し、算出した周波数において当該信号強度を足し合せていく処理を行う畳込み効果判定手段とを備えることを特徴とする周波数分析装置。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、受信した信号の周波数分布を短時間フーリエ変換によって分析する周波数分析装置に関するものである。
【０００２】
【従来の技術】
レーダなどの無線信号を受信し、その信号の周波数成分を短時間フーリエ変換によって分析する周波数分析装置においては、レーダ送信波の周波数ホッピング等に対応するため、分析する周波数帯域が広いことが望まれる。
【０００３】
ところで、一般的に時間間隔τでＡＤ変換した信号は、元々の信号の周波数ｆのほかに、ｆ＋ｒ／τ（ｒは整数）の信号成分を持つことが知られている。この成分をエイリアジングと呼ぶ。このエイリアジングの影響を避けながら周波数帯域を広帯域化するためには、帯域の逆数に比例してデータを細かい時間間隔で変換した信号を扱う必要がある。これは即ち、扱うデータ量が周波数帯域に比例して多くのデータに対する処理を行う必要があることを表す。例えば、周波数帯域を１６倍にする場合、ＡＤ変換される信号の時間間隔も１／１６細かくする必要があり、このため扱うデータ数は１６倍になる。短時間フーリエ変換を効率的に行うＦＦＴアルゴリズムにおいては、データ数がＺであるとすると、そのデータ数Ｚに対して計算量はＺｌｏｇ２Ｚに比例して増大することが知られており、データ数ＺをＹ倍にすると、その計算量は（ＺＹｌｏｇ２ＺＹ）／（ＺｌｏｇＺ）＝Ｙ×（１＋（ｌｏｇ２Ｙ／ｌｏｇ２Ｚ））となる。一例として、Ｚ＝１２８、Ｙ＝１６の場合、約２５倍の計算量が要求され、この処理負荷の増大が広帯域信号の分析を困難にしている。
【０００４】
また、短時間フーリエ変換には、窓関数による畳込み演算によって真の周波数成分が拡散してしまい精度良く周波数の算出ができない等の問題がある。
【０００５】
以下に、短時間フーリエ変換を用いた従来の周波数分析装置について説明し、まずエイリアジングの発生による問題点を明らかにする。図１９は従来の周波数分析装置の構成を示している。図において、１は受信される電波信号から、監視する対象の周波数成分のみを選択する帯域フィルタ、２はアナログ信号を一定の時間間隔でサンプリングし、デジタル信号に変換するＡＤ変換器、３はデジタル信号に変換された入力信号を順次記録するメモリ、４はメモリ３に記録された信号を一定の長さずつ読み出し、窓関数を乗じて重み付けを行う重み付け手段、５は重み付けされた一定の長さの信号に順次フーリエ変換を行い、周波数分析を行うフーリエ変換手段、６は分析結果を出力する表示装置である。
【０００６】
次に、動作について説明する。図１９に示した周波数分析装置は、短時間フーリエ変換の原理を使用して受信する信号の周波数分布を算出するものである。短時間フーリエ変換については、貴家仁志「マルチレート信号処理」（昭晃堂）ｐｐ．１４９〜１５３等の文献で説明されているが、次のようなものである。
【０００７】
図２０に示したように、短時間フーリエ変換は、信号に有限な長さの窓関数を乗じて切り出し、それをフーリエ変換するという処理技法である。図２０（ａ）はサンプリングにより得られた実際の信号列であり、（ｂ）が窓関数の形状、（ｃ）が信号列と窓関数との畳込み演算により得られる重み付け後の信号であり、（ｄ）はそれぞれをフーリエ変換して周波数成分に直した図であり、（ｅ）は各時間における周波数成分を並べた時間変化を表す図である。窓関数を時間に従いシフトし、処理を繰り返すことで時間と周波数の情報を調べる。
【０００８】
短時間フーリエ変換は算術的には以下のように表現できる。
クロック周期τでサンプリングされた信号をｓ（ｎ）、窓関数をｇ（ｎ）（ｇ（ｎ）はｎ＝０，τ，２τ，…，（Ｎ−１）τで０でない）とする。ｓ（ｎ）は図２０（ａ）に、ｇ（ｎ）は図２０（ｂ）に示したような波形であるとする。ここで、ｇ（ｎ）をＮτずつシフトした関数ｇ（ｎ−ｋＮ）を考え、このｇ（ｎ−ｋＮ）とｓ（ｎ）を乗ずることで式（１）のようにｋブロック目の信号ｓｋ（ｎ）を抽出する。
ｓｋ（ｎ）＝ｇ（ｎ−ｋＮ）ｓ（ｎ）                  （１）
である。図２０（ｃ）にｓｋ（ｎ）を示す。
【０００９】
（１）式のｓｋ（ｎ）をフーリエ変換する処理が、短時間フーリエ変換である。従って、ｓＮ（ｎ）の短時間フーリエ変換は、
【数１】
0003646702.tif 000002
である。ただし、ω＝２πｐ／Ｎτである。これより、
【数２】
0003646702.tif 000003
である。（３）式は入力信号ｓ（ｎ）について、周波数ｐ／Ｎτにおける周波数成分を表すものである。
【００１０】
ｇ（ｎ）は窓関数と呼ばれるもので、例えば、次のようなハニング窓関数などが一般的に用いられる。
ｇ（ｎ）＝１／２｛１＋ｃｏｓ２π／Ｎ（ｎ−Ｎ／２）｝
ただし、ｎ＝０〜Ｎ−１、それ以外のｎではｇ（ｎ）＝０
【００１１】
図１９は、以上の原理を実現する短時間フーリエ変換を用いた周波数分析装置の構成である。例えば、無線信号が受信されると、適当な中間周波数信号に変換され帯域フィルタ１に入力される。帯域フィルタ１では入力信号について、分析対象とする周波数帯の信号のみを取り出し、Ａ／Ｄ変換器２ではこの取り出した信号について適当な時間間隔τでサンプリングして、デジタル信号に変換する。このデジタル変換した信号列が上述のｓ（ｎ）であり、ｓ（ｎ）は順次メモリ３に記憶される。なお、以下、Ａ／Ｄ変換については、ＩおよびＱチャネルを使用するなどにより、信号が複素数として取得されることを想定するが、これは、信号が実数としてサンプリングされる場合に比べると、サンプリング間隔が２倍で良い他は本質的な違いはない。
【００１２】
重み付け手段４では、メモリ２に記録された信号を順次読み出し、信号列ｓ（ｎ）に窓関数ｇ（ｎ）を掛け合わせ（１）式に表されたような演算を行い、高速フーリエ変換手段５では、重み付けの行われた信号ｓｋ（ｎ）に対し、順次（３）式の処理を行い、各時点での信号に含まれる周波数成分を分析し、その結果は表示装置６に表示する。
【００１３】
以上のように、従来の周波数分析装置は、短時間フーリエ変換を使用して受信される電波信号の周波数分布を求めている。ここで、短時間フーリエ変換の性質を考えると、（３）式の性質から以下の特徴が導き出せる。
【００１４】
短時間フーリエ変換を行った結果、周波数ｆ＝ｐ／Ｎτ（ｐは整数）となり、離散的な値を持つ。その内、最も高い周波数はｐ＝Ｎ−１のときの（Ｎ−１）／Ｎτである。即ち、サンプリング定理により、エイリアジングの発生しない周波数範囲は１／τである。即ち、分析しようとする最大周波数をｆｗとすると、ｆｗ＜１／τならば、周波数の特定が可能であることを意味する。周波数１／τはナイキスト周波数ともいう。
【００１５】
また、（３）式の右辺は、
【数３】
0003646702.tif 000004
のような関係が成り立つ。ここでｒは整数である。これは、真の周波数ｐ／Ｎτからｒ／τだけずれた周波数に同一の値を持った成分が現れることを意味する。即ち、図２１に示したように、本来の信号の周波数の他に、そこから１／τごとにエイリアジングと呼ばれる周波数成分が発生する。
【００１６】
帯域フィルタ１によって分析する周波数成分の範囲を１／τ以下に制限すればこのようなエイリアジングの発生を防ぐことができる。しかし、上述のように周波数分析装置には分析周波数の広帯域化が望まれるが、これにより、帯域フィルタの通過させる周波数範囲を１／τ（τはサンプリング間隔）より広帯域化すれば、それだけ多くのエイリアジングが生じ、信号の真の周波数と分離することができない。以上が、エイリアジング発生による問題点である。
【００１７】
次に、上述したように、有限の時間範囲で窓関数によって畳込み演算を行うことで原理的にスペクトルが拡散する誤差について説明する。即ち、ｇ（ｎ）のフーリエ変換をＧ（２πｐ／Ｎτ）とすると、
0003646702.tif 000005
また、（１）式および（２）式から、ｓ（ｎ）のフーリエ変換をＳ（ω）とすると、
【数４】
0003646702.tif 000006
ここで、便宜上ｍ＝ｎ−ｋＮのように置換すると、
【数５】
0003646702.tif 000007
従って、Ｇ（ω）が複数の値を持つことから、短時間フーリエ変換で計算される周波数分布は、真の周波数の周囲に、窓関数の性質にしたがって拡散する。この拡散の様子を示したものを図２２に示す。図２２（ａ）は信号の元の周波数分布であり、（ｂ）は窓関数によって拡散した周波数分布を表す。拡散した各成分の強度の和は、（８）式の畳込み演算の性質から元の信号の強度に一致する。また、その拡散の程度は窓関数ｇ（ｎ）の性質によるが、性質の良いハニング窓関数であっても、真の周波数の高低に３つの拡散成分が発生することは免れない。
【００１８】
この拡散成分を抑え真のスペクトルを算出する周波数分析装置の一例として、例えば、特開平６−１６０４４５号公報に示された周波数分析装置がある。図２３に特開平６−１６０４４５号公報の周波数分析装置を示す。図中、１はＡＤ変換器、２はメモリ、３は窓関数重み付け手段、４は高速フーリエ変換手段、５は高速フーリエ変換器４の出力側に設けられたスペクトル補間判定手段、６は表示器、７は補間判定手段５の判定結果によりスペクトル補間を実行するスペクトル補間手段、８はスペクトル補間手段７によって算出したスペクトルの周波数、振幅、位相データを時間軸データに逆変換する逆フーリエ変換器、９は逆フーリエ変換器８で逆フーリエ変換した時間軸データを累積する累積器、１１は累積器９によって累積された累積結果をメモリ２から読み出される時間データから減算し、この減算結果を再び窓関数重み付け手段３に入力する減算器、１２はスペクトル補間手段７によって算出した周波数領域データを累積する累積器、１３はこの累積器１２に累積した周波数領域データを減算器１１で減算した残差データを再度高速フーリエ変換して得られた周波数領域データＳ（ｔ）に加算する加算器である。
【００１９】
この周波数分析装置では、フーリエ変換して得られた離散周波数スペクトルの中で３本の周波数スペクトルが互いに隣接して存在し、しかもその上位と下側のスペクトルの位相がほぼ同じであったら、この３本のスペクトルは真のスペクトルが分散されているものと判定し、補間処理によりその真のスペクトルの周波数と振幅、位相を求める。
【００２０】
この真のスペクトルデータを累積器１２に累積するとともに、逆フーリエ変換して時間軸データに逆変換し、この時間軸データを累積器９で累積し、その累積データを入力データ列から減算するから、補間処理により真のスペクトル線を算出したあとでは残差エネルギーが低下し、線スペクトルは広がりを持つことはなく、線スペクトルとして表示される。この結果、スペクトル分解能が向上する。また、補間処理により振幅、位相も真値近くに近似して算出するから振幅値および位相も正確に測定することができる。
【００２１】
しかし、この場合、フィードバックループにより累積したデータを繰り返し減算して真の周波数成分を取り出していく構成であるため、値が収束するまでに何回もサンプリングを行う必要があり、時間がかかり、また、フィードバックループを無くした場合は、微少信号を抽出できない。
【００２２】
また、上記の例では真の周波数の周辺の周波数を正しく３つ選択できることが必要だが、真の周波数が、離散周波数値のほぼ中間に存在した場合、例えば、ｎ／Ｎτと（ｎ＋１）／Ｎτの中間の場合、真の周波数に最も近い周波数を（ｎ＋１）／Ｎτと誤判定しやすく、真の周波数の周辺の３つの周波数を間違えてｎ／Ｎτ、（ｎ＋１）／Ｎτ、（ｎ＋２）／Ｎτのように選んでしまうことが起こりやすいため、実際の周波数算出は不安定になりやすい。
【００２３】
さらに、この場合、真の周波数の算出方法は特開平６−１６０４４５号公報の（１４）式などのように引き算が主体であり、ｓｉｎやｃｏｓが０の付近では計算結果が不安定であり、周波数算出に誤差が出やすい。
【００２４】
【発明が解決しようとする課題】
本発明は、以上のような問題点を解決するためになされたもので、短時間フーリエ変換の際に、受信信号を広帯域化するに伴う計算量の増大を抑え、周波数分布を効率よく正確に算出することを目的とする。
【００２５】
【課題を解決するための手段】
本発明に係わる周波数分析装置は、周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、上記帯域制限されたアナログ入力信号を第一の時間間隔で順次デジタルデータ列に変換するＡＤ変換器と、上記ＡＤ変換器によって変換された入力データ列を順次記憶するメモリと、上記メモリに記憶された入力データ列から上記第一の時間間隔よりも長い第二の時間間隔を持った第一の分配データ列、およびこの第一の分配データ列から上記第一の時間間隔ずつ遅れた第二の分配データ列とを抽出するデータ分配手段と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、上記第一または第二の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、信号が検出された場合、上記第一および第二の周波数領域データの位相を比較して真の周波数を算出するエイリアジング判定手段とを備えるものである。
【００２６】
また、周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、上記帯域制限されたアナログ入力信号を２つに分配する信号分配器と、上記分配された一方のアナログ入力信号を一定の時間間隔でデジタル信号に変換し、第一の分配データ列を生成する第一のＡ／Ｄ変換器と、上記分配された他方のアナログ入力信号を、上記第一のＡ／Ｄ変換器のデジタル変換より上記一定の時間間隔より短い時間間隔だけ遅れて、一定の時間間隔でデジタル信号に変換し、第二の分配データ列を生成する第二のＡ／Ｄ変換器と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、上記第一または第二の周波数領域データを各周波数ごとに信号が存在するかを検出する信号検出手段と、信号が検出された場合、上記第一および第二の周波数領域データを比較することで、エイリアジングの影響を排除し、真の周波数を算出するエイリアジング判定手段とを備えるものである。
【００２７】
また、周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、上記帯域制限されたアナログ入力信号を第一の時間間隔で順次デジタルデータ列に変換するＡＤ変換器と、上記ＡＤ変換器によって変換された入力データ列を順次記憶するメモリと、上記メモリに記憶された入力データ列から上記第一の時間間隔よりも長い第二の時間間隔を持った第一の分配データ列と、この第一の分配データ列から上記第一の時間間隔ずつ遅れた第二の分配データ列と、この第一の分配データ列から上記第一および第二の時間間隔とは異なる第三の時間間隔ずつ遅れた第三の分配データ列とを抽出するデータ分配手段と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、上記第三の分配データ列に所定の窓関数を乗算し、重み付けを行う第三の重み付け手段と、この重み付けされた第三の分配データ列をフーリエ変換し、第三の周波数領域データに変換する第三の高速フーリエ変換手段と、上記第一または第二または第三の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、信号が検出された場合、上記第一、第二、第三の周波数領域データのいずれか２つの単位時間当たりの位相差を比較して、信号が検出された周波数における信号成分が単一か否かを判定し、信号成分が単一の場合のみ真の周波数を算出するエイリアジング判定手段とを備えるものである。
【００２８】
また、周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、上記帯域制限されたアナログ入力信号を第一の時間間隔で順次デジタルデータ列に変換するＡＤ変換器と、上記ＡＤ変換器によって変換された入力データ列を順次記憶するメモリと、上記メモリに記憶された入力データ列から上記第一の時間間隔よりも長い第二の時間間隔を持った第一の分配データ列と、この第一の分配データ列から上記第一の時間間隔ずつ遅れた第二の分配データ列と、この第一の分配データ列から上記第一および第二の時間間隔とは異なる第三の時間間隔ずつ遅れた第三の分配データ列とを抽出するデータ分配手段と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、上記第三の分配データ列に所定の窓関数を乗算し、重み付けを行う第三の重み付け手段と、この重み付けされた第三の分配データ列をフーリエ変換し、第三の周波数領域データに変換する第三の高速フーリエ変換手段と、上記第一または第二または第三の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、信号が検出された場合、上記第一、第二、第三の各周波数における周波数領域データを要素とする行列に、当該信号が検出された周波数から上記第二の時間間隔の逆数の整数倍ずつずれた複数の周波数を要素とする行列の逆行列を乗算することで真の周波数を算出するエイリアジング判定手段とを備えるものである。
【００２９】
また、周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、上記帯域制限されたアナログ入力信号を第一の時間間隔であらかじめ決められた一定回数だけ順次デジタルデータに変換するＡＤ変換器と、上記ＡＤ変換器によって変換された入力データを順次記憶するメモリと、上記メモリに記憶された入力データから上記第一の時間間隔を持ったデータを抽出してなる第一の分配データ列、および上記第一の時間間隔に上記デジタル変換した一定の回数を掛けた数の二分の一より短い第二の時間間隔ずつ上記第一の分配データ列より遅れた第二の分配データ列とを抽出するデータ分配手段と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、上記第一または第二の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、信号が検出された場合、上記第一および第二の周波数領域データの位相を比較して、この位相が同じの場合、真の周波数を算出し、算出した周波数において当該信号強度を足し合せていく処理を行う畳込み効果判定手段とを備えるものである。
【００３０】
また、周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、上記帯域制限されたアナログ入力信号を第一の時間間隔であらかじめ定められた一定回数だけ順次デジタルデータ列に変換するＡＤ変換器と、上記ＡＤ変換器によって変換された入力データ列を順次記憶するメモリと、上記メモリに記憶された入力データから上記第一の時間間隔を持ったデータを抽出してなる第一の分配データ列、および上記第一の時間間隔に上記デジタル変換した一定の回数を掛けた数の二分の一より短い第二の時間間隔ずつ上記第一の分配データ列より遅れた第二の分配データ列と、上記第一の時間間隔に上記デジタル変換した一定の回数を掛けた数の二分の一より短かく、かつ上記第一および第二の時間間隔とは異なる第三の時間間隔ずつ上記第一の分配データ列より遅れた第三の分配データ列とを抽出するデータ分配手段と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、上記第三の分配データ列に所定の窓関数を乗算し、重み付けを行う第三の重み付け手段と、この重み付けされた第三の分配データ列をフーリエ変換し、第三の周波数領域データに変換する第三の高速フーリエ変換手段と、上記第一または第二または第三の周波数領域データいずれかの各周波数ごとに信号が存在するかを検出する信号検出手段と、信号が検出された場合、上記第一、第二、第三の周波数領域データのいずれか２つの単位時間当たりの位相差を比較して、信号が検出された周波数における信号成分が単一か否かを判定し、信号成分が単一の場合のみ上記第一および第二の周波数領域データの位相を比較して、この位相が同じの場合、真の周波数を算出し、算出した周波数において当該信号強度を足し合せていく処理を行う畳込み効果判定手段とを備えるものである。
【００３１】
【発明の実施の形態】
実施の形態１．
本発明の実施の形態１について説明する。本実施の形態１に係わる周波数分析装置は、Ａ／Ｄ変換された信号から微少タイミングずらしてサンプルして得られた２系列の信号列をそれぞれ短時間フーリエ変換処理して局所的な周波数分布を求めるものである。また、その際生じるエイリアジングを微少タイミングずらしたことにより生じる信号列間の位相差を検出することで真の周波数を特定するものである。
【００３２】
以下、詳細について説明する。図１は本実施の形態１に係わる周波数分析装置の構成を表すブロック図である。図１において、図１９と同じ構成要素には同じ符号を付す。７はメモリ３内に記録した信号のデータからサンプリングタイミングのわずかに異なる２つのデータ列を抽出、分配するデータ分配手段、８はフーリエ変換結果から各周波数に関して信号の存在を判定する信号検出手段、９は信号検出手段８で判定された信号の存在する周波数に関して、２つの短時間フーリエ変換の結果を比較してエイリアジングの影響を排除し、正しい周波数を算出するエイリアジング判定手段である。
【００３３】
図２は本実施の形態１に係わる周波数分析装置のサンプリング方法を表す図である。図２に示した通り、時間間隔Δτでサンプリングしたデータの中から、適当な時間間隔τの信号列ｓ（ｎ）と、当該ｓ（ｎ）の各信号からΔτだけ遅れた信号列ｓΔ（ｎ）を考える。信号の周波数をｆ（＝ｐ／Ｎτ，ｐは整数）とすると、信号列ｓ（ｎ）とｓΔ（ｎ）とは、位相遅れ２πｆΔτが発生するため、
【数６】
0003646702.tif 000008
の関係が成立する。ここで、係数ｅｘｐ（２πｊｆΔτ）は定数であるため、重み付け、フーリエ変換とあった短時間フーリエ変換の各プロセスによって変化することはなく、信号列ｓΔ（ｎ）は、ｓ（ｎ）に対して一様な周波数ｆによる位相差を持つ。このため、ｓΔ（ｎ）のフーリエ変換をＳＴＦＴΔ（２πｐ／Ｎτ）とすると、
【数７】
0003646702.tif 000009
（１０）式より、ＳＴＦＴ（２πｐ／Ｎτ）≠０となる周波数ｆに対して、元のデータ列ｓ（ｎ）と、微少時間Δτだけずらしてサンプリングを行ったデータ列ｓΔ（ｎ）のそれぞれの短時間フーリエ変換を比較すると、係数ｅｘｐ（２πｊｆΔτ）の違いとなる。従って、（１０）式を変形すると、
【数８】
0003646702.tif 000010
であり、（１１）式右辺の各ＦＦＴ成分については、それぞれ間隔τのデータを使用しているため、ｐに対する周波数ｆは、
ｆ＝ｐ／Ｎτ，（ｐ＋Ｎ）／Ｎτ，（ｐ＋２Ｎ）／Ｎτ…
のいずれかである。従って、ｒを整数としてｆ＝（ｐ＋ｒＮ）／Ｎτのように表現でき、（１１）式は、
【数９】
0003646702.tif 000011
である。その際、−２πｊｒＮΔτ／Ｎτの項によって位相が一回転しない場合は、ｒを一意に求めることができる。即ち、
【数１０】
0003646702.tif 000012
を満たす場合である。（１２）式を変形すると、
１／Δτ＞ｒＭ／τ                                        （１４）の関係が得られる。１／Δτは本実施の形態１の周波数分析装置の特定可能な周波数の範囲である。ｒＭが余りに大きい場合には（１４）式を満たすことができなくなるため、帯域フィルタ１で適当な周波数帯域を選択する。
【００３４】
図３に（１３）式の関係を表すグラフを示す。図３は縦軸を振幅、横軸を周波数とした周波数分析結果の図である。周波数の値は離散値でその最小単位は１／Ｎτである。また、１／Δτがこの周波数分析装置で分析できる最大の周波数であり、１／Δτより上の周波数が入力されないように帯域フィルタ１で帯域制限する。一方、サンプリング間隔τの信号列ｓ（ｎ）を短時間フーリエ変換するのであるから、エイリアジングの発生しない周波数範囲は１／τである関係上、真の信号周波数から１／τずつずれたところにエイリアジングが発生する。
【００３５】
図では例えば、周波数（ｐ＋２Ｎ）／Ｎτに真の信号成分がある場合である。信号列ｓ（ｎ）およびｓΔ（ｎ）から見るとエイリアジングの発生しない周波数範囲は１／τであるのに対し、帯域フィルタ１から入力される受信信号はその周波数範囲を超えているので、（１１）式の右辺を実測によって求め、そこから左辺のｐおよびｒを求める。
【００３６】
以上が、本実施の形態１に係わる周波数分析装置の原理である。図１の周波数分析装置は以上の原理に基いて短時間フーリエ変換を行うものである。受信した電波を適当な中間周波数に変換した受信信号が、帯域フィルタ１に入力され、帯域制限される。Ａ／Ｄ変換器では、選択された信号をナイキスト周波数の逆数以下である適当な短い時間間隔Δτでサンプリングし、メモリ３に順次蓄積する。データ分配手段７では、メモリに蓄積した信号からそれぞれ適当な間隔τの２系列の信号列ｓ（ｎ）およびｓΔ（ｎ）を選択し、それぞれ重み付け手段４ａおよび４ｂに分配する。重み付け手段４ａおよび４ｂにおいて窓関数による重み付けがなされた後、フーリエ変換手段５ａおよび５ｂに入力される。フーリエ変換手段５ａおよび５ｂでは、重み付けされたｓ（ｎ）およびｓΔ（ｎ）のフーリエ変換を行い、（２）式のように、ＳＴＦＴ（ω）およびＳＴＦＴΔ（ω）を求める。
【００３７】
信号検出手段８ではフーリエ変換手段５ａの出力ＳＴＦＴ（ω）の結果を調べ、信号の存在する周波数を判定し、それらの周波数をエイリアジング判定手段９に通知する。エイリアジング判定手段９では、ＳＴＦＴ（ω）およびＳＴＦＴΔ（ω）を比較し、（１０）式の左辺ｅｘｐ（２πｊｆΔτ）を算出する。
【００３８】
信号検出手段８およびエイリアジング判定手段９で行う処理を図４に示したフローチャートにしたがって説明する。まず、信号検出手段８で離散値の各周波数ｆ＝（ｐ＋ｒＮ）／Ｎτをｐ＝０（ステップ１００）から順番に信号が存在するかの判定が行われる（ステップ１０１）。そして、信号が未検出ならば、そのときのｐ（今の場合、ｐ＝０）に対して、全てのｒについて受信信号の周波数成分が０である（ステップ１０２）。即ち、ｐ＝０で信号が未検出であれば、そのエイリアジングの可能性がある周波数ｆ＝Ｎ／Ｎτ，２Ｎ／Ｎτ，３Ｎ／Ｎτ，…における振幅＝０ということである。
【００３９】
一方、ｐ＝０で信号が検出された場合は、ＳＴＦＴ（ω）およびＳＴＦＴΔ（ω）から（１１）式を満足するｒを算出する（ステップ１０３）。そして、算出された値以外のｒについての受信信号の周波数成分は０である（ステップ１０４）。例えば、ｒ＝３と算出された場合は、それ以外の周波数ｆ＝（３＋Ｎ）／Ｎτ，（３＋２Ｎ）／Ｎτ，（３＋４Ｎ）／Ｎτ，…における振幅＝０ということである。
【００４０】
そして、ｐ＝０のとき、信号が検出、または未検出いずれかの場合の処理が終了すると、次の周波数、即ちｐ＝１（ステップ１０５）の場合について、同様の処理を行う。その際、ｐ≧Ｎであれば、演算する部分が重なってしまうため、ｐ≧Ｎであれば分析を終了し、ｐ＜Ｎの場合のみステップ１０１に戻る（ステップ１０６）。
【００４１】
このように、ｐ＝０からＮまで、即ち離散値の周波数ｆを０から１／τまで信号の有無を検出するだけで、実際は図３の０からｒＭ／τまでの信号帯域を分析することになるので、全ての周波数を分析する必要がなく、計算量が大幅に軽減することができる。
【００４２】
例えば、サンプル数Ｎ＝１２８の場合に、周波数帯域を１６倍に広げることを考えると、
（１）従来の周波数分析装置で周波数帯域を１６倍に広げる場合
フーリエ変換を行うデータ数が１６倍であるため、計算量は（１２８×１６）ｌｏｇ２（１２８×１６）／１６ｌｏｇ２１２８倍、即ち２５倍になる。
（２）本実施の形態１の周波数分析装置で周波数帯域を１６倍に広げる場合
フーリエ変換を行うデータ数自体は変化がない。そして、２つの信号列ｓ（ｎ）とｓΔ（ｎ）をフーリエ変換するため計算を２回行うためほぼ２倍である。ただし、（１１）式の計算も行う必要があるが、これは高々Ｎ回（今の場合１２８回）であり、フーリエ変換の計算量と比べれば１／１０程度であるため、計算量の増加にはほとんど影響しない。
このように、本実施の形態１の周波数分析装置を用いれば、約１２倍の計算負荷が軽減できる。
【００４３】
以上のように、本実施の形態１では、時間間隔τでサンプリングした信号列ｓ（ｎ）と、この信号列ｓ（ｎ）から微少時間Δτだけずらして時間間隔τでサンプリングした信号列ｓΔ（ｎ）とをそれぞれ短時間フーリエ変換し、比較することでΔτだけずらしたことにより生じる位相差ｅｘｐ（２πｊｆΔτ）を算出する。この位相差によって、受信信号を広帯域化に伴って生じるエイリアジングを利用して（１１）式によりｒを求めることで、帯域フィルタ１から出力される周波数帯域全てを走査しなくても受信信号の存在する周波数を算出することができるため、広帯域化しても計算量の増大を抑え、エイリアジングの効果を排除し、正確な周波数分析結果を得ることができる。
【００４４】
実施の形態２．
上記実施の形態１では、２つの信号列ｓ（ｎ）、ｓΔ（ｎ）を得るために、一律にΔτでサンプリングしたデータ群の中から適当な時間間隔τとなるデータを抽出してデータを分配する例を示したが、例えば、図５のように、タイミングのΔτだけ異なるクロックでサンプリングを行う２つのＡ／Ｄ変換器を使用して実現してもよい。
【００４５】
図６に本実施の形態２に係わる周波数分析装置の装置構成を表すブロック図を示す。図１と同じ構成要素には同じ符号を付す。１０は信号分配器である。本実施の形態２の周波数分析装置は図のように２つのＡ／Ｄ変換器２ａ、２ｂを使用するが、無駄な信号をサンプリングすることがなくなるため、メモリ３ａ、３ｂの容量を小さくて済む。
【００４６】
以下、動作について説明する。帯域フィルタ１は入力された受信信号を帯域制限して分析対象の周波数信号のみ選択する。信号分配器１０では、信号を２つのＡ／Ｄ変換器２ａ、２ｂに分配する。各Ａ／Ｄ変換器２ａ、２ｂではそれぞれ時間間隔τでサンプリングを行うが、例えば、Ａ／Ｄ変換器２ｂの方がＡ／Ｄ変換器２ａよりもΔτだけ遅れたタイミングでサンプリングを行う。これによって、実施の形態１で示したような２つの信号列ｓ（ｎ）およびｓΔ（ｎ）が生成され、順次メモリ３ａ、３ｂに記憶される。以下の動作は上記実施の形態１の周波数分析装置の動作と同様である。それぞれの信号列の短時間フーリエ変換ＳＴＦＴ（ω）、ＳＴＦＴΔ（ω）を比較し、（１１）式からｅｘｐ（２πｊｆΔτ）を算出する。
【００４７】
以上のように、本実施の形態２では、帯域制限した信号を予め分配してＡ／Ｄ変換器に入力するため、無駄な信号をサンプリングすることがなく、メモリの容量が小さくて済む。
【００４８】
実施の形態３．
上記（１１）式はそれぞれのｐについてＳＴＦＴ（ω）に含まれる信号が単一である場合のみ成立する。従って、上記実施の形態１および２では、受信信号の周波数分布が単一もしくは比較的疎な場合のみ分析が可能であった。本実施の形態３では、上記実施の形態１または２の周波数分析装置を用いた場合に互いにエイリアジング成分となるような複数の周波数成分が受信信号に含れていた場合、信号が複数のエイリアジング成分を含むか含まないかの判定を行い、含まないと判定された場合にのみ分析を行う周波数分析装置について説明する。
【００４９】
図７は本実施の形態３の周波数分析装置が対象とする信号の周波数スペクトル図である。点線は実際に測定されるスペクトル、実線は信号の真の周波数成分を表す。図７のように上記実施の形態１および２の周波数分析装置を用いた場合に、互いにエイリアジング成分となるような独立した信号成分ｗ１、ｗ３が受信信号に含まれていた場合、その周波数をｆ１、ｆ３とすると、
【数１１】
0003646702.tif 000013
また、
【数１２】
0003646702.tif 000014
これらの式をそれぞれフーリエ変換すると、
【数１３】
0003646702.tif 000015
また、
【数１４】
0003646702.tif 000016
である。
【００５０】
実際の信号から求まるのはＳＴＦＴΔ（２πｐ／Ｎτ）およびＳＴＦＴ（２πｐ／Ｎτ）の値のみである。従って、上記実施の形態１または２のようにＳＴＦＴΔ（２πｐ／Ｎτ）／ＳＴＦＴ（２πｐ／Ｎτ）を求めても、ｗ１、ｗ３が混ざっているため正しいｆ１、ｆ３の値を求めることができない。
【００５１】
そのために、本実施の形態３では、図８に示したように周波数分析装置を、時間間隔τを持った一つの信号系列ｓ（ｎ）と、そのｓ（ｎ）に対してそれぞれ２つの異なる微少な時間差Δτ１、Δτ２を持った信号系列ｓΔ１（ｎ）、ｓΔ２（ｎ）を生成する構成とし、それぞれ短時間フーリエ変換を行う。図において９ａは信号が単一か否かを判定しながら真の周波数を算出するエイリアジング判定手段である。
【００５２】
ここで、エイリアジングの周波数をｆ１＝ｐ／Ｎτ、ｆ２＝（ｐ＋Ｎ）／Ｎτ、ｆ３＝（ｐ＋２Ｎ）／Ｎτ、…ｆｒ＋１＝（ｐ＋ｒＮ）／Ｎτ、また各エイリアジング周波数における各信号強度をｗ１、ｗ２、ｗ３、…ｗｒとする。すると、信号が複数含まれている場合を一般化して表現すると、
【数１５】
0003646702.tif 000017
【数１６】
0003646702.tif 000018
である。ｓΔ２（ｎ）も（１６）式のΔτ１を、それぞれΔτ２と置き換えたものと同様である。
【００５３】
ここで、例えば、周波数ｆ３に真の信号成分が１つのみあるとすると、ｗ３のみが０以外の値を持ち、かつｗ３以外のｗ１、ｗ２、ｗ４…ｗｒは全て０である。
【００５４】
従って、その場合にｓ（ｎ）、ｓΔ１（ｎ）、ｓΔ２（ｎ）をフーリエ変換すると、それぞれ、以下のようになる。
【数１７】
0003646702.tif 000019
【００５５】
逆に、信号成分が単一でなく、例えば、周波数ｆ３の他に周波数ｆ１に信号成分ｗ１が存在していた場合は、そのフーリエ変換は（１７）、（１８）、（１９）式の代わりに、
【数１８】
0003646702.tif 000020
【００５６】
ここで、各信号列間の単位時間当たりの位相差として以下の量を考える。
【数１９】
0003646702.tif 000021
ただし、（２０）式において∠は位相を表す。
例えば、信号成分が周波数ｆ３に一つの場合、（１７）〜（１９）式の右辺を（２１）式に代入すると、全て２πｆ３が得られるが、信号成分が周波数ｆ３およびｆ１にあるような場合、（１７ａ）〜（１９ａ）式の右辺を（２０）式に代入しても一定の値は得られない。そのため、この（２０）式を用いて信号成分が複数あるか否かを判定する。
【００５７】
図８は本実施の形態３に係わる周波数分析装置の構成を表すブロック図である。図１と同じ構成要素には同じ符号を付す。１つの信号系列ｓ（ｎ）およびそのｓ（ｎ）から時間Δτ１およびΔτ２だけずれた２つの信号系列ｓΔ１（ｎ）およびｓΔ２（ｎ）の計３つの信号系列をサンプリングし、それぞれ短時間フーリエ変換する。そこから、信号検出手段８ｂおよびエイリアジング判定手段９ｂによって、上記（２０）式のような単位時間当たりの位相差を求め、信号が一つか否かを判定しながら真の周波数を求める。
【００５８】
以下、図９に示したフローチャートに従って、信号検出手段８ｂおよびエイリアジング判定手段９ｂの動作を説明する。まず、信号検出手段８ｂで離散値の各周波数ｆ＝（ｐ＋ｒＮ）／Ｎτをｐ＝０（ステップ２００）から順番に信号が存在するかの判定が行われる（ステップ２０１）。そして、信号が未検出ならば、そのときのｐ（今の場合、ｐ＝０）に対して、全てのｒについて受信信号の周波数成分が０である（ステップ２０２）。即ち、上記実施の形態１の場合と同様に、ｐ＝０で信号が未検出であれば、そのエイリアジングの可能性がある周波数ｆ＝Ｎ／Ｎτ，２Ｎ／Ｎτ，３Ｎ／Ｎτ，…についての周波数成分も０ということである。
【００５９】
一方、ｐ＝０で信号が検出された場合は、ＳＴＦＴ（ω）、ＳＴＦＴΔ１（ω）、ＳＴＦＴΔ２（ω）、ＳＴＦＴΔ３（ω）から（２０）式を計算し（ステップ２０４）、全て同じ値を算出するか判定する（ステップ２０５）。もし、（２１）式の計算結果が一致しないならば、上述したように信号成分が複数存在することになるので、当該周波数の成分は不定である（ステップ２０３）。一方、（２０）式の結果が同じ値ならば信号は１つであり、周波数の分析が可能であるので実施の形態１と同様に（１１）式を満足するｒＴを算出する（ステップ２０６）。そして、算出したｒＴ以外のｒに対し、周波数（ｐ＋ｒＮ）／Ｎτの成分が０である（ステップ２０７）。
【００６０】
以上、いずれかの処理が終わったら周波数（ｐ＋ｒＮ）／Ｎτのｐをｐ＝ｐ＋１として、次の周波数成分について同様の処理を行う（ステップ２０８）。その際、ｐ≧Ｎであれば、演算する部分が重なってしまうため、ｐ≧Ｎであれば分析を終了し、ｐ＜Ｎの場合のみステップ２０１に戻る（ステップ２０９）。
【００６１】
このように、ｐ＝０からＮ−１まで、即ち離散値の周波数ｆを０から１／τまで信号の有無を検出するだけで、全ての周波数を分析する必要がなく、計算量が大幅に軽減することができる。
【００６２】
また、信号がいくつあるか判定しながら周波数分析を行うため、信号成分が複数ある場合でも、誤った周波数を算出しない。
【００６３】
実施の形態４．
本実施の形態４では、それぞれ異なる時間だけずれた２つ以上の信号系列を短時間フーリエ変換して、その離散値の各周波数に含まれる周波数成分が１つか否かを判定し、複数の信号が含まれていてもエイリアジングの影響を除くことができる周波数分析装置について説明する。
【００６４】
ここでは上記実施の形態３の図７のように互いのエイリアジング成分となるような周波数に独立した複数の信号がある場合に、３つ信号系列を短時間フーリエ変換する例を示す。
【００６５】
図１０は本実施の形態４の周波数分析装置の構成を表すブロック図である。図において、図８と同じ構成要素には同じ符号を付す。９ｂは本実施の形態３に係わるエイリアジング判定手段である。前記実施の形態３と同様に、信号系列ｓ（ｎ）に対し複数の時間差Δτ１、Δτ２、Δτ３を持った信号系列ｓΔ１（ｎ）、ｓΔ２（ｎ）、ｓΔ３（ｎ）をサンプリングして生成し、短時間フーリエ変換を行う。ただし、ここで、Δτ１＜Δτ２＜Δτ３であって、図７のｒ／τ＜１／Δτ１であるとする。これら４つの信号系列の短時間フーリエ変換後のスペクトルＳＴＦＴ（２πｐ／Ｎτ）、ＳＴＦＴΔ１（２πｐ／Ｎτ）、ＳＴＦＴΔ２（２πｐ／Ｎτ）、ＳＴＦＴΔ３（２πｐ／Ｎτ）それぞれについて（１１）式のような位相差を考える。
【００６６】
ここで、エイリアジングの性質として、真の周波数から１／τごとにエイリアジング成分が生じることから、エイリアジングの生じる可能性のある周波数は、
【数２０】
0003646702.tif 000022
のように表現でき、この各周波数における振幅をｗ１、ｗ２、…ｗｒ＋１とすれば、４つの信号系列は以下のような行列にて表現できる。
【数２１】
0003646702.tif 000023
として、
Ｘ＝Ｓ・Ａ                                            （２２）
である。
【００６７】
ここで、行列Ｘの各要素については４つの信号系列を実際に短時間フーリエ変換することで、その値が得られ、また、行列Ｓの各要素についても、上述のようにエイリアジング成分の性質から各周波数が定められるため、（２１）各式中のｐを固定すればその値が定まる。従って、Ｓの一般逆行列を用いて、
Ａ＝Ｓ−１Ｘ＝（Ｓ＊Ｓ）−１・Ｓ＊・Ｘ                （２２ａ）
を計算することで、行列Ａの各要素を求めることができる。
【００６８】
エイリアジング判定手段９ｂでは、以上のような原理に基づいて周波数０から１／τまでのｐについて信号の有無を判定し、信号がある場合にはそのｐについて上記（２３ａ）式によって行列Ａの各要素を求める。
【００６９】
この過程をフローチャートにしたものが図１１である。以下、図に従って説明する。まず、（２１）式のｐ＝０から分析を開始する（ステップ３００）。そして、信号の有無を確認するが（ステップ３０１）、例えば、信号がｐ＝０のとき発見されなければ、（２１）式の各周波数において振幅＝０、即ち、信号はない（ステップ３０２）として、ｐを１だけ増やして（ステップ３０４）再びステップ３０１に戻って、信号の有無を確認する。一方、信号が発見された場合には、上述したように、（２３ａ）式を用いて３つの信号系列から行列Ａの各要素、即ち、各周波数における振幅を計算することで、信号成分が複数あっても正しく算出することができる。そして、算出できたらｐを１だけ増やして（ステップ３０４）再び信号の分析を行う。そして、ｐ＝Ｎ−１まで分析したら、周波数０からｒＭ／τまで分析が終わったことになるので、分析を終了する。
【００７０】
このように、周波数を０から１／τまで分析するだけで、周波数０からｒＭ／τまで分析を行うことになるため、信号を広帯域化しても効率のよい分析を行うことができる。
【００７１】
また、エイリアジング成分と複数の真の周波数成分とが重なりあっている場合でも、互いに異なる微少時間だけずれた信号系列を３つ生成し、位相差を比較することで、複数の周波数成分を分離して正しい周波数を算出することができるため、正確な分析を行うことができる。
【００７２】
以上の説明では、全ての信号成分を（２２）式を用いて算出する例について示したが、（２２）式は（１１）式に比べて計算量が多くなるため、実施の形態３に示したように、信号を単一か否かを判定するステップを挿入し、信号が単一の場合には（１１）式を用いて周波数を算出するようにしてもよい。
【００７３】
図１２に信号成分が単一か否かを判定しながら、信号成分を算出する場合のフローチャートを示す。図１１と異なる点は、ステップ４０４にて実施の形態３で行ったように（２０）式を用いて信号が単一か否かを判定する点である。信号成分が単一ならば実施の形態１と同様に算出し、複数ある場合は本実施の形態４で上述したように（２２）式を用いて全ての信号成分を算出する。
【００７４】
なお、本実施の形態４では、生成する信号系列を３つとしたため、周波数成分が２つ以内の場合に対応できるが、重み付け手段と、高速フーリエ変換手段を更に増やせば、さらに多くの周波数成分がエイリアジング成分と重なっても正確に分離することができる。
【００７５】
以上のように、本実施の形態４では、生成した３つの信号系列を短時間フーリエ変換し、エイリアジング成分の性質を利用して行列の要素を決定し、位相差を比較して信号成分を算出するため、複数の信号が同じ周波数に重なっている場合でも、互いに分離して算出することができる。
【００７６】
実施の形態５．
以上、実施の形態１〜４では、広帯域化に伴って発生するエイリアジング成分を排除しつつ、計算量の増大を抑え効率のよい周波数分析装置について説明したが、本実施の形態５では、窓関数の効果による周波数成分の拡散による誤差を排除しつつ、短時間のサンプリングで正確な周波数分析が可能な周波数分析装置について説明する。
【００７７】
図１３は本実施の形態５に係わる周波数分析装置の装置構成を表すブロック図である。図において図１と同じ構成要素には同じ符号を付す。１１は畳込み効果判定手段である。
【００７８】
以下、動作について説明する。帯域フィルタ１で信号の周波数帯を適当な帯域に制限し、Ａ／Ｄ変換器２で適当な時間間隔τでサンプリングする。サンプリングしたデジタル信号はメモリ３に記憶し、データ分配手段７によって、データを重み付け手段４ａ、４ｂに分配し、一方の信号系列ｓ（ｎ）に対して、他方ｓ（ｎ）に対してΔτだけずれたｓΔ（ｎ）が生成するようにする。それぞれの重み付け手段４ａ、４ｂでは窓関数による重み付けを行い、高速フーリエ変換手段５ａ、５ｂにおいてフーリエ変換を行う。信号検出手段８ではフーリエ変換後の各周波数をチェックして信号の有無を確認し、そして、畳込み効果判定手段１１では、このフーリエ変換後の周波数スペクトル同士の位相を比較することで、真の周波数を特定する。
【００７９】
このように構成することで、短時間のサンプリングで畳込みによる周波数成分の拡散を抑えることができ、連続波のみならず、断続波の場合でも、正確な周波数スペクトルを得る。
【００８０】
ここで、窓関数による畳込みの効果には以下に示すような特徴がある。即ち、（Ａ）畳込みにより拡散する各周波数の信号強度の和は元の信号強度にほぼ一致する。（Ｂ）１つの信号成分（周波数ｆ）から畳込みにより拡散した各周波数成分について、（１１）式は同じ値を持つ。（Ｃ）ハニング窓関数を用いることで、１つの信号の畳込みによる拡散の影響は、真の周波数周辺の３つの周波数値に限定することができる。
【００８１】
（Ａ）は、図２３および、（８）式の説明で示したように、畳込み演算の性質によるものである。この場合、周波数は離散値を取るため、近似的に一致する。また、（Ｂ）については、ＳＴＦＴ、ＳＴＦＴΔとも同じ元の周波数成分を窓関数ｇ（ｎ）により拡散したものであるから相互の位相の関係は同じになり、即ち、位相を表す（１０）式が異なる値を取ることはないからである。
【００８２】
従って、これらから畳込みによる周波数分布の拡散を排除して正しい周波数分布を求めるには、検出された周波数の信号成分が畳込みによる拡散の結果であることを（１０）式により判定し、畳込みによる拡散によると判定された場合には、拡散した各信号成分の信号強度を足し合せ、（１０）式のより求められた周波数に足し合せた信号強度の信号成分があるとすればよい。
【００８３】
なお、（Ｃ）で示したように、周波数の拡散はハニング窓関数を用いることで、真の周波数周辺の３つの周波数値に限定することができる。そのため、Δτについて（１０）式のｆが一意の値を持つためには、３つの周波数の幅の内で、（１０）式左辺の位相が１回転しないような値に限定する必要がある。即ち、ｅｘｐ（２πｊｆΔτ）が、（ｐ−１）／Ｎτ≦ｆ＜（ｐ＋１）／Ｎτで一意となるには、
Δτ＜Ｎτ／２                                  （２３）
であることが必要である。従って、本実施の形態５では、（２３）式のようなΔτとなるよう、データ分配手段７を動作させる。
【００８４】
図１４に本実施の形態５の信号のサンプリング方法を表す。また、フーリエ変換後の周波数スペクトルを図１５に示す。図１５のように、本実施の形態５はエイリアジングの影響を排除するための前記実施の形態１から４までと異なり、Δτはτより長い時間間隔となっている。（２３）式の条件からすれば、Δτはτより短くても問題はないが、その場合、あまり短く取りすぎると、雑音の影響が出てくるため、Δτの間隔は適用した各装置の特性に応じて定められる。
【００８５】
図１６に本実施の形態５の信号検出手段８と畳込み判定手段１１の動作のフローチャートを示す。まず、ｐ＝０から分析を開始し（ステップ５００）、信号の有無を検出する（ステップ５０１）。信号が検出された場合は、（１０）式から左辺のｆを算出する（ステップ５０２）。このｆにおける、信号強度に｜ＳＴＦＴ（２πｐ／Ｎτ）｜２を加算する（ステップ５０３）。これが周波数ｐ／Ｎτにおける１サイクルであり、以上が終了したらｐの値を１増やして（ステップ５０４）、再び信号の有無を確認してから（１０）式を用いてｆを計算する。
【００８６】
上述したように、畳込みによって１の周波数が３つの周波数に拡散するが、（１０）式のｆは同じ値が出てくるため、同じｆならば畳込みによる拡散で生じたスペクトルであるとして、ステップ５０２および５０３の手順によって畳込みの影響を除去する。このように従来例のように真の周波数周辺の周波数を３つ選択するのではなく、（１０）式によってあらかじめ真の周波数を求めて、この周波数における信号強度を足し合せていくため、真の周波数が離散周波数値のほぼ中間にあっても、精度よく真の周波数を算出できる。これをｐ＝０〜Ｎまで行って（ステップ５０５）、分析可能範囲をすべて分析する。
【００８７】
以上のように、本実施の形態５では、サンプリングした信号から同時に２つの信号列を生成し、それらの位相を比較することで、畳込みによってスペクトルが拡散する効果を抑え、また、加算を主体にして信号強度を算出するため、正確な周波数分析を行う周波数分析装置を得る。
【００８８】
実施の形態６．
上記実施の形態５では、２つの信号列を得るために、時間間隔τでサンプリングを行った結果を記録し、その中から、適当なタイミングのデータをそれぞれ抽出し、２つの信号列の時間差がΔτとなるようにしたものを示したが、本実施の形態６では、図１７のように、タイミングのΔτだけ異なるクロックでサンプリングを行う２つのＡ／Ｄ変換器を使用して実現してもよい。
【００８９】
図１７に本実施の形態６に係わる周波数分析装置の装置構成を表すブロック図で示す。図６と同じ構成要素には同じ符号を付す。１１は上記実施の形態５と同様の畳込み効果判定手段である。本実施の形態６の周波数分析装置は図のように２つのＡ／Ｄ変換器２ａ、２ｂを使用するが、無駄な信号をサンプリングすることがなくなるため、メモリ３ａ、３ｂの容量を小さくて済む。
【００９０】
以下、動作について説明する。帯域フィルタ１は入力された受信信号を帯域制限して分析対象の周波数信号のみ選択する。信号分配器１０では、信号を２つのＡ／Ｄ変換器２ａ、２ｂに分配する。各Ａ／Ｄ変換器２ａ、２ｂではそれぞれ時間間隔τでサンプリングを行うが、例えば、Ａ／Ｄ変換器２ｂの方がＡ／Ｄ変換器２ａよりもΔτだけ遅れたタイミングでサンプリングを行う。これによって、実施の形態５で示したような２つの信号列ｓ（ｎ）およびｓΔ（ｎ）が生成され、順次メモリ３ａ、３ｂに記憶される。以下の動作は上記実施の形態５の周波数分析装置の動作と同様である。
【００９１】
以上のように、本実施の形態６では、帯域制限した信号を予め分配してＡ／Ｄ変換器に入力するため、無駄な信号をサンプリングすることがなく、メモリの容量が小さくて済む。
【００９２】
実施の形態７．
上記実施の形態５および６では、受信信号の周波数分布が単一もしくは比較的疎な場合のみ分析が可能であったが、本実施の形態７では同じ周波数に畳込みによって拡散した異なる信号成分が重なりあった場合でも、信号成分が単一か否かを判定し、信号成分が単一の場合のみ分析を行う周波数分析装置について説明する。
【００９３】
図１８は本実施の形態７の周波数分析装置の装置構成を表すブロック図である。本実施の形態７では、信号系列を３つ用いてるため、実施の形態３の図８と同じ構成要素には同じ符号を付す。１１ａは本実施の形態７の特徴である、各周波数における信号成分が単一か否かを判定しながら畳込みの効果を抑えて正しい周波数を算出する畳込み効果判定手段である。
【００９４】
次に、動作について説明する。１つの信号系列ｓ（ｎ）およびそのｓ（ｎ）から時間Δτ１およびΔτ２だけずれた２つの信号系列ｓΔ１（ｎ）およびｓΔ２（ｎ）の計３つの信号系列をサンプリングし、それぞれ短時間フーリエ変換する。そこから、信号検出手段８および畳込み効果判定手段１１ａによって、上記実施の形態３と同様に（２１）式のような単位時間当たりの位相差を求め、信号が一つか否かを判定しながら真の周波数を求める。
【００９５】
以下、信号検出手段８および畳込み効果判定手段１１ａの動作については、図１９に示したフローチャートにしたがって説明する。まず、信号検出手段８ｂで離散値の各周波数ｆ＝ｐ／Ｎτをｐ＝０（ステップ６００）から順番に信号が存在するかの判定が行われる（ステップ６０１）。そして、信号が未検出ならば、ｐ＝ｐ＋１（ステップ６０６）とし、次のｐについて再び信号が存在するか判定を行う。
【００９６】
一方、ｐ＝０で信号が検出された場合は、ＳＴＦＴ（ω）、ＳＴＦＴΔ１（ω）、ＳＴＦＴΔ２（ω）、ＳＴＦＴΔ３（ω）から（２０）式を計算し（ステップ６０２）、全て同じ値を算出するか判定する（ステップ６０３）。もし、（２１）式の計算結果が一致しないならば、上述したように信号成分が複数存在することになるので、当該周波数の成分は不定であるとして、次のｐの値について信号の存在の有無を判定する。一方、（２０）式がすべて同じ値ならば信号は１つであり、周波数の分析が可能であるので実施の形態５と同様に（１１）式を満足する周波数ｆを算出する（ステップ６０４）。そして、算出した周波数ｆにおける信号強度に補正を加える（ステップ６０５）。
【００９７】
以上、いずれかの処理が終わったら周波数ｐ／Ｎτのｐをｐ＝ｐ＋１として、次の周波数成分について同様の処理を行う（ステップ６０６）。その際、ｐ≧Ｎであれば、演算する部分が重なってしまうため、ｐ≧Ｎであれば分析を終了し、ｐ＜Ｎの場合のみステップ５０１に戻る（ステップ６０７）。
【００９８】
以上のように、本実施の形態７では、信号がいくつあるか判定しながら周波数分析を行うため、信号の周波数成分が複数ある場合でも、誤った周波数を算出せず、短時間のサンプリングで周波数分析が終了し、また、加算を主体にして信号強度を算出するため、正確な周波数分析を行うことが可能な周波数分析装置を得る。
【００９９】
【発明の効果】
以上のように、本発明の周波数分析装置は、サンプリングしたデジタルデータから同時に２つの信号系列を生成し、この２つの信号系列を短時間フーリエ変換する。そして、それぞれの位相差を算出し、エイリアジングを利用して、周波数帯域全てを走査しなくても受信信号の存在する周波数を算出することができるため、分析すべき信号の周波数帯が広帯域化しても計算量の増大を抑え、エイリアジングの効果を排除し、正確な周波数分析結果を得ることができる。
【０１００】
また、帯域制限した信号を予め分配してＡ／Ｄ変換器に入力するため、無駄な信号をサンプリングすることがなく、メモリの容量が小さくて済む。
【０１０１】
また、生成した３つの信号系列を短時間フーリエ変換し、各信号系列の単位時間当たりの位相差を算出して各周波数において信号成分が単一か否かを判定しながら分析を行うため、信号成分が複数ある場合でも、誤った周波数を算出しない。
【０１０２】
また、生成した３つの信号系列を短時間フーリエ変換し、エイリアジング成分の性質を利用して行列の要素を決定し、位相差を比較して信号成分を算出するため、複数の信号が同じ周波数に重なっている場合でも、互いに分離して算出することができる。
【０１０３】
また、加算によって信号強度を算出する畳込み効果判定手段を設けたため、窓関数の効果による周波数成分の拡散による誤差を排除し、正確な周波数分析結果を得る。
【図面の簡単な説明】
【図１】  本発明の実施の形態１に係わる周波数分析装置の構成を表すブロック図である。
【図２】  本発明の実施の形態１に係わる周波数分析装置のサンプリング方法を表す図である。
【図３】  本発明の実施の形態１に係わる周波数分析装置の周波数分析方法を表す図である。
【図４】  本発明の実施の形態１に係わる周波数分析装置の周波数算出方法のフローチャートである。
【図５】  本発明の実施の形態２に係わる周波数分析装置のサンプリング方法を表す図である。
【図６】  本発明の実施の形態２に係わる周波数分析装置の装置構成を表すブロック図である。
【図７】  本発明の実施の形態３が対象とする周波数分布を表す図である。
【図８】  本発明の実施の形態３に係わる周波数分析装置の装置構成を表すブロック図である。
【図９】  本発明の実施の形態３に係わる周波数分析装置の周波数算出方法のフローチャートである。
【図１０】  本発明の実施の形態４に係わる周波数分析装置の装置構成を表すブロック図である。
【図１１】  本発明の実施の形態４に係わる周波数分析装置の周波数算出方法のフローチャートである。
【図１２】  本発明の実施の形態４に係わる周波数分析装置の別の周波数算出方法のフローチャートである。
【図１３】  本発明の実施の形態５に係わる周波数分析装置の装置構成を表すブロック図である。
【図１４】  本発明の実施の形態５に係わる周波数分析装置のサンプリング方法を表す図である。
【図１５】  本発明の実施の形態５に係わる周波数分析結果の一例である。
【図１６】  本発明の実施の形態５に係わる周波数分析装置の周波数算出方法のフローチャートである。
【図１７】  本発明の実施の形態６に係わる周波数分析装置の装置構成を表すブロック図である。
【図１８】  本発明の実施の形態７に係わる周波数分析装置の装置構成を表すブロック図である。
【図１９】  本発明の実施の形態７に係わる周波数分析装置の周波数算出方法のフローチャートである。
【図２０】  従来の周波数分析装置の装置構成を表すブロック図である。
【図２１】  短時間フーリエ変換の原理を表す図である。（ａ）サンプリングにより生成された信号列である。（ｂ）窓関数の形状を表す図である。（ｃ）上記生成された信号列と窓関数とを畳込み演算によって重み付けした後の図である。（ｄ）それぞれの時間における周波数分布を表す図である。（ｅ）各時間における周波数分布を並べて時間変化を見る図である。
【図２２】  エイリアジングの効果を表す図である。
【図２３】  窓関数との畳込みの効果を表す図である。
【図２４】  従来の周波数分析装置を表すブロック図である。
【符号の説明】
１  帯域フィルタ、    ２、２ａ、２ｂ  Ａ／Ｄ変換器、
３、３ａ、３ｂ  メモリ、    ４ａ、４ｂ  重み付け手段、
５ａ、５ｂ  高速フーリエ変換手段、    ６  表示器、
７  データ分配手段、    ８  信号検出手段、
９、９ａ  エイリアジング判定手段、    １０  信号分配器、
１１、１１ａ  畳込み効果判定手段。
【図１】
0003646702.tif 000024
【図２】
0003646702.tif 000025
【図３】
0003646702.tif 000026
【図４】
0003646702.tif 000027
【図５】
0003646702.tif 000028
【図６】
0003646702.tif 000029
【図７】
0003646702.tif 000030
【図８】
0003646702.tif 000031
【図９】
0003646702.tif 000032
【図１０】
0003646702.tif 000033
【図１１】
0003646702.tif 000034
【図１２】
0003646702.tif 000035
【図１３】
0003646702.tif 000036
【図１４】
0003646702.tif 000037
【図１５】
0003646702.tif 000038
【図１６】
0003646702.tif 000039
【図１７】
0003646702.tif 000040
【図１８】
0003646702.tif 000041
【図１９】
0003646702.tif 000042
【図２０】
0003646702.tif 000043
【図２１】
0003646702.tif 000044
【図２２】
0003646702.tif 000045
【図２３】
0003646702.tif 000046
【図２４】
0003646702.tif 000047