JP2005181168A - 波数推定装置および波数推定方法 - Google Patents

Abstract

【課題】 到来電波の波数を高精度に推定する。
【解決手段】 本発明は、波数推定装置は、到来電波を受信する複数のアンテナ素子１と、複数のA/D変換部２と、第１の相関行列計算部３と、フーリエ変換部４と、信号周波数帯域抑圧部５と、逆フーリエ変換部６と、第２の相関行列計算部７と、コレスキー分解計算部８と、白色化計算部９と、固有値計算部１０と、到来波数推定部１１と、到来方向推定部１２とを有する。信号周波数帯域を抑圧した信号成分に基づいて白色化を行って固有値を計算するため、到来波分布が不確定な空電雑音のようなアンテナ素子１間で有相関な雑音が存在する環境においても、高精度に到来波数を推定することができる。
【選択図】 図１

Description

本発明は、到来波数を推定する波数推定装置および波数推定方法に関する。

アレーアンテナに到来する電波の方向を推定する到来方向推定装置が知られており、到来方向を推定するための種々のアルゴリズムが提案されている。これらのアルゴリズムの中には、到来波数の情報を用いるアルゴリズムが数多く存在する。到来波数の情報を用いるアルゴリズムを採用すると、到来波数が実際の到来数と異なる場合に到来方向推定結果の精度が劣化してしまう。従って、到来波数を高精度に推定することが重要となる。

到来波数を推定するアルゴリズムのうち代表的な方法として、AIC法とMDL法が知られている。これらのアルゴリズムは、各アンテナ素子の受信信号間で無相関な熱雑音を前提としている。熱雑音は、回路内部で発生すると考えられるので、各アンテナ素子の受信信号間で無相関となる。

ここで、到来電波の周波数が約20MHz以下の場合には、空電雑音が熱雑音よりも強くなる。空電雑音とは、外来雑音の一種であり、雷を波源とする電磁波ノイズである。空電雑音は広帯域に存在しているが、周波数によって強度が異なる。周波数が約20MHz以下では、空電雑音が熱雑音よりも強くなり、雑音として空電雑音を考慮に入れる必要がある。

空電雑音が存在する場合には、各アンテナ素子の受信信号間の雑音は有相関となる。その理由は、遠方から到来する空電雑音は全てのアンテナで共通して受信され、各アンテナ素子の受信信号の中に共通の空電雑音が含まれるためである。従って、従来のように、アンテナ素子間で無相関な熱雑音を前提とした到来波数推定方法は適用できず、無理に適用すると、推定精度が大幅に劣化するという問題がある。

この問題を解決する従来の手法として、白色化を行う手法が提案されている（非特許文献１）。白色化とは、アンテナ素子間の受信信号の有相関雑音を無相関にする演算である。ここで、白色化を行うためには、空電雑音のみの受信信号から計算される相関行列が必要である。
Junhao Xie, Yeshu Yuan, and Yongtan Liu, "Super-resolution processing for HF surface wave radar based on pre-whitened MUSIC", IEEE Journal of Oceanic Engineering, vol. 23, no. 4, pp. 313-321, October 1998.

しかしながら、空電雑音と方向推定或いは波数推定を行いたい到来信号は共に受信信号の中に含まれているので、空電雑音のみの受信信号を取得するのは難しい。非特許文献１では、空電雑音の受信信号を解析的に導出している。導出するにあたっては、空電雑音の到来波分布を仮定して、アンテナパターンの利得を考慮して計算している。ところが、この方法では、空電雑音の到来波分布が既知でない場合には適用できない。または、仮定した空電雑音の到来波分布が実際の到来波分布と異なる場合には精度が劣化してしまう。空電雑音は雷を波源としているので、発生頻度は安定していない。従って、精度良く空電雑音の到来波分布を予測することは困難である。従って、従来の方法では、安定した白色化を行うことができず、到来波数推定精度が劣化し、その結果、方向推定結果も劣化してしまうという問題点がある。

