JP4808984B2 - 到来波の方向推定装置 - Google Patents

到来波の方向推定装置 Download PDF

Info

Publication number
JP4808984B2
JP4808984B2 JP2005100314A JP2005100314A JP4808984B2 JP 4808984 B2 JP4808984 B2 JP 4808984B2 JP 2005100314 A JP2005100314 A JP 2005100314A JP 2005100314 A JP2005100314 A JP 2005100314A JP 4808984 B2 JP4808984 B2 JP 4808984B2
Authority
JP
Japan
Prior art keywords
calculation
matrix
householder
conversion calculation
conversion
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.)
Active
Application number
JP2005100314A
Other languages
English (en)
Other versions
JP2006284180A (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.)
Denso Corp
Denso IT Laboratory Inc
Original Assignee
Denso Corp
Denso IT Laboratory Inc
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 Denso Corp, Denso IT Laboratory Inc filed Critical Denso Corp
Priority to JP2005100314A priority Critical patent/JP4808984B2/ja
Priority to DE102006013448.6A priority patent/DE102006013448B4/de
Priority to US11/390,866 priority patent/US7372404B2/en
Publication of JP2006284180A publication Critical patent/JP2006284180A/ja
Application granted granted Critical
Publication of JP4808984B2 publication Critical patent/JP4808984B2/ja
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S3/00Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
    • G01S3/02Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using radio waves
    • G01S3/74Multi-channel systems specially adapted for direction-finding, i.e. having a single antenna system capable of giving simultaneous indications of the directions of different signals
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S3/00Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received
    • G01S3/80Direction-finders for determining the direction from which infrasonic, sonic, ultrasonic, or electromagnetic waves, or particle emission, not having a directional significance, are being received using ultrasonic, sonic or infrasonic waves
    • G01S3/8006Multi-channel systems specially adapted for direction-finding, i.e. having a single aerial system capable of giving simultaneous indications of the directions of different signals

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Radar Systems Or Details Thereof (AREA)
  • Radio Transmission System (AREA)

Description

