(57)【要約】【課題】本発明は、マルチパス環境下においても、ＭＵＳＩＣ法よりも多くの到来波の電波到来方向を推定でき、また、到来波が強い相関を持つ場合においても、ＭＯＤＥよりも電波到来方向の高い推定精度を維持することができる電波到来方向推定装置を提供する。【解決手段】本発明に係る電波到来方向推定装置１７０は、アレイアンテナ１０２の受信信号から得られる共分散行列を固有値分解する固有値分解部１０６を備える。また、前記固有値から所定の固有値を抽出する固有値抽出部１０７を備える。また、所定の固有値に対応する固有ベクトルを用いて構成される固有空間の基底と、固有値分解部１０６で求められた固有値の個数以上の数の信号成分で構成される信号空間との、直交成分が最小となるという条件を用いて、到来角度を推定する最適化部１０８を備える。【選択図】  図１
2010121991.tif 000002
【特許請求の範囲】
【請求項１】
電波を受信するアレイアンテナと、
前記アレイアンテナにより受信された電波を検波する受信部と、
前記受信部から送信されるアナログ信号をデジタル信号に変換するＡ／Ｄ変換部と、
前記Ａ／Ｄ変換部から送信されるデジタル信号から、共分散行列を求める共分散行列計算部と、
前記共分散行列の固有値を求める固有値分解部と、
前記固有値分解部によって得られた前記固有値のうち、所定の前記固有値を抽出する固有値抽出部と、
前記アレイアンテナに対する前記電波の到来角度を推定する最適化部とを、備えており、
前記最適化部は、
前記所定の固有値に対応する固有ベクトルを用いて構成される固有空間の基底と、前記固有値分解部で求められた前記固有値の個数以上の数の信号成分で構成される信号空間との、直交成分が最小となるという条件を用いて、前記到来角度を推定する、
ことを特徴とする電波到来方向推定装置。
【請求項２】
前記信号空間は、
前記最適化部で予め設定されている（前記アレイアンテナの素子数−前記電波の到来数）行×（前記アレイのアンテナ数）列から成る行列、で規定される空間である、
ことを特徴とする電波到来方向推定装置。
【請求項３】
前記所定の固有値は、
前記固有値分解部で求められた前記固有値のうち、最大の値を有する最大固有値である、
ことを特徴とする請求項１に記載の電波到来方向推定装置。
【請求項４】
前記最適化部は、
（Ｍ−Ｎ）×Ｍ行列Ｂbの複素共役転置行列を、式（Ｄ）として定義しており、
（ａ）式（Ａ）で表される行列のフロベニウスノルムが最小化となる、前記行列Ｂbの前記行列要素を求め、
（ｂ）前記行列Ｂbの行列要素ｂ0の偏角γを求め、
（ｃ）前記行列Ｂbの全行列要素を係数として構成される多項式である式（Ｂ）を用いて、前記多項式の解の偏角φnを求め、
（ｄ）前記偏角φnと式（Ｃ）とを用いて、前記到来角度を推定し、
さらに、前記最適化部は、
前記工程（ａ）および前記工程（ｂ）を、当該順に繰り返し実施し、
前記繰り返し実施の度に求まる前記行列Ｂbの全行列要素の変化が収束したとき、前記繰り返しを終了し、
前記繰り返し実施の度に求まる前記偏角γの値に応じて、式（Ｅ）および式（Ｆ）の条件を切替えて用い、
前記収束時点で得られた前記行列Ｂbの全行列要素を用いて、前記工程（ｃ）および前記工程（ｄ）を実施する、
ここで、前記式（Ａ）、（Ｂ）、（Ｃ）、（Ｄ）、（Ｅ）および（Ｆ）は、
【数１】
2010121991.tif 000033
【数２】
2010121991.tif 000034
【数３】
2010121991.tif 000035
【数４】
2010121991.tif 000036
【数５】
2010121991.tif 000037
であり、
Ｍは、前記アレイアンテナの素子数であり、
ｓVは、前記基底であり、
Ｎは、到来電波数であり、
θnは、前記到来角度であり、
ｃは、光速度であり、
πは、円周率であり、
ｆは、前記電波の搬送周波数であり、
ｄは、前記アンテナアレイの素子間隔である、
ことを特徴とする請求項３に記載の電波到来方向推定装置。
【請求項５】
前記所定の固有値は、
前記固有値分解部で求められた前記固有値のうち、最大値から予め規定されている所定の番目までの大きさの値を有する、複数の固有値である、
ことを特徴とする請求項１に記載の電波到来方向推定装置。
【請求項６】
前記所定の固有値は、
前記固有値分解部で求められた前記固有値のうち、予め規定されている所定のしきい値を用いて決定されるしきい値固有値より大きな値を有する、複数の固有値である、
ことを特徴とする請求項１に記載の電波到来方向推定装置。
【請求項７】
前記最適化部は、
（Ｍ−Ｎ）×Ｍ行列Ｂbの複素共役転置行列を、式（Ｄ）として定義しており、
（ａ）式（Ａ）で表される行列のフロベニウスノルムが最小化となる、前記行列Ｂbの全行列要素を求め、
（ｂ）前記行列Ｂbの行列要素ｂ0の偏角γを求め、
（ｃ）前記行列Ｂbの全行列要素を係数として構成される多項式である式（Ｂ）を用いて、前記多項式の解の偏角φnを求め、
（ｄ）前記偏角φnと式（Ｃ）とを用いて、前記到来角度を推定し、
さらに、前記最適化部は、
前記工程（ａ）および前記工程（ｂ）を、当該順に繰り返し実施し、
前記繰り返し実施の度に求まる前記行列Ｂbの全行列要素の変化が収束したとき、前記繰り返しを終了し、
前記繰り返し実施の度に求まる前記偏角γの値に応じて、式（Ｅ）および式（Ｆ）の条件を切替えて用い、
前記収束時点で得られた前記行列Ｂbの全行列要素を用いて、前記工程（ｃ）および前記工程（ｄ）を実施する、
ここで、前記式（Ａ）、（Ｂ）、（Ｃ）、（Ｄ）、（Ｅ）および（Ｆ）は、
【数６】
2010121991.tif 000038
【数７】
2010121991.tif 000039
【数８】
2010121991.tif 000040
【数９】
2010121991.tif 000041
【数１０】
2010121991.tif 000042
であり、
Ｍは、前記アレイアンテナの素子数であり、
ｓVは、前記基底であり、
Ｎは、到来電波数であり、
θnは、前記到来角度であり、
ｃは、光速度であり、
πは、円周率であり、
ｆは、前記電波の搬送周波数であり、
ｄは、前記アンテナアレイの素子間隔である、
ことを特徴とする請求項４または請求項５に記載の電波到来方向推定装置。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、電波到来方向推定装置に係る発明であり、たとえば、マルチパス環境下（特に、到来経路長の差が小さいマルチパスが存在する環境下）で、電波到来方向を推定することができる電波到来方向推定装置に適用することができる。
【背景技術】
【０００２】
基地局装置が受信する信号は、様々な方向から到来する信号による干渉を受ける。当該干渉により、受信信号の品質劣化が生ずる。当該干渉を抑制し、所望の方向から到来する信号（電波）のみを強く受信する方法として、アレイアンテナを用いる受信方法が知られている。
【０００３】
アレイアンテナは、指向性アンテナの一種であり、小さなアンテナを平面状に多数配列してなる。当該小さなアンテナの数は、多いときには数千個にも上り、個々のアンテナの出力は微弱である。しかしながら、実際には複数のアンテナからの出力が合成されるので、アレイアンテナ全体としては結果として大出力が得られる。
【０００４】
それぞれのアンテナから同時に電波が出力される場合は、合成波はアンテナ面と直角な方向に向かって進むが、アンテナ毎に電波の位相をずらすと、合成波はアンテナ面から斜め方向に出力される。これとは逆に電波受信の場合には、アンテナ毎の位相の差を調べると、電波の入射方向を調べることができる。
【０００５】
当該アレイアンテナを用いる受信方法では、各アンテナ素子での受信信号に乗算する重み付け係数を調整し、アンテナ素子毎の受信信号に対して与える振幅と位相を調整する。これにより、所望の方向から到来する信号のみを強く受信することができる。
【０００６】
電波到来方向推定を行う従来技術として、非特許文献１に示されているＭＵＳＩＣ（Ｍｕｌｔｉｐｌｅ  Ｓｉｇｎａｌ  Ｃｌａｓｓｉｆｉｃａｔｉｏｎ）法が知られている。また、非特許文献２には、空間平均法をＭＵＳＩＣ法とともに用いる手法が示されている。さらに、電波到来方向推定を行う技術として、非特許文献３には、ＭＯＤＥ（Ｍｅｔｈｏｄ  ｏｆ  Ｄｉｒｅｃｔｉｏｎ  Ｅｓｔｉｍａｔｉｏｎ）に関する記述がある。
【０００７】
さらに、電波到来方向推定を行う技術として、他に、ビームフォーマ（ｂｅａｍ ｆｏｒｍｅｒ）法、Ｃａｐｏｎ法、線型予測法等の様々な方法が存在する。
【０００８】
【非特許文献１】Ｒ．Ｏ．Ｓｃｈｍｉｄｔ，“Ｍｕｌｔｉｐｌｅ  ｅｍｉｔｔｅｒ  ｌｏｃａｔｉｏｎ  ａｎｄ  ｓｉｇｎａｌ  ｐａｒａｍｅｔｅｒ  ｅｓｔｉｍａｔｉｏｎ，”ＩＥＥＥ．Ｔｒａｎｓ．Ａｎｔｅｎｎａｓ ａｎｄ Ｐｒｏｐａｇａｔｉｏｎ，ｖｏｌ．３４  ｐｐ．２７６−２８０  １９８６
【非特許文献２】Ｋ．Ｔａｋａｏ，Ｎ．Ｋｉｋｕｍａ，“Ａｎ  ａｄａｐｔｉｖｅａｒｒａｙ  ｕｔｉｌｉｚｉｎｇ  ａｎ  ａｄａｐｔｉｖｅ  ｓｐａｔｉａｌ  ａｖｅｒａｇｉｎｇ  ｔｅｃｈｎｉｑｕｅ  ｆｏｒ  ｍｕｌｔｉｐａｔｈ  ｅｎｖｉｒｏｎｍｅｎｔｓ，” ＩＥＥＥ．Ｔｒａｎｓ．Ｏｎ  Ａｎｔｅｎｎａｓ  ａｎｄ  Ｐｒｏｐａｇａｔｉｏｎ，ｖｏｌ．３５，Ｎｏ．１２，１９８７
【非特許文献３】Ｐ．Ｓｔｏｉｃａ  ａｎｄ  Ｋ．Ｃ．Ｓｈａｒｍａｎ，“Ｎｏｖｅｌ  ｅｉｇｅｎａｎａｌｙｓｉｓ  ｍｅｔｈｏｄ  ｆｏｒ  ｄｉｒｅｃｔｉｏｎ  ｅｓｔｉｍａｔｉｏｎ,，”ＩＥＥＥ．Ｐｒｏｃ．Ｒａｄａｒ  ａｎｄ  Ｓｉｇｎａｌ  Ｐｒｏｃ． ｖｏｌ．１３７．ｐｐ． １９−２６  １９９０
【発明の開示】
【発明が解決しようとする課題】
【０００９】
到来経路長の差が小さいマルチパスにおいては、反射波は直接波と強い相関を持っている。この場合、ＭＵＳＩＣ法では正しい方向を推定するこができない。到来波が互いに強い相関を持つ場合の電波到来方向推定方法として、非特許文献２に示されている、空間平均法をＭＵＳＩＣ法とともに用いる手法が存在する。
【００１０】
しかしながら、当該非特許文献２に示されている方法では、同じ素子数のアレイアンテナに対し推定可能な到来波数が少なくなり、また推定精度が理論的な限界に漸近しない、という問題点がある。
【００１１】
また、上記ＭＯＤＥでは、到来波数が少なくなる問題点は無く、また相関のある到来波にも対応している。しかしながら、到来波が非常に強い相関を持つ場合に、推定精度が悪化する。
【００１２】
そこで、本発明は、到来経路長の差が小さいマルチパスが発生する場合においても、ＭＵＳＩＣ法よりも多くの到来波の電波到来方向を推定でき、また、到来波が強い相関を持つ場合においても、ＭＯＤＥよりも電波到来方向の高い推定精度を維持することができる電波到来方向推定装置を提供することを目的とする。
【課題を解決するための手段】
【００１３】
上記の目的を達成するために、本発明に係る請求項１に記載の電波到来方向推定装置は、電波を受信するアレイアンテナと、前記アレイアンテナにより受信された電波を検波する受信部と、前記受信部から送信されるアナログ信号をデジタル信号に変換するＡ／Ｄ変換部と、前記Ａ／Ｄ変換部から送信されるデジタル信号から、共分散行列を求める共分散行列計算部と、前記共分散行列の固有値を求める固有値分解部と、前記固有値分解部によって得られた前記固有値のうち、所定の前記固有値を抽出する固有値抽出部と、前記アレイアンテナに対する前記電波の到来角度を推定する最適化部とを、備えており、前記最適化部は、前記所定の固有値に対応する固有ベクトルを用いて構成される固有空間の基底と、前記固有値分解部で求められた前記固有値の個数以上の数の信号成分で構成される信号空間との、直交成分が最小となるという条件を用いて、前記到来角度を推定する。
【発明の効果】
【００１４】
本発明の請求項１に記載の電波到来方向推定装置は、上記構成の固有値抽出部および上記構成の最適化部を含んでいる。したがって、マルチパス環境下、特に到来距離差の小さなマルチパスが存在する環境下においても、ＭＵＳＩＣ法よりも多い、「アレイアンテナの素子数−１」個の電波の到来方向を精度良く推定することができる。
【００１５】
また、当該電波到来方向推定装置は、上記構成の固有値抽出部を含んでいるので、到来波が相互に強い相関を持つ場合においても、ＭＯＤＥよりも高い電波到来方向の推定精度が保たれる。
【発明を実施するための最良の形態】
【００１６】
以下、この発明をその実施の形態を示す図面に基づいて具体的に説明する。
【００１７】
＜実施の形態１＞
図１は、本実施の形態に係る電波到来方向推定装置１７０の概略構成を示すブロック図である。
【００１８】
図１に示すように、電波到来方向推定装置１７０は、アレイアンテナ１０２、受信器１０３、Ａ／Ｄ変換部１０４、共分散行列計算部１０５、固有値分解部１０６、固有値抽出部１０７、最適化部１０８および表示部１０９を、少なくとも備えている。
【００１９】
アレイアンテナ１０２は、指向性アンテナの一種であり、小さなアンテナを平面状に多数配列してなる。各小さなアンテナ（アレイアンテナの素子と把握できる）は、異なる方向から到来する複数の到来波（電波と把握できる）１０１を受信する。各小さなアンテナには、所定の到来角度θ1，θ2で、各到来波１０１が入射される。ここで、本願では、到来角度は、アレイアンテナ１０２の各素子（各小さなアンテナ）の配列方向に直交する方向を０ラジアンとして、規定されることとする。
【００２０】
各小さなアンテナには、受信部１０３が各々接続されている。受信部１０３は、アレイアンテナ１０２により受信された到来波（電波）１０１を検波する。つまり、受信部１０３は、到来波１０１を検波し、当該到来波１０１を、Ｉ（同相信号）Ｑ（直交信号）信号に変換する。
【００２１】
各受信部１０３には、Ａ／Ｄ変換部１０４が各々接続されている。Ａ／Ｄ変換部１０４は、受信部１０３から送信される検波後のアナログ信号であるＩＱ信号を、複素デジタル信号ｙ（ｔ）に変換する。各Ａ／Ｄ変換部１０４からの各複素デジタル信号ｙ（ｔ）は、一つの共分散行列計算部１０５へと送信される。
【００２２】
複素デジタル信号ｙ（ｔ）は、次の通りである。
【００２３】
素子数Ｍ（本願では、Ｍはアレイアンテナ１０２の素子の数を示すものとする）のリニアアレイアンテナ１０２に、N個（本願では、Ｎはアレイアンテナ１０２に到来してくる到来電波数を示すものとする）の電波が各々異なる方向から到来しているとする。第ｎ到来波の到来方向をθn（ラジアン）とおく。このとき、第ｍ素子で受信される信号ｙm（ｔ）は、式（１）で表される。
【００２４】
【数１】
2010121991.tif 000003
【００２５】
ここで、ｔは、時間であり、ｊは、虚数単位である。ｗn（ｔ）は、第ｎ到来波の信号であり、ｘm（ｔ）は、第ｍ素子で受信されるノイズである。また、φnは、式（２）で定義される角度（偏角であると把握できる）である。
【００２６】
【数２】
2010121991.tif 000004
【００２７】
ここで、ｃは、光速（ｍ／ｓ）であり、πは、円周率であり、ｆは、搬送波周波数（Ｈｚ）であり、ｄは、アレイアンテナ１０２の素子間隔（ｍ）である。
【００２８】
上記式（１）を行列でまとめて表すと、式（３）となる。
【００２９】
【数３】
2010121991.tif 000005
【００３０】
ただし、
【００３１】
【数４】
2010121991.tif 000006
【００３２】
【数５】
2010121991.tif 000007
【００３３】
【数６】
2010121991.tif 000008
【００３４】
であり、行列Ａφは、式（４）で定義される。
【００３５】
【数７】
2010121991.tif 000009
【００３６】
共分散行列計算部１０５では、Ａ／Ｄ変換部１０４から送信される複素デジタル信号ｙ（ｔ）から、共分散行列を求める。具体的に、共分散行列計算部１０５は、式（５）で定義される共分散行列Ｒを計算する。
【００３７】
【数８】
2010121991.tif 000010
【００３８】
ここで、Ｔは、信号のサンプル数である。
【００３９】
当該計算された共分散行列Ｒは、共分散行列計算部１０５から固有値分解部１０６へと送信される。そして、固有値分解部１０６では、共分散行列Ｒを固有値分解する。つまり、固有値分解部１０６では、共分散行列Ｒの固有値を求める。
【００４０】
その後、本実施の形態に係る固有値抽出部１０７では、固有値分解部１０６によって得られた固有値のうち、最大の固有値（最大固有値であり、本実施の形態における所定の固有値と把握できる）を抽出する。そして、固有値抽出部１０７は、当該最大の固有値に対応する固有ベクトルｓVを、最適化部１０８へ送信する。
【００４１】
最適化部１０８は、アレイアンテナ１０２に対する電波の到来角度を推定する到来角度推定部（図２の符号２１２参照）を備えている。最適化部１０８では、上記最大の固有値に対応する固有ベクトルｓVを用いて構成される固有空間の基底と、固有値分解部１０６で求められた固有値の個数以上の数の信号成分で構成される信号空間との、直交成分が最小となるという条件を用いて、前記到来角度を推定する。
【００４２】
具体的に、上記最適化部１０８では、下式（６）の条件の下、固有空間の基底（本実施の形態では、基底は固有ベクトルｓVである）の行列Ｂで張られる空間への射影であるフロベニウスノルム（すなわち、下式（７）で表される目的関数ＦD（ｂ））を、最小化するｂ（行列Ｂbを構成する全行列要素ｂ0からｂN。以下同様。）を求める。
【００４３】
【数９】
2010121991.tif 000011
【００４４】
【数１０】
2010121991.tif 000012
【００４５】
ここで、式（６）および式（７）の「ｂ」は、（Ｍ−Ｎ）×Ｍ行列Ｂbの要素である。また、行列Ｂbは、最適化部１０８に予め設定されている。「（Ｂb*Ｂb）-1/2Ｂb*」は、固有値分解部１０６で求められた固有値の個数以上の数の信号成分で構成される信号空間であると把握でき、行列Ｂbの複素共役転置行列Ｂb*は、下式（８）で定義される。換言すれば、当該信号空間は、（Ｍ−Ｎ）×Ｍの行列で規定される空間である。また、式（６）の「ｂN-1*」は、「ｂN-1」の複素共役である。また、式（６）の「Ｒｅｂ0」は、「ｂ0」の実数部であり、「Ｉｍｂ0」は、「ｂ0」の虚数部である。
【００４６】
【数１１】
2010121991.tif 000013
【００４７】
さらに、上記最適化部１０８では、上記で求められた行列Ｂの要素ｂから到来角度θ＝［θ1・・・θN］を計算する。このようにして得られた到来角度θは、表示または他の目的使用のため、表示部１０９やその他の回路ブロックに送信される。
【００４８】
次に、最適化部１０８の内部機能構成を示す図２のブロック図を用いて、上記到来角度θの計算方法を詳細に説明する。
【００４９】
固有値抽出部１０７から送信された最大の固有値に対応する固有ベクトルｓVは、２次最適化手段２０２に送信される。２次最適化手段２０２は、拘束条件式「ｂi＝ｂN-i*  （ｉ＝０，１，・・・，Ｎ）、および、Ｒｅｂ0＝１」を用いて、下式（９）を最小化するｂを求める。
【００５０】
【数１２】
2010121991.tif 000014
【００５１】
ただし、２次最適化手段２０２における上記固有ベクトルｓV受信後の最初の計算では、式（９）におけるＷは、単位行列とする。
【００５２】
式（９）は、拘束条件と合わせて上記ｂに関する２次最適化問題となっているため、ＱＲ分解などを用いて最適点を求めることができる。２次最適化手段２０２で求めたｂは拘束条件判定手段２０３へ送信される。
【００５３】
拘束条件判定手段２０３では、γをｂ0の偏角とし、ｔａｎ（γ）＜１ の場合は、上記ｂをＷ計算手段２０４へ送信する。これに対して、ｔａｎ（γ）＞１の場合には、上記ｂを制約切替器２０６へ送信する。Ｗ計算手段２０４は、式（８）を用いて行列Ｂbを構成し、Ｗを下式（１０）で定義する。
【００５４】
【数１３】
2010121991.tif 000015
【００５５】
上記式（１０）を定義した後、収束判定手段２０５は、行列Ｂbの要素ｂの変化を追跡し、収束を判定する。当該ｂが収束していないと収束判定手段２０５が判断すれば、式（１０）で定義されたＷを２次最適化手段２０２へ送信する。そして、２次最適化手段２０２は、拘束条件式「ｂi＝ｂN-i*  （ｉ＝０，１，・・・，Ｎ）、および、Ｒｅｂ0＝１」を用いて、下式（９）を最小化するｂを再度計算する。ここで、当該再度の計算の際に、式（９）において行列Ｗは上記ｂを用いて構成されたものであるが、最適化の際は当該行列Ｗは最適化の変数には含めない。
【００５６】
ところで、制約切替器２０６は、式（８）および式（１０）を用いて行列Ｗを構成し、当該構成した行列Ｗを２次最適化手段２０７へ送信する。
【００５７】
２次最適化手段２０７では、拘束条件式「ｂi＝ｂN-i*  （ｉ＝０，１，・・・，Ｎ）、および、Ｉｍｂ0＝１」を用いて、式（９）を最小化する新たなｂを求める。その後、拘束条件判定手段２０８は、γをｂ0の偏角として、ｔａｎ（γ）＞１ の場合は、上記新たなｂをＷ計算手段２０９へ送信する。これに対して、ｔａｎ（γ）＜１の場合には、上記新たなｂを制約切替器２１１へ送信する。
【００５８】
なお、Ｗ計算手段２０９の動作は、Ｗ計算手段２０４の動作と同様である。また、収束判定手段２１０の動作は、収束判定手段２０５の動作と同様である。また、制約切替器２１１は、式（８）および式（１０）を用いて行列Ｗを構成し、当該構成した行列Ｗを２次最適化手段２０２へ送信する。
【００５９】
収束判定手段２０５または収束判定手段２１０によって上述したｂの収束が検知されると、その時点でｂを到来角計算手段２１２へ送信する。到来角計算手段２１２では、上述したｂを用いて、下式（１１）で表される多項式の零点ｚ1，・・・，ｚNをすべて求める。
【００６０】
【数１４】
2010121991.tif 000016
【００６１】
式（１１）は、送信された行列Ｂbの要素ｂを係数として構成される多項式である。つまり、ｚを不定項、ｂkを係数とする多項式である。当該多項式の零点とは、すなわち、下式（１２）の解ｚを表す。
【００６２】
【数１５】
2010121991.tif 000017
【００６３】
到来角度計算手段２１２は、当該式（１２）を用いて、上記多項式の解の偏角φnを、まず求める。具体的に、式（１２）は、複素係数のＮ＋１次多項式なので、式（１２）の解はＮ＋１個の複素数となる。当該Ｎ＋１個の複素数をそれぞれｚnと表したとき、
【００６４】
【数１６】
2010121991.tif 000018
【００６５】
のように極座標表示されたときのφn（偏角）｛＝φ1，・・・，φN｝を、到来角度計算手段２１２は求める。ここで、ｊは虚数単位である。多項式の求解は、一般的なアルゴリズムが多数存在する。簡単なものでは、行列の固有値問題に帰着し、ＱＲ分解などによって解くことができる。
【００６６】
上記零点の偏角φn｛φ1，・・・，φN｝導出後、到来角度計算手段２１２は、当該偏角φnと式（１３）とを用いて、到来角度θnを推定する。
【００６７】
【数１７】
2010121991.tif 000019
【００６８】
ここで、上述したように、ｃは、光速（ｍ／ｓ）であり、πは、円周率であり、ｆは、搬送波周波数（Ｈｚ）であり、ｄはアレイアンテナ１０２の素子間隔（ｍ）である。
【００６９】
最適化部１０８（より具体的には、到来角度計算手段２１２）は、式（１３）により導出された到来角度を、表示部１０９などに出力する。
【００７０】
ここで、２次最適化手段２０２の内部機能構成を示す図３のブロック図を用いて、上記ｂの出力までの動作をより詳細に説明する。
【００７１】
Ｓ計算手段３０４は、上記最大の固有値に対応する上記固有ベクトルｓV（下式（１４））を入力し、行列Ｓを下式（１５）で定義する。
【００７２】
【数１８】
2010121991.tif 000020
【００７３】
【数１９】
2010121991.tif 000021
【００７４】
すると 、
【００７５】
【数２０】
2010121991.tif 000022
【００７６】
となる。したがって、式（１５）で定義されるＳを用いると、式（９）は、式（１６）に変形される。
【００７７】
【数２１】
2010121991.tif 000023
【００７８】
なお、Ｓ計算手段３０４は、到来角度の推定に際し、上記行列Ｓの計算は一度だけ行えばよい。
【００７９】
コレスキー分解手段３０３は、入力された行列Ｗをコレスキー分解する。
【００８０】
そして、Ｈ計算手段３０５は、コレスキー分解手段３０３から当該Ｈ計算手段３０５に向けて出力された、行列Ｗのコレスキー分解である行列Ｗ1/2を用いて、行列Ｈ＝Ｗ1/2Ｓを、計算する。Ｈ計算手段３０５は、当該行列Ｈの計算結果を、Ｆ計算手段３０６へ送信する。
【００８１】
信号数（到来電波数）Ｎが奇数の場合には、行列Ｈの第１列〜（ｎ＋１）／２列を「Ｈ1」とし、行列Ｈの第（ｎ＋１）／２＋１列〜（ｎ＋１）列を「Ｈ2」として、Ｆ計算手段３０６は、行列Fを下式（１７）のように計算する。
【００８２】
【数２２】
2010121991.tif 000024
【００８３】
上記に対して、信号数（到来電波数）Ｎが偶数の場合には、行列Ｈの第１列〜ｎ／２列を「Ｈ1」とし、行列Ｈの第ｎ／２＋１列を「ｈ」とし、行列Ｈの第（ｎ＋１）／２＋１列〜（ｎ＋１）列を「Ｈ2」として、Ｆ計算手段３０６は、行列Fを下式（１８）のように計算する。
【００８４】
【数２３】
2010121991.tif 000025
【００８５】
ここで、式（１７），（１８）において、「Ｒｅ」は、複素数の実数部を示し、「Ｉｍ」は、複素数の虚数部を示す。また、式（１７），（１８）において、「ｕｄ」は、行列の上下を反転させる操作を意味する。
【００８６】
ＱＲ分解手段３０７は、Ｆ計算手段３０６から出力される式（１７），（１８）の計算結果を受信する。そして、ＱＲ分解手段３０７は、行列Ｆの第１列を「ｆ」、行列Ｆの第２〜ｎ＋１列を「Ｆ2」とおき、当該Ｆ2をＱＲ分解し、行列ＦQおよび行列ＦRを求める。そして、ＱＲ分解手段３０７は、当該行列ＦQおよび行列ＦRを方程式求解手段３０８へ送信する。
【００８７】
方程式求解手段３０８は、受信した行列ＦQおよび行列ＦRを用いて、連立１次方程式ＦQTｆ＝ＦRμを解き、解μ（μ：列ベクトル）を求める。なお、行列ＦRは上三角行列であるから、解μの求解は容易に行うことが可能である。ここで、「ＦQT」は、行列ＦQの転置行列である。
【００８８】
β計算手段３０９は、方程式求解手段３０８から出力された列ベクトルμを受信する。そして、β計算手段３０９は、行ベクトルβ＝［１  μT］を作成し、ｂ計算手段３１０へ出力する。ここで、「μT（μT：行ベクトル）」は、列ベクトルμの転置ベクトルである。
【００８９】
信号数（到来電波数）Ｎが奇数の場合には、ｂ計算手段３１０は、受信した行ベクトルβを用いて、行ベクトルβの第１〜（ｎ＋１）／２の要素を「βr」とし、行ベクトルβの第（ｎ＋１）／２＋１〜ｎ＋１の要素を「βi」として、下式（１９）により列ベクトルｂを計算する。
【００９０】
【数２４】
2010121991.tif 000026
【００９１】
これに対して、信号数Ｎが偶数の場合には、ｂ計算手段３１０は、受信した行ベクトルβを用いて、行ベクトルβの第１〜ｎ／２の要素を「βr」とし、行ベクトルβの第ｎ／２＋１〜ｎの要素を「βi」とし、行ベクトルβの第ｎ＋１の要素を「ν」として、下式（２０）により列ベクトルｂを計算する。
【００９２】
【数２５】
2010121991.tif 000027
【００９３】
ここで、式（１９），（２０）の左辺は、列ベクトルｂの転置ベクトルを示している。また、式（１９），（２０）において、「ｊ」は、虚数単位である。また、式（１９），（２０）において、ｌｒ（βr＋ｊβi）は、行ベクトル（βr＋ｊβi）の左右反転操作（ｆｌｉｐｌｒ）を示す。また、式（１９），（２０）において、（βr＋ｊβi）*は、（βr＋ｊβi）の複素共役であることを示す。
【００９４】
式（１９），（２０）を用いたベクトルｂの計算後、ｂ計算手段３１０は、当該計算結果であるベクトルｂを出力する。
【００９５】
上記までが、２次最適化手段２０２の構成・動作である。２次最適化手段２０７の機能構成は、２次最適化手段２０２の機能構成（図３に示すブロック構成）と同じである。しかしながら、２次最適化手段２０７の動作と２次最適化手段２０２の動作とは、下記の点において相違する。
【００９６】
つまり、２次最適化手段２０７では、ＱＲ分解手段３０７は、信号数Ｎが奇数のときに、式（１７）または式（１８）の行列Ｆの第１〜（ｎ＋１）／２列をＦ1、当該行列Ｆの第（ｎ＋１）／２＋１列をｆ、当該行列Ｆの第（ｎ＋１）／２＋２〜ｎ＋１列をＦ2と、決定する。
【００９７】
これに対して、ＱＲ分解手段３０７は、信号数Ｎが偶数のときに、式（１７）または式（１８）の行列Ｆの第１〜ｎ／２列をＦ1、当該行列Ｆの第ｎ／２＋１列をｆ、当該行列Ｆの第ｎ／２＋２〜ｎ＋１列をＦ2と、決定する。
【００９８】
さらに、当該ＱＲ分解手段３０７は、上記で決定したＦ1，Ｆ2を用いて行列［Ｆ1  Ｆ2］を構成する。そして、当該ＱＲ分解手段３０７は、行列［Ｆ1  Ｆ2］をＱＲ分解し、行列ＦQおよび行列ＦRを求める。
【００９９】
また、２次最適化手段２０７では、β計算手段３０９は、列ベクトルμをＦ1，Ｆ2の列数に合わせて分割し、それぞれμ1，μ2とする。そして、β計算手段３０９は、行ベクトルβ＝［μ1T  １  μ2T］を作成し、ｂ計算手段３１０へ出力する。ここで、μ1およびμ2の上付きの「Ｔ」は、転置操作を示す。
【０１００】
たとえば、行列Ｆ1が３行２列の行列であり、行列Ｆ2が３行１列の行列であるとする。この場合、行列Ｆ1は列数は２であり、行列Ｆ2の列数は１であるので、列ベクトルμは、要素数３（＝２列＋１列）の列ベクトルとなる。そして、前記場合では、行列Ｆ1の列数は「２」なので、当該列数に合わせて、列ベクトルμの第一の要素および第二要素から成る列ベクトルが「μ1」となる。さらに、行列Ｆ2の列数は「１」ので、当該列数に合わせて、列ベクトルμの第三要素が「μ2」となる。なお、この場合、行ベクトルβは、要素数が「４」の行ベクトルとなる。
【０１０１】
なお、上記ＱＲ分解手段３０７およびβ計算手段３０９以外の、２次最適化手段２０７内の各機能構成ブロック３０３〜３０６，３０８，３１０の動作は、２次最適化手段２０２内の各機能構成ブロック３０３〜３０６，３０８，３１０の動作と同じである。
【０１０２】
以上のように、本実施の形態に係る電波到来方向推定装置１７０は、ＭＵＳＩＣ 法と異なり、上記動作を行う固有値抽出部１０７および上記動作を行う最適化部１０８を含んでいる。
【０１０３】
したがって、マルチパス環境下、特に到来距離差の小さなマルチパスが存在する環境下においても、ＭＵＳＩＣ法よりも多い、「アレイアンテナ１０２の素子数−1」個の電波の到来方向を精度良く推定することができる。
【０１０４】
さらに、本実施の形態に係る電波到来方向推定装置１７０は、ＭＯＤＥと異なり、上記動作を行う固有値抽出部１０７を含んでいるので、到来波が相互に強い相関を持つ場合においても、ＭＯＤＥよりも高い電波到来方向の推定精度が保たれる。
【０１０５】
以下において、本実施の形態に係る発明を用いた場合の、測角精度を計算機シミュレーションによって求め、公知手法と比較することにより発明の効果を示す。
【０１０６】
シミュレーション条件は、アレイアンテナ１０２の素子数Ｍ＝４、到来電波数Ｎ＝２、サンプリング周波数ｆ＝２０ＭＨｚ、サンプル数Ｔ＝８５４０とした。２つの到来波はスカラー倍の関係、すなわち信号間の相関は１とした。また、第２の波は第１の波と比較して、−１０ｄＢの強度をもつとした。第１の波はθ1＝０度方向から、第２の波はθ2＝５度方向から到来しているとした。
【０１０７】
図4は、ノイズ強度（横軸：単位ｄＢ）を変化させたときの本発明（○印）、ＭＵＳＩＣ法（＊印）およびＭＯＤＥ（△印）での第1の波の到来方向の推定誤差のＲＭＳＥ（縦軸：単位ｄｅｇ）を、それぞれ図示したものである。また、図４中実線で示したものは、Ｃｒａｍｅｒ−Ｒａｏの下限であり、理論的な推定精度の限界を示す。
【０１０８】
信号間の相関が非常に大きい場合、図４に示すように、ＭＵＳＩＣ法はＣｒａｍｅｒ−Ｒａｏの下限に漸近せず、ＭＯＤＥにおいてもノイズ強度１５ｄＢ〜２５ｄＢにかけてＣｒａｍｅｒ−Ｒａｏの下限と比較して大きく精度が悪化している。
【０１０９】
これに対して、本発明では、図４に示すように、推定精度がＣｒａｍｅｒ−Ｒａｏの下限に漸近している。これは、信号間の相関が非常に大きい場合に、本発明の手法が公知の手法（ＭＵＳＩＣ法やＭＯＤＥ）に比べ、高い推定精度を実現することを示している。
【０１１０】
また、本実施の形態に係る電波到来方向推定装置１７０では、最適化部１０８は、行列Ｂbの要素ｂの収束するまでの計算において、当該要素ｂの偏角φnの値に応じて、第一条件（ｂi＝ｂN-i*  （ｉ＝０，１，・・・，Ｎ）、および、Ｉｍｂ0＝１）および第二の条件（ｂi＝ｂN-i*  （ｉ＝０，１，・・・，Ｎ）、および、Ｒｅｂ0＝１）を切替えて用いている。
【０１１１】
したがって、電波到来方向を推定できない特異な条件を回避し、また前記特異な条件に近い条件で電波到来方向推定精度が悪化することを防ぐことができる。
【０１１２】
＜実施の形態２＞
本実施の形態では、固有値抽出部１０７は、固有値分解部１０６で求められた固有値のうち、最大値から予め規定されている所定の番目までの大きさの値を有する、複数の固有値を、所定の固有値として抽出する。たとえば、上記所定の番目が３番目と予め固有値抽出部１０７に規定されているとする。このとき、固有値抽出部１０７は、固有値分解部１０６で求められた固有値のうち、１番大きな値を有する固有値、２番目に大きな値を有する固有値、および３番目に大きな値を有する固有値を抽出する。
【０１１３】
さらに、本実施の形態では、固有値抽出部１０７は、上記で抽出した複数の固有値を対角に並べた行列Λaを構成する。さらに、当該固有値抽出部１０７は、上記で抽出した複数の固有値に対応する固有ベクトルを列に並べた列ベクトルｓaを構成する。
【０１１４】
そして、本実施の形態では、固有値抽出部１０７は、Λaに含まれない固有値の平均λmを用いて、（Λa−λmＩ）を計算し、当該計算結果をコレスキー分解する。ここで、「Ｉ」は単位行列である。さらに、当該固有値抽出部１０７は、前記コレスキー分解後の行列に対して、左側から上記列ベクトルｓaを乗じたものを、固有空間の基底ｓVとして最適化部１０８へ送信する。当該基底ｓVは、上記で抽出した複数の固有値に対応する固有ベクトルを用いて構成される固有空間の基底であると把握できる。
【０１１５】
たとえば、固有値分解部１０６で、固有値「１．２」、「０．６」、「０．５」および「０．３」が求められたとする。そして、上記所定の番目が３番目であると、予め固有値抽出部１０７に規定されているとする。この場合には、行列Λaは、式（２１）となる。
【０１１６】
【数２６】
2010121991.tif 000028
【０１１７】
ここで、固有値「１．２」に対応する固有ベクトルをｖ1とおき、固有値「０．６」に対応する固有ベクトルをｖ2とおき、固有値「０．５」に対応する固有ベクトルをｖ3とおく。すると、列ベクトルｓaは、［ｖ1  ｖ2  ｖ3］となる。
【０１１８】
また、上記場合では、Λaに含まれない固有値は「０．３」のみである。したがって、λmは、当該固有値「０．３」のみの平均（もちろん、平均は「０．３」である）となる。なお、上記場合では、固有値抽出部１０７で抽出された固有値の数は「３」個であるので、単位行列Ｉは、３行３列の単位行列となる。
【０１１９】
さて、最適化部１０８では、上記固有空間の基底ｓVと、上記で抽出した複数の固有値の個数以上の数の信号成分で構成される信号空間との、直交成分が最小となるという条件を用いて、前記到来角度を推定する。
【０１２０】
具体的に、上記最適化部１０８では、式（６）の条件の下、下式（２２）で表される目的関数（フロベニウスノルム）ＦD（ｂ）を、最小化するｂを求める。
【０１２１】
【数２７】
2010121991.tif 000029
【０１２２】
ここで、式（６）および式（２２）の「ｂ」は、実施の形態１と同様、（Ｍ−Ｎ）×Ｍ行列Ｂbの要素である（式（８）参照）。また、行列Ｂbは、最適化部１０８に予め設定されている。「Ｂb（Ｂb*Ｂb）-1Ｂb*」は、固有値分解部１０６で求められた固有値の個数以上の数の信号成分で構成される信号空間であると把握できる。換言すれば、当該信号空間は、（Ｍ−Ｎ）×Ｍの行列で規定される空間である。なお、式（２２）において、行列ｓV*は、行列（上記基底）ｓVの複素共役転置行列である。
【０１２３】
また、本実施の形態では、Ｓ計算手段３０４は、下式（２３）で定義される基底ｓVを入力し、行列Ｓ1〜ＳLを下式（２４）で定義する。
【０１２４】
【数２８】
2010121991.tif 000030
【０１２５】
【数２９】
2010121991.tif 000031
【０１２６】
ここで、上記「Ｌ」は、固有値抽出部１０７で抽出した固有値の数である。
【０１２７】
また、本実施の形態では、Ｈ計算手段３０５は、コレスキー分解手段３０３から当該Ｈ計算手段３０５に向けて出力された、行列Ｗのコレスキー分解である行列Ｗ1/2を用いて、行列Ｈを下式（２５）を、計算する。Ｈ計算手段３０５は、当該行列Ｈの計算結果を、Ｆ計算手段３０６へ送信する。
【０１２８】
【数３０】
2010121991.tif 000032
【０１２９】
なお、上記以外の構成・動作は、実施の形態１と同様であるので、ここでの説明は省略する。
【０１３０】
本実施の形態に係る電波到来方向推定装置は上記のように構成されているので、固有値抽出部１０７が所定の固有値として複数の固有値を抽出する場合においても、実施の形態１で説明した効果と同様な効果を有することができる。
【０１３１】
＜実施の形態３＞
本実施の形態では、固有値抽出部１０７は、固有値分解部１０６で求められた複数の固有値のうち、所定の固有値として、次の固有値を抽出する。
【０１３２】
固有値抽出部１０７には、しきい値γ（0＜γ＜１）が予め規定・設定されている。当該固有値抽出部１０７は、しきい値固有値（＝γ×（最大固有値））より大きな固有値を、所定の固有値として複数抽出する。
【０１３３】
たとえば、固有値抽出部１０７には、しきい値γとして「０．４」が予め設定されているとする。また、固有値分解部１０６では、固有値「１．５」、「１．２」、「０．９」、「０．７」、「０．５」、「０．３」が求められたとする。この場合、最大固有値は「１．５」であるので、しきい値固有値は、０．４×１．５＝０．６となる。したがって、固有値抽出部１０７は、所定の固有値として、当該「０．６」より大きな固有値である「１．５」、「１．２」、「０．９」、「０．７」を抽出する。
【０１３４】
上記所定の固有値の抽出方法以外は、実施の形態２と同様の構成・動作である。
【０１３５】
本実施の形態に係る電波到来方向推定装置は上記のように構成されているので、固有値抽出部１０７が所定の固有値として複数の固有値を抽出する場合においても、実施の形態１で説明した効果と同様な効果を有することができる。
【０１３６】
なお、実施の形態１乃至３では、固有値抽出部１０７における所定の固有値の抽出方法を各々言及した。しかしながら、当該所定の固有値の抽出方法は、これに限定されるものではなく、他の方法を採用しても良い。
【図面の簡単な説明】
【０１３７】
【図１】本発明に係る電波到来方向推定装置の構成を示すブロック図である。
【図２】最適化部の機能構成を示すブロック図である。
【図３】２次最適化手段の機能構成を示すブロック図である。
【図４】本発明に係る電波到来方向推定装置が奏する効果を説明するための図である。
【符号の説明】
【０１３８】
１０１  到来波（電波）、１０２  アレイアンテナ、１０３  受信機、１０４  Ａ／Ｄ変換器、１０５  共分散行列計算部、１０６  固有値分解部、１０７  固有値抽出部、１０８  最適化部、１０９  表示部、１７０  電波到来方向推定装置、２０２，２０７  ２次最適化手段、２０３，２０８  拘束条件判定手段、２０４，２０９  Ｗ計算手段、２０５，２１０  収束判定手段、２０６，２１１  制約切替器、２１２  到来角計算手段、３０３  コレスキー分解手段、３０４  Ｓ計算手段、３０５  Ｈ計算手段、３０６  Ｆ計算手段、３０７  ＱＲ分解手段、３０８  方程式求解手段、３０９  β計算手段、３１０  ｂ計算手段。
【図１】
2010121991.tif 000043
【図２】
2010121991.tif 000044
【図３】
2010121991.tif 000045
【図４】
2010121991.tif 000046