このように、非特許文献１等の従来技術では、空電雑音の到来波分布を仮定できないと適用できないために、任意の空電雑音に対して白色化を適用できないという問題があった。また、厳密に空電雑音の到来波分布を推定することは困難であり、白色化の効果が小さく、到来波数推定の精度劣化と到来方向推定の精度が劣化するという問題があった。

その他、各アンテナ素子のキャリブレーション、波数推定および方向推定に関係する演算の高速化、装置の小型化、装置の低コスト化などが課題となっている。

本発明は、このような問題点に鑑みてなされたものであり、その目的は、到来信号の波数を高精度に推定可能な波数推定装置および波数推定方法を提供することにある。

上記の目的を達成するため、本発明は、到来信号を受信する複数のアンテナ素子と、前記複数のアンテナ素子の受信信号間の相関度を表す第１の相関行列を計算する第１の相関行列計算手段と、前記複数のアンテナ素子の受信信号に含まれる所定の周波数帯域成分の振幅を小さくする帯域抑圧手段と、前記帯域抑圧手段の出力信号に基づいて、前記第１の相関行列を白色化した白色化行列を計算する白色化計算手段と、前記白色化行列に基づいて、到来信号の波数を推定する波数推定手段と、を備える。

本発明によれば、受信信号に含まれる所定の周波数帯域成分を抑圧した上で白色化を行って到来信号の波数を推定するため、アンテナ素子の受信信号間で有相関な空電雑音を無相関化でき、到来波数を高精度に推定できる。

以下、図面を参照しながら本発明の一実施形態について説明する。

（第１の実施形態）
図１は本発明に係る波数推定装置の概略構成の一例を示すブロック図である。図１の波数推定装置は、到来電波を受信する複数のアンテナ素子１と、複数のA/D変換部２と、第１の相関行列計算部３と、複数のフーリエ変換部４と、複数の信号周波数帯域抑圧部５と、複数の逆フーリエ変換部６と、第２の相関行列計算部７と、コレスキー分解計算部８と、白色化計算部９と、固有値計算部１０と、到来波数推定部１１と、到来方向推定部１２とを備えている。

複数のA/D変換部２は、複数のアンテナ素子１の受信信号をデジタル信号に変換する。第１の相関行列計算部３は、複数のA/D変換部２の出力信号から第１の相関行列を計算する。複数のフーリエ変換部４は、複数のA/D変換部２の出力信号を周波数領域の信号に変換する。複数の信号周波数帯域抑圧部５は、複数のフーリエ変換部４の出力の中で所定のしきい値以上の振幅をもつ周波数帯域の振幅をゼロにする。複数の逆フーリエ変換部６は、複数の信号周波数帯域抑圧部５の出力を時間領域の信号に変換する。第２の相関行列計算部７は、複数の逆フーリエ変換部６の出力に基づいて第２の相関行列を計算する。コレスキー分解計算部８は、第２の相関行列をコレスキー分解する。白色化計算部９は、コレスキー分解の結果に基づいて、第１の相関行列を白色化して新たな行列（白色化行列）を計算する。固有値計算部１０は、白色化行列の固有値を計算する。到来波数推定部１１は、この固有値に基づいて、到来信号の波数を推定する。到来方向推定部１２は、推定された到来信号の波数を用いて到来方向推定を行う。

複数のアンテナ素子１は、必要とする周波数で到来する電波を受信するように設計されている。ダイポールアンテナ、モノポールアンテナおよびパッチアンテナ等のいずれのアンテナでもよい。複数のアンテナ素子１は、方向推定および波数推定を行いたい到来信号と雷を波源とする空電雑音の両方を受信可能である。すなわち、所望の周波数で到来してくる電磁波のすべてを受信する。

図２は図１の波数推定装置の処理手順の一例を示すフローチャートであり、以下、このフローチャートに基づいて、本実施形態の処理手順を説明する。アンテナ素子１で受信した受信信号は、A/D変換部２でアナログ信号からデジタル信号に変換される。受信信号をA/D変換部２に入力する前に、不図示の増幅器、周波数変換器およびフィルタを用いて、受信信号の増幅、周波数変換および周波数帯域の制限をそれぞれ行ってもよい。