本発明は、複数のアンテナ素子に到来する、電波や音波などの到来波の方向推定装置に関する。
従来、複数のアンテナ素子で到来する同一周波数の一つ以上の電波や音波などの波を捕捉し、その到来方向を推定演算する技術が知られている。こうした方法として、一般にMUSIC(Multiple Signal Classification)法が知られている(例えば、特許文献1参照)。
MUSIC法は、到来波の信号データベクトルから、それら信号データベクトルの相関行列を求め、当該求められた相関行列を、固有ベクトルと固有値とに分解する固有値分解演算を行い、求められた固有ベクトルから、到来波の方位を推定するものである。上記した到来方位を求める際には、事前に到来波の数を推定しておく必要がある。到来波数の推定は、一般に固有値を用いて行われ、AIC、MDL、しきい値法などがある。
特開平2000−121716号公報
この、固有値分解演算には、ハウスホルダや、ヤコビ法を利用した方法等が知られているが、到来波のSN比が高い、即ち、到来波がノイズに比して非常に大きな場合や、複数の到来波の到来方位差が小さな場合に、固有値の最大値と最小値の比が大きくなるために、特に小さい方の固有値の演算中に桁落ちが生じ、演算された固有値が負の値を取ってしまったり、信頼度の低い値をとってしまったりする場合がある。特に、演算に不十分な語長の固定小数点方式を用いていると、上述した不都合が現れやすい。
固有値の値が信頼性の低いものとなってしまうと、当該固有値から演算される到来波の数は不正確なものとなり、正しい到来波の数の推定が困難となる不都合が生じるばかりか、結果的に信頼性の低い固有値を求めるために、固有値分解演算において無駄な演算動作が繰り返され、演算処理装置に無用な演算負荷を与えることとなる。
本発明は、上記した事情に鑑み、無駄な固有値分解演算を極力なくすことが出来、信頼性の高い到来波の数の判定が可能な、到来波の方向推定装置の提供を可能とするものである。
請求項1の発明は、複数のアンテナで受信した複数の到来波の受信データから、前記到来波の到来方向を推定する、到来波の方向推定装置において、
複数のアンテナ(15)で受信した、複数の到来波の受信データから、その相関行列(R)を演算する相関行列演算手段(6)、
前記相関行列演算手段で演算された前記相関行列に対して、ハウスホルダ変換演算(図3のステップS5からステップS8など)を実施して、前記相関行列を、固有ベクトル(Q)と固有値(λ)とに分解処理する固有値分解演算手段(7)、
前記固有値分解演算手段に設けられ、前記相関行列について、固有値を求めるための行列(QT)と固有ベクトルを求めるための行列(A)を生成する行列生成手段(7,図3のステップS1、ステップS2、ステップS10など)、
前記固有値分解演算手段に設けられ、前記固有値を求める行列において、前記ハウスホルダ変換演算を行うたび毎に、前記ハウスホルダ変換演算の対象となる行列に対して、各列の成分の列ノルムを演算する列ノルム演算手段(7,ステップS3など)、
前記固有値分解演算手段に設けられ、前記ハウスホルダ変換演算を行うたび毎に、前記列ノルム演算手段により演算されたハウスホルダ変換演算の対象となる行列の列ノルムの最大値が、所定の第1の閾値(Th)以下であるか否かを判定する、列ノルム最大値判定手段(7,図3のステップS4など)、
前記固有値分解演算手段に設けられ、前記列ノルム最大値判定手段により、前記列ノルム演算手段により演算されたハウスホルダ変換演算の対象となる行列の列ノルムの最大値が、前記所定の第1の閾値以下でないものと判定された場合に限り、前記固有値を求めるための行列に対してハウスホルダ変換演算を、繰り返し実行するハウスホルダ変換演算実行手段(7,図3のステップS3乃至ステップS8)、
前記固有値分解演算手段に設けられ、前記列ノルム最大値判定手段により、前記列ノルム演算手段により演算されたハウスホルダ変換演算の対象となる行列の列ノルムの最大値が、前記所定の第1の閾値以下であるものと判定された場合に、前記ハウスホルダ変換演算実行手段によるハウスホルダ変換演算の実行を中断すると共に、前記ハウスホルダ変換演算実行手段によるハウスホルダ変換演算の実行回数を所定のメモリに格納する変換演算実行制御手段(7,図3のステップS4など)、
前記固有値分解演算手段に設けられ、前記変換演算実行制御手段によりハウスホルダ変換演算の実行が中断された前記固有値を求めるための行列に対して、前記ハウスホルダ変換演算実行手段及び変換演算実行制御手段によるハウスホルダ変換演算を、所定回数に達するまで更に実行させて、その結果から、前記相関行列の前記固有ベクトルと固有値を求める、収束演算実行手段(7.図3のステップS11、S12,S14など)、及び、
前記求められた固有値を構成する値と所定の第2の閾値を比較し、当該第2の閾値よりも大きな値については到来波の信号の固有値とし、小さなものは雑音の固有値と判定して、到来波の数を推定演算すると共に、該推定演算された到来波の数と前記メモリ手段に格納されたハウスホルダ変換演算の実行回数を比較し、その少ない方の値を到来波の数と判定する、信号ランク判定手段(9,10)、
から構成される。
請求項2の発明は、前記第1の閾値(Th)として、0又は0に近い値格納するメモリ手段を有し、
前記変換演算実行制御手段は、前記列ノルムの最大値が負の場合に、前記ハウスホルダ変換演算の更なる実行を中断させる、ことを特徴として構成される。
請求項3の発明は、複数のアンテナで受信した複数の到来波の受信データから、前記到来波の到来方向を推定する、到来波の方向推定装置において、
複数のアンテナ(15)で受信した、複数の到来波の受信データから、その相関行列(R)を演算する相関行列演算手段(6)、
前記相関行列演算手段で演算された前記相関行列に対して、ハウスホルダ変換演算(図3のステップS5からステップS8など)を実施して、前記相関行列を、固有ベクトル(Q)と固有値(λ)とに分解処理する固有値分解演算手段(7)、
前記固有値分解演算手段に設けられ、前記相関行列について、固有値を求めるための行列(QT)と固有ベクトルを求めるための行列(A)を生成する行列生成手段(7,図3のステップS1、ステップS2、ステップS10など)、
前記固有値分解演算手段に設けられ、前記固有値を求める行列において、前記ハウスホルダ変換演算を行うたび毎に、前記ハウスホルダ変換演算の対象となる行列に対して、各列の成分の列ノルムを演算する列ノルム演算手段(7,ステップS3など)、
所定の第1の閾値を、前記受信データにおいて雑音のみが捕捉された際の、前記固有値を求めるための行列の列ノルムよりも大きく、前記受信データの到来波の信号についての、前記固有値を求めるための行列の列ノルムよりも小さな値となるように設定する閾値設定手段、
前記固有値分解演算手段に設けられ、前記ハウスホルダ変換演算を行うたび毎に、前記列ノルム演算手段により演算されたハウスホルダ変換演算の対象となる行列の列ノルムの最大値が、前記所定の第1の閾値(Th)以下であるか否かを判定する、列ノルム最大値判定手段(7,図3のステップS4など)、
前記固有値分解演算手段に設けられ、前記列ノルム最大値判定手段により、前記列ノルム演算手段により演算されたハウスホルダ変換演算の対象となる行列の列ノルムの最大値が、前記所定の第1の閾値以下でないものと判定された場合に限り、前記固有値を求めるための行列に対してハウスホルダ変換演算を、繰り返し実行するハウスホルダ変換演算実行手段(7,図3のステップS3乃至ステップS8)、
前記固有値分解演算手段に設けられ、前記列ノルム最大値判定手段により、前記列ノルム演算手段により演算されたハウスホルダ変換演算の対象となる行列の列ノルムの最大値が、前記所定の第1の閾値以下であるものと判定された場合に、前記ハウスホルダ変換演算実行手段によるハウスホルダ変換演算の実行を中断すると共に、前記ハウスホルダ変換演算実行手段によるハウスホルダ変換演算の実行回数を所定のメモリに格納する変換演算実行制御手段(7,図3のステップS4など)、
前記固有値分解演算手段に設けられ、前記変換演算実行制御手段によりハウスホルダ変換演算の実行が中断された前記固有値を求めるための行列に対して、前記ハウスホルダ変換演算実行手段及び変換演算実行制御手段によるハウスホルダ変換演算を、所定回数に達するまで更に実行させて、その結果から、前記相関行列の前記固有ベクトルと固有値を求める、収束演算実行手段(7.図3のステップS11、S12,S14など)、及び、
前記メモリ手段に格納されたハウスホルダ変換演算の実行回数を到来波の数と判定する、信号ランク判定手段、
を有することを特徴として構成される。
請求項の発明は、請求項3の発明において、前記閾値設定手段は、前記所定の第1の閾値を、周囲の環境や観測距離に対応させて変動的に設定することを特徴として構成される。
請求項の発明は、請求項1又は3の発明において、前記ハウスホルダ変換演算実行手段によるハウスホルダ変換演算の実行回数の最大数を設定する演算回数設定手段を有し、
前記ハウスホルダ変換演算実行手段によるハウスホルダ変換演算を、前記演算回数設定手段に設定されたハウスホルダ変換演算の実行回数の最大数以内で行うことを特徴として構成される。
請求項の発明は、請求項の発明において、前記演算回数設定手段は、前記ハウスホルダ変換演算の実行回数の最大数を、周囲の環境や観測距離に対応させて変動的に設定することを特徴として構成される。
請求項の発明は、請求項1又は3の発明において、前記収束演算実行手段による、前記ハウスホルダ変換演算実行手段及び変換演算実行制御手段による前記ハウスホルダ変換演算の、各回の実行回数を制御する変換演算回数制御手段を設け、
前記変換演算回数制御手段は、N回目の前記ハウスホルダ変換演算実行手段及び変換演算実行制御手段によるハウスホルダ変換演算の実行回数を、それ以前の回におけるハウスホルダ変換演算の実行回数以下に制御することを特徴として構成される。
請求項1の発明によれば、ハウスホルダ変換演算は、固有値を求める行列から演算される列ノルムの最大値が、所定の第1の閾値以下の場合に、それ以上のハウスホルダ変換演算を行わないので、第1の閾値を適宜設定することにより、信号に対応せず、かつ演算結果の信頼性の低い列についてハウスホルダ変換演算が実行されることを未然に防止することが出来る。これにより、無駄な固有値分解演算の発生を極力防止することが出来る。
また、固有値を求める行列から演算される列ノルムの最大値が、所定の第1の閾値以下の場合には、それ以降のハウスホルダ変換演算を行わないので、その時点までのハウスホルダ変換演算の実行回数が、信号数(ランク)と対応することとなり、より信頼性の高い到来波のランクの判定が可能となる。これにより、たとえ、固有値から求められた到来波の数が、誤って多く演算されたとしても、ハウスホルダ変換演算の実行回数に基づいて修正される形でその到来波の数が判定されるので、より正確な判定が可能となる。
請求項2の発明によると、第1の閾値は、0又は0に近い値であることから、固有値が負となる列に対する無駄なハウスホルダ変換演算の実行を排除することが可能となる。
請求項3の発明によると、第1の閾値は、前記受信データにおいて雑音のみが捕捉された際の、前記固有値を求めるための行列の列ノルムよりも大きくなるように予め設定されているので、雑音に対するハウスホルダ変換演算の実行を排除し、演算負荷を軽減することができる。
更に、第1の閾値は、前記受信データにおいて雑音のみが捕捉された際の、前記固有値を求めるための行列の列ノルムよりも大きく、前記受信データの到来波の信号についての、前記固有値を求めるための行列の列ノルムよりも小さな値となるように予め設定されているので、雑音に対するハウスホルダ変換演算の実行を排除し、演算負荷を軽減することができばかりか、ハウスホルダ変換演算回数を直ちに到来波の信号の数とすることが出来、固有値に基づく到来波の数の判定処理をなくすことが出来る。
請求項の発明によると、例えば、カーナビゲーションシステムなどを用いて周囲の環境を検出し、その結果、検出されたトンネル内などにおいては、レーダから発射した電波が、多くの場所(車両、壁など)に反射されて帰ってくることが予想される。こうした到来波の数が多くなることが予想され場合に、第1の閾値(Th)を低く設定して、より多くの到来波を的確に捕捉し、それ以外の到来波の数が少ないことが予想される開放された路上などにおいては、第1の閾値(Th)を高めに設定して、無駄な雑音列ノルムに係わる信号の演算動作を行わないようにし、効率的な演算動作を可能とすることが出来る。また、観測距離が遠くなる場合も、多くの場所から反射してくることが予想されるため、上述と同様のことが可能である。
請求項の発明によると、到来波の最大数が予め判っているような場合には、ハウスホルダ変換演算実行手段によるハウスホルダ変換演算の実行回数を予め当該最大数に設定しておき、無駄な演算動作が発生することを防止するように構成することも可能である。
請求項の発明によると、例えば、カーナビゲーションシステムなどを用いて周囲の環境を検出し、その結果、検出されたトンネル内などにおいては、レーダから発射した電波が、多くの場所(車両、壁など)に反射されて帰ってくることが予想される。こうした到来波の数が多くなることが予想され場合に、ハウスホルダ変換演算の実行回数を多く設定して、より多くの到来波についての固有値を的確に演算し、それ以外の到来波の数が少ないことが予想される開放された路上などにおいては、ハウスホルダ変換演算の実行回数を少なく設定して、無駄な雑音列ノルムに係わる信号の固有値演算動作を行わないようにし、効率的な演算動作を可能とすることが出来る。また、観測距離が遠くなる場合も、多くの場所から反射してくることが予想されるため、上述と同様のことが可能である。
請求項の発明によると、前の回の演算サイクルにおいて桁落ちなどが生じてハウスホルダ変換演算が行われていない行列部分に対してハウスホルダ変換演算が行なわれることが防止され、演算効率を高めることが出来る。
なお、括弧内の番号等は、図面における対応する要素を示す便宜的なものであり、従って、本記述は図面上の記載に限定拘束されるものではない。
以下、図面に基づき、本発明の実施例を説明する。
図1は、到来波方位測定装置の一例を示すブロック図、図2は、信号処理の一例を示すチャート、図3は、固有値分解演算処理の一例を示すフローチャート、図4乃至図10は、図3のフローチャートにおける、具体的な演算例の一例を示す図、図11は、列ノルムの演算方法の一例を示す図、図12は、到来波数の判定処理の別の例を示すチャート、図13は、演算式の一例を示す図、図14は、固有値分解演算の基本概念の一例を示すフローチャートである。
到来波方位測定装置1は、図1に示すように、主制御部2を有しており、主制御部2には、バス線3を介して、信号受信部5、相関行列演算部6、固有値分解演算部7、信号ランク推定演算部9、信号ランク判定部10、スペクトル演算部11、到来波方位演算部12及び出力部13などが接続している。信号受信部5には、K個(複数)のアンテナ15が設けられている。
なお、図1に示した到来波方位測定装置1は一例であり、図面にブロックで示された各構成部分は、コンピュータを構成するCPU(中央演算処理装置)及び適宜なメモリ、さらに、当該コンピュータを活性化させるプログラムにより、置換することも当然可能である。この場合、コンピュータがプログラムの各ステップを実行することにより、図1に示した到来波方位測定装置1を構成する各ブロックとして機能することとなる。
以下、図1に示す、到来波方位測定装置1の信号受信部5,相関行列演算部6,固有値分解演算部7,信号ランク推定演算部9,信号ランク判定部10の機能について、図14と共に説明する。
[信号受信部5]
サンプリング時刻mΔT1(ΔT1はサンプリング間隔、mは自然数)において、各アンテナ素子15-1から15-K(Kはアンテナ素子数)で得られる複素データをx1(m),x2(m),…,xK(m)とすると、信号受信部では、式(1)で表されるサイズKのデータベクトルx(m)を得る。上付きのTは転置を表す。
Figure 0004808984
{相関行列演算部6}
相関行列演算部では式(2)を用いて、時刻nΔT2(ΔT2は相関行列を算出する時間間隔)における相関行列Rxx(n)を求める。Rxxは(K×K)の半正定値あるいは正定値エルミート行列である。あとでも触れるが、半正定値エルミート行列の場合、固有値はすべて0以上の実数となる(正定値の場合、固有値は0より大きい実数となる)。上付きのHは複素共役転置を表す。SSNはスナップショット数で、相関行列演算に用いるデータベクトルの数を表す。βは忘却係数で0〜1の範囲の値を取る。Rxx(n-1)は時刻(n-1)ΔT2における相関行列である。
Figure 0004808984
なお、アンテナ素子が中心対称に配置されている場合、式(3)に示すような変換をすることで実数対称な相関行列Ryyを得られる。QKはK×Kのユニタリ行列である。詳細は文献[2][4]を参照。realは実数部のみを取り出す演算を表す。
Figure 0004808984
実数相関行列Ryyを用いれば、後の固有値分解演算量を大きく削減することができる。
こうして求められた相関行列Ryy(あるいはRxx)は、固有値分解演算部に出力され、固有ベクトルQと固有値λに分解される。
{固有値分解演算部7}
固有値分解演算は、ハウスホルダ変換演算を利用したQR分解を繰りかえし利用することで行う。1回の固有値分解演算は複数回のQR分解演算を含み、1回のQR分解演算は複数回のハウスホルダ変演算換を含む。ここでは、1回の固有値分解演算中に行うQR分解演算の回数をNQR[回]とする。また、1回のQR分解演算においては、最大(K-1)回(Kは素子数)のハウスホルダ変換演算が実行される。このハウスホルダ変換演算の実行回数は適宜なメモリ手段に格納される。
フローチャート例(図14)を用いて処理の流れを説明する。最初は、列ピボット選択が行われる(図14-S100)。列ピボット選択演算は、演算対象行列(上記の場合RxxやRyy)のランクが不足している場合にでも、得られる固有ベクトルが直交し、演算精度が向上する利点がある。具体的には、演算対象となる行列の中で「最大のノルムを持つ列」を一番左の列に移動(交換)する操作を行う。列ピボット選択を用いたハウスホルダ変換演算については文献[3]p248-p250に記述がある。
列ピボット選択後、ハウスホルダ変換処理を実行するか否かを判定する(図14-S200)。この判定は、列ノルムの最大値が、予め設定されたしきい値(Th)より大きいか否かによって決定する(以下)。
列ノルムの最大値>=しきい値 ハウスホルダ変換を実行
列ノルムの最大値< しきい値 ハウスホルダ変換を中止
通常、しきい値Thは列ノルムの最大値が負であることをチェックすることのできる値、即ち、「0」または「0」に近い値が予め設定されている。先にも述べたが,Ryy(あるいはRxx)が半正定値のエルミート行列であるので、得られる固有値は必ず0以上の実数となる(正定値の場合0より大きな実数)。固有値と列ノルムは相互に強く関連している。列ノルムが負または0に近い値となった場合には、桁落ちなどの演算誤差が生じていることが考えられ、以後、更なるハウスホルダ変換演算を行っても求められる固有値の信頼性が低下してしまい、演算負荷ばかりが増加してしまうことから、ハウスホルダ変換を中止しステップ図14-S500へ入る。
一方ハウスホルダ変換を実行と判定された場合、ハウスホルダ変換を実行する(図14-S300)。
例としてK=4の場合のj(=1,2,3)回目のハウスホルダ変換終了後に得られる行列A(j)を示す。ここで×印は不特定の値が入る。Hは(K×K)のハウスホルダ行列である。A(0)は初期値であり、最初は相関行列Ryyが入る。
Figure 0004808984
Figure 0004808984
Figure 0004808984
Figure 0004808984
このように、j回目のハウスホルダ変換により、行列A(j-1)の第j列の第j+1行目以降が0となる。このような変換を行うハウスホルダ行列Hの算出法については文献[2]p209-p210参照。ハウスホルダ変換演算を最大(K-1)実行すれば、行列A(K-1)は右上三角行列(対角項より下側成分がすべて0の行列)となる。実際の処理では、2回目以降のハウスホルダ変換はしきい値判定により中断することがあるため、完全な上三角行列にはならない場合がある。この場合、ハウスホルダ変換演算が実行されなかった部分の固有値は正しい値に収束しない。ただし、最終的に得られる固有ベクトル行列Qは各列ベクトルが正規直交するユニタリ行列となるため、最低でも信号ランク数(到来波の数)だけハウスホルダ変換が実行されれば(=信号の固有値・固有ベクトルが正しく得られていれば)、雑音固有ベクトルの張る空間(雑音部分空間)は正しく得られる。これは式(8)の関係から導かれる。ここでEsEsHは信号部分空間を表すK×K行列、EnEnHは雑音部分空間を表すK×K行列である。Es=[e(1) e(2) … e(L)]は信号固有ベクトルをL個横並べたK×L行列、En=[e(K-L) … e(K)]は雑音固有ベクトルを(K-L)個横に並べたK×(K-L)行列、IkはK×Kの単位行列である。このように、信号固有ベクトルEsが正しくもとまれば正しい雑音固有空間が一意に決定される。
Figure 0004808984
つまりは、MUSICスペクトル算出する場合に必要なのは雑音固有空間であるので、正確なスペクトルを得ることができる。
図14-S200において、ハウスホルダ変換が中断するか、図14-S400において、ハウスホルダ変換が規定回数(K-1)回に達したと判定された場合、図14-S500に入り式(9)(10)を計算する。
Figure 0004808984
Figure 0004808984
式(10)において、Qの初期値はK×Kの単位行列を入れておく。ここでNh(ite)(=<K-1)はite回目のQR分解演算において実行されたハウスホルダ変換実行回数であり、後の信号ランクの推定に利用するため適宜メモリ手段に記憶しておく。次に図14-S600においてQR分解演算が規定回数(NQR[回])に達したかの判定を行う。QR分解演算が規定回数に達していない場合、A(0) は次回QR分解演算への入力として使用される。QR分解演算が規定回数に達した場合、固有値分解演算を終了する。このときA(0)の対角項が固有値λ=[λ1 λ2 … λK]Tであり、行列QH=[e(1) e(2) … e(K)]Hの各列が各固有値λ1〜λKに対応した固有ベクトルである。
{信号ランク推定演算部9}
こうして固有値分解演算部7において、相関行列Ryy(あるいはRxx)の固有値λが求められたところで、当該求められた固有値λを、図1の信号ランク推定演算部9に出力し、信号ランク推定演算部9は、固有値λから、到来波数(信号ランク)を推定する。推定法としては、AICやMDL、しきい値法などがある。しきい値法では、例えば予め設定さえた所定しきい値より大きな固有値については到来波の信号の固有値とし、小さなものは雑音の固有値と判定して、信号ランクLt、即ち到来波の数Ltを推定する。
{信号ランク判定部10}
信号ランク判定部10は、信号ランク推定演算部9で推定された信号ランク、即ち到電波の数Ltと、固有値分解演算部7で実行されたハウスホルダ変換実行回数Nh(i)(i=1,2,..,NQR)のうち、一番小さい値を信号ランク、即ち到来波の数Lとする演算を行い、求められた到来波(信号)の数Lを、スペクトル演算部11に対して出力する。
{スペクトル演算部11}
こうして、信号ランク推定部10により、ハウスホルダ変換演算回数Nhが考慮された形で、信号ランクLが判定され、その値Lが、図1に示す、スペクトル演算部11に出力されると、スペクトル演算部11は、公知の手法により、図2に示すようにMUSICスペクトラムSPEを演算して求め、当該求められたMUSICスペクトラムSPEを、到来波方位演算部12に出力する。
{到来波方位演算部12}
到来波方位演算部12は、当該MUSICスペクトラムSPEから、公知の手法で、到来波の到来する方位θ1、θ2を演算することが出来るが、MUSICスペクトラムSPEは、既に述べたように、信号の数Lが、相関行列Ryy(あるいはRxx)の固有値分解演算に際したハウスホルダ変換演算回数Nhに基づいて補正され、信頼性の低いデータが排除されているので、実態のない所謂「偽ピーク」を排除した形で、正確な到来波の方位を算出することができる。ただし「偽ピーク」を完全に排除できるわけではないので、より正確な推定には電力推定演算などを利用する必要がある(関連参考特許文献[5])。
文献[2]菊間信良、‘アダプティブアンテナ技術’、オーム社、ISBN4-274-03611-1
文献[3]Gene H Golub et al,‘MATRIX COMPUTATION’(3rd edition),Jhon Hopkins Univ.Press,ISBN0-8018-5414-8
参考特許文献[4] 特願平10−152629
参考特許文献[5] 特願平10−290766
以下、より具体的に説明すると、到来波方位測定装置1は、発信位置(レーダの場合は反射位置)の異なる1個以上の、同一周波数の複数の到来波16が方位θ1、θ2、……などでK個のアンテナ15にそれぞれ入射すると、信号受信部5は、それら各アンテナ15で受信した、一個以上の到来波16を含む受信データxを、受信データベクトルXとして、相関行列演算部6に出力する。受信データベクトルXは、例えば、図13の(a)で示すように表示される。図中xi(iは個々のアンテナを示す序数で、この場合、アンテナは例えば、3個有り、順にi=1,2,3の番号が付されているものとする。なお、アンテナの数は任意である)は、各アンテナ15で受信した信号の受信データである。以後、MUSIC法に基づいて、受信データに対する処理を実行する。
K個のアンテナ素で受信された受信データxi(i=1,2……K)からなる受信データベクトルXは、相関行列演算部6に出力され、相関行列演算部6は、当該受信データベクトルXに基づいて、図13(b)に示すように、K次の正方行列からなる受信データベクトルXの相関行列Rを、演算生成する。こうして求められた相関行列Rは、固有値分解演算部7に出力され、ここで、固有ベクトルQと固有値λとに分解演算される。
この固有値分解演算部7による固有値分解演算は、図3に示すような手順で行われる。これらの手順は、コンピュータによりプログラマブルに実行することも出来る。即ち、固有値分解演算部7は、ステップS1に示すように、相関行列演算部6から出力された相関行列Rを、固有値を求めるための行列Aに代入し、到来波方位測定装置1に設けられたメモリであるRAM(ランダムアクセスメモリ)17に、図3及び図4(1)及び(2)に示すように、格納する、初期化処理を実行する。なお、図4に示す演算例は、説明を簡略化するために、アンテナの数が3個の場合で説明するが、既に述べたように、アンテナの数は任意である。
また、初期化に際しては、図4(2)に示すように、固有ベクトルを求めるための行列QTとして、素子数K(図4(2)の場合、素子数は3)の単位行列を代入し、更に、固有値を求める収束演算のサイクル演算回数を表示するパラメータとしてIte=1を代入し、図3(2)に示すように、RAM17に格納する。
次に、図3のステップS2に入り、図4(3)に示すように、ハウスホルダ変換演算回数Nh=1としてRAM17に格納すると共に、行列Qiとして、素子数Kに対応したK次の単位行列(図4(2)の場合、素子数は3なので、3次の単位行列)を代入し、同様にRAM17に格納する。
次に、図3のステップS3に入り、以下ステップS3からステップS8から構成される、1回以上のハウスホルダ変換演算を含むQR分解演算を、図3のステップS11で示す所定のサイクル回数NQRだけ実行する、固有値分解演算を行う。
即ち、まず行列Aの列ノルムccを演算する。図4(2)に示す行列Aの列ノルムccは、図4(4)に示すように求められる。いまだハウスホルダ変換演算を実行していない状態の行列Aの列ノルムccの演算は、図11(1)に示すように、行列Aの全ての成分について行われる。
ここで、列ノルムccの最大値を示す列を抽出し、ind_kに代入して、その値をRAM17に格納する。図4(4)の場合、2列目の列ノルムccが最大値:2.9260を取るので、ind_k=2を代入する。
次に、図3のステップS4に入り、列ノルムccの最大値が、予め設定され、RAM17などのメモリに格納された閾値Thより大きいか否かを判定すると共に、前述のハウスホルダ変換演算回数Nhが、アンテナの素子数K未満であるか否かを判定する。通常、閾値Thは列ノルムccの最大値が負であることをチェックすることの出来る値、即ち、「0」又は「0」に近い値(後述する、列ノルムccの最大値が負となる場合をチェックすることが出来れば、必ずしも0でなくても良いという意味である)が、固有値分解演算部7により、予め設定されている(プログラムで設定することも当然可能である)ので、列ノルムccが負の値となった状態で、後述するハウスホルダ演算が行われる事態を回避することが出来る。通常、列ノルムccは負の値を取ることはないので、列ノルムccが負又は「0」に近い値となった場合には、桁落ちなどの演算誤差が生じていることが考えられ、以後、更なるハウスホルダ演算を行っても求められる固有値の信頼性が低下してしまい、演算負荷ばかりが増加してしまうことから、後述するステップS5からステップS8のハウスホルダ演算は行なわず、ステップS9(詳細は後述する)に入る。
また、信号受信部5で捕捉し、解析可能な到来波の数は、通常、(素子数K−1)であることから、固有値を求めるのに必要な、K次の相関行列Rに対するハウスホルダ変換演算回数Nhも、最大(素子数K−1)となる。従って、ステップS4では、ハウスホルダ変換演算回数Nhが、素子数K未満となるように制御し、無駄な変換演算の発生を防止する。
ステップS4で、列ノルムccの最大値が閾値Th以上であり、かつハウスホルダ変換演算回数Nhが、素子数K未満と判定された場合には、図3のステップS5に入り、ステップS8までの間で、行列Aに対するハウスホルダ変換演算を実行する(図5の(1))。
まず、ステップS5では、ステップS3の列ノルムccが最大値を示した列と、ハウスホルダ変換演算により対角要素以外の成分が0に未だ変換されていない列の内、最左方の列を入れ替える、ピポット選択を行う。即ち、図4(2)の行列Aの場合、列ノルムccの最大値を示す列は、図4(4)のind_k=2より、2列目であり、行列Aに対するハウスホルダ変換演算は未だ実行されていないので、行列Aの最左方である1列目の列と、2列目の列を入れ替えるピポット選択を演算実行し、図5(2)の行列Aを求め、RAM17に更新格納する。
更に、ステップS6で、ハウスホルダ演算回数Nhに1を加算して、当該加算された値をNh=2を、新たなハウスホルダ演算回数Nhとして、RAM17に更新格納する。
次に、ステップS7に入り、図5(3)に示すように、行列Aに対するハウスホルダ行列Hを演算して求め、更にステップS8で、A=H*A、Qi=H*Qiなる行列演算を行い、求められた、図5(4)に示す、行列A及びQiを、RAM17にそれまで格納されていた行列A及びQiを更新する形で格納する。この場合、行列Aの第1列目の対角要素より下の二つの要素が0に変換される。
次に、図3のステップS3に戻り、行列Aの、次にハウスホルダ変換演算の対象となる行列に対して、再度列ノルムccを演算し、図5(5)に示す値を得る。即ち、一度ハウスホルダ変換演算が行われた行列Aについての列ノルムccの演算は、図11(2)に示すように、既に「0」にされた1列目を除外した、次にハウスホルダ変換演算を行って、対角成分より下の部分を「0」とする対象となる行列部分、従って、2列目及び3列目の対角成分a′22、a′23、a′32、a′33を含む、図中四角で囲んだ正方行列に対して求める。その演算方法は、(A)又は(B)の二通りがある。得られた列ノルムccから、既に1回目のハウスホルダ変換演算に際して最大値が採用された一列目の列ノルムccを除いた2列目以降の列ノルムccを比較し、最大値を取る列を求める。図5(5)の場合、第3列目の列ノルムccが最大値となるので、ind_k=3とする演算を行う。
同様に、ステップS4で、列ノルムccの最大値が、予め設定された閾値Thより大きいか否かを判定すると共に、前述のハウスホルダ変換演算回数Nhが、アンテナの素子数K未満であるか否かを判定する。このばあい、列ノルムcc=0.6555で、「0」より大きく、また、ハウスホルダ演算回数Nhも2であり、素子数K=3よりも少ないので、再度ステップS5に入り、ハウスホルダ変換演算を再度実行する。
ステップS6では、図6(2)に示すように、列ノルムccの最大値を取る列と、ハウスホルダ変換演算により既に「0」変換された列の右方に隣接する列を入れ替えるピポット選択を実行する。この場合、第2列と第3列を入れ替え、図6(2)に示す行列Aを得る。また、ステップS6で、ハウスホルダ変換演算回数Nhを、1加算し、Nh=3とする。
この状態で、ステップS7を実行して、図6(3)に示すように、図6(2)の行列Aに対するハウスホルダ行列Hを演算して求め、更にステップS8で、A=H*A、Qi=H*Qiなる行列演算を行い、求められた、図6(4)に示す、行列A及びQiを、RAM17にそれまで格納されていた行列A及びQiを更新する形で格納する。この演算により、行列Aの第2列目の対角要素より下の1つの要素が0に変換される。これにより、行列Aは、第1列目及び第2列目の対角要素より下の要素がすべて「0」に変換されたこととなる。
次に、図3のステップS3に戻り、再度行列Aの列ノルムを求め、ステップS4に入るが、この際ハウスホルダ変換演算回数Nhが既に3となっており、素子数Kと等しくなっているので、これ以上のハウスホルダ変換演算は不要なものと判定して、固有値を求める収束演算の最初の演算サイクルを抜けてステップS9に入る。
ステップS9では、ステップS5でピポット選択により入れ替えた列順を戻す。即ち、図6(4)の、ハウスホルダ変換演算を2回行った行列Aは、相関行列Rの1列、2列、3列が、2列,3列、1列と入れ替えられた形となっており、その入れ替え状態を示すために、図6(5)に示すように、列入れ替え情報pivがpiv=231と設定され、RAM17に格納される。
列入れ替え情報pivに基づいて、図7(1)に示すように、図6(4)の行列Aの各列を、相関行列Rの列に対応する形で入れ替える演算をおこなう。即ち、図6(4)の行列Aの1列を2列に、2列を3列に、3列を1列に入れ替える演算処理を行い、その結果を行列Aiに代入する。
次に、ステップS10に入り、図7(2)に示すように、行列A=Ai*Qi、QT=Qi*QTの演算を行い、求められた行列A及びQTをRAM17に格納する。次に、ステップS11に入り、RAM17内の、固有値を求めるための収束演算のサイクル回数を表示するパラメータIte=1が、所定のサイクル演算回数NQRよりも大きいか否かを判定し、演算回数Ite=1が所定のサイクル演算回数NQRよりも小さな場合には、ステップS12に入り、サイクル回数を表示するパラメータIteに1を加算してIte=2としてRAM17の内容を更新する。以後、固有値分解演算の収束演算のサイクルは2周目入り、図7(2)で求めた行列Aに対して、再度ハウスホルダ変換演算を実行して、行列Aの対角要素を固有値に収束させてゆく演算処理を行う。
即ち、収束演算の演算サイクルの2周目に入った固有値分解演算は、図3のステップS2で、ハウスホルダ変換演算回数NhをNh=1に戻してRAM17の値を更新すると共に、行列Qiとして、素子数Kに対応したKの単位行列を再代入し、RAM17の値を更新する。
以後、図7(5)から図9(3)により、演算サイクルの2周目の固有値分解演算において、図7(2)の行列Aに対して、2度のハウスホルダ変換演算を更に、実行し、図3のステップS4で、Nh(=3)がアンテナの数K(=3)となって、ステップS9に入る。
ステップS9では、図8(1)及び図9(1)で示すように、ステップS5でピポット選択により入れ替えた列順を戻して(図8(1)及び図9(1)の場合、列ノルムccの最大値は、それぞれ1列目(図8(1)のind_k=1参照)と2列目(図8(4)のind_k=2参照)なので、行列Aのピポット選択は、実際には行われていない)、行列Aiとしてその値を、RAM17に格納する。
なお、固有値分解演算部7は、演算サイクルの2周目以降の固有値分解演算においては、1周目で行った行列A(相関行列R)に対して行ったハウスホルダ変換演算の回数を超えて、1周目で求められた行列Aに対するハウスホルダ変換演算を行うことはしないように制御する。これは、1周目の固有値分解演算において、ステップS4で、列ノルムccの最大値が所定の閾値Th以下と判定され、それ以上ハウスホルダ変換演算を行っても、演算結果の信頼性が低いものと判断された場合、2週目以降の収束演算において、1周目のハウスホルダ変換演算回数を超えた回数の変換演算を行っても、演算結果の信頼性は同様に低くなる。従って、2週目以降の固有値分解演算におけるハウスホルダ変換演算回数は、1周目で行った行列A(相関行列R)に対して行ったハウスホルダ変換演算の回数を超えないように制御し、その分の演算処理に要する負担をなくすようにする。即ち、各演算サイクルで実行するハウスホルダ変換演算の回数は、最初の演算サイクルで実行されたハウスホルダ変換演算の回数が最大値となる。また、ハウスホルダ変換演算回数Nhは、MUSIC法の原理に基づき、最大でも(アンテナの数K−1)回である。
次に、ステップS10に入り、図10(1)に示すように、行列A=Ai*Qi、QT=Qi*QTの演算を行い、求められた行列A及びQTをRAM17に格納する。次に、ステップS11に入り、RAM17内の演算回数を表示するパラメータIte=2が、所定の演算回数NQRよりも大きいか否かを判定し、演算回数Ite=2が所定の演算回数NQRよりも等しいか大きな場合には、行列Aの対角要素は固有値に収束したものと判定して、ステップS14に入る。
ステップS14では、図10(1)で求められた行列Aから、対角成分を取り出して(diag(A))、図10(3)に示すように、固有値λを求め、更に固有ベクトルQを、図10(1)の行列QTの共役転置行列として求め、RAM17などの適宜なメモリに格納する。
こうして、固有値分解演算部7において、相関行列Rの固有値λが求められたところで、当該求められた固有値λを、図1の信号ランク推定演算部9に出力し、信号ランク推定演算部9は、固有値λから、当該固有値を構成する値と所定の閾値を比較して、当該閾値よりも大きな値については到来波の信号の固有値とし、小さなものは雑音の固有値と判定して、信号ランクLt、即ち到来波の数Ltを推定演算する。例えば、図10(3)の固有値λの場合、閾値を1.0とすると、当該閾値をこえる固有値の値は、2.2192だけなので、信号ランクLtを「1」と推定して、信号ランク判定部10に当該推定された到来波の(ランク)数Ltを出力する。
信号ランク判定部10は、信号ランク推定演算部9で推定された信号ランク、即ち到来波の数Ltと、固有値分解演算部7で実行されたハウスホルダ変換演算回数Nhを比較し、その少ない方の値を、信号ランク、即ち到来波の数Lとする演算を行い、求められた到来波(信号)の数Lを、スペクトル演算部11に対して出力する。
ハウスホルダ変換演算回数Nhは、既に述べたように、列ノルムccの最大値が所定の閾値Th以下となった場合には、それ以上のハウスホルダ変換演算は、演算負荷が増大する反面、計算結果の信頼性が低くなることから実行しないによう制御されるので、各固有値分解演算のサイクルにおいてハウスホルダ変換演算が実行された回数Nhが信号ランク推定演算部9で推定された信号の数Ltよりも少ない場合には、推定された信号の数Ltには、ハウスホルダ変換演算Nhが行われていない列についての演算結果に基づいて信号の数Ltが判定されている部分が有るものと判断して、ハウスホルダ変換演算が行われた回数に等しい数を、信号ランク(到来波の数)として判定する。
特に、図3のステップS4における閾値Thを「0」ないしは「0」に近い値とした場合には、列ノルムccの最大値が「0」に近い値以下となることは、通常の固有値分解演算においてはあり得ないので、そうしたデータに基づいて判定された到来波の数を修正することが出来、信頼性の高い信号ランクの推定が可能となる。更に、閾値Thを「0」ないしは「0」に近い値とした場合、ハウスホルダ変換演算回数Nhは、固有値から推定される到来波の数よりも大きくなることはないので、演算された固有値から信号ランクが過大に判定されることを防止することが出来る。
こうして、信号ランク判定部10により、ハウスホルダ変換演算回数Nhが考慮された形で、信号ランクLが判定され、その値Lが、図1に示す、スペクトル演算部11に出力されると、スペクトル演算部11は、公知の手法により、図2に示すように、MUICスペクトラムSPEを演算して求め、当該求められたMUSICスペクトラムSPEを、到来波方位演算部12に出力する。
到来波方位演算部12は、当該MUSICスペクトラムSPEから、公知の手法で、到来波の到来する方位θ1、θ2を演算することが出来るが、MUSICスペクトラムSPEは、既に述べたように、信号の数Lが、相関行列Rの固有値分解演算に際したハウスホルダ変換演算回数Nhに基づいて補正され、信頼性の低いデータが排除されているので、実体のない所謂「偽電波」を排除した形で、正確な到来波の方位を算出することが出来る。
なお、固有値分解演算の演算処理の、図3のステップS4(又は図14のステップS200)における閾値Thとしては別の値を用いることが出来る。即ち、到来波において、各アンテナ15で雑音のみが捕捉された際の相関行列Rを求め、当該相関行列Rについて、固有値を求めるための行列Aについての雑音列ノルムを、実験などにより求めて演算しておく。次に、到来波の信号について、固有値を求めるための行列Aの信号列ノルムccを実験などにより求めて演算する。そして雑音列ノルムの値(平均値、最大値など)よりも大きく、信号列ノルムの値(平均値、最大値など)よりも小さな値を信号閾値Th′として設定することにより、ハウスホルダ変換演算回数Nhも、捕捉した本来の信号の数に応じた回数だけ、実行するように構成することも出来る。この場合、雑音に対応する列ノルムccの最大値(<Th)が、図3のステップS3で得られた時点で、ステップS4でハウスホルダ変換演算を中止することが出来、余分なハウスホルダ変換演算を行わずに済み、演算負荷を少なくすることが出来る。
また、この場合、図12に示すように、信号の数(到来波数)Lをハウスホルダ変換演算回数Nhと一致する形で判定することが出来るので、信号ランク推定部9及び信号ランク判定部10による信号の数の推定及び判定演算が不要となり、直ちにスペクトル演算部11によるMUSICスペクトラムの演算に入ることが出来、高速な処理動作が可能となる。
また、閾値Thを、信号列ノルムの値とは無関係に、雑音列ノルムの値のみを基準にして設定することも可能である。これにより、設定動作が簡単になる。更に、閾値Thを、到来波を受信する状況に応じて、変動させるようにすることも可能である。例えば、カーナビゲーションシステムなどを用いて周囲の環境を検出し、その結果、検出されたトンネル内などにおいては、レーダから発射した電波が、多くの車両や壁などに反射されて帰ってくることが予想される。こうした到来波の数が多くなることが予想され場合に、閾値Thを低く設定して、より多くの到来波を的確に捕捉し、それ以外の到来波の数が少ないことが予想される開放された路上などにおいては、閾値Thを高めに設定して、無駄な雑音列ノルムに係わる信号の演算動作を行わないようにし、効率的な演算動作を可能とすることが出来る。また、観測距離が遠い場合にも、到来波が多くの場所から反射してくるので、閾値を低くすることが効果的である。
なお、到来波の数、あるいは最大数が予め判っているような場合には、ハウスホルダ変換演算実行手段によるハウスホルダ変換演算の実行回数を予め当該数、あるいは当該最大数に設定しておき、無駄な演算動作が発生することを防止するように構成することも可能である。更に、各演算サイクル(図3のステップS2〜S8、又は図14におけるステップS100からS400の、「QR分解演算」と称された、点線で囲まれた演算サイクル)におけるハウスホルダ変換演算の実行回数を、周囲の環境に対応させて変動させることも可能である。例えば、カーナビゲーションシステムなどを用いて周囲の環境を検出し、その結果、検出されたトンネル内などにおいては、レーダから発射した電波が、多くの車両に反射されて帰ってくることが予想される。こうした到来波の数が多くなることが予想され場合に、ハウスホルダ変換演算の実行回数を多く設定して、より多くの到来波についての固有値を的確に演算し、それ以外の到来波の数が少ないことが予想される開放された路上などにおいては、ハウスホルダ変換演算の実行回数を少なく設定して、無駄な雑音列ノルムに係わる信号の固有値演算動作を行わないようにし、効率的な演算動作を可能とすることが出来る。
更に、固有値分解演算部などの収束演算実行手段は、N−1回目の所定回数のハウスホルダ変換演算からなる演算サイクル(QR分解演算)を行った後、N回目の演算サイクル(QR分解演算)を実行し、これにより、図3のステップS11に規定されている所定回数(NQR)の演算サイクル数だけ、QR分解演算を実行する。この際、N回目の演算サイクル(QR分解演算)におけるハウスホルダ変換演算の実行回数は、それ以前に実行された演算サイクル(QR分解演算)におけるハウスホルダ変換演算の実行回数以下に制御する手段を設けることが、望ましい。直前に行われた演算サイクル(QR分解演算)のハウスホルダ変換演算の実行回数を上回るハウスホルダ変換演算を、次の演算サイクル(QR分解演算)で実行しても、直前にハウスホルダ変換演算が行われていない行列部分(桁落ちや雑音部分など)に対してハウスホルダ変換演算を行うこととなり、演算の意味がないからである。
本発明は、車載用のレーダなど、簡易な構成で複数のアンテナ素子に到来する到来波の数、即ちランクを判定することの出来る、到来波の方向推定装置として、利用することが出来る。
図1は、到来電波方位測定装置の一例を示すブロック図。 図2は、信号処理の一例を示すチャート。 図3は、固有値分解演算処理の一例を示すフローチャート。 図4は、図3のフローチャートにおける、具体的な演算例の一例を示す図。 図5は、図3のフローチャートにおける、具体的な演算例の一例を示す図。 図6は、図3のフローチャートにおける、具体的な演算例の一例を示す図。 図7は、図3のフローチャートにおける、具体的な演算例の一例を示す図。 図8は、図3のフローチャートにおける、具体的な演算例の一例を示す図。 図9は、図3のフローチャートにおける、具体的な演算例の一例を示す図。 図10は、図3のフローチャートにおける、具体的な演算例の一例を示す図。 図11は、列ノルムの演算方法の一例を示す図。 図12は、到来波数の判定処理の別の例を示すチャート。 図13は、演算式の一例を示す図。 図14は、固有値分解演算の基本概念の一例を示すフローチャートである。
符号の説明
1……到来波方位測定装置
6……相関行列演算部
7……固有値分解演算部
9……信号ランク推定演算部
10……信号ランク判定部
15……アンテナ
Th……第1の閾値

