JP3646702B2 - 周波数分析装置 - Google Patents

周波数分析装置 Download PDF

Info

Publication number
JP3646702B2
JP3646702B2 JP2002023783A JP2002023783A JP3646702B2 JP 3646702 B2 JP3646702 B2 JP 3646702B2 JP 2002023783 A JP2002023783 A JP 2002023783A JP 2002023783 A JP2002023783 A JP 2002023783A JP 3646702 B2 JP3646702 B2 JP 3646702B2
Authority
JP
Japan
Prior art keywords
frequency
signal
data string
distribution data
time interval
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
JP2002023783A
Other languages
English (en)
Other versions
JP2003222646A (ja
Inventor
充 田沼
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Mitsubishi Electric Corp
Original Assignee
Mitsubishi Electric Corp
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Mitsubishi Electric Corp filed Critical Mitsubishi Electric Corp
Priority to JP2002023783A priority Critical patent/JP3646702B2/ja
Publication of JP2003222646A publication Critical patent/JP2003222646A/ja
Application granted granted Critical
Publication of JP3646702B2 publication Critical patent/JP3646702B2/ja
Anticipated expiration legal-status Critical
Expired - Fee Related legal-status Critical Current

Links

Images

Landscapes

  • Complex Calculations (AREA)

Description

【0001】
【発明の属する技術分野】
本発明は、受信した信号の周波数分布を短時間フーリエ変換によって分析する周波数分析装置に関するものである。
【0002】
【従来の技術】
レーダなどの無線信号を受信し、その信号の周波数成分を短時間フーリエ変換によって分析する周波数分析装置においては、レーダ送信波の周波数ホッピング等に対応するため、分析する周波数帯域が広いことが望まれる。
【0003】
ところで、一般的に時間間隔τでAD変換した信号は、元々の信号の周波数fのほかに、f+r/τ(rは整数)の信号成分を持つことが知られている。この成分をエイリアジングと呼ぶ。このエイリアジングの影響を避けながら周波数帯域を広帯域化するためには、帯域の逆数に比例してデータを細かい時間間隔で変換した信号を扱う必要がある。これは即ち、扱うデータ量が周波数帯域に比例して多くのデータに対する処理を行う必要があることを表す。例えば、周波数帯域を16倍にする場合、AD変換される信号の時間間隔も1/16細かくする必要があり、このため扱うデータ数は16倍になる。短時間フーリエ変換を効率的に行うFFTアルゴリズムにおいては、データ数がZであるとすると、そのデータ数Zに対して計算量はZlog2Zに比例して増大することが知られており、データ数ZをY倍にすると、その計算量は(ZYlog2ZY)/(ZlogZ)=Y×(1+(log2Y/log2Z))となる。一例として、Z=128、Y=16の場合、約25倍の計算量が要求され、この処理負荷の増大が広帯域信号の分析を困難にしている。
【0004】
また、短時間フーリエ変換には、窓関数による畳込み演算によって真の周波数成分が拡散してしまい精度良く周波数の算出ができない等の問題がある。
【0005】
以下に、短時間フーリエ変換を用いた従来の周波数分析装置について説明し、まずエイリアジングの発生による問題点を明らかにする。図19は従来の周波数分析装置の構成を示している。図において、1は受信される電波信号から、監視する対象の周波数成分のみを選択する帯域フィルタ、2はアナログ信号を一定の時間間隔でサンプリングし、デジタル信号に変換するAD変換器、3はデジタル信号に変換された入力信号を順次記録するメモリ、4はメモリ3に記録された信号を一定の長さずつ読み出し、窓関数を乗じて重み付けを行う重み付け手段、5は重み付けされた一定の長さの信号に順次フーリエ変換を行い、周波数分析を行うフーリエ変換手段、6は分析結果を出力する表示装置である。
【0006】
次に、動作について説明する。図19に示した周波数分析装置は、短時間フーリエ変換の原理を使用して受信する信号の周波数分布を算出するものである。短時間フーリエ変換については、貴家仁志「マルチレート信号処理」(昭晃堂)pp.149〜153等の文献で説明されているが、次のようなものである。
【0007】
図20に示したように、短時間フーリエ変換は、信号に有限な長さの窓関数を乗じて切り出し、それをフーリエ変換するという処理技法である。図20(a)はサンプリングにより得られた実際の信号列であり、(b)が窓関数の形状、(c)が信号列と窓関数との畳込み演算により得られる重み付け後の信号であり、(d)はそれぞれをフーリエ変換して周波数成分に直した図であり、(e)は各時間における周波数成分を並べた時間変化を表す図である。窓関数を時間に従いシフトし、処理を繰り返すことで時間と周波数の情報を調べる。
【0008】
短時間フーリエ変換は算術的には以下のように表現できる。
クロック周期τでサンプリングされた信号をs(n)、窓関数をg(n)(g(n)はn=0,τ,2τ,…,(N−1)τで0でない)とする。s(n)は図20(a)に、g(n)は図20(b)に示したような波形であるとする。ここで、g(n)をNτずつシフトした関数g(n−kN)を考え、このg(n−kN)とs(n)を乗ずることで式(1)のようにkブロック目の信号sk(n)を抽出する。
sk(n)=g(n−kN)s(n) (1)
である。図20(c)にsk(n)を示す。
【0009】
(1)式のsk(n)をフーリエ変換する処理が、短時間フーリエ変換である。従って、sN(n)の短時間フーリエ変換は、
【数1】
Figure 0003646702
である。ただし、ω=2πp/Nτである。これより、
【数2】
Figure 0003646702
である。(3)式は入力信号s(n)について、周波数p/Nτにおける周波数成分を表すものである。
【0010】
g(n)は窓関数と呼ばれるもので、例えば、次のようなハニング窓関数などが一般的に用いられる。
g(n)=1/2{1+cos2π/N(n−N/2)}
ただし、n=0〜N−1、それ以外のnではg(n)=0
【0011】
図19は、以上の原理を実現する短時間フーリエ変換を用いた周波数分析装置の構成である。例えば、無線信号が受信されると、適当な中間周波数信号に変換され帯域フィルタ1に入力される。帯域フィルタ1では入力信号について、分析対象とする周波数帯の信号のみを取り出し、A/D変換器2ではこの取り出した信号について適当な時間間隔τでサンプリングして、デジタル信号に変換する。このデジタル変換した信号列が上述のs(n)であり、s(n)は順次メモリ3に記憶される。なお、以下、A/D変換については、IおよびQチャネルを使用するなどにより、信号が複素数として取得されることを想定するが、これは、信号が実数としてサンプリングされる場合に比べると、サンプリング間隔が2倍で良い他は本質的な違いはない。
【0012】
重み付け手段4では、メモリ2に記録された信号を順次読み出し、信号列s(n)に窓関数g(n)を掛け合わせ(1)式に表されたような演算を行い、高速フーリエ変換手段5では、重み付けの行われた信号sk(n)に対し、順次(3)式の処理を行い、各時点での信号に含まれる周波数成分を分析し、その結果は表示装置6に表示する。
【0013】
以上のように、従来の周波数分析装置は、短時間フーリエ変換を使用して受信される電波信号の周波数分布を求めている。ここで、短時間フーリエ変換の性質を考えると、(3)式の性質から以下の特徴が導き出せる。
【0014】
短時間フーリエ変換を行った結果、周波数f=p/Nτ(pは整数)となり、離散的な値を持つ。その内、最も高い周波数はp=N−1のときの(N−1)/Nτである。即ち、サンプリング定理により、エイリアジングの発生しない周波数範囲は1/τである。即ち、分析しようとする最大周波数をfwとすると、fw<1/τならば、周波数の特定が可能であることを意味する。周波数1/τはナイキスト周波数ともいう。
【0015】
また、(3)式の右辺は、
【数3】
Figure 0003646702
のような関係が成り立つ。ここでrは整数である。これは、真の周波数p/Nτからr/τだけずれた周波数に同一の値を持った成分が現れることを意味する。即ち、図21に示したように、本来の信号の周波数の他に、そこから1/τごとにエイリアジングと呼ばれる周波数成分が発生する。
【0016】
帯域フィルタ1によって分析する周波数成分の範囲を1/τ以下に制限すればこのようなエイリアジングの発生を防ぐことができる。しかし、上述のように周波数分析装置には分析周波数の広帯域化が望まれるが、これにより、帯域フィルタの通過させる周波数範囲を1/τ(τはサンプリング間隔)より広帯域化すれば、それだけ多くのエイリアジングが生じ、信号の真の周波数と分離することができない。以上が、エイリアジング発生による問題点である。
【0017】
次に、上述したように、有限の時間範囲で窓関数によって畳込み演算を行うことで原理的にスペクトルが拡散する誤差について説明する。即ち、g(n)のフーリエ変換をG(2πp/Nτ)とすると、
Figure 0003646702
また、(1)式および(2)式から、s(n)のフーリエ変換をS(ω)とすると、
【数4】
Figure 0003646702
ここで、便宜上m=n−kNのように置換すると、
【数5】
Figure 0003646702
従って、G(ω)が複数の値を持つことから、短時間フーリエ変換で計算される周波数分布は、真の周波数の周囲に、窓関数の性質にしたがって拡散する。この拡散の様子を示したものを図22に示す。図22(a)は信号の元の周波数分布であり、(b)は窓関数によって拡散した周波数分布を表す。拡散した各成分の強度の和は、(8)式の畳込み演算の性質から元の信号の強度に一致する。また、その拡散の程度は窓関数g(n)の性質によるが、性質の良いハニング窓関数であっても、真の周波数の高低に3つの拡散成分が発生することは免れない。
【0018】
この拡散成分を抑え真のスペクトルを算出する周波数分析装置の一例として、例えば、特開平6−160445号公報に示された周波数分析装置がある。図23に特開平6−160445号公報の周波数分析装置を示す。図中、1はAD変換器、2はメモリ、3は窓関数重み付け手段、4は高速フーリエ変換手段、5は高速フーリエ変換器4の出力側に設けられたスペクトル補間判定手段、6は表示器、7は補間判定手段5の判定結果によりスペクトル補間を実行するスペクトル補間手段、8はスペクトル補間手段7によって算出したスペクトルの周波数、振幅、位相データを時間軸データに逆変換する逆フーリエ変換器、9は逆フーリエ変換器8で逆フーリエ変換した時間軸データを累積する累積器、11は累積器9によって累積された累積結果をメモリ2から読み出される時間データから減算し、この減算結果を再び窓関数重み付け手段3に入力する減算器、12はスペクトル補間手段7によって算出した周波数領域データを累積する累積器、13はこの累積器12に累積した周波数領域データを減算器11で減算した残差データを再度高速フーリエ変換して得られた周波数領域データS(t)に加算する加算器である。
【0019】
この周波数分析装置では、フーリエ変換して得られた離散周波数スペクトルの中で3本の周波数スペクトルが互いに隣接して存在し、しかもその上位と下側のスペクトルの位相がほぼ同じであったら、この3本のスペクトルは真のスペクトルが分散されているものと判定し、補間処理によりその真のスペクトルの周波数と振幅、位相を求める。
【0020】
この真のスペクトルデータを累積器12に累積するとともに、逆フーリエ変換して時間軸データに逆変換し、この時間軸データを累積器9で累積し、その累積データを入力データ列から減算するから、補間処理により真のスペクトル線を算出したあとでは残差エネルギーが低下し、線スペクトルは広がりを持つことはなく、線スペクトルとして表示される。この結果、スペクトル分解能が向上する。また、補間処理により振幅、位相も真値近くに近似して算出するから振幅値および位相も正確に測定することができる。
【0021】
しかし、この場合、フィードバックループにより累積したデータを繰り返し減算して真の周波数成分を取り出していく構成であるため、値が収束するまでに何回もサンプリングを行う必要があり、時間がかかり、また、フィードバックループを無くした場合は、微少信号を抽出できない。
【0022】
また、上記の例では真の周波数の周辺の周波数を正しく3つ選択できることが必要だが、真の周波数が、離散周波数値のほぼ中間に存在した場合、例えば、n/Nτと(n+1)/Nτの中間の場合、真の周波数に最も近い周波数を(n+1)/Nτと誤判定しやすく、真の周波数の周辺の3つの周波数を間違えてn/Nτ、(n+1)/Nτ、(n+2)/Nτのように選んでしまうことが起こりやすいため、実際の周波数算出は不安定になりやすい。
【0023】
さらに、この場合、真の周波数の算出方法は特開平6−160445号公報の(14)式などのように引き算が主体であり、sinやcosが0の付近では計算結果が不安定であり、周波数算出に誤差が出やすい。
【0024】
【発明が解決しようとする課題】
本発明は、以上のような問題点を解決するためになされたもので、短時間フーリエ変換の際に、受信信号を広帯域化するに伴う計算量の増大を抑え、周波数分布を効率よく正確に算出することを目的とする。
【0025】
【課題を解決するための手段】
本発明に係わる周波数分析装置は、周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、上記帯域制限されたアナログ入力信号を第一の時間間隔で順次デジタルデータ列に変換するAD変換器と、上記AD変換器によって変換された入力データ列を順次記憶するメモリと、上記メモリに記憶された入力データ列から上記第一の時間間隔よりも長い第二の時間間隔を持った第一の分配データ列、およびこの第一の分配データ列から上記第一の時間間隔ずつ遅れた第二の分配データ列とを抽出するデータ分配手段と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、上記第一または第二の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、信号が検出された場合、上記第一および第二の周波数領域データの位相を比較して真の周波数を算出するエイリアジング判定手段とを備えるものである。
【0026】
また、周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、上記帯域制限されたアナログ入力信号を2つに分配する信号分配器と、上記分配された一方のアナログ入力信号を一定の時間間隔でデジタル信号に変換し、第一の分配データ列を生成する第一のA/D変換器と、上記分配された他方のアナログ入力信号を、上記第一のA/D変換器のデジタル変換より上記一定の時間間隔より短い時間間隔だけ遅れて、一定の時間間隔でデジタル信号に変換し、第二の分配データ列を生成する第二のA/D変換器と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、上記第一または第二の周波数領域データを各周波数ごとに信号が存在するかを検出する信号検出手段と、信号が検出された場合、上記第一および第二の周波数領域データを比較することで、エイリアジングの影響を排除し、真の周波数を算出するエイリアジング判定手段とを備えるものである。
【0027】
また、周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、上記帯域制限されたアナログ入力信号を第一の時間間隔で順次デジタルデータ列に変換するAD変換器と、上記AD変換器によって変換された入力データ列を順次記憶するメモリと、上記メモリに記憶された入力データ列から上記第一の時間間隔よりも長い第二の時間間隔を持った第一の分配データ列と、この第一の分配データ列から上記第一の時間間隔ずつ遅れた第二の分配データ列と、この第一の分配データ列から上記第一および第二の時間間隔とは異なる第三の時間間隔ずつ遅れた第三の分配データ列とを抽出するデータ分配手段と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、上記第三の分配データ列に所定の窓関数を乗算し、重み付けを行う第三の重み付け手段と、この重み付けされた第三の分配データ列をフーリエ変換し、第三の周波数領域データに変換する第三の高速フーリエ変換手段と、上記第一または第二または第三の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、信号が検出された場合、上記第一、第二、第三の周波数領域データのいずれか2つの単位時間当たりの位相差を比較して、信号が検出された周波数における信号成分が単一か否かを判定し、信号成分が単一の場合のみ真の周波数を算出するエイリアジング判定手段とを備えるものである。
【0028】
また、周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、上記帯域制限されたアナログ入力信号を第一の時間間隔で順次デジタルデータ列に変換するAD変換器と、上記AD変換器によって変換された入力データ列を順次記憶するメモリと、上記メモリに記憶された入力データ列から上記第一の時間間隔よりも長い第二の時間間隔を持った第一の分配データ列と、この第一の分配データ列から上記第一の時間間隔ずつ遅れた第二の分配データ列と、この第一の分配データ列から上記第一および第二の時間間隔とは異なる第三の時間間隔ずつ遅れた第三の分配データ列とを抽出するデータ分配手段と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、上記第三の分配データ列に所定の窓関数を乗算し、重み付けを行う第三の重み付け手段と、この重み付けされた第三の分配データ列をフーリエ変換し、第三の周波数領域データに変換する第三の高速フーリエ変換手段と、上記第一または第二または第三の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、信号が検出された場合、上記第一、第二、第三の各周波数における周波数領域データを要素とする行列に、当該信号が検出された周波数から上記第二の時間間隔の逆数の整数倍ずつずれた複数の周波数を要素とする行列の逆行列を乗算することで真の周波数を算出するエイリアジング判定手段とを備えるものである。
【0029】
また、周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、上記帯域制限されたアナログ入力信号を第一の時間間隔であらかじめ決められた一定回数だけ順次デジタルデータに変換するAD変換器と、上記AD変換器によって変換された入力データを順次記憶するメモリと、上記メモリに記憶された入力データから上記第一の時間間隔を持ったデータを抽出してなる第一の分配データ列、および上記第一の時間間隔に上記デジタル変換した一定の回数を掛けた数の二分の一より短い第二の時間間隔ずつ上記第一の分配データ列より遅れた第二の分配データ列とを抽出するデータ分配手段と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、上記第一または第二の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、信号が検出された場合、上記第一および第二の周波数領域データの位相を比較して、この位相が同じの場合、真の周波数を算出し、算出した周波数において当該信号強度を足し合せていく処理を行う畳込み効果判定手段とを備えるものである。
【0030】
また、周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、上記帯域制限されたアナログ入力信号を第一の時間間隔であらかじめ定められた一定回数だけ順次デジタルデータ列に変換するAD変換器と、上記AD変換器によって変換された入力データ列を順次記憶するメモリと、上記メモリに記憶された入力データから上記第一の時間間隔を持ったデータを抽出してなる第一の分配データ列、および上記第一の時間間隔に上記デジタル変換した一定の回数を掛けた数の二分の一より短い第二の時間間隔ずつ上記第一の分配データ列より遅れた第二の分配データ列と、上記第一の時間間隔に上記デジタル変換した一定の回数を掛けた数の二分の一より短かく、かつ上記第一および第二の時間間隔とは異なる第三の時間間隔ずつ上記第一の分配データ列より遅れた第三の分配データ列とを抽出するデータ分配手段と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、上記第三の分配データ列に所定の窓関数を乗算し、重み付けを行う第三の重み付け手段と、この重み付けされた第三の分配データ列をフーリエ変換し、第三の周波数領域データに変換する第三の高速フーリエ変換手段と、上記第一または第二または第三の周波数領域データいずれかの各周波数ごとに信号が存在するかを検出する信号検出手段と、信号が検出された場合、上記第一、第二、第三の周波数領域データのいずれか2つの単位時間当たりの位相差を比較して、信号が検出された周波数における信号成分が単一か否かを判定し、信号成分が単一の場合のみ上記第一および第二の周波数領域データの位相を比較して、この位相が同じの場合、真の周波数を算出し、算出した周波数において当該信号強度を足し合せていく処理を行う畳込み効果判定手段とを備えるものである。
【0031】
【発明の実施の形態】
実施の形態1.
本発明の実施の形態1について説明する。本実施の形態1に係わる周波数分析装置は、A/D変換された信号から微少タイミングずらしてサンプルして得られた2系列の信号列をそれぞれ短時間フーリエ変換処理して局所的な周波数分布を求めるものである。また、その際生じるエイリアジングを微少タイミングずらしたことにより生じる信号列間の位相差を検出することで真の周波数を特定するものである。
【0032】
以下、詳細について説明する。図1は本実施の形態1に係わる周波数分析装置の構成を表すブロック図である。図1において、図19と同じ構成要素には同じ符号を付す。7はメモリ3内に記録した信号のデータからサンプリングタイミングのわずかに異なる2つのデータ列を抽出、分配するデータ分配手段、8はフーリエ変換結果から各周波数に関して信号の存在を判定する信号検出手段、9は信号検出手段8で判定された信号の存在する周波数に関して、2つの短時間フーリエ変換の結果を比較してエイリアジングの影響を排除し、正しい周波数を算出するエイリアジング判定手段である。
【0033】
図2は本実施の形態1に係わる周波数分析装置のサンプリング方法を表す図である。図2に示した通り、時間間隔Δτでサンプリングしたデータの中から、適当な時間間隔τの信号列s(n)と、当該s(n)の各信号からΔτだけ遅れた信号列sΔ(n)を考える。信号の周波数をf(=p/Nτ,pは整数)とすると、信号列s(n)とsΔ(n)とは、位相遅れ2πfΔτが発生するため、
【数6】
Figure 0003646702
の関係が成立する。ここで、係数exp(2πjfΔτ)は定数であるため、重み付け、フーリエ変換とあった短時間フーリエ変換の各プロセスによって変化することはなく、信号列sΔ(n)は、s(n)に対して一様な周波数fによる位相差を持つ。このため、sΔ(n)のフーリエ変換をSTFTΔ(2πp/Nτ)とすると、
【数7】
Figure 0003646702
(10)式より、STFT(2πp/Nτ)≠0となる周波数fに対して、元のデータ列s(n)と、微少時間Δτだけずらしてサンプリングを行ったデータ列sΔ(n)のそれぞれの短時間フーリエ変換を比較すると、係数exp(2πjfΔτ)の違いとなる。従って、(10)式を変形すると、
【数8】
Figure 0003646702
であり、(11)式右辺の各FFT成分については、それぞれ間隔τのデータを使用しているため、pに対する周波数fは、
f=p/Nτ,(p+N)/Nτ,(p+2N)/Nτ…
のいずれかである。従って、rを整数としてf=(p+rN)/Nτのように表現でき、(11)式は、
【数9】
Figure 0003646702
である。その際、−2πjrNΔτ/Nτの項によって位相が一回転しない場合は、rを一意に求めることができる。即ち、
【数10】
Figure 0003646702
を満たす場合である。(12)式を変形すると、
1/Δτ>rM/τ (14)の関係が得られる。1/Δτは本実施の形態1の周波数分析装置の特定可能な周波数の範囲である。rMが余りに大きい場合には(14)式を満たすことができなくなるため、帯域フィルタ1で適当な周波数帯域を選択する。
【0034】
図3に(13)式の関係を表すグラフを示す。図3は縦軸を振幅、横軸を周波数とした周波数分析結果の図である。周波数の値は離散値でその最小単位は1/Nτである。また、1/Δτがこの周波数分析装置で分析できる最大の周波数であり、1/Δτより上の周波数が入力されないように帯域フィルタ1で帯域制限する。一方、サンプリング間隔τの信号列s(n)を短時間フーリエ変換するのであるから、エイリアジングの発生しない周波数範囲は1/τである関係上、真の信号周波数から1/τずつずれたところにエイリアジングが発生する。
【0035】
図では例えば、周波数(p+2N)/Nτに真の信号成分がある場合である。信号列s(n)およびsΔ(n)から見るとエイリアジングの発生しない周波数範囲は1/τであるのに対し、帯域フィルタ1から入力される受信信号はその周波数範囲を超えているので、(11)式の右辺を実測によって求め、そこから左辺のpおよびrを求める。
【0036】
以上が、本実施の形態1に係わる周波数分析装置の原理である。図1の周波数分析装置は以上の原理に基いて短時間フーリエ変換を行うものである。受信した電波を適当な中間周波数に変換した受信信号が、帯域フィルタ1に入力され、帯域制限される。A/D変換器では、選択された信号をナイキスト周波数の逆数以下である適当な短い時間間隔Δτでサンプリングし、メモリ3に順次蓄積する。データ分配手段7では、メモリに蓄積した信号からそれぞれ適当な間隔τの2系列の信号列s(n)およびsΔ(n)を選択し、それぞれ重み付け手段4aおよび4bに分配する。重み付け手段4aおよび4bにおいて窓関数による重み付けがなされた後、フーリエ変換手段5aおよび5bに入力される。フーリエ変換手段5aおよび5bでは、重み付けされたs(n)およびsΔ(n)のフーリエ変換を行い、(2)式のように、STFT(ω)およびSTFTΔ(ω)を求める。
【0037】
信号検出手段8ではフーリエ変換手段5aの出力STFT(ω)の結果を調べ、信号の存在する周波数を判定し、それらの周波数をエイリアジング判定手段9に通知する。エイリアジング判定手段9では、STFT(ω)およびSTFTΔ(ω)を比較し、(10)式の左辺exp(2πjfΔτ)を算出する。
【0038】
信号検出手段8およびエイリアジング判定手段9で行う処理を図4に示したフローチャートにしたがって説明する。まず、信号検出手段8で離散値の各周波数f=(p+rN)/Nτをp=0(ステップ100)から順番に信号が存在するかの判定が行われる(ステップ101)。そして、信号が未検出ならば、そのときのp(今の場合、p=0)に対して、全てのrについて受信信号の周波数成分が0である(ステップ102)。即ち、p=0で信号が未検出であれば、そのエイリアジングの可能性がある周波数f=N/Nτ,2N/Nτ,3N/Nτ,…における振幅=0ということである。
【0039】
一方、p=0で信号が検出された場合は、STFT(ω)およびSTFTΔ(ω)から(11)式を満足するrを算出する(ステップ103)。そして、算出された値以外のrについての受信信号の周波数成分は0である(ステップ104)。例えば、r=3と算出された場合は、それ以外の周波数f=(3+N)/Nτ,(3+2N)/Nτ,(3+4N)/Nτ,…における振幅=0ということである。
【0040】
そして、p=0のとき、信号が検出、または未検出いずれかの場合の処理が終了すると、次の周波数、即ちp=1(ステップ105)の場合について、同様の処理を行う。その際、p≧Nであれば、演算する部分が重なってしまうため、p≧Nであれば分析を終了し、p<Nの場合のみステップ101に戻る(ステップ106)。
【0041】
このように、p=0からNまで、即ち離散値の周波数fを0から1/τまで信号の有無を検出するだけで、実際は図3の0からrM/τまでの信号帯域を分析することになるので、全ての周波数を分析する必要がなく、計算量が大幅に軽減することができる。
【0042】
例えば、サンプル数N=128の場合に、周波数帯域を16倍に広げることを考えると、
(1)従来の周波数分析装置で周波数帯域を16倍に広げる場合
フーリエ変換を行うデータ数が16倍であるため、計算量は(128×16)log2(128×16)/16log2128倍、即ち25倍になる。
(2)本実施の形態1の周波数分析装置で周波数帯域を16倍に広げる場合
フーリエ変換を行うデータ数自体は変化がない。そして、2つの信号列s(n)とsΔ(n)をフーリエ変換するため計算を2回行うためほぼ2倍である。ただし、(11)式の計算も行う必要があるが、これは高々N回(今の場合128回)であり、フーリエ変換の計算量と比べれば1/10程度であるため、計算量の増加にはほとんど影響しない。
このように、本実施の形態1の周波数分析装置を用いれば、約12倍の計算負荷が軽減できる。
【0043】
以上のように、本実施の形態1では、時間間隔τでサンプリングした信号列s(n)と、この信号列s(n)から微少時間Δτだけずらして時間間隔τでサンプリングした信号列sΔ(n)とをそれぞれ短時間フーリエ変換し、比較することでΔτだけずらしたことにより生じる位相差exp(2πjfΔτ)を算出する。この位相差によって、受信信号を広帯域化に伴って生じるエイリアジングを利用して(11)式によりrを求めることで、帯域フィルタ1から出力される周波数帯域全てを走査しなくても受信信号の存在する周波数を算出することができるため、広帯域化しても計算量の増大を抑え、エイリアジングの効果を排除し、正確な周波数分析結果を得ることができる。
【0044】
実施の形態2.
上記実施の形態1では、2つの信号列s(n)、sΔ(n)を得るために、一律にΔτでサンプリングしたデータ群の中から適当な時間間隔τとなるデータを抽出してデータを分配する例を示したが、例えば、図5のように、タイミングのΔτだけ異なるクロックでサンプリングを行う2つのA/D変換器を使用して実現してもよい。
【0045】
図6に本実施の形態2に係わる周波数分析装置の装置構成を表すブロック図を示す。図1と同じ構成要素には同じ符号を付す。10は信号分配器である。本実施の形態2の周波数分析装置は図のように2つのA/D変換器2a、2bを使用するが、無駄な信号をサンプリングすることがなくなるため、メモリ3a、3bの容量を小さくて済む。
【0046】
以下、動作について説明する。帯域フィルタ1は入力された受信信号を帯域制限して分析対象の周波数信号のみ選択する。信号分配器10では、信号を2つのA/D変換器2a、2bに分配する。各A/D変換器2a、2bではそれぞれ時間間隔τでサンプリングを行うが、例えば、A/D変換器2bの方がA/D変換器2aよりもΔτだけ遅れたタイミングでサンプリングを行う。これによって、実施の形態1で示したような2つの信号列s(n)およびsΔ(n)が生成され、順次メモリ3a、3bに記憶される。以下の動作は上記実施の形態1の周波数分析装置の動作と同様である。それぞれの信号列の短時間フーリエ変換STFT(ω)、STFTΔ(ω)を比較し、(11)式からexp(2πjfΔτ)を算出する。
【0047】
以上のように、本実施の形態2では、帯域制限した信号を予め分配してA/D変換器に入力するため、無駄な信号をサンプリングすることがなく、メモリの容量が小さくて済む。
【0048】
実施の形態3.
上記(11)式はそれぞれのpについてSTFT(ω)に含まれる信号が単一である場合のみ成立する。従って、上記実施の形態1および2では、受信信号の周波数分布が単一もしくは比較的疎な場合のみ分析が可能であった。本実施の形態3では、上記実施の形態1または2の周波数分析装置を用いた場合に互いにエイリアジング成分となるような複数の周波数成分が受信信号に含れていた場合、信号が複数のエイリアジング成分を含むか含まないかの判定を行い、含まないと判定された場合にのみ分析を行う周波数分析装置について説明する。
【0049】
図7は本実施の形態3の周波数分析装置が対象とする信号の周波数スペクトル図である。点線は実際に測定されるスペクトル、実線は信号の真の周波数成分を表す。図7のように上記実施の形態1および2の周波数分析装置を用いた場合に、互いにエイリアジング成分となるような独立した信号成分w1、w3が受信信号に含まれていた場合、その周波数をf1、f3とすると、
【数11】
Figure 0003646702
また、
【数12】
Figure 0003646702
これらの式をそれぞれフーリエ変換すると、
【数13】
Figure 0003646702
また、
【数14】
Figure 0003646702
である。
【0050】
実際の信号から求まるのはSTFTΔ(2πp/Nτ)およびSTFT(2πp/Nτ)の値のみである。従って、上記実施の形態1または2のようにSTFTΔ(2πp/Nτ)/STFT(2πp/Nτ)を求めても、w1、w3が混ざっているため正しいf1、f3の値を求めることができない。
【0051】
そのために、本実施の形態3では、図8に示したように周波数分析装置を、時間間隔τを持った一つの信号系列s(n)と、そのs(n)に対してそれぞれ2つの異なる微少な時間差Δτ1、Δτ2を持った信号系列sΔ1(n)、sΔ2(n)を生成する構成とし、それぞれ短時間フーリエ変換を行う。図において9aは信号が単一か否かを判定しながら真の周波数を算出するエイリアジング判定手段である。
【0052】
ここで、エイリアジングの周波数をf1=p/Nτ、f2=(p+N)/Nτ、f3=(p+2N)/Nτ、…fr+1=(p+rN)/Nτ、また各エイリアジング周波数における各信号強度をw1、w2、w3、…wrとする。すると、信号が複数含まれている場合を一般化して表現すると、
【数15】
Figure 0003646702
【数16】
Figure 0003646702
である。sΔ2(n)も(16)式のΔτ1を、それぞれΔτ2と置き換えたものと同様である。
【0053】
ここで、例えば、周波数f3に真の信号成分が1つのみあるとすると、w3のみが0以外の値を持ち、かつw3以外のw1、w2、w4…wrは全て0である。
【0054】
従って、その場合にs(n)、sΔ1(n)、sΔ2(n)をフーリエ変換すると、それぞれ、以下のようになる。
【数17】
Figure 0003646702
【0055】
逆に、信号成分が単一でなく、例えば、周波数f3の他に周波数f1に信号成分w1が存在していた場合は、そのフーリエ変換は(17)、(18)、(19)式の代わりに、
【数18】
Figure 0003646702
【0056】
ここで、各信号列間の単位時間当たりの位相差として以下の量を考える。
【数19】
Figure 0003646702
ただし、(20)式において∠は位相を表す。
例えば、信号成分が周波数f3に一つの場合、(17)〜(19)式の右辺を(21)式に代入すると、全て2πf3が得られるが、信号成分が周波数f3およびf1にあるような場合、(17a)〜(19a)式の右辺を(20)式に代入しても一定の値は得られない。そのため、この(20)式を用いて信号成分が複数あるか否かを判定する。
【0057】
図8は本実施の形態3に係わる周波数分析装置の構成を表すブロック図である。図1と同じ構成要素には同じ符号を付す。1つの信号系列s(n)およびそのs(n)から時間Δτ1およびΔτ2だけずれた2つの信号系列sΔ1(n)およびsΔ2(n)の計3つの信号系列をサンプリングし、それぞれ短時間フーリエ変換する。そこから、信号検出手段8bおよびエイリアジング判定手段9bによって、上記(20)式のような単位時間当たりの位相差を求め、信号が一つか否かを判定しながら真の周波数を求める。
【0058】
以下、図9に示したフローチャートに従って、信号検出手段8bおよびエイリアジング判定手段9bの動作を説明する。まず、信号検出手段8bで離散値の各周波数f=(p+rN)/Nτをp=0(ステップ200)から順番に信号が存在するかの判定が行われる(ステップ201)。そして、信号が未検出ならば、そのときのp(今の場合、p=0)に対して、全てのrについて受信信号の周波数成分が0である(ステップ202)。即ち、上記実施の形態1の場合と同様に、p=0で信号が未検出であれば、そのエイリアジングの可能性がある周波数f=N/Nτ,2N/Nτ,3N/Nτ,…についての周波数成分も0ということである。
【0059】
一方、p=0で信号が検出された場合は、STFT(ω)、STFTΔ1(ω)、STFTΔ2(ω)、STFTΔ3(ω)から(20)式を計算し(ステップ204)、全て同じ値を算出するか判定する(ステップ205)。もし、(21)式の計算結果が一致しないならば、上述したように信号成分が複数存在することになるので、当該周波数の成分は不定である(ステップ203)。一方、(20)式の結果が同じ値ならば信号は1つであり、周波数の分析が可能であるので実施の形態1と同様に(11)式を満足するrTを算出する(ステップ206)。そして、算出したrT以外のrに対し、周波数(p+rN)/Nτの成分が0である(ステップ207)。
【0060】
以上、いずれかの処理が終わったら周波数(p+rN)/Nτのpをp=p+1として、次の周波数成分について同様の処理を行う(ステップ208)。その際、p≧Nであれば、演算する部分が重なってしまうため、p≧Nであれば分析を終了し、p<Nの場合のみステップ201に戻る(ステップ209)。
【0061】
このように、p=0からN−1まで、即ち離散値の周波数fを0から1/τまで信号の有無を検出するだけで、全ての周波数を分析する必要がなく、計算量が大幅に軽減することができる。
【0062】
また、信号がいくつあるか判定しながら周波数分析を行うため、信号成分が複数ある場合でも、誤った周波数を算出しない。
【0063】
実施の形態4.
本実施の形態4では、それぞれ異なる時間だけずれた2つ以上の信号系列を短時間フーリエ変換して、その離散値の各周波数に含まれる周波数成分が1つか否かを判定し、複数の信号が含まれていてもエイリアジングの影響を除くことができる周波数分析装置について説明する。
【0064】
ここでは上記実施の形態3の図7のように互いのエイリアジング成分となるような周波数に独立した複数の信号がある場合に、3つ信号系列を短時間フーリエ変換する例を示す。
【0065】
図10は本実施の形態4の周波数分析装置の構成を表すブロック図である。図において、図8と同じ構成要素には同じ符号を付す。9bは本実施の形態3に係わるエイリアジング判定手段である。前記実施の形態3と同様に、信号系列s(n)に対し複数の時間差Δτ1、Δτ2、Δτ3を持った信号系列sΔ1(n)、sΔ2(n)、sΔ3(n)をサンプリングして生成し、短時間フーリエ変換を行う。ただし、ここで、Δτ1<Δτ2<Δτ3であって、図7のr/τ<1/Δτ1であるとする。これら4つの信号系列の短時間フーリエ変換後のスペクトルSTFT(2πp/Nτ)、STFTΔ1(2πp/Nτ)、STFTΔ2(2πp/Nτ)、STFTΔ3(2πp/Nτ)それぞれについて(11)式のような位相差を考える。
【0066】
ここで、エイリアジングの性質として、真の周波数から1/τごとにエイリアジング成分が生じることから、エイリアジングの生じる可能性のある周波数は、
【数20】
Figure 0003646702
のように表現でき、この各周波数における振幅をw1、w2、…wr+1とすれば、4つの信号系列は以下のような行列にて表現できる。
【数21】
Figure 0003646702
として、
X=S・A (22)
である。
【0067】
ここで、行列Xの各要素については4つの信号系列を実際に短時間フーリエ変換することで、その値が得られ、また、行列Sの各要素についても、上述のようにエイリアジング成分の性質から各周波数が定められるため、(21)各式中のpを固定すればその値が定まる。従って、Sの一般逆行列を用いて、
A=S−1X=(S*S)−1・S*・X (22a)
を計算することで、行列Aの各要素を求めることができる。
【0068】
エイリアジング判定手段9bでは、以上のような原理に基づいて周波数0から1/τまでのpについて信号の有無を判定し、信号がある場合にはそのpについて上記(23a)式によって行列Aの各要素を求める。
【0069】
この過程をフローチャートにしたものが図11である。以下、図に従って説明する。まず、(21)式のp=0から分析を開始する(ステップ300)。そして、信号の有無を確認するが(ステップ301)、例えば、信号がp=0のとき発見されなければ、(21)式の各周波数において振幅=0、即ち、信号はない(ステップ302)として、pを1だけ増やして(ステップ304)再びステップ301に戻って、信号の有無を確認する。一方、信号が発見された場合には、上述したように、(23a)式を用いて3つの信号系列から行列Aの各要素、即ち、各周波数における振幅を計算することで、信号成分が複数あっても正しく算出することができる。そして、算出できたらpを1だけ増やして(ステップ304)再び信号の分析を行う。そして、p=N−1まで分析したら、周波数0からrM/τまで分析が終わったことになるので、分析を終了する。
【0070】
このように、周波数を0から1/τまで分析するだけで、周波数0からrM/τまで分析を行うことになるため、信号を広帯域化しても効率のよい分析を行うことができる。
【0071】
また、エイリアジング成分と複数の真の周波数成分とが重なりあっている場合でも、互いに異なる微少時間だけずれた信号系列を3つ生成し、位相差を比較することで、複数の周波数成分を分離して正しい周波数を算出することができるため、正確な分析を行うことができる。
【0072】
以上の説明では、全ての信号成分を(22)式を用いて算出する例について示したが、(22)式は(11)式に比べて計算量が多くなるため、実施の形態3に示したように、信号を単一か否かを判定するステップを挿入し、信号が単一の場合には(11)式を用いて周波数を算出するようにしてもよい。
【0073】
図12に信号成分が単一か否かを判定しながら、信号成分を算出する場合のフローチャートを示す。図11と異なる点は、ステップ404にて実施の形態3で行ったように(20)式を用いて信号が単一か否かを判定する点である。信号成分が単一ならば実施の形態1と同様に算出し、複数ある場合は本実施の形態4で上述したように(22)式を用いて全ての信号成分を算出する。
【0074】
なお、本実施の形態4では、生成する信号系列を3つとしたため、周波数成分が2つ以内の場合に対応できるが、重み付け手段と、高速フーリエ変換手段を更に増やせば、さらに多くの周波数成分がエイリアジング成分と重なっても正確に分離することができる。
【0075】
以上のように、本実施の形態4では、生成した3つの信号系列を短時間フーリエ変換し、エイリアジング成分の性質を利用して行列の要素を決定し、位相差を比較して信号成分を算出するため、複数の信号が同じ周波数に重なっている場合でも、互いに分離して算出することができる。
【0076】
実施の形態5.
以上、実施の形態1〜4では、広帯域化に伴って発生するエイリアジング成分を排除しつつ、計算量の増大を抑え効率のよい周波数分析装置について説明したが、本実施の形態5では、窓関数の効果による周波数成分の拡散による誤差を排除しつつ、短時間のサンプリングで正確な周波数分析が可能な周波数分析装置について説明する。
【0077】
図13は本実施の形態5に係わる周波数分析装置の装置構成を表すブロック図である。図において図1と同じ構成要素には同じ符号を付す。11は畳込み効果判定手段である。
【0078】
以下、動作について説明する。帯域フィルタ1で信号の周波数帯を適当な帯域に制限し、A/D変換器2で適当な時間間隔τでサンプリングする。サンプリングしたデジタル信号はメモリ3に記憶し、データ分配手段7によって、データを重み付け手段4a、4bに分配し、一方の信号系列s(n)に対して、他方s(n)に対してΔτだけずれたsΔ(n)が生成するようにする。それぞれの重み付け手段4a、4bでは窓関数による重み付けを行い、高速フーリエ変換手段5a、5bにおいてフーリエ変換を行う。信号検出手段8ではフーリエ変換後の各周波数をチェックして信号の有無を確認し、そして、畳込み効果判定手段11では、このフーリエ変換後の周波数スペクトル同士の位相を比較することで、真の周波数を特定する。
【0079】
このように構成することで、短時間のサンプリングで畳込みによる周波数成分の拡散を抑えることができ、連続波のみならず、断続波の場合でも、正確な周波数スペクトルを得る。
【0080】
ここで、窓関数による畳込みの効果には以下に示すような特徴がある。即ち、(A)畳込みにより拡散する各周波数の信号強度の和は元の信号強度にほぼ一致する。(B)1つの信号成分(周波数f)から畳込みにより拡散した各周波数成分について、(11)式は同じ値を持つ。(C)ハニング窓関数を用いることで、1つの信号の畳込みによる拡散の影響は、真の周波数周辺の3つの周波数値に限定することができる。
【0081】
(A)は、図23および、(8)式の説明で示したように、畳込み演算の性質によるものである。この場合、周波数は離散値を取るため、近似的に一致する。また、(B)については、STFT、STFTΔとも同じ元の周波数成分を窓関数g(n)により拡散したものであるから相互の位相の関係は同じになり、即ち、位相を表す(10)式が異なる値を取ることはないからである。
【0082】
従って、これらから畳込みによる周波数分布の拡散を排除して正しい周波数分布を求めるには、検出された周波数の信号成分が畳込みによる拡散の結果であることを(10)式により判定し、畳込みによる拡散によると判定された場合には、拡散した各信号成分の信号強度を足し合せ、(10)式のより求められた周波数に足し合せた信号強度の信号成分があるとすればよい。
【0083】
なお、(C)で示したように、周波数の拡散はハニング窓関数を用いることで、真の周波数周辺の3つの周波数値に限定することができる。そのため、Δτについて(10)式のfが一意の値を持つためには、3つの周波数の幅の内で、(10)式左辺の位相が1回転しないような値に限定する必要がある。即ち、exp(2πjfΔτ)が、(p−1)/Nτ≦f<(p+1)/Nτで一意となるには、
Δτ<Nτ/2 (23)
であることが必要である。従って、本実施の形態5では、(23)式のようなΔτとなるよう、データ分配手段7を動作させる。
【0084】
図14に本実施の形態5の信号のサンプリング方法を表す。また、フーリエ変換後の周波数スペクトルを図15に示す。図15のように、本実施の形態5はエイリアジングの影響を排除するための前記実施の形態1から4までと異なり、Δτはτより長い時間間隔となっている。(23)式の条件からすれば、Δτはτより短くても問題はないが、その場合、あまり短く取りすぎると、雑音の影響が出てくるため、Δτの間隔は適用した各装置の特性に応じて定められる。
【0085】
図16に本実施の形態5の信号検出手段8と畳込み判定手段11の動作のフローチャートを示す。まず、p=0から分析を開始し(ステップ500)、信号の有無を検出する(ステップ501)。信号が検出された場合は、(10)式から左辺のfを算出する(ステップ502)。このfにおける、信号強度に|STFT(2πp/Nτ)|2を加算する(ステップ503)。これが周波数p/Nτにおける1サイクルであり、以上が終了したらpの値を1増やして(ステップ504)、再び信号の有無を確認してから(10)式を用いてfを計算する。
【0086】
上述したように、畳込みによって1の周波数が3つの周波数に拡散するが、(10)式のfは同じ値が出てくるため、同じfならば畳込みによる拡散で生じたスペクトルであるとして、ステップ502および503の手順によって畳込みの影響を除去する。このように従来例のように真の周波数周辺の周波数を3つ選択するのではなく、(10)式によってあらかじめ真の周波数を求めて、この周波数における信号強度を足し合せていくため、真の周波数が離散周波数値のほぼ中間にあっても、精度よく真の周波数を算出できる。これをp=0〜Nまで行って(ステップ505)、分析可能範囲をすべて分析する。
【0087】
以上のように、本実施の形態5では、サンプリングした信号から同時に2つの信号列を生成し、それらの位相を比較することで、畳込みによってスペクトルが拡散する効果を抑え、また、加算を主体にして信号強度を算出するため、正確な周波数分析を行う周波数分析装置を得る。
【0088】
実施の形態6.
上記実施の形態5では、2つの信号列を得るために、時間間隔τでサンプリングを行った結果を記録し、その中から、適当なタイミングのデータをそれぞれ抽出し、2つの信号列の時間差がΔτとなるようにしたものを示したが、本実施の形態6では、図17のように、タイミングのΔτだけ異なるクロックでサンプリングを行う2つのA/D変換器を使用して実現してもよい。
【0089】
図17に本実施の形態6に係わる周波数分析装置の装置構成を表すブロック図で示す。図6と同じ構成要素には同じ符号を付す。11は上記実施の形態5と同様の畳込み効果判定手段である。本実施の形態6の周波数分析装置は図のように2つのA/D変換器2a、2bを使用するが、無駄な信号をサンプリングすることがなくなるため、メモリ3a、3bの容量を小さくて済む。
【0090】
以下、動作について説明する。帯域フィルタ1は入力された受信信号を帯域制限して分析対象の周波数信号のみ選択する。信号分配器10では、信号を2つのA/D変換器2a、2bに分配する。各A/D変換器2a、2bではそれぞれ時間間隔τでサンプリングを行うが、例えば、A/D変換器2bの方がA/D変換器2aよりもΔτだけ遅れたタイミングでサンプリングを行う。これによって、実施の形態5で示したような2つの信号列s(n)およびsΔ(n)が生成され、順次メモリ3a、3bに記憶される。以下の動作は上記実施の形態5の周波数分析装置の動作と同様である。
【0091】
以上のように、本実施の形態6では、帯域制限した信号を予め分配してA/D変換器に入力するため、無駄な信号をサンプリングすることがなく、メモリの容量が小さくて済む。
【0092】
実施の形態7.
上記実施の形態5および6では、受信信号の周波数分布が単一もしくは比較的疎な場合のみ分析が可能であったが、本実施の形態7では同じ周波数に畳込みによって拡散した異なる信号成分が重なりあった場合でも、信号成分が単一か否かを判定し、信号成分が単一の場合のみ分析を行う周波数分析装置について説明する。
【0093】
図18は本実施の形態7の周波数分析装置の装置構成を表すブロック図である。本実施の形態7では、信号系列を3つ用いてるため、実施の形態3の図8と同じ構成要素には同じ符号を付す。11aは本実施の形態7の特徴である、各周波数における信号成分が単一か否かを判定しながら畳込みの効果を抑えて正しい周波数を算出する畳込み効果判定手段である。
【0094】
次に、動作について説明する。1つの信号系列s(n)およびそのs(n)から時間Δτ1およびΔτ2だけずれた2つの信号系列sΔ1(n)およびsΔ2(n)の計3つの信号系列をサンプリングし、それぞれ短時間フーリエ変換する。そこから、信号検出手段8および畳込み効果判定手段11aによって、上記実施の形態3と同様に(21)式のような単位時間当たりの位相差を求め、信号が一つか否かを判定しながら真の周波数を求める。
【0095】
以下、信号検出手段8および畳込み効果判定手段11aの動作については、図19に示したフローチャートにしたがって説明する。まず、信号検出手段8bで離散値の各周波数f=p/Nτをp=0(ステップ600)から順番に信号が存在するかの判定が行われる(ステップ601)。そして、信号が未検出ならば、p=p+1(ステップ606)とし、次のpについて再び信号が存在するか判定を行う。
【0096】
一方、p=0で信号が検出された場合は、STFT(ω)、STFTΔ1(ω)、STFTΔ2(ω)、STFTΔ3(ω)から(20)式を計算し(ステップ602)、全て同じ値を算出するか判定する(ステップ603)。もし、(21)式の計算結果が一致しないならば、上述したように信号成分が複数存在することになるので、当該周波数の成分は不定であるとして、次のpの値について信号の存在の有無を判定する。一方、(20)式がすべて同じ値ならば信号は1つであり、周波数の分析が可能であるので実施の形態5と同様に(11)式を満足する周波数fを算出する(ステップ604)。そして、算出した周波数fにおける信号強度に補正を加える(ステップ605)。
【0097】
以上、いずれかの処理が終わったら周波数p/Nτのpをp=p+1として、次の周波数成分について同様の処理を行う(ステップ606)。その際、p≧Nであれば、演算する部分が重なってしまうため、p≧Nであれば分析を終了し、p<Nの場合のみステップ501に戻る(ステップ607)。
【0098】
以上のように、本実施の形態7では、信号がいくつあるか判定しながら周波数分析を行うため、信号の周波数成分が複数ある場合でも、誤った周波数を算出せず、短時間のサンプリングで周波数分析が終了し、また、加算を主体にして信号強度を算出するため、正確な周波数分析を行うことが可能な周波数分析装置を得る。
【0099】
【発明の効果】
以上のように、本発明の周波数分析装置は、サンプリングしたデジタルデータから同時に2つの信号系列を生成し、この2つの信号系列を短時間フーリエ変換する。そして、それぞれの位相差を算出し、エイリアジングを利用して、周波数帯域全てを走査しなくても受信信号の存在する周波数を算出することができるため、分析すべき信号の周波数帯が広帯域化しても計算量の増大を抑え、エイリアジングの効果を排除し、正確な周波数分析結果を得ることができる。
【0100】
また、帯域制限した信号を予め分配してA/D変換器に入力するため、無駄な信号をサンプリングすることがなく、メモリの容量が小さくて済む。
【0101】
また、生成した3つの信号系列を短時間フーリエ変換し、各信号系列の単位時間当たりの位相差を算出して各周波数において信号成分が単一か否かを判定しながら分析を行うため、信号成分が複数ある場合でも、誤った周波数を算出しない。
【0102】
また、生成した3つの信号系列を短時間フーリエ変換し、エイリアジング成分の性質を利用して行列の要素を決定し、位相差を比較して信号成分を算出するため、複数の信号が同じ周波数に重なっている場合でも、互いに分離して算出することができる。
【0103】
また、加算によって信号強度を算出する畳込み効果判定手段を設けたため、窓関数の効果による周波数成分の拡散による誤差を排除し、正確な周波数分析結果を得る。
【図面の簡単な説明】
【図1】 本発明の実施の形態1に係わる周波数分析装置の構成を表すブロック図である。
【図2】 本発明の実施の形態1に係わる周波数分析装置のサンプリング方法を表す図である。
【図3】 本発明の実施の形態1に係わる周波数分析装置の周波数分析方法を表す図である。
【図4】 本発明の実施の形態1に係わる周波数分析装置の周波数算出方法のフローチャートである。
【図5】 本発明の実施の形態2に係わる周波数分析装置のサンプリング方法を表す図である。
【図6】 本発明の実施の形態2に係わる周波数分析装置の装置構成を表すブロック図である。
【図7】 本発明の実施の形態3が対象とする周波数分布を表す図である。
【図8】 本発明の実施の形態3に係わる周波数分析装置の装置構成を表すブロック図である。
【図9】 本発明の実施の形態3に係わる周波数分析装置の周波数算出方法のフローチャートである。
【図10】 本発明の実施の形態4に係わる周波数分析装置の装置構成を表すブロック図である。
【図11】 本発明の実施の形態4に係わる周波数分析装置の周波数算出方法のフローチャートである。
【図12】 本発明の実施の形態4に係わる周波数分析装置の別の周波数算出方法のフローチャートである。
【図13】 本発明の実施の形態5に係わる周波数分析装置の装置構成を表すブロック図である。
【図14】 本発明の実施の形態5に係わる周波数分析装置のサンプリング方法を表す図である。
【図15】 本発明の実施の形態5に係わる周波数分析結果の一例である。
【図16】 本発明の実施の形態5に係わる周波数分析装置の周波数算出方法のフローチャートである。
【図17】 本発明の実施の形態6に係わる周波数分析装置の装置構成を表すブロック図である。
【図18】 本発明の実施の形態7に係わる周波数分析装置の装置構成を表すブロック図である。
【図19】 本発明の実施の形態7に係わる周波数分析装置の周波数算出方法のフローチャートである。
【図20】 従来の周波数分析装置の装置構成を表すブロック図である。
【図21】 短時間フーリエ変換の原理を表す図である。(a)サンプリングにより生成された信号列である。(b)窓関数の形状を表す図である。(c)上記生成された信号列と窓関数とを畳込み演算によって重み付けした後の図である。(d)それぞれの時間における周波数分布を表す図である。(e)各時間における周波数分布を並べて時間変化を見る図である。
【図22】 エイリアジングの効果を表す図である。
【図23】 窓関数との畳込みの効果を表す図である。
【図24】 従来の周波数分析装置を表すブロック図である。
【符号の説明】
1 帯域フィルタ、 2、2a、2b A/D変換器、
3、3a、3b メモリ、 4a、4b 重み付け手段、
5a、5b 高速フーリエ変換手段、 6 表示器、
7 データ分配手段、 8 信号検出手段、
9、9a エイリアジング判定手段、 10 信号分配器、
11、11a 畳込み効果判定手段。

Claims (6)

  1. 周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、
    上記帯域制限されたアナログ入力信号を第一の時間間隔で順次デジタルデータ列に変換するAD変換器と、
    上記AD変換器によって変換された入力データ列を順次記憶するメモリと、
    上記メモリに記憶された入力データ列から上記第一の時間間隔よりも長い第二の時間間隔を持った第一の分配データ列、およびこの第一の分配データ列から上記第一の時間間隔ずつ遅れた第二の分配データ列とを抽出するデータ分配手段と、
    上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、
    この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、
    上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、
    この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、
    上記第一または第二の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、
    信号が検出された場合、上記第一および第二の周波数領域データの位相を比較して真の周波数を算出するエイリアジング判定手段とを備えることを特徴とする周波数分析装置。
  2. 周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、
    上記帯域制限されたアナログ入力信号を2つに分配する信号分配器と、
    上記分配された一方のアナログ入力信号を一定の時間間隔でデジタル信号に変換し、第一の分配データ列を生成する第一のA/D変換器と、
    上記分配された他方のアナログ入力信号を、上記第一のA/D変換器のデジタル変換より上記一定の時間間隔より短い時間間隔だけ遅れて、一定の時間間隔でデジタル信号に変換し、第二の分配データ列を生成する第二のA/D変換器と、上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、
    この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、
    上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、
    この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、
    上記第一または第二の周波数領域データを各周波数ごとに信号が存在するかを検出する信号検出手段と、
    信号が検出された場合、上記第一および第二の周波数領域データを比較することで、エイリアジングの影響を排除し、真の周波数を算出するエイリアジング判定手段とを備えることを特徴とする周波数分析装置。
  3. 周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、
    上記帯域制限されたアナログ入力信号を第一の時間間隔で順次デジタルデータ列に変換するAD変換器と、
    上記AD変換器によって変換された入力データ列を順次記憶するメモリと、
    上記メモリに記憶された入力データ列から上記第一の時間間隔よりも長い第二の時間間隔を持った第一の分配データ列と、この第一の分配データ列から上記第一の時間間隔ずつ遅れた第二の分配データ列と、この第一の分配データ列から上記第一および第二の時間間隔とは異なる第三の時間間隔ずつ遅れた第三の分配データ列とを抽出するデータ分配手段と、
    上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、
    この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、
    上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、
    この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、
    上記第三の分配データ列に所定の窓関数を乗算し、重み付けを行う第三の重み付け手段と、
    この重み付けされた第三の分配データ列をフーリエ変換し、第三の周波数領域データに変換する第三の高速フーリエ変換手段と、
    上記第一または第二または第三の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、
    信号が検出された場合、上記第一、第二、第三の周波数領域データのいずれか2つの単位時間当たりの位相差を比較して、信号が検出された周波数における信号成分が単一か否かを判定し、信号成分が単一の場合のみ真の周波数を算出するエイリアジング判定手段とを備えることを特徴とする周波数分析装置。
  4. 周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、
    上記帯域制限されたアナログ入力信号を第一の時間間隔で順次デジタルデータ列に変換するAD変換器と、
    上記AD変換器によって変換された入力データ列を順次記憶するメモリと、
    上記メモリに記憶された入力データ列から上記第一の時間間隔よりも長い第二の時間間隔を持った第一の分配データ列と、この第一の分配データ列から上記第一の時間間隔ずつ遅れた第二の分配データ列と、この第一の分配データ列から上記第一および第二の時間間隔とは異なる第三の時間間隔ずつ遅れた第三の分配データ列とを抽出するデータ分配手段と、
    上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、
    この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、
    上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、
    この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、
    上記第三の分配データ列に所定の窓関数を乗算し、重み付けを行う第三の重み付け手段と、
    この重み付けされた第三の分配データ列をフーリエ変換し、第三の周波数領域データに変換する第三の高速フーリエ変換手段と、
    上記第一または第二または第三の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、
    信号が検出された場合、上記第一、第二、第三の各周波数における周波数領域データを要素とする行列に、当該信号が検出された周波数から上記第二の時間間隔の逆数の整数倍ずつずれた複数の周波数を要素とする行列の逆行列を乗算することで真の周波数を算出するエイリアジング判定手段とを備えることを特徴とする周波数分析装置。
  5. 周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、
    上記帯域制限されたアナログ入力信号を第一の時間間隔であらかじめ決められた一定回数だけ順次デジタルデータに変換するAD変換器と、
    上記AD変換器によって変換された入力データを順次記憶するメモリと、
    上記メモリに記憶された入力データから上記第一の時間間隔を持ったデータを抽出してなる第一の分配データ列、および上記第一の時間間隔に上記デジタル変換した一定の回数を掛けた数の二分の一より短い第二の時間間隔ずつ上記第一の分配データ列より遅れた第二の分配データ列とを抽出するデータ分配手段と、
    上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、
    この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、
    上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、
    この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、
    上記第一または第二の周波数領域データの各周波数ごとに信号が存在するかを検出する信号検出手段と、
    信号が検出された場合、上記第一および第二の周波数領域データの位相を比較して、真の周波数を算出し、算出した周波数において当該信号強度を足し合せていく処理を行う畳込み効果判定手段とを備えることを特徴とする周波数分析装置。
  6. 周波数分析すべきアナログ入力信号を一定の周波数帯域に制限する帯域フィルタと、
    上記帯域制限されたアナログ入力信号を第一の時間間隔であらかじめ定められた一定回数だけ順次デジタルデータ列に変換するAD変換器と、
    上記AD変換器によって変換された入力データ列を順次記憶するメモリと、
    上記メモリに記憶された入力データから上記第一の時間間隔を持ったデータを抽出してなる第一の分配データ列、および上記第一の時間間隔に上記デジタル変換した一定の回数を掛けた数の二分の一より短い第二の時間間隔ずつ上記第一の分配データ列より遅れた第二の分配データ列と、上記第一の時間間隔に上記デジタル変換した一定の回数を掛けた数の二分の一より短かく、かつ上記第一および第二の時間間隔とは異なる第三の時間間隔ずつ上記第一の分配データ列より遅れた第三の分配データ列とを抽出するデータ分配手段と、
    上記第一の分配データ列に所定の窓関数を乗算し、重み付けを行う第一の重み付け手段と、
    この重み付けされた第一の分配データ列をフーリエ変換し、第一の周波数領域データに変換する第一の高速フーリエ変換手段と、
    上記第二の分配データ列に所定の窓関数を乗算し、重み付けを行う第二の重み付け手段と、
    この重み付けされた第二の分配データ列をフーリエ変換し、第二の周波数領域データに変換する第二の高速フーリエ変換手段と、
    上記第三の分配データ列に所定の窓関数を乗算し、重み付けを行う第三の重み付け手段と、
    この重み付けされた第三の分配データ列をフーリエ変換し、第三の周波数領域データに変換する第三の高速フーリエ変換手段と、
    上記第一または第二または第三の周波数領域データいずれかの各周波数ごとに信号が存在するかを検出する信号検出手段と、
    信号が検出された場合、上記第一、第二、第三の周波数領域データのいずれか2つの単位時間当たりの位相差を比較して、信号が検出された周波数における信号成分が単一か否かを判定し、信号成分が単一の場合のみ上記第一および第二の周波数領域データの位相を比較して、真の周波数を算出し、算出した周波数において当該信号強度を足し合せていく処理を行う畳込み効果判定手段とを備えることを特徴とする周波数分析装置。
JP2002023783A 2002-01-31 2002-01-31 周波数分析装置 Expired - Fee Related JP3646702B2 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2002023783A JP3646702B2 (ja) 2002-01-31 2002-01-31 周波数分析装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2002023783A JP3646702B2 (ja) 2002-01-31 2002-01-31 周波数分析装置

Publications (2)

Publication Number Publication Date
JP2003222646A JP2003222646A (ja) 2003-08-08
JP3646702B2 true JP3646702B2 (ja) 2005-05-11

Family

ID=27746394

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2002023783A Expired - Fee Related JP3646702B2 (ja) 2002-01-31 2002-01-31 周波数分析装置

Country Status (1)

Country Link
JP (1) JP3646702B2 (ja)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2006064549A (ja) * 2004-08-27 2006-03-09 Nippon Telegr & Teleph Corp <Ntt> スペクトル解析方法、スペクトル解析装置、およびスペクトル解析プログラム
CN114157322B (zh) * 2021-11-16 2023-03-21 山东轻工职业学院 基于加权分数阶傅里叶变换的低截获信号产生方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2904425B2 (ja) * 1991-07-17 1999-06-14 株式会社アドバンテスト スペクトラム・アナライザ
JP3163573B2 (ja) * 1992-11-20 2001-05-08 株式会社アドバンテスト 高分解能周波数分析装置及びこの装置を用いたホログラム観測装置、ベクトルスペクトル解析装置
JP3478300B2 (ja) * 1994-01-31 2003-12-15 株式会社アドバンテスト ジッタ周波数成分の検出方法
JP3305496B2 (ja) * 1994-04-28 2002-07-22 株式会社アドバンテスト 高分解能周波数分析装置を用いたドップラ補償装置
JPH10213613A (ja) * 1996-11-29 1998-08-11 Anritsu Corp 周波数測定装置
JP2000055949A (ja) * 1998-08-10 2000-02-25 Hitachi Building Systems Co Ltd 周波数分析方法及び周波数分析装置

Also Published As

Publication number Publication date
JP2003222646A (ja) 2003-08-08

Similar Documents

Publication Publication Date Title
KR102065603B1 (ko) 우세 신호 검출 방법 및 장치
KR101376556B1 (ko) 사이클로스테이션너리 툴박스를 이용하여 잡음에 삽입된텔레비전 신호의 존재 검출
US5576978A (en) High resolution frequency analyzer and vector spectrum analyzer
CN110784222A (zh) Adc输出曲线的生成方法、装置、设备及介质
CN101558318B (zh) 信号分析器和频域数据产生方法
JP5448452B2 (ja) スペクトル・トレースを発生するデータ圧縮
WO2002073222A1 (fr) Procede d&#39;analyse de frequence, appareil d&#39;analyse de frequence et analyseur de spectre
CN111812404B (zh) 一种信号处理方法以及处理装置
CN107478883A (zh) 一种实现任意n倍等效采样的方法和装置
JPH09243679A (ja) 任意区間波形を用いた非調和的周波数分析法
CN118688777B (zh) 一种应用于干扰信道场景的多载波相位测距方法
US7994959B2 (en) System and method of signal sensing, sampling and processing through the exploitation of channel mismatch effects
Wu et al. Time domain averaging based on fractional delay filter
US20080222228A1 (en) Bank of cascadable digital filters, and reception circuit including such a bank of cascaded filters
JP3163573B2 (ja) 高分解能周波数分析装置及びこの装置を用いたホログラム観測装置、ベクトルスペクトル解析装置
JP5035815B2 (ja) 周波数測定装置
JP2003222646A (ja) 周波数分析装置
JP2000055949A (ja) 周波数分析方法及び周波数分析装置
Greitans Time-frequency representation based chirp-like signal analysis using multiple level crossings
JP3599994B2 (ja) 電波諸元測定装置
JP5550203B2 (ja) 波形分析装置
JP3650767B2 (ja) ジッタ測定装置、ジッタ測定方法、及び試験装置
JPH11251969A (ja) 周波数ホッピングスペクトラム拡散方式の受信装置
JP4344356B2 (ja) 検波装置、方法、プログラム、記録媒体
US12546807B2 (en) Spectrum analyzer, system and method for outputting data from a spectrum analyzer

Legal Events

Date Code Title Description
RD01 Notification of change of attorney

Free format text: JAPANESE INTERMEDIATE CODE: A7421

Effective date: 20040709

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20041222

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

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20050118

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20050131

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20080218

Year of fee payment: 3

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20090218

Year of fee payment: 4

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20100218

Year of fee payment: 5

LAPS Cancellation because of no payment of annual fees