まず、第１の相関行列計算部３は、デジタル化された各アンテナ素子１の受信信号を用いて第１の相関行列Ｒ_xxを計算する（ステップＳ１）。第１の相関行列Ｒ_xxは、（１）式により計算される。

Ｒ_xx＝Ｅ［ｘ（ｔ）ｘ^H（ｔ）］ …（１）
ここで、アンテナ素子数をｐ、アンテナ素子１の受信信号をｘi（ｔ）（ｉ＝１、２、…、ｐ）とすると、ｘ（ｔ）＝［ｘ1（ｔ）、ｘ2（ｔ）、…、ｘp（ｔ）］は受信信号ｘi（ｔ）（ｉ＝１、２、…、ｐ）のベクトルである。Ｅ［］はアンサンブル平均であり、有限の時間サンプル数Ｎの平均である。Ｈは複素共役配置、Ｔは転置である。

次に、複数のフーリエ変換部４は、各アンテナ素子１で受信した受信信号のデジタル信号を周波数領域に変換する（ステップＳ２）。ここで、フーリエ変換を行う時に用いる受信信号は、第１の相関行列計算部３に入力される受信信号と同じでよい。複数のフーリエ変換部４でFFT（Fast Fourier Transformer）を行えば、より高速に計算することができる。

図３は受信信号をフーリエ変換して得られる周波数スペクトラム図であり、中心周波数に対して±５kHzの範囲を示している。中心周波数付近に到来信号が存在し、空電雑音が広帯域に存在していることがわかる。

図２のステップＳ２の処理が終了すると、次に、複数の信号周波数帯域抑圧部５は、フーリエ変換部４の変換結果から、ある閾値以上の振幅を有する帯域の振幅を小さくする計算を行う。本実施形態では、例えば振幅をゼロにする計算を行う（ステップＳ３）。一般に、空電雑音に比べて到来信号の強度は強く、フーリエ変換した結果の最大振幅は到来信号と考えることができる。従って、この計算を行うことで、受信信号の中から到来信号を取り除くことができる。そして、空電雑音のデータのみを取り出すことができる。

なお、閾値は、最大振幅を基準として、その10分の1や100分の1のように設定してもよい。この場合、到来信号の電力が既知である場合に有効である。あるいは、空電雑音の推定電力から決めることもできる。現在、周波数ごとの空電雑音の電力分布は調べられており、空電雑音電力を閾値として、それ以上を到来信号成分と判別することもできる。

図４は信号周波数帯域抑圧部５で信号成分を抑圧した周波数スペクトラム図である。図４は、振幅の最大値の130分の１をしきい値として信号成分を抑圧する例を示している。

図２のステップＳ３の処理が終了すると、次に、逆フーリエ変換部６は、到来信号が抑圧された周波数領域の受信信号を時間領域に変換する（ステップＳ４）。逆FFTを行うことで、高速に演算を行うこともできる。

次に、第２の相関行列計算部７は、時間領域に変換された複数の受信信号を用いて、第２の相関行列を計算する（ステップＳ５）。計算方法は、第１の相関行列と同じ（１）式を用いる。この結果、空電雑音のみの相関行列Ｒ_aaを計算できる。

次に、コレスキー分解計算部８で、相関行列を（２）式のようにコレスキー分解する（ステップＳ６）。

Ｒ_aa＝Ｃ^HＣ …（２）
白色化計算部９は、第１の相関行列計算部３で得られた相関行列とコレスキー分解計算部８で計算された行列を用いて、以下の（３）式にて白色化計算を行って、白色化行列Ｒ_whiteを計算する（ステップＳ７）。

Ｒ_white＝Ｃ^-HＲ_xxＣ^-1 …（３）
この計算を行うことで、各アンテナ素子１の受信信号間で有相関な空電雑音を無相関化できる。