Claims (7)

  1. 複数のアンテナで受信した複数の到来波の受信データから、前記到来波の到来方向を推定する、到来波の方向推定装置において、
    前記複数のアンテナで受信した、複数の到来波の受信データから、その相関行列を演算する相関行列演算手段、
    前記相関行列演算手段で演算された前記相関行列に対して、ハウスホルダ変換演算を実施して、前記相関行列を、固有ベクトルと固有値とに分解処理する固有値分解演算手段、
    前記固有値分解演算手段に設けられ、前記相関行列について、固有値を求めるための行列と固有ベクトルを求めるための行列を生成する行列生成手段、
    前記固有値分解演算手段に設けられ、前記固有値を求めるための行列において、前記ハウスホルダ変換演算を行うたび毎に、前記ハウスホルダ変換演算の対象となる行列に対して、各列の成分の列ノルムを演算する列ノルム演算手段、
    前記固有値分解演算手段に設けられ、前記ハウスホルダ変換演算を行うたび毎に、前記列ノルム演算手段により演算されたハウスホルダ変換演算の対象となる行列の列ノルムの最大値が、所定の第1の閾値以下であるか否かを判定する、列ノルム最大値判定手段、
    前記固有値分解演算手段に設けられ、前記列ノルム最大値判定手段により、前記列ノルム演算手段により演算されたハウスホルダ変換演算の対象となる行列の列ノルムの最大値が、前記所定の第1の閾値以下でないものと判定された場合に限り、前記固有値を求めるための行列に対してハウスホルダ変換演算を、繰り返し実行するハウスホルダ変換演算実行手段、
    前記固有値分解演算手段に設けられ、前記列ノルム最大値判定手段により、前記列ノルム演算手段により演算されたハウスホルダ変換演算の対象となる行列の列ノルムの最大値が、前記所定の第1の閾値以下であるものと判定された場合に、前記ハウスホルダ変換演算実行手段によるハウスホルダ変換演算の実行を中断すると共に、前記ハウスホルダ変換演算実行手段によるハウスホルダ変換演算の実行回数を所定のメモリに格納する変換演算実行制御手段、
    前記固有値分解演算手段に設けられ、前記変換演算実行制御手段によりハウスホルダ変換演算の実行が中断された前記固有値を求めるための行列に対して、前記ハウスホルダ変換演算実行手段及び変換演算実行制御手段によるハウスホルダ変換演算を、所定回数に達するまで更に実行させて、その結果から、前記相関行列の前記固有ベクトルと固有値を求める、収束演算実行手段、及び、
    前記求められた固有値を構成する値と所定の第2の閾値を比較し、当該第2の閾値よりも大きな値については到来波の信号の固有値とし、小さなものは雑音の固有値と判定して、到来波の数を推定演算すると共に、該推定演算された到来波の数と前記メモリ手段に格納されたハウスホルダ変換演算の実行回数を比較し、その少ない方の値を到来波の数と判定する、信号ランク判定手段、
    を有する、到来波の方向推定装置。
  2. 前記第1の閾値として、0又は0に近い値を格納するメモリ手段を有し、
    前記変換演算実行制御手段は、前記列ノルムの最大値が負の場合に、前記ハウスホルダ変換演算の更なる実行を中断させる、ことを特徴とする、請求項1記載の到来波の方向推定装置。
  3. 複数のアンテナで受信した複数の到来波の受信データから、前記到来波の到来方向を推定する、到来波の方向推定装置において、
    前記複数のアンテナで受信した、複数の到来波の受信データから、その相関行列を演算する相関行列演算手段、
    前記相関行列演算手段で演算された前記相関行列に対して、ハウスホルダ変換演算を実施して、前記相関行列を、固有ベクトルと固有値とに分解処理する固有値分解演算手段、
    前記固有値分解演算手段に設けられ、前記相関行列について、固有値を求めるための行列と固有ベクトルを求めるための行列を生成する行列生成手段、
    前記固有値分解演算手段に設けられ、前記固有値を求めるための行列において、前記ハウスホルダ変換演算を行うたび毎に、前記ハウスホルダ変換演算の対象となる行列に対して、各列の成分の列ノルムを演算する列ノルム演算手段、
    所定の第1の閾値を、前記受信データにおいて雑音のみが捕捉された際の、前記固有値を求めるための行列の列ノルムよりも大きく、前記受信データの到来波の信号についての、前記固有値を求めるための行列の列ノルムよりも小さな値となるように設定する閾値設定手段、
    前記固有値分解演算手段に設けられ、前記ハウスホルダ変換演算を行うたび毎に、前記列ノルム演算手段により演算されたハウスホルダ変換演算の対象となる行列の列ノルムの最大値が、前記所定の第1の閾値以下であるか否かを判定する、列ノルム最大値判定手段、
    前記固有値分解演算手段に設けられ、前記列ノルム最大値判定手段により、前記列ノルム演算手段により演算されたハウスホルダ変換演算の対象となる行列の列ノルムの最大値が、前記所定の第1の閾値以下でないものと判定された場合に限り、前記固有値を求めるための行列に対してハウスホルダ変換演算を、繰り返し実行するハウスホルダ変換演算実行手段、
    前記固有値分解演算手段に設けられ、前記列ノルム最大値判定手段により、前記列ノルム演算手段により演算されたハウスホルダ変換演算の対象となる行列の列ノルムの最大値が、前記所定の第1の閾値以下であるものと判定された場合に、前記ハウスホルダ変換演算実行手段によるハウスホルダ変換演算の実行を中断すると共に、前記ハウスホルダ変換演算実行手段によるハウスホルダ変換演算の実行回数を所定のメモリに格納する変換演算実行制御手段、
    前記固有値分解演算手段に設けられ、前記変換演算実行制御手段によりハウスホルダ変換演算の実行が中断された前記固有値を求めるための行列に対して、前記ハウスホルダ変換演算実行手段及び変換演算実行制御手段によるハウスホルダ変換演算を、所定回数に達するまで更に実行させて、その結果から、前記相関行列の前記固有ベクトルと固有値を求める、収束演算実行手段、及び、
    前記メモリ手段に格納されたハウスホルダ変換演算の実行回数を到来波の数と判定する、信号ランク判定手段、
    を有することを特徴とする、到来波の方向推定装置。
  4. 前記閾値設定手段は、前記所定の第1の閾値を、周囲の環境や観測距離に対応させて変動的に設定することを特徴とする、請求項3記載の到来波の方向推定装置。
  5. 前記ハウスホルダ変換演算実行手段によるハウスホルダ変換演算の実行回数の最大数を設定する演算回数設定手段を有し、
    前記ハウスホルダ変換演算実行手段によるハウスホルダ変換演算を、前記演算回数設定手段に設定されたハウスホルダ変換演算の実行回数の最大数以内で行うことを特徴とする、請求項1又は3記載の到来波の方向推定装置。
  6. 前記演算回数設定手段は、前記ハウスホルダ変換演算の実行回数の最大数を、周囲の環境や観測距離に対応させて変動的に設定することを特徴とする、請求項5記載の到来波の方向推定装置。
  7. 前記収束演算実行手段による、前記ハウスホルダ変換演算実行手段及び変換演算実行制御手段による前記ハウスホルダ変換演算の、各回の実行回数を制御する変換演算回数制御手段を設け、
    前記変換演算回数制御手段は、N回目の前記ハウスホルダ変換演算実行手段及び変換演算実行制御手段によるハウスホルダ変換演算の実行回数を、それ以前の回におけるハウスホルダ変換演算の実行回数以下に制御することを特徴とする、請求項1または3記載の到来波の方向推定装置。