次に、固有値計算部１０は、相関行列Ｒ_xxの固有値を計算する（ステップＳ８）。アンテナ素子１数がｐの場合、相関行列はｐ×ｐの行列であり、ｐ個の固有値が計算される。計算された固有値を大きい順に、λ1≧λ2≧…≧λpとする。固有値は、電力に比例する値であり、ｎ（ｎは１以上の整数）波の電波が到来すると、固有値にはｎ個のピークが現れる。このように、固有値を検出することで、到来波数を精度よく推定できる。

図５は白色化を行わずに第１の相関行列から固有値を計算した結果を示す図である。図５は、４つの固有値のうち、１つだけが大きく、３つが小さい例を示している。図５の固有値を用いてAIC法とMDL法で波数推定した結果は２波となった。実際には、１波で測定を行ったため、精度の悪い結果となっている。図５の４つの固有値ｃ１，ｃ２，ｃ３，ｃ４の具体的な値はそれぞれ、1.457E-11、2.818E-16、1.401E-16、1.349E-16である。

図６は図２の手順で白色化を行った後に固有値を計算した結果を示す図である。図６の４つの固有値ｃ１，ｃ２，ｃ３，ｃ４の具体的な値はそれぞれ、1.519E0、2.159E-3、2.100E-3、2.041E-3である。

図６の固有値が図５の固有値と異なるのは、白色化を行ったためである。図５と図６のいずれも、相対的に大きな固有値が１つで、相対的に小さな固有値が３つである。しかしながら、小さな３つの固有値同士の大きさを比較すると、図５と図６で異なる。白色化しない場合には、これら３つの固有値間に大きさのばらつきがあるのに対し、白色化を行うと、大きさのばらつきがなくなる。

これら３つの固有値は、雑音電力に比例する値である。熱雑音のように、アンテナ素子１間で無相関な雑音の場合は、これら３つの固有値は同じ値になる性質がある。また、空電雑音のように、アンテナ素子１間で有相関の場合には、これら３つの固有値は異なる値になる性質がある。したがって、図６によれば、白色化を行った結果、これら３つの固有値にばらつきがなくなり、雑音が無相関化されたことがわかる。

図２のステップＳ８の処理が終了すると、次に、到来波数推定部１１は、計算されたｐ個の固有値を用いて到来波推定を行う（ステップＳ９）。波数推定は、AIC法やMDL法を用いることができる（Mati Wax、 and Thomas Kailath、 “Detection of signals by information theoretic criteria”, IEEE Trans. Acoustic, Speech, and Signal Processing, vol. ASSP-33, no. 2, pp. 387-392, April 1985参照）。

AIC法とMDL法では、それぞれ（４）式および（５）式の評価関数を計算する。

評価関数は、到来波数ｋの関数となっていて、その関数の最小値から到来波数ｋを推定する。

上述した図６の結果を用いて、AIC法やMDL法にて波数推定を行うと、精度よく１波と推定された。

図２のステップＳ９の処理が終了すると、次に、到来方向推定部１２は、得られた到来波数情報を用いて、その波数個の到来方向を推定する（ステップＳ１０）。高精度に波数が推定されていると、到来方向も高精度に推定することができる。

このように、第１の実施形態によれば、信号周波数帯域を抑圧した信号成分に基づいて白色化を行って固有値を計算するため、到来波分布が不確定な空電雑音のようなアンテナ素子間で有相関な雑音が存在する環境においても、高精度に到来波数を推定することができる。

また、到来方向を推定する際にも到来波数が必要となり、特に、MUSIC法、ESPRIT法、MODE法および最尤推定法などの高分解能の到来方向推定アルゴリズムでは、波数が既知であることが高精度化の条件となるが、本実施形態の手法で到来波数を推定すれば、結果として、到来方向も高精度に推定することができる。

上述した実施形態では、波数推定アルゴリズムにAIC法とMDL法を用いる例を説明したが、他のアルゴリズムを用いても、同様に高精度に波数推定を行うことができる。

（第２の実施形態）
第２の実施形態は、到来信号の最大振幅を有する周波数成分のみを信号周波数帯域抑圧部５で抑圧するものである。

第２の実施形態による波数推定装置は、図１と同様のブロック構成を備えているため、構成の詳細な説明を省略する。第２の実施形態では、広帯域の到来信号が受信される場合を想定しており、信号周波数帯域抑圧部５の処理動作が第１の実施形態とは異なっている。

図７は本実施形態で想定する到来信号の周波数スペクトラム図である。図７の到来信号の周波数帯域は、図７の矢印ｙ１で示した範囲である。仮に、この範囲を信号周波数帯域抑圧部５で抑圧すると、図８に示すように、到来信号だけでなく、空電雑音成分も大幅に抑圧されてしまう。この結果、白色化の効果が小さくなる。

そこで、本実施形態では、図７のような広帯域の到来信号が受信された場合には、図９に示すように、到来信号の最大振幅を取る周波数成分のみを抑圧する。この結果、ある一部分の周波数帯域以外の受信信号と空電雑音との合成信号になる。したがって、この合成信号を用いて白色化を行うと、ある一部分の周波数帯域以外の空電雑音が白色化され、同時に、ある一部分の周波数帯域以外の受信信号も白色化される。白色化されると、各受信信号は、アンテナ素子１間で無相関な熱雑音のように変換される。したがって、ある一部分の周波数帯域のみの到来信号を受信できるとともに、その帯域以外では熱雑音が存在しているかのように、受信信号が変換される。

このように、第２の実施形態によれば、到来信号の最大振幅の周波数成分のみを信号周波数帯域抑圧部５で抑圧するため、広帯域な到来信号が受信された場合でも、空電雑音を的確に白色化することができ、到来波数を高精度に推定できる。また、到来波数を用いる到来方向推定アルゴリズムも高精度化される。

なお、本実施形態では、最大振幅を有する周波数成分のみを抑圧しているが、信号電力が小さい場合には、雑音電力よりも信号電力の方が小さくなり、波数推定の精度が劣化する場合がある。この場合には、信号周波数帯域抑圧部５で抑圧する周波数に微小な周波数幅を設定することで、到来波数推定精度の改善が可能になる。

（第３の実施形態）
第３の実施形態では、到来信号中に複数のピークがある場合に、１つのピークの周波数成分のみを抑圧するものである。

第３の実施形態による波数推定装置は、図１と同様のブロック構成を備えているため、詳細な説明を省略する。第３の実施形態では、信号周波数帯域抑圧部５の処理動作が第１の実施形態とは異なっている。

第３の実施形態では、図１０に示すように、到来信号が複数のスペクトラルで構成されている例を考える。なお、説明を簡略化するために、２つのスペクトラルがあり、かつアンテナ素子１数は４とする。

第１の実施形態の手法で信号周波数帯域を抑圧すると、２つのスペクトラル内に存在する到来波の波数を推定することになる。このとき、アレーアンテナで推定可能な波数の上限は（素子数−１）波である。したがって、アンテナ素子１数が４の場合は、最大３波までしか推定できない。

ここで、２つのスペクトル内に存在する到来波数の合計が４であり、各スペクトラル内に２波ずつある場合を想定する。この場合、白色化を実施して空電雑音が無相関化されたとしても、アレーアンテナの素子数から来る制限のために、波数推定精度は劣化してしまう。

本実施形態では、２つのスペクトラムが存在する場合には、両方のスペクトラムを抑圧せずに、１つのスペクトラムのみを抑圧する。その結果、抑圧されたスペクトラムの到来信号のみが残り、他の到来信号は白色化され、熱雑音のように変換される。したがって、あるスペクトラムに対応する到来信号のみの到来波数の推定を行うことになる。

このように、第３の実施形態では、周波数スペクトラムに複数のピークが存在する場合には、いずれか１つのピークに対応する周波数帯域のみを抑圧して白色化することにより、各ピーク内の到来波数を順番に推定することができる。したがって、到来信号数が（アンテナ素子数−１）より多い場合でも、到来波数を精度よく推定することができる。