JP2005100314A 2005-03-31 2005-03-31 到来波の方向推定装置 Active JP4808984B2 (ja)

Priority Applications (3)

Application Number Priority Date Filing Date Title
JP2005100314A JP4808984B2 (ja) 2005-03-31 2005-03-31 到来波の方向推定装置
DE102006013448.6A DE102006013448B4 (de) 2005-03-31 2006-03-20 Vorrichtung zur Schätzung der Richtung einer ankommenden Welle
US11/390,866 US7372404B2 (en) 2005-03-31 2006-03-27 Apparatus for estimating direction of arrival wave

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2005100314A JP4808984B2 (ja) 2005-03-31 2005-03-31 到来波の方向推定装置

Publications (2)

Publication Number Publication Date
JP2006284180A JP2006284180A (ja) 2006-10-19
JP4808984B2 true JP4808984B2 (ja) 2011-11-02

Family

ID=37071873

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2005100314A Active JP4808984B2 (ja) 2005-03-31 2005-03-31 到来波の方向推定装置

Country Status (3)

Country Link
US (1) US7372404B2 (ja)
JP (1) JP4808984B2 (ja)
DE (1) DE102006013448B4 (ja)

Families Citing this family (22)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1903348B1 (en) * 2005-07-11 2015-12-30 Fujitsu Ltd. Number-of-arriving-waves estimating method, number-of-arriving-waves estimating device, and radio device
US7839835B2 (en) * 2006-08-22 2010-11-23 Nec Laboratories America, Inc. Quantized precoding over a set of parallel channels
JP4531738B2 (ja) * 2006-11-21 2010-08-25 日本電信電話株式会社 行列における数値分解方法
KR100932789B1 (ko) * 2007-12-15 2009-12-21 한국전자통신연구원 다중입력 다중출력 시스템에서 qr 분해 장치 및 그 방법
JP5600866B2 (ja) * 2008-03-04 2014-10-08 富士通株式会社 探知測距装置および探知測距方法
JP5519132B2 (ja) * 2008-07-28 2014-06-11 株式会社デンソー レーダ装置
JP5102165B2 (ja) * 2008-09-22 2012-12-19 株式会社デンソー レーダ装置
JP5511479B2 (ja) * 2010-04-16 2014-06-04 三菱電機株式会社 測定装置
US8224627B2 (en) 2010-04-28 2012-07-17 Mitov Iliya P Technique for determination of the signal subspace dimension
US20110286306A1 (en) * 2010-05-21 2011-11-24 Leo Eisner Determining origin and mechanism of microseismic events in the earth's subsurface by deviatoric moment inversion
US9294310B2 (en) * 2010-10-08 2016-03-22 Blackberry Limited Method and apparatus for LTE channel state information estimation
KR20120071851A (ko) * 2010-12-23 2012-07-03 한국전자통신연구원 통신 시스템에서 신호 도래 방향 추정 장치 및 방법
JP5628732B2 (ja) * 2011-04-04 2014-11-19 富士通テン株式会社 レーダ装置用の演算装置、レーダ装置、レーダ装置用の演算方法およびプログラム
JP6195490B2 (ja) * 2013-08-23 2017-09-13 日本放送協会 Mimo受信装置
NL2014562A (en) * 2014-04-04 2015-10-13 Asml Netherlands Bv Control system, positioning system, lithographic apparatus, control method, device manufacturing method and control program.
CN104297594A (zh) * 2014-10-13 2015-01-21 电子科技大学 一种智能家居的家电设备运行状态监测方法
RU2604871C2 (ru) * 2015-04-15 2016-12-20 федеральное государственное автономное образовательное учреждение высшего профессионального образования "Южный федеральный университет" (Южный федеральный университет") Способ определения местоположения объекта навигации
US11067659B2 (en) 2016-12-15 2021-07-20 Bae Systems Information And Electronic Systems Integration Inc. System and method for rank estimation of electromagnetic emitters
US10473746B1 (en) * 2016-12-15 2019-11-12 Bae Systems Information And Electronic Systems Integration Inc. System and method for rank estimation of electromagnetic emitters
US11255955B2 (en) * 2018-12-28 2022-02-22 Panasonic Intellectual Property Management Co., Ltd. Estimation method, estimation device, and recording medium
CN111381228A (zh) * 2018-12-28 2020-07-07 松下知识产权经营株式会社 推测方法以及推测装置
US10986509B1 (en) * 2020-05-14 2021-04-20 At&T Intellectual Property I, L.P. Placement of antennas for fifth generation (5G) or other next generation networks

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH09212489A (ja) * 1996-01-31 1997-08-15 Fujitsu Ltd 対称行列の固有値問題を解く並列処理装置および方法
JP3081522B2 (ja) * 1996-02-14 2000-08-28 株式会社エイ・ティ・アール光電波通信研究所 受信信号処理装置
JPH11344517A (ja) 1998-06-02 1999-12-14 Nec Corp 電波環境分析装置
JP2000121716A (ja) * 1998-10-13 2000-04-28 Anritsu Corp 電波到来方向推定装置
JP2005176020A (ja) * 2003-12-12 2005-06-30 Rikogaku Shinkokai 復号方法および復号装置

Also Published As

Publication number Publication date
JP2006284180A (ja) 2006-10-19
DE102006013448B4 (de) 2022-06-09
DE102006013448A1 (de) 2007-01-18
US7372404B2 (en) 2008-05-13
US20060224655A1 (en) 2006-10-05

Similar Documents

Publication Publication Date Title
JP4808984B2 (ja) 到来波の方向推定装置
Bazzi et al. Detection of the number of superimposed signals using modified MDL criterion: A random matrix approach
Athley Threshold region performance of maximum likelihood direction of arrival estimators
JP4138825B2 (ja) ウェイト算出方法、ウェイト算出装置、アダプティブアレーアンテナ、及びレーダ装置
US6567034B1 (en) Digital beamforming radar system and method with super-resolution multiple jammer location
US7576682B1 (en) Method and system for radar target detection and angle estimation in the presence of jamming
EP1425607A1 (en) Radar system and method including superresolution raid counting
CN107870315B (zh) 一种利用迭代相位补偿技术估计任意阵列波达方向方法
WO2006067869A1 (ja) 到来方向推定装置及びプログラム
JP2008096137A (ja) レーダ装置及び測角装置
CN109471063B (zh) 基于延迟快拍的均匀线列阵高分辨波达方向估计方法
CN109901103B (zh) 基于非正交波形的mimo雷达doa估算方法及设备
CN116881385B (zh) 轨迹平滑方法、装置、电子设备及可读存储介质
CN117092585B (zh) 单比特量化DoA估计方法、***和智能终端
Yan et al. Computationally efficient direction of arrival estimation with unknown number of signals
WO2010066306A1 (en) Apparatus and method for constructing a sensor array used for direction of arrival (doa) estimation
CN106844886B (zh) 基于主分量分析的目标波达方向获取方法
JP5152949B2 (ja) ウェイト算出方法、ウェイト算出装置、アダプティブアレーアンテナ、及びレーダ装置
CN115421098A (zh) 嵌套面阵下降维求根music的二维doa估计方法
JP2008032437A (ja) ウェイト算出方法、ウェイト算出装置、アダプティブアレーアンテナ、及びレーダ装置
CN108051773A (zh) 基于盖式圆盘准则估计信源数目的epuma方法
CN109683128B (zh) 冲击噪声环境下的单快拍测向方法
Roemer et al. Higher order SVD based subspace estimation to improve multi-dimensional parameter estimation algorithms
CN113281698A (zh) 一种嵌套阵中基于级联的非高斯信源测向方法
CN111880167A (zh) 一种基于先随机后优化的波达方向估计方法

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20070614

A977 Report on retrieval

Free format text: JAPANESE INTERMEDIATE CODE: A971007

Effective date: 20100402

A131 Notification of reasons for refusal

Free format text: JAPANESE INTERMEDIATE CODE: A131

Effective date: 20101026

A521 Request for written amendment filed

Free format text: JAPANESE INTERMEDIATE CODE: A523

Effective date: 20101214

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: 20110726

A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20110818

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

Free format text: PAYMENT UNTIL: 20140826

Year of fee payment: 3

R150 Certificate of patent or registration of utility model

Ref document number: 4808984

Country of ref document: JP

Free format text: JAPANESE INTERMEDIATE CODE: R150

Free format text: JAPANESE INTERMEDIATE CODE: R150

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

Free format text: PAYMENT UNTIL: 20140826

Year of fee payment: 3

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250