本発明に係る波数推定装置の概略構成の一例を示すブロック図。 図１の波数推定装置の処理手順の一例を示すフローチャート。 受信信号をフーリエ変換して得られる周波数スペクトラム図。 信号周波数帯域抑圧部５で信号成分を抑圧した周波数スペクトラム図。 白色化を行わずに第１の相関行列から固有値を計算した結果を示す図。 図２の手順で白色化を行った後に固有値を計算した結果を示す図。 本実施形態で想定する到来信号の周波数スペクトラム図。 図７から所定の周波数帯域を抑圧した周波数スペクトラム図。 図７から最大振幅周波数成分のみを抑圧した周波数スペクトラム図。 到来信号中に複数のピークがある場合の周波数スペクトラム図。

符号の説明

１ アンテナ素子
２ A/D変換部
３ 第１の相関行列計算部
４ フーリエ変換部
５ 信号周波数帯域抑圧部
６ 逆フーリエ変換部
７ 第２の相関行列計算部
８ コレスキー分解計算部
９ 白色化計算部
１０ 固有値計算部
１１ 到来波数推定部
１２ 到来方向推定部

Claims (7)

1. 到来信号を受信する複数のアンテナ素子と、
   前記複数のアンテナ素子の受信信号間の相関度を表す第１の相関行列を計算する第１の相関行列計算手段と、
   前記複数のアンテナ素子の受信信号に含まれる所定の周波数帯域成分の振幅を小さくする帯域抑圧手段と、
   前記帯域抑圧手段の出力信号に基づいて、前記第１の相関行列を白色化した白色化行列を計算する白色化計算手段と、
   前記白色化行列に基づいて、到来信号の波数を推定する波数推定手段と、を備えることを特徴とする波数推定装置。
2. 前記複数のアンテナ素子の受信信号をデジタル信号に変換する複数のＡ／Ｄ変換手段と、
   前記デジタル信号を周波数領域の信号に変換するフーリエ変換手段と、を備え、
   前記第１の相関行列計算手段は、前記デジタル信号に基づいて前記第１の相関行列を計算し、
   前記帯域抑圧手段は、前記フーリエ変換手段の出力信号に含まれる、所定のしきい値以上の振幅をもつ前記所定の周波数帯域成分の振幅をゼロにすることを特徴とする請求項１に記載の波数推定装置。
3. 前記帯域抑圧手段の出力信号を時間領域の信号に変換する逆フーリエ変換手段と、
   前記逆フーリエ変換手段の出力信号間の相関度を表す第２の相関行列を計算する第２の相関行列計算手段と、
   前記第２の相関行列をコレスキー分解するコレスキー分解手段と、を備え、
   前記白色化計算手段は、前記コレスキー分解手段の出力に基づいて、前記第１の相関行列を白色化した前記白色化行列を計算することを特徴とする請求項２に記載の波数推定装置。
4. 前記白色化行列の固有値を計算する固有値計算手段を備え、
   前記波数推定手段は、前記固有値に基づいて到来信号の波数を推定することを特徴とする請求項１乃至３のいずれかに記載の波数推定装置。
5. 前記帯域抑圧手段は、前記複数のアンテナ素子の受信信号のうち最大振幅を持つ周波数成分のみの振幅をゼロにすることを特徴とする請求項１乃至４のいずれかに記載の波数推定装置。
6. 前記帯域抑圧手段は、前記複数のアンテナ素子の受信信号に複数のピークが存在する場合には、いずれか１つのピークの振幅をゼロにすることを特徴とする請求項１乃至４のいずれかに記載の波数推定装置。
7. 到来信号を受信する複数のアンテナ素子の受信信号間の相関度を表す第１の相関行列を計算し、
   前記複数のアンテナ素子の受信信号に含まれる所定の周波数帯域成分の振幅を小さくする帯域抑圧を行い、
   帯域抑圧後の信号に基づいて、前記第１の相関行列を白色化した白色化行列を計算し、
   前記白色化行列に基づいて、到来信号の波数を推定することを特徴とする波数推定方法。

Images