JP2004112508A - Receiver - Google Patents

Abstract

<P>PROBLEM TO BE SOLVED: To provide a receiver which can perform predetermined signal processing only based on an effective observation signal to improve processing speed and accuracy of the processing results. <P>SOLUTION: A number-of-signal estimating and processing part 20 generates a first covariance matrix from time series data of the observation signal and the number of data by a first eigen value analyzing part 22. The first covariance matrix is analyzed with the eigen value to obtain a first eigen value. A second eigen value analyzing part 23 generates a second covariance matrix from time series data and the number of data of a noise jamming signal without the observation signal. The second covariance matrix is analyzed with the eigen value to obtain a second eigen value. The difference between the first and the second eigen values is obtained by a subtractor 24. A number-of-signal determination part 25 outputs selection permit signals by the number of the difference larger than the reference value. A target signal selection part 26 selects the observation signals by the number of the signals according to the selection permit signals. <P>COPYRIGHT: (C)2004,JPO

Description

とかける。
【００７１】
結合分布と周辺分布の距離Ｋｕｌｌｂａｃｋ−Ｌｅｉｂｌｅｒ　Ｄｉｖｅｒｇｅｎｃｅ（Ｋｕｌｌｂａｃｋ−Ｌｅｉｂｌｅｒ情報量）は、
【数１２】

となる。但し、Ｈ（Ｙ；Ｂ）は同時分布ｐ（ｙ）のエントロピー、Ｈ（Ｙ_ｉ；Ｂ）は周辺分布ｐ（ｙ_ｉ）のエントロピーである。
【００７２】
これは｛Ｙ_ｉ｝（ｉ＝１，・・・，ｎ）の相互情報量である。信号源が正規分布でないという仮定からＫＬ（Ｂ）はｐ（ｙ_ｉ）が互いに独立な場合に限り「０」となる。これらはｐ（ｘ）とＢによって定まるものである。
【００７３】
ここで、ｐ（ｙ）ｄｙ＝ｐ（ｘ）ｄｘ、ｐ（ｙ）＝ｐ（ｘ）／｜Ｂ｜（｜Ｂ｜はＢの行列式）であることに注意すると、
【数１３】

となる。
【００７４】
一方、周辺分布のエントロピーは、
【数１４】

である。よって、
【数１５】

となり、
【数１６】

のようにすれば最急降下法として正しいＢを求めることができる。上記式（１６）の中で問題となるのは逆行列（Ｂ^Ｔ）^−１を計算している点である。
【００７５】
収束性に関しては、これにいかなる正定値行列を掛けても構わないことから、Ｂ^ＴＢを掛ければ（これは正則な行列の多様体上でのリーマン（Ｒｉｅｍａｎ）計量に対応している）、
【数１７】
ΔＢ∝（Ｉ−Ｅ_ｘ［ψ（ｙ）ｙ^Ｔ］Ｂ　　…式（１７）
が新たな学習則となる。定常性の仮定よりｐ（ｓ）、ｐ（ｘ）とｐ（ｙ）は時間的に独立である。この仮定のもと、アンサンブル平均を時間平均に置き換えることができる。
【００７６】
【数１８】

したがって、ηを正の定数とし、データが観測される毎に下記式（１９）に従ってパラメータを更新すればＢ_ｔが得られる。
【００７７】
【数１９】
Ｂ_ｔ＋１＝Ｂ_ｔ＋η（Ｉ−ψ（ｙ）ｙ^Ｔ）Ｂ_ｔ　　…式（１９）
ここで、当然問題になるのは、上記式（１５）のｐ（ｙ_ｉ）或いはψ（ｙ）をいかに定義するかである。通常、これはパラメトリックな非線形関数や統計的な展開法が用いられる。
【００７８】
大雑把な考え方を示すと、もしｐ（ｙ_ｉ）が正規分布ならはψ（ｙ）は線形関数となる。一方、正規分布より裾が”重い”場合（ｓｕｂ−Ｇａｕｓｓｉａｎ）多項式などで近似するのがよく、正規分布より裾が”軽い”場合（ｓｕｐｅｒ−Ｇａｕｓｓｉａｎ）シグモイド（Ｓｉｇｍｏｉｄ）関数などで近似するのがよいとされている。音声信号などは裾が”軽い”ので、シグモイド関数などがうまく働く。
【００７９】
ブラインド処理部３０で分離・抽出された信号Ｙ_{１，・・・，４}（ｔ）は分離信号出力処理部４０に送られる。分離信号出力処理部４０では、ブラインド処理部３０からのデジタル信号として送られてくる信号Ｙ_{１，・・・，４}（ｔ）をアナログ信号に変換し、分離信号Ｏ_１〜Ｏ_４として外部にそれぞれ送出する。
【００８０】
なお、上記のＩＣＡの解法は、一例を示すものであり、これに限るものではない。
【００８１】
一方、ブラインド処理部３０で求められた行列Ａは、到来方位演算部５０に送られる。到来方位演算部５０は、受け取った行列Ａに基づいて、以下に示すセンサ間位相差演算を行う。
【００８２】
すなわち、先ず、行列Ａから対応する信号に対するベクトルを抽出する。ここで、ベクトルを下記式（２０）のように表す。
【００８３】
【数２０】
ｓｉｇ１＝［ａ１１＋ｊｂ１１，ａ１２＋ｊｂ１２，ａ１３＋ｊｂ１３，ａ１４＋ｊｂ１４］…式（２０）
ｓｉｇ２＝［ａ２１＋ｊｂ２１，ａ２２＋ｊｂ２２，ａ２３＋ｊｂ２３，ａ２４＋ｊｂ２４］
ｓｉｇ３＝［ａ３１＋ｊｂ３１，ａ３２＋ｊｂ３２，ａ３３＋ｊｂ３３，ａ３４＋ｊｂ３４］
ｓｉｇ４＝［ａ４１＋ｊｂ４１，ａ４２＋ｊｂ４２，ａ４３＋ｊｂ４３，ａ４４＋ｊｂ４４］
上記式（２０）において、Ｓｉｇ１は、入射信号Ｓ_１のみが存在する時の状態を表しており、複素ベクトル「ａ１１＋ｊｂ１１」は第１アンテナ１１で観測された信号の振幅と位相を、「ａ１２＋ｊｂ１２」は第２アンテナ１２で観測された信号の振幅と位相を、「ａ１３＋ｊｂ１３」は第３アンテナ１３で観測された信号の振幅と位相を、「ａ１４＋ｊｂ１４」は第４アンテナ１４で観測された信号の振幅と位相をそれぞれ表す。Ｓｉｇ２、Ｓｉｇ３及びＳｉｇ４についても同様である。
【００８４】
この式（２０）から、行列Ａは、各アンテナ１〜７で観測された信号の位相を反映しており、この位相の分布から入射信号Ｓ_１〜Ｓ_７の到来方位を推定できることが理解できる。後述するインタフェロメータ方式は、この位相の分布から信号の到来方位を求めるものである。
【００８５】
このベクトルからセンサ間（ビーム間）の位相差を求める。上記式（２０）からは、例えば「位相角（ａ１１＋ｊｂ１１）−位相角（ａ１２＋ｊｂ１２）」から第１アンテナ１と第２アンテナ２の位相差を求めることができ、「位相角（ａ１３＋ｊｂ１３）−位相角（ａ１４＋ｊｂ１４）」から第３アンテナ３と第４アンテナ４の位相差を求めることができる。
【００８６】
次に、到来方位演算部５０は、方位演算を行う。すなわち、上記で求められた位相差の何れかとアンテナの物理的位置関係から通常のインタフェロメータ方式の手法により入射信号Ｓ_１の到来方位を求める。
【００８７】
以下、入射信号Ｓ_１の到来方位の求め方を詳細に説明する。今、図４に示すように、第１アンテナ１と第２アンテナ２との距離をｄとし、第１アンテナ１と第２アンテナ２とを結ぶ線に直交する線に角度θで入射信号Ｓ_１が入射するものとする。この場合、入射信号Ｓ_１が第１アンテナ１に到着してから第２アンテナ２に至るまでの距離ｙは「ｄｓｉｎθ」である。これを位相に直すと、第１アンテナ１と第２アンテナ２の位相差φは、
【数２１】
φ＝（２π／λ）・ｄｓｉｎθ・・・式（２１）
で表される。ここで、λは波長である。この位相差φが観測されるデータである。
【００８８】
従って、上記式（２１）からθを求めることにより、入射信号Ｓ_１の到来方位を求めることができる。すなわち、
【数２２】
φ＝（２π／λ）・ｄｓｉｎθ・・・式（２２）
よって、θは、下記式（２３）から求めることができる。
【００８９】
【数２３】
θ＝ｓｉｎ^−１（φ・１／（（２π／λ）・ｄ））・・・式（２３）
以下、順次、入射信号Ｓ_２、入射信号Ｓ_３及び入射信号Ｓ_４に対応する到来方位を求める。このようにして求められた到来方位を示す信号Ｚ_{１，・・・，４}（ｔ）は、到来方位出力処理部６０に送られる。
【００９０】
到来方位出力処理部６０では、到来方位演算部５０からデジタル信号として送られてくる信号Ｚ_{１，・・・，４}（ｔ）をアナログ信号に変換し、方位信号Ｄ_２〜Ｏ_４として外部にそれぞれ送出する。
【００９１】
これにより、分離された４つの信号が、到来方位を表す信号とともに、時間軸波形として出力される。このようにして、観測信号Ｘ_{１，・・・，４}（ｔ）の所定の区間のサンプリングデータに対する処理が終了すると、次の区間のサンプリングデータに処理が実行され、以下同様の処理が繰り返される。
【００９２】
なお、第１の実施の形態では、センサとして、アンテナを用いた例を説明したが本発明はこれに限定されるものではない。例えば、センサとして、音声を検知するマイクロフォン、生体の種々の状態を検知するセンサ等を用いることができる。
【００９３】
以上詳述したように、この発明の第１の実施の形態に係る受信装置によれば、複数の入射信号が同時に入射される複数のセンサからの観測信号に基づいて、信号数推定処理部２０により入射信号の信号数を推定して複数の観測信号の中から当該信号数分の観測信号を選択して出力するので、有効な観測信号のみに基づいて所定の信号処理を行うことができ、受信装置の処理速度と処理結果の精度を向上することができる。
【００９４】
また、信号数推定処理部２０では、第１固有値解析部２２により観測信号の時系列データとデータ数から第１の共分散行列を生成し、この第１の共分散行列を固有値解析して第１の固有値を求め、第２固有値解析部２３により観測信号がない状態での雑音妨害信号の時系列データとデータ数から第２の共分散行列を生成し、この第２の共分散行列を固有値解析して第２の固有値を求め、減算器２４−１〜２４−７では求まった第１および第２の固有値からその差を求め、信号数決定部２５では求められた差が基準値よりも大きくなる信号数分の選択許可信号を出力し、対象信号選択部２６によりこの選択許可信号に応じて複数の観測信号の中から信号数分の観測信号を選択して出力するので、観測信号と雑音との境界を明確に区別して有効な観測信号のみに基づいて所定の信号処理を行うことができ、受信装置の処理速度と処理結果の精度を向上することができる。
【００９５】
さらに、ブラインド処理部３０（ベクトル演算部）では、選択された信号数分の観測信号に基づいて、複数の入射信号をブラインド処理により分離し、複数の入射信号の各々について、複数のセンサにおける振幅及び位相分布の相対関係を示す複素ベクトル量を求め、到来方位演算部５０により複数のセンサの物理的位置と複素ベクトル量とに基づいて、前記複数の入射信号の到来方位を求め、分離された信号と到来方位を表す信号とを関連付けて出力することで、入射信号の到来方位を正確に求めることができる。
【００９６】
また、信号分離部７０では、信号数推定処理部２０により選択された信号数分の観測信号に基づいて固有値解析を行った際に、同時に求まる固有ベクトルを信号固有ベクトルと雑音固有ベクトルに分離し、この信号固有ベクトルにより捕捉信号の信号成分のみを写像し、写像したデータの独立成分分析により同時に存在する複数の信号を分離することで、入射信号を正確に分離することができる。
【００９７】
また、到来方位演算部８０では、信号数推定処理部２０により選択された信号数分の観測信号に基づいて固有値解析を行った際に、同時に求まる固有ベクトルから信号固有ベクトルと雑音固有ベクトルに分離し、この固有ベクトルと前記センサの応答データとの演算により信号の到来方位を探知することで、入射信号の到来方位を正確に求めることができる。
【００９８】
（第２の実施の形態）
図５は、本発明の第２の実施の形態に係る受信装置の構成を概略的に示すブロック図である。
【００９９】
本実施の形態における特徴は、信号分離部７０のアルゴリズムとして、確率分布の独立性に基づく分離法を採用している。
【０１００】
この分離手法では、観測信号の各成分ｓ_ｉ（ｔ）が正規分布ではない確率分布に従っていると仮定する。ただし、その確率分布は未知である。
【０１０１】
今、ｙ（ｔ）の同時確率密度変数を、
【数２４】
ｙ＝（ｙ_１，…，ｙ_ｎ）^Ｔ　　　…式（２４）
と置く、行列Ｗによって観測信号ｘ（ｔ）が正しく分離できたならば、ｙの各成分ｙ_ｉは独立になる。ｙの各成分ｙ_ｉが独立と言うことは、ｙ_ｉの値がわかってもｙ_ｉ以外のｙ_ｉ，…，ｙ_ｎの値が何かわからないことであり、数式で表すと、
【数２５】

同時密度関数　＝　周辺密度関数の積
となる。
【０１０２】
ここで、基準の１つとして復元信号から推定される結合分布と周辺分布の距離が考えられる。分布間の統計的距離には様々なものが考えられるが、独立成分分析においては、Ｋｕｌｌｂａｃｋ−Ｌｅｉｂｌｅｒ　ｄｉｖｅｒｇｅｎｃｅが良く用いられる。
【０１０３】
Ｋｕｌｌｂａｃｋ−Ｌｅｉｂｌｅｒ　ｄｉｖｅｒｇｅｎｃｅの定義は、次のとおりである。
【０１０４】
【数２６】

ここで、ｚの分布をｘの分布で置き換えたときの増加情報量を表しており、値が大きい場合は、分布が似ていないことを表し、値が小さい場合は、分布が似ていることを表す。
【０１０５】
そこで、　　（１２）式　　　より、
【数２７】

として、　　（１３）式　　　に代入すると、
【数２８】

となり、ｙ_ｉが独立であれば（９）式はゼロとなり、すなわちｙ（ｔ）が独立の場合に限り最小となる。行列Ｗを求めるには、Ｉ（ｙ｜Ｗ）のＷに関する勾配を求め、最急降下法により求まる。
【０１０６】
【数２９】

ただし、
【数３０】

と更新していくことでＷを求められる。計算上（１１）式中の逆行列（Ｗ^Ｔ）^−１が問題となる。収束性に関しては正定値行列を掛けてもかまわないことから（Ｗ^Ｔ）^Ｗを掛けて、
【数３１】

とした方が、計算量も少なく収束性が速い。
【０１０７】
この更新則はψ（ｙ）を含んでおり、密度関数の形がわからなければ計算できないが、適当なパラメトリックな非線形関数や統計的な展開法によって近似しても解は正しく求まる。
【０１０８】
（変形例１）
変形例１では、図５に示す信号分離部７０のアルゴリズムとして、勾配法を採用している。
【０１０９】
勾配法は、（９）式に示した相互情報量を最小とするＷを最急降下法で求める手法であり、更新則（１７）式により求める手法である。
【０１１０】
具体的には、不動点法がある。
【０１１１】
不動点法では、ベクトルｗを引数とする関数ψ（ｗ）をノルムが‖ｗ‖という条件下の基で最大あるいは最小化する条件付最適化問題を考える。
【０１１２】
この問題は、Ｌａｇｒａｎｇｅの未定係数法を用いて、
【数３２】
ψ（ｗ）＋λ（‖ｗ‖−１）→ｍｉｎ　　…式（３２）
と書き換えられる。最小ではこの微分が０となることから、
【数３３】

これを第２項について解くと、
【数３４】

となる。ｗは、‖ｗ‖＝１で正規化されているので、λを求めるのには正規化によればよい。実際の計算は、
【数３５】

という計算を収束するまで繰り返すこととなる。
【０１１３】
（変形例２）
変形例２では、図５に示す信号分離部７０のアルゴリズムとして、ＦａｓｔＩＣＡ法を採用している。
【０１１４】
Ｈｙｖａｒｉｎｅｎらは、不動点法の収束速度の条件に関して論議しており、特に尖度（ｋｕｒｔｏｓｉｓ）を用いた彼らの方法（ＦａｓｔＩＣＡ）は、３次収束するとしている。
【０１１５】
φ（ｗ）として、
【数３６】
φ（ｗ）＝Ｅ（（ｗ・ｘ）^４）−３Ｅ（（ｗ・ｘ）^２）^２　　　…式（３６）
と言う評価関数を使い、これの最大値や最小値を不動点法で求めて、
によって独立な成分を抽出するアルゴリズムである。
【０１１６】
今、原信号の平均は０、分散は１、混合行列Ａは直交行列にして、
【数３７】
ｖ（ｔ）Ｔ＝ｗ（ｔ）ＴＡ　　‖ｖ‖＝１　　…式（３７）
というものの振る舞いを考える。
【０１１７】
（Ａを直交行列に限るのは、ｗｈｉｔｅｎｉｎｇを行い、一旦信号を無相関にする。）
この書き換えにより、ｗが１つの成分を捉えているのならば、ｖはどれかひとつの成分が１で、あとは０のベクトルに近づいていく。
【０１１８】
評価関数を原信号ｓとｖを使って書くと、
【数３８】
Ｅ（（ｖ・ｓ）^４）−３Ｅ（（ｖ・ｓ）^２）^２　＝Ｅ（（ｖ・ｓ））^４−３‖ｖ‖^４　　式（３８）
となる。
【０１１９】
ｓの成分は平均０、分散１であること、および互いに独立であることを用い第ｉ成分を考えると、
【数３９】

となる。Ｅ（ｓ_ｉ ^４）＝１であるから最後の括弧の中は４次のキュムラントになっている。これを正規化したものが新しいベクトルになるが、比を考えると、
【数４０】

と評価され、ｖ_ｉが１近くてｖ_ｊが０に近いと仮定すれば、上式はｖ_ｊへの０への近づく速さがだいたい３乗のオーダーであることを意味している。
【０１２０】
（変形例３）
変形例３では、図５に示す信号分離部７０のアルゴリズムとして、相関関数行列の同時対角化法を採用する。
【０１２１】
Ｊａｃｏｂｉ法は、実対称行列を対角化する。すなわち、固有値を求める古典的な方法である。その基本方針は、Ｇｉｖｅｎｓ変換による相似変換を繰り返すことにより直交変換を求める方法である。
【０１２２】
提案されている方法としては、相関関数行列と縮約された尖度（ｋｕｒｔｏｓｉｓ）の同時対角化算法がこれにあたる。
【０１２３】
ここで、相関関数行列の同時対角化法として、時間構造に基づく分離法について説明する。
【０１２４】
相関関数行列による算法は、時間差のある相関信号を考える。観測信号が、ｘ（ｔ）＝Ａｓ（ｔ）と書けるとき、原信号の独立性を仮定すればその相関関数行列は、
【数４１】

となる。ただし、ここでは複素数値時系列を考え＊は複素共役転置を、＜・＞は平均操作を、またＲ_ｓｉ（τ）は原信号ｓ_ｉ（ｔ）の自己相関関数を表すものとする。
【０１２５】
これも対角行列を同じ行列で相似変換した対称行列だから、いろいろな時間差における行列を同時に対角化してやるＷを探せばよい。
【０１２６】
具体的には、観測信号ｘ（ｔ）をτ＝０における相関行列を使って無相関化しておき、τ≠０における相関行列の非対角化成分を小さくするように、Ｊａｃｏｂｉ法により算出する。
【０１２７】
（変形例４）
変形例４は、図５に示す信号分離部７０のアルゴリズムとして、Ｋｕｒｔｏｓｉｓの同時対角化による算法であるＪＡＤＥ（Ｊｏｉｎｔ　Ａｐｐｒｏｘｉｍａｔｅ　Ｄｉａｇｏｎａｌｉｚａｔｉｏｎ　ｏｆ　Ｅｉｇｅｎ−ｍａｔｒｉｃｅｓ）法を採用する。
【０１２８】
縮約された尖度（ｋｕｒｔｏｓｉｓ）の同時対角化算法は、４次のｃｒｏｓｓ−ｋｕｒｔｏｓｉｓとして、
【数４２】

の縮約を考える。
【０１２９】
観測信号が平均０でｗｈｉｔｅｎｉｎｇにより無相関化されているとすると、原信号と観測信号はある直交変換Ｕ＝（ｕ_{１， … ，ｎ}）により、Ｘ＝ＵＳという関係で結ばれている。独立性の仮定から、
【数４３】

として、行列Ｍ＝（ｍ_ｉｊ）により縮約されたｋｕｒｔｏｓｉｓの行列として、
【数４４】

を考えると、
【数４５】

さらに、
【数４６】

として、
【数４７】

となる。従って、適当な組の｛Ｎ_ｒ｝を取って｛Ｃ（Ｎ_ｒ）｝を計算し、これを同時対角化する直交行列を求めればよいことがわかる。この事情は
行列でも同様である。
【０１３０】
ここで、Ｃａｒｄｏｓｏ　ａｎｄ　Ａ．　Ｓｏｕｌｏｕｍｉａｃは｛Ｎ_ｒ｝として、
【数４８】

ただし、ｅ_ｋはｋ成分だけが１の単位ベクトル、あるいは尖度（ｋｕｒｔｏｓｉｓ）そのものを固有値行列展開しておき、固有値の大きいものから適当な数だけ固有行列を使う方法を提案している。
【０１３１】
（変形例５）
変形例５は、図５に示す信号分離部７０のアルゴリズムとして、ノイズを含むモデルを採用する。
【０１３２】
ＩＣＡのモデルとして（３）式を示したが、実際の測定において外来雑音が無視できない場合がある。外来雑音が大きい場合、信号の分離・復元に影響を及ぼす。ここで、観測雑音を、
【数４９】
ε（ｔ）＝［Ｅ_１（ｔ），…，Ｅ_ｎ（ｔ）］^Ｔ　　　…式（４９）
とすると、センサの観測信号は、
【数５０】
ｘ（ｔ）＝Ａｓ（ｔ）＋ε（ｔ）　　　…式（５０）
となる。計算をする場合、アルゴリズムは信号とノイズを区別しないため、
【数５１】

のように書き直せば、信号の次元をｍとするとＡ’は、ｎ×（ｍ＋ｎ）の行列となり、信号源の数が観測信号の次元より多くなり線形な行列の解は存在しなくなる。また、ｗｈｉｔｅｎｉｎｇによる無相関化もノイズを含めた信号を無相関とするため意味がなくなってしまう。
【０１３３】
変形例５は、図５に示す信号分離部７０のアルゴリズムとして、ノイズを低減する手法としてＩＣＡの前処理として因子分析法を採用する。
【０１３４】
この因子分析法では、因子分析混合・復元モデルを採用する。
【０１３５】
因子分析法は、観測ベクトルを要素間にまたがって存在する共通因子と要素間に独立因子の２つの部分に分けて考えるものである。
【０１３６】
すなわち、観測値ｘ∈Ｒ_ｎを、
【数５２】
ｘ＝ａ_１ｓ_１＋ａ_２ｓ_２＋…ａ_ｍｓ_ｍ＋ξ＝Ａｓ＋ξ　　…式（５２）
のように共通因子ｓ∈Ｒ^ｍと独立因子ξ∈Ｒ^ｎを用いて分解できるとする。ここで、Ａは因子負荷量と呼ばれ、ｓ，ξは互いに独立で、ξは各成分が独立とする。
【０１３７】
ここで、原信号ｓの共分散行列がｍ次元の単位行列Ｉ_ｍ、ξが正規分布に従うとして、共分散がｎ次元の対角行列Ψと仮定し、
【数５３】
Σ＝Ｅ（ｘｘ^Ｔ）＝ＡＡ^Ｔ＋Ψ　　　…式（５３）
とし、観測値から求められる共分散は、
【数５４】

に基づいて推定を行う。もっとも単純な方法が最小二乗推定であり、Σとｒとの誤差を最小化するものである。
【０１３８】
すなわち、損失関数Ｌ（Ａ，Ψ）をｒ＝ＡＡ^Ｔ＋Ψにより、
【数５５】
Ｌ（Ａ，Ψ）＝ｔｒ（ｒ−ＡＡ^Ｔ−Ψ）^Ｔ（ｒ−ＡＡ^Ｔ−Ψ）　　…式（５５）
として最小化する。（４２）式をＡ，Ψにより偏微分することで最小二乗法によるＡ，Ψは、ＡとΨに関する次の解が得られる。
【０１３９】
【数５６】
（ｒ−Ψ）Ａ＝Ａ（Ａ^ＴＡ）　　　…式（５６）
【数５７】
Ψ＝ｄｉａｇ（ｒ−ＡＡ^Ｔ）　　　　…式（５７）
次に、自由度を考えると、Σの自由度はｎ（ｎ＋１）／２である。一方、ＡとΣのパラメータの数は（ｍ＋１）ｎ。ただし、直交行列分の任意性を差し引かないといけないのでｍ（ｍ−１）／２を引くと、この自由パラメータがΣの自由度より少ないということから、
【数５８】

ここで、ｍ＜ｎと言う条件を加えてこれを解くと、
【数５９】

となる。これが解が存在するための必要十分条件である。
【０１４０】
以上よりノイズを含むモデルについてセンサによる観測信号ｘと因子分析後のセンサによる信号χ＾との関係は、一般化逆行列として、
【数６０】
Ｑ＝［ＡΨ^−１Ａ^ＴΨ^−１］^−１　　　…式（６０）
を用いて
【数６１】

より得られる。
【０１４１】
（変形例６）
変形例６は、図５に示す信号分離部７０のアルゴリズムとして、情報量最大化（Ｉｎｆｏｍａｘ）する手法を採用する。
【０１４２】
今、観測信号の確率変数をＸとすると通信路Ｇに通して確率変数Ｙを得たとする。この間の関係を、
【数６２】
Ｙ＝Ｇ（Ｘ）＋Ｎ　　　…式（６２）
Ｎ：雑音
とし、このとき出力Ｙが持っている入力Ｘに関する情報量ができる限り大きくなるように通信路を設計するというのが情報化最大量（Ｉｎｆｏｍａｘ）の考えかたの基本となる。
【０１４３】
具体的には、ＸとＹとの間の相互情報量をある制約のもとで最大化する問題ということになる。
【０１４４】
独立成分分析の場合は線形演算と非線形変換の組み合わせで、
【数６３】

と考え、その相互情報量は、
【数６４】

となる。
【０１４５】
この場合ＸとＹの間には雑音成分の曖昧さしかないから、条件つきエントロピーは雑音のエントロピーに等しくなる。情報化最大量（Ｉｎｆｏｍａｘ）の考え方では、この部分が大切な役割を果たし、結局ここから先は雑音を考えないで討議することになる。
【０１４６】
最急降下法によって相互情報量の微分はＹのエントロピーの微分は、
【数６５】

となる。この場合非線形変換の取り扱いに注意して、ＸとＹの分布の関係はｙ＝Ｇ（Ｗｘ）の対応関係が一対一だとすれば、
【数６６】

と書ける。ただし、Ｊはこの対応関係のＪａｃｏｂｉａｎで、
【数６７】

で定義され、具体的に計算するとＧの微分をＧ’として、
【数６８】

【数６９】

となる。ＹのエントロピーをＸの分布で表すと、
【数７０】

となる。これは、
【数７１】

となることを使って、
【数７２】

と書けば、［１］で示したようにψ（ｙ_ｉ）を更新していくことで、Ｗを求めることができる。
【０１４７】
（変形例７）
変形例７は、図５に示す到来方位演算部８０のアルゴリズムとして、ＭＵＳＩＣ　（ＭＵｌｔｉｐｌｅ　ＳＩｇｎａｌ　Ｃｌａｓｓｉｆｉｃａｔｉｏｎ）法を採用する。
【０１４８】
ＭＵＳＩＣ法は、Ｓｕｐｅｒ　Ｒｅｓｏｌｕｔｉｏｎアルゴリズムの中でも、最も代表的なものである。共分散行列を固有値分解して信号と雑音の部分空間に分け、信号が雑音の部分空間に直交する性質を用いて到来方向を求める方法で、Ｎｏｉｓｅ　Ｓｕｂｓｐａｃｅ　Ｆｉｔｔｉｎｇ　（ＮＳＦ）の一種である。空間の直交性を利用しているため、１回の方向サーチで複数の到来方向を求めることができる。しかしながら、相関性のある入射信号は分離できない。
【０１４９】
以下、ＭＵＳＩＣ法に関するアルゴリズムについて説明する。
【０１５０】
（１）　Ｒを固有値解析し、観測信号と雑音の部分空間に分解する。なお、上述した到来波数推定が必要になる。
【０１５１】
（２）　サーチする方向θのＭＵＳＩＣスペクトラムは、
【数７３】

となる。観測信号と雑音の部分空間は互いに直交する補空間であることから、Ｐ_ＭＵは入射信号の到来方向に鋭いピークを持つ。方向サーチして、大きい順にＬ個のピークを抽出する。
【０１５２】
（変形例８）
変形例８は、図５に示す到来方位演算部８０のアルゴリズムとして、ｃｕｍｕｌａｎｔ−ＭＵＳＩＣ法を採用する。
【０１５３】
その特徴は、アルゴリズムはＭＵＳＩＣ法と同じであるが、共分散行列の計算に４次のｃｕｍｕｌａｎｔを導入している。４次のｃｕｍｕｌａｎｔはガウス雑音を抑圧すると言われているため、ＭＵＳＩＣ法よりも雑音に強いことが期待される。
【０１５４】
アルゴリズムは、基本的にＭＵＳＩＣ法と同じであり、２次の共分散行列Ｒの代わりに４次のｃｕｍｕｌａｎｔによる行列Ｃを計算する。ここで、行列Ｃは、
【数７４】
Ｃ　＝　｛ｃ_ｉｊ｝　（ｉ，ｊ＝１．．Ｋ）　　　　　　…式（７４）
ｃ_ｉｊ　＝　Σ　ｃｕｍ（Ｘ_ｉ ^＊，　Ｘ_ｊ，　Ｘ_ｋ ^＊，　Ｘ_ｋ）
となる。
【０１５５】
（変形例９）
変形例９は、図５に示す到来方位演算部８０のアルゴリズムとして、ＷＳＦ　（Ｗｅｉｇｈｔｅｄ　Ｓｕｂｓｐａｃｅ　Ｆｉｔｔｉｎｇ）法を採用する。
【０１５６】
その特徴は、ＭＵＳＩＣ法が雑音の部分空間を用いて到来方向を求めるのに対して、ＷＳＦ法では信号の部分空間を用いて到来方向を求めており、Ｓｉｇｎａｌ　Ｓｕｂｓｐａｃｅ　Ｆｉｔｔｉｎｇ　（ＳＳＦ）　の一種である。
【０１５７】
相関性のある観測信号も分離可能である。Ｌ波の到来方向を多次元サーチしなければならず、計算量が多くなるため、予めＭＵＳＩＣ等で初期値を決めておき、サーチする範囲を限定するのが良い。
【０１５８】
アルゴリズムは、（１）　Ｒを固有値解析し、観測信号と雑音の部分空間に分解する。なお、前もって到来波数の推定を行う必要がある。（２）　以下のＶ_ＷＳＦを最小にするθ＝　［θ_１，…，　θ_Ｌ　］をサーチして求める。
【０１５９】
【数７５】
Ｖ_ＷＳＦ　（θ）　＝　Ｔｒ　（　Ｐ_Ａ ^⊥Ｅ_ｓ　Ｗ_ｓ　Ｅ_ｓ ^Ｈ　）　　　…式（７５）
Ｗ_ｓ　＝　（Λ_ｓ　−　σ^２Ｉ　）^２Λ_ｓ ^−１
（変形例１０）
変形例１０は、図５に示す到来方位演算部８０のアルゴリズムとして、ＵＳＦ　（Ｕｎｗｅｉｇｈｔｅｄ　Ｓｕｂｓｐａｃｅ　Ｆｉｔｔｉｎｇ）　法を採用する。
【０１６０】
ＵＳＦ法の特徴は、ＷＳＦ法と同じくＳＳＦ法の一種であり、ＷＳＦ法のウエイトＷ_ｓを単位行列Ｉに置き換えたものであり、（ｕｎｗｅｉｇｈｔｅｄ）　ＭＤ−ＭＵＳＩＣ　（Ｍｕｌｔｉｄｉｍｅｎｔｉｏｎａｌ　ＭＵＳＩＣ）法も呼ばれる。
【０１６１】
Ｌ波の到来方向の多次元サーチが必要で、予めＭＵＳＩＣ法で初期値を決めておき、サーチする範囲を限定するのが良く、繰り返し計算により収束させるアルゴリズムもよい。
【０１６２】
アルゴリズムは、（１）　Ｒを固有値解析し、観測信号と雑音の部分空間に分解する。この時、到来波数を推定しておく必要がある。（２）　以下のＶ_ＵＳＦを最小にするθ＝　［θ_１，…，　θ_Ｌ　］をサーチ、または収束アルゴリズムにより求める。
【０１６３】
【数７６】
Ｖ_ＵＳＦ　（θ）　＝　Ｔｒ　（　Ｐ_Ａ ^⊥Ｅ_ｓ　Ｅ_ｓ ^Ｈ　）　　　　…式（７６）
（変形例１１）
変形例１１は、図５に示す到来方位演算部８０のアルゴリズムとして、ＭＬ　（Ｍａｘｉｍｕｍ　Ｌｉｋｅｌｉｈｏｏｄ）　法を採用する。
【０１６４】
ＭＬ法の特徴は、観測量から隠れた未知のパラメータを推定する一般的な方法である。観測された量の確立密度関数（ＰＤＦ）の積で表されるｌｉｋｅｌｉｈｏｏｄ　ｆｕｎｃｔｉｏｎが最大になるようにパラメータを決定する。
【０１６５】
ＭＬ法には　ｕｎｃｏｎｄｉｔｉｏｎａｌ　ＭＬ　（ＵＭＬ）（またはｓｔｏｃｈａｓｔｉｃ　ＭＬ）と　ｃｏｎｄｉｔｉｏｎａｌ　ＭＬ　（ＣＭＬ）　（またはｄｅｔｅｒｍｉｎｉｓｔｉｃ　ＭＬ）の２つがあるが、いずれもＬ波の到来方向の多次元サーチが必要である。予めＭＵＳＩＣ法等で初期値を決めておき、サーチする範囲を限定するのが良く、相関性のある観測信号にも有効である。
【０１６６】
アルゴリズムは、（１）　到来波数推定アルゴリズムにより到来波数を推定する。
（２）　以下のＶ_ＵＭＬ　又はＶ_ＣＭＬを最小にするθ＝　［θ_１，…，　θ_Ｌ　］をサーチ、または収束アルゴリズムにより求める。
【０１６７】
【数７７】
Ｖ_ＵＭＬ　（θ）　＝　ｌｏｇ　（　ｄｅｔ　［Ｐ_Ａ　ＲＰ_Ａ　−　Ｔｒ（Ｐ_Ａ ^⊥）　Ｐ_Ａ ^⊥　／　（Ｋ−Ｌ）］　）　　…式（７７）
Ｖ_ＣＭＬ　（θ）　＝　Ｔｒ（Ｐ_Ａ ^⊥Ｒ）
（変形例１２）
変形例１２は、図５に示す到来方位演算部８０のアルゴリズムとして、最小分散法を採用する。
【０１６８】
その特徴は、Ｃａｐｏｎ法とも呼ばれ、メインローブをある到来波方向に向け、他の到来波方向にヌルを向けており、ＭＵＳＩＣ法に比べて分解能は悪く、相関性のある観測信号は分離できない。
【０１６９】
アルゴリズムは、（１）　サーチする方向θのスペクトラムＰ_ＣＰを計算する。
【０１７０】
【数７８】

方向サーチして、大きい順にＬ個のピークを抽出する。なお、アルゴリズムとして直接的には到来波数の推定を行う必要はないが、Ｌ個のピークをサーチするために必要となる。
【０１７１】
（変形例１３）
変形例１３は、図５に示す到来方位演算部８０のアルゴリズムとして、ＢＦ法　（ｂｅａｍｆｏｒｍｅｒ）を採用する。
【０１７２】
その特徴は、メインローブを走査してアレイの出力電力が大きくなる方向をサーチする方法である。但し、角度分解能が悪く、複数の到来波がある場合分離が困難で方向も不正確になる。
【０１７３】
アルゴリズムは、（１）　サーチする方向θのスペクトラムＰ_ＢＦを計算する。
【０１７４】
【数７９】

方向サーチして、大きい順にピークを抽出する。
【０１７５】
（変形例１４）
変形例１４は、図５に示す到来方位演算部８０のアルゴリズムとして、ＤＯＳＥ　（Ｄｉｒｅｃｔｉｏｎ　ｏｆ　Ｓｏｕｒｃｅ　Ｅｌｉｍｉｎａｔｉｏｎ）法　を採用する。
【０１７６】
その特徴は、到来方向にヌルを作ってスペクトラムをサーチし、別な到来波のピークがあればその方向にもヌルを作り、残余電力が十分小さくなるまで段階的にヌルを作って到来波を一つ一つ特定していく方法である。
【０１７７】
他の到来波数推定アルゴリズムとは、独立に到来波数を推定できる。比較的スナップショット数が小さい範囲でＭＵＳＩＣ法と同等以上の精度が出る。また、固有値解析を行わないため、到来波数が小さい場合はＭＵＳＩＣ法よりも演算時間が短いこともある。また、固有値解析をしないため、相関性のある信号でも分離できると考えられるが、スペクトラムは分解能が悪く（ヌルを作らない場合ＢＦ法と同じ）、近接した到来波の分離は難しいと考えられ、到来方向推定の初段のサーチに向いている可能性がある。
【０１７８】
アルゴリズムは、（１）　最初はヌルを指定せずに、サーチする方向θについてスペクトラムＰ_Ｄを計算し、最大になる方向θ_１を求める。
【０１７９】
【数８０】

（２）　θ_１方向にヌルをつくってスペクトラムＰ_Ｄを計算する。
【０１８０】
【数８１】

但し、
【数８２】
Ｑ　＝　Ｉ　−　Ａ（Ａ^ＨＡ）^−１Ａ^Ｈ
Ａ　：　ヌルを作る方向のステアリング行列　　　…式（８２）
Ｙ　＝　ＱＸ　　Ｒ_Ｃ　＝　ＹＹ^Ｈ　／　Ｎ
（３）　サーチしてＰ_Ｄの最大値を求める。最大値がスレッショルド値より大きければ、その方向θ_２を求め、θ_１，　θ_２の方向にヌルをつくるようにして（２）の手順に戻る。
【０１８１】
（４）　（２），（３）の手順を、Ｐ_Ｄの最大値がスレッショルド値を下回るまで繰り返す。Ｌ回繰り返したとすると、θ＝　［θ_１，…，　θ_Ｌ　］が到来方向である。
【０１８２】
（変形例１５）
変形例１５は、図５に示す到来方位演算部８０のアルゴリズムとして、Ｂｌｉｎｄ−ＭＡＸＣＯＲ　（ＭＡＸｉｍｕｍ　ｓｐａｔｉａｌ　ＣＯＲｒｅｌａｔｉｏｎ　ｍｅｔｈｏｄ　ａｆｔｅｒ　Ｂｌｉｎｄ　ｉｄｅｎｔｉｆｉｃａｔｉｏｎ　ｏｆｔｈｅ　ｓｏｕｒｃｅｓ　ｓｔｅｅｒｉｎｇ　ｖｅｃｔｏｒ）　法　を採用する。
【０１８３】
Ｂｌｉｎｄ−ＭＡＸＣＯＲ法の特徴は、ブラインド分離アルゴリズムにより計算で求められた信号のステアリングベクトルと、サーチする方向θのステアリングベクトルの相関をとり、相関係数が最大になる方向をサーチすることにより到来方向を求める方法である。
【０１８４】
アルゴリズムは、（１）　ブラインド分離アルゴリズムにより信号分離を行い、信号のステアリングベクトルを求める。得られたｉ番目の信号のステアリングベクトルをａ’_ｉ（ｉ＝１．．Ｌ）とする。（２）　１〜Ｌの到来波について、サーチする方向θの範囲で相関係数α_ｉ（ｉ＝１．．Ｌ）を計算し、最大になる方位を求める。
【０１８５】
【数８３】

（変形例１６）
変形例１６は、図５に示す到来方位演算部８０のアルゴリズムとして、Ｂｌｉｎｄ−ＭＵＳＩＣ法を採用する。
【０１８６】
Ｂｌｉｎｄ−ＭＵＳＩＣ法の特徴は、ブラインド分離アルゴリズムにより計算で求められた信号のステアリングベクトルと、サーチする方向θのステアリングベクトルでスペクトラムを計算し、最大になる方向をサーチすることにより到来方向を求める方法である。
【０１８７】
アルゴリズムは、（１）　ブラインド分離アルゴリズムにより信号分離を行い、信号のステアリングベクトルを求める。得られたｉ番目の信号のステアリングベクトルをａ’_ｉ（ｉ＝１．．Ｌ）とし、ステアリング行列をＡ’＝［ａ’_１，…，　ａ’_Ｌ］とする。
【０１８８】
（２）　サーチする方向θの範囲でＢｌｉｎｄ−ＭＵＳＩＣスペクトラムＰ_ＢＭＵを計算し、大きい順にＬ個のピークとその方位を求める。
【０１８９】
【数８４】

（変形例１７）
変形例１７は、図５に示す到来方位演算部８０のアルゴリズムとして、　ＥＳＰＲＩＴ　（Ｅｓｔｉｍａｔｉｏｎ　ｏｆ　Ｓｉｇｎａｌ　Ｐａｒａｍｅｔｅｒｓ　ｖｉａ　Ｒｏｔａｔｉｏｎａｌ　Ｉｎｖａｒｉａｎｃｅ　Ｔｅｃｈｎｉｑｕｅｓ）法を採用する。
【０１９０】
ＥＳＰＲＩＴの特徴は、干渉計の原理と同じく、２つの同一のサブアレイを含むアレイに適用される。２つのサブアレイは既知のベクトルΔだけ平行移動した関係にあるため、それぞれの到来波についてサブアレイ間で受信する信号は位相が異なっているだけである。到来方向θ_ｌは、サブアレイ間の位相差から求まる。
【０１９１】
他のアルゴリズムと異なり、スペクトラムサーチをしない。また、２つの同一パターンをもったサブアレイという制約はあるが、ステアリングベクトルは未知でも良いという特徴がある。但し、円形アレイには適さないため、直接ＥＳＰＲＩＴの適用はしない。
【０１９２】
アルゴリズムは、（１）　Ｒを固有値解析し、信号部分空間Ｅ_Ｓを求める。なお、到来波数を推定しておく必要がある。
【０１９３】
（２）　Ｅ_Ｓから２つのサブアレイの部分Ｅ_ＸとＥ_Ｙを取り出す。（全体のアレイから２つのサブアレイを取り出す行列をＪ_１，Ｊ_２で表すと、Ｅ_Ｘ　＝　Ｊ_１Ｅ_Ｓ，　Ｅ_Ｙ　＝　Ｊ_２Ｅ_Ｓ）
（３）　Ｅ_Ｙ　＝　Ｅ_ＸΨをΨについて解く。
【０１９４】
▲１▼　ＬＳ−ＥＳＰＲＩＴ　（Ｌｅａｓｔ　Ｓｑｕａｒｅｓ　ＥＳＰＲＩＴ）
【数８５】
Ψ　＝　（Ｅ_Ｘ ^ＨＥ_Ｘ）^−１Ｅ_Ｘ ^ＨＥ_Ｙ　　　　…式（８５）
▲２▼　ＴＬＳ−ＥＳＰＲＩＴ　（Ｔｏｔａｌ　Ｌｅａｓｔ　Ｓｑｕａｒｅｓ　ＥＳＰＲＩＴ）
Ｅ_Ｈを計算する。
【０１９５】
【数８６】
Ｅ_ＸＹ　＝　［Ｅ_Ｘ　Ｅ_Ｙ］　　　　…式（８６）
Ｅ_Ｈ　＝　Ｅ_ＸＹ ^Ｈ　Ｅ_ＸＹ
Ｅ_Ｈを固有値分解し、固有値の大きい順に固有ベクトルを並べる。Ｖ　＝　［ｖ_１，…，　ｖ_２Ｌ］　　　　　　　：　ｖ_１，…，　ｖ_２Ｌ　：　Ｒの固有ベクトル（固有値の大きい順）
【数８７】
Ｖ_Ｎ　＝　　　　＝　［ｖ_Ｌ＋１，…，　ｖ_２Ｌ］：Ｇ_Ｘ：　Ｖ_Ｎの上半分、Ｇ_Ｙ：　Ｖ_Ｎの下半分
Ψ　＝　−Ｇ_Ｘ Ｇ_Ｙ ^−１　　　　　…式（８７）
（４）　Ψを固有値解析し、固有値φ_ｌ　（　ｌ　＝　１，．．，Ｌ　）からθ_ｌ（　ｌ　＝　１，．．，Ｌ　）を求める。
【０１９６】
【数８８】
θ_ｌ＝　ｓｉｎ^−１　　　　　ａｒｇ（φ_ｌ）　　　　…式（８８）
（変形例１８）
変形例１８は、図５に示す到来方位演算部８０のアルゴリズムとして、　ＶＥＳＰＡ　（Ｖｉｒｔｕａｌ　ＥＳＰｒｉｔ　Ａｌｇｏｒｉｔｈｍ）　法を採用する。
【０１９７】
ＥＳＰＲＩＴ法には、２つの同一パターンを持つサブアレイを含むアレイという制約があったが、仮想素子を作ることによって、２つの同一パターンを持つアンテナを含むアレイにＥＳＰＲＩＴを適用できるように拡張したアルゴリズムがＶＥＳＰＡである。仮想素子はｃｕｍｕｌａｎｔを用いて作る。
【０１９８】
ＥＳＰＲＩＴと同様、スペクトラムサーチ不要で、同一パターンのアンテナ２本があればステアリングベクトルは未知でも良いという特徴がある。
【０１９９】
アルゴリズムは、（１）　４次のｃｕｍｕｌａｎｔにより共分散行列Ｃの各要素を計算する。（１番目とＫ番目のアンテナパターンが等しいとする。）
【数８９】

（２）　共分散行列Ｒの代わりにＣについてＥＳＰＲＩＴアルゴリズムを適用する。２つのサブアレイはアンテナ１〜ＫとアンテナＫ〜２Ｋ−１（Ｋ＋１〜２Ｋ−１は仮想アンテナ）である。
【０２００】
以上詳述したように、信号分離部７０は、確率分布の独立性に基づく分離法、勾配法、不動点法、ＦａｓｔＩＣＡ法、相関関数行列の同時対角化法、ＪＡＤＥ法、因子分析法、情報量を最大化する手法のいずれか１つを処理アルゴリズムとして採用することで、信号数推定処理部２０から出力された観測信号から同時に存在する入射信号をそれぞれ分離することができる。
【０２０１】
また、到来方位探知部８０は、ＭＵＳＩＣ法、ｃｕｍｕｌａｎｔ−ＭＵＳＩＣ法、ＷＳＦ法、ＵＳＦ法、ＭＬ　法、最小分散法、ＢＦ法、ＤＯＳＥ法、Ｂｌｉｎｄ−ＭＡＸＣＯＲ法、Ｂｌｉｎｄ−ＭＵＳＩＣ法、ＥＳＰＲＩＴ法、ＶＥＳＰＡ法のいずれか１つを処理アルゴリズムとして採用することで、信号数推定処理部２０から出力された観測信号から同時に存在する入射信号の到来方位を分離することができる。
【０２０２】
【発明の効果】
請求項１記載の本発明によれば、複数の入射信号が同時に入射される複数のセンサからの観測信号に基づいて、入射信号の信号数を推定して複数の観測信号の中から当該信号数分の観測信号を選択して出力するので、有効な観測信号のみに基づいて所定の信号処理を行うことができ、受信装置の処理速度と処理結果の精度を向上することができる。
【図面の簡単な説明】
【図１】本発明の第１の実施の形態に係る受信装置の構成を概略的に示すブロック図である。
【図２】信号数推定処理部２０の内部構成を示す図である。
【図３】観測信号の固有値を示すグラフ（ａ）、観測信号がないときの固有値を示すグラフ（ｂ）、グラフ（ａ），（ｂ）間の差を示す図である。
【図４】入射信号の到来方位に関する測定方法を説明するための図である。
【図５】本発明の第２の実施の形態に係る受信装置の構成を概略的に示すブロック図である。
【符号の説明】
１〜７　第１〜第７アンテナ
１０　サンプリング部
２０　信号数推定処理部
２１　入力バッファ部
２２　第１固有値解析部
２３　第２固有値解析部
２４−１〜２４−７　減算器
２５　信号数決定部
２６　対象信号選択部
３０　ブラインド処理部
４０　分離信号出力処理部
５０　到来方位演算部
６０　到来方位出力処理部
７０　信号分離部
８０　到来方位探知部[0001]
TECHNICAL FIELD OF THE INVENTION
The present invention relates to a receiving device that measures a direction of arrival of a signal by arranging a plurality of sensors, and separates and receives a plurality of simultaneously existing signals.
[0002]
[Prior art]
Generally, in direction detection by super-resolution (super-resolution) and signal separation by independent component analysis, first, a covariance matrix is generated based on time-series data of an observed signal, and eigenvalue analysis is performed.
[0003]
It is necessary to estimate the number of signals from the distribution of the eigenvalues. Generally, in estimating the number of signals, a technique such as AIC (Akaike Information Criteria) or MDL (Minimum Description Length) is adopted.
[0004]
[Problems to be solved by the invention]
In general, among eigenvalues equal to the number of sensors, it is easy to make it difficult to distinguish which is a signal component and which is a noise component.
[0005]
Here, even if the method such as AIC or MDL described above is adopted, if an error occurs when estimating the number of signals, the performance is degraded in subsequent direction detection processing and signal separation processing.
[0006]
Therefore, when estimating the number of signals simultaneously included in observation signals from a plurality of sensors from the eigenvalue distribution, there has been a demand to clarify a boundary between a signal and noise.
[0007]
The present invention has been made in view of the above, and an object of the present invention is to perform predetermined signal processing based on only valid observation signals, and to improve processing speed and accuracy of processing results. It is to provide a device.
[0008]
[Means for Solving the Problems]
According to an embodiment of the present invention, there is provided a receiving apparatus for performing predetermined signal processing based on a received observation signal, wherein a plurality of sensors to which a plurality of incident signals are simultaneously incident, and the plurality of sensors are provided. And a signal number estimation processing section for estimating the number of incident signals based on the observation signals from the plurality of observation signals and selecting observation signals corresponding to the number of observation signals from a plurality of observation signals.
[0009]
According to a second aspect of the present invention, in order to solve the above-mentioned problem, the number estimation processing unit generates a first covariance matrix from time-series data and the number of data of the observation signal, and generates the first covariance matrix. A first eigenvalue analyzing unit for eigenvalue analyzing the matrix to obtain a first eigenvalue, and a second covariance matrix generated from time-series data and the number of data of the noise interference signal in the absence of the observation signal. A second eigenvalue analysis unit for eigenvalue analysis of the covariance matrix of Eq. 2 to obtain a second eigenvalue, and a difference between the first and second eigenvalues obtained from the first eigenvalue analysis unit and the second eigenvalue analysis unit. A subtractor to be obtained, a signal number determining unit that outputs a selection permission signal for the number of signals whose difference obtained by the subtractor is larger than a reference value, and a selection permission signal output from the signal number determination unit. Out of the plurality of observation signals And gist that a target signal selecting section for selecting and outputting measurement signals.
[0010]
According to a third aspect of the present invention, in order to solve the above-described problem, a separation unit that separates the plurality of incident signals by blind processing based on observation signals corresponding to the number of signals selected by the signal number estimation processing unit. For each of the plurality of incident signals, a vector operation unit that obtains a complex vector amount indicating a relative relationship between amplitude and phase distribution in the plurality of sensors, and a physical position of the plurality of sensors and the vector operation unit. Based on the complex vector amount, an arrival direction calculation unit that determines the arrival directions of the plurality of incident signals, and a signal that is separated by the separation unit and a signal that indicates the arrival direction obtained by the arrival direction calculation unit. The gist of the present invention is to provide an output unit for outputting in association with the output.
[0011]
In order to solve the above problem, the invention according to claim 4 performs eigenvalue analysis based on the number of observation signals corresponding to the number of signals selected by the signal number estimation processing unit, and sets eigenvectors obtained simultaneously as signal eigenvectors. A gist of the present invention is to provide a signal separation unit that separates a signal component of a captured signal into noise eigenvectors by using the signal eigenvector, and separates a plurality of simultaneously existing signals by independent component analysis of the mapped data.
[0012]
In order to solve the above problem, the invention according to claim 5 performs eigenvalue analysis based on the number of observation signals corresponding to the number of signals selected by the number-of-signals estimation processing unit, and calculates a signal eigenvector from an eigenvector obtained simultaneously. The gist of the invention is to provide a direction-of-arrival detecting unit that separates the eigenvector into noise eigenvectors and detects the direction of arrival of a signal by calculating the eigenvectors and the response data of the sensor.
[0013]
According to a sixth aspect of the present invention, in order to solve the above-mentioned problem, the signal separation unit includes a separation method based on independence of a probability distribution, a gradient method, a fixed point method, a FastICA method, and a simultaneous diagonalization method of a correlation function matrix. , The JADE method, the factor analysis method, or the method for maximizing the amount of information is adopted as the processing algorithm.
[0014]
According to a seventh aspect of the present invention, in order to solve the above-mentioned problem, the direction-of-arrival detecting unit includes a MUSIC method, a cumulant-MUSIC method, a WSF method, a USF method, a ML method, a minimum dispersion method, a BF method, a DOSE method, and a Blind method. -The gist is to adopt any one of the MAXCOR method, the Blind-MUSIC method, the ESPRIT method, and the VESPA method as the processing algorithm.
[0015]
BEST MODE FOR CARRYING OUT THE INVENTION
Hereinafter, embodiments of the present invention will be described with reference to the drawings.
[0016]
(First Embodiment)
In the receiving device according to the first embodiment of the present invention, an antenna is used as a sensor. In the following description, the number n of antennas is set to “7” and the number m of incident signals is set to “4”. However, the number of antennas and the number of incident signals are not limited to these and are arbitrary.
[0017]
The receiving apparatus according to the first embodiment separates and outputs an original signal from a mixed signal of a plurality of incident signals incident on a plurality of antennas, and two-dimensionally measures an arrival direction of each incident signal. And output.
[0018]
FIG. 1 is a block diagram schematically showing a configuration of a receiving apparatus according to the first embodiment of the present invention.
[0019]
The receiving apparatus includes first to seventh antennas 1 to 7, a sampling unit 10, a signal number estimation processing unit 20, a blind processing unit 30, a separated signal output processing unit 40, an arrival direction calculation unit 50, and an arrival direction output processing unit 60. It is composed of
[0020]
The first to seventh antennas 1 to 7 correspond to a plurality of sensors of the present invention. As these first to seventh antennas 1 to 7, antennas having directivity such as a parabolic antenna and a Yagi antenna are used. The intervals and heights at which the first to seventh antennas 1 to 7 are installed are arbitrary.
[0021]
The first antenna 1 has a plurality of incident signals (radio waves) S from the air.₁~ S₇Are sent to the sampling unit 10. Similarly, the second to seventh antennas 2 to 7 output a plurality of incident signals S from the air.₁~ S₇Are sent to the sampling unit 10.
[0022]
Although not shown, the sampling unit 10 includes a local oscillator, an intermediate frequency converter, an oscillator, and an A / D converter.
[0023]
The local oscillator generates a signal having an oscillation frequency necessary for converting the received high frequency to an intermediate frequency. The intermediate frequency converter amplifies the high frequency in the reception frequency band from the first to seventh antennas 1 to 7 and mixes the signal with the signal from the local oscillator to generate a sum or difference frequency between the signals. The signal is converted to a signal and further amplified. The amplified intermediate frequency signal is supplied to the A / D converter.
[0024]
The oscillator generates a sampling clock for sampling an analog signal. The A / D converter uses the sampling clock from the oscillator to sample the analog signal from the intermediate frequency converter and converts it into a digital signal. The digital signal converted by the A / D converter is supplied to the signal number estimation processing unit 20.
[0025]
As shown in FIG. 2, the signal number estimation processing unit 20 includes an input buffer unit 21, a first eigenvalue analysis unit 22, a second eigenvalue analysis unit 23, subtracters 24-1 to 24-7, a signal number determination unit 25, The target signal selector 26 is provided.
[0026]
The signal number estimation processing unit 20 estimates the number of incident signals based on the observation signals from the first to seventh antennas 1 to 7 and divides the number of observation signals from the plurality of observation signals into observation signals. Select and supply to the blind processing unit 30.
[0027]
Here, the input buffer unit 21 outputs the time-series data of the observation signal input from the sampling unit 10 to the first eigenvalue analysis unit 22, the second eigenvalue analysis unit 23, and the target signal selection unit 26.
[0028]
The first eigenvalue analysis unit 22 generates a covariance matrix Rxx from the time series data and the number of data of the observation signal input from the sampling unit 10 via the input buffer unit 21, and performs eigenvalue analysis on the covariance matrix Rxx. The eigenvalues Ex1, Ex2,..., Exn are obtained.
[0029]
The second eigenvalue analysis unit 23 calculates a covariance matrix Ryy from time-series data and the number of data of the observation signal Y (noise interfering signal) without the observation signal input from the sampling unit 10 via the input buffer unit 21. , And the eigenvalues Ey1, Ey2,..., Eyn are obtained by eigenvalue analysis of the covariance matrix Ryy.
[0030]
The subtracters 24-1 to 24-7 determine the difference ΔR from the two eigenvalue groups Rxx and Ryy obtained from the first eigenvalue analysis unit 22 and the second eigenvalue analysis unit 23.
[0031]
The signal number determination unit 25 counts +1 when each eigenvalue indicating the difference ΔR, that is, the value of the eigenvalue Ex1-Ey1, Ex2-Ey2,..., Exn-Eyn is larger than the reference value ΔRref. By setting it to 0, the number of signals k is determined from each of ΔR from 1 to n, and the selection permission signals EN1 to EN7 are output to the target signal selection unit 26.
[0032]
The target signal selection unit 26 uses the selection permission signals EN1 to EN7 output from the signal number determination unit 25 to generate n observation signals x₁(T)-x_nK observation signals selected and permitted from (t) are represented by observation signals X_{1, ..., k}It is sent to the blind processing unit 30 as (t).
[0033]
The blind processing unit 30 corresponds to the separation unit and the vector operation unit according to the present invention. The blind processing unit 30 uses an original signal (incident signal S₁~ S₇Are separated and extracted from the observation signal obtained by mixing any one of the two. The signals separated and extracted by the blind processing unit 30 are sent to the separated signal output processing unit 40. Further, the blind processing unit 30 calculates the above-described matrix A (see Expression (6) described below) and sends the matrix A to the arrival direction calculation unit 50.
[0034]
The separation signal output processing unit 40 corresponds to a part of the output unit of the present invention. The separation signal output processing unit 40 is configured by a D / A converter. The D / A converter converts the digital signal from the blind processing unit 30 into an analog signal and outputs the separated signal O₁~ O₄Is output as
[0035]
The arrival direction calculation unit 50 corresponds to the arrival direction calculation unit of the present invention. The arrival direction calculation unit 50 performs a phase difference calculation between sensors and a direction calculation by an interferometer method based on the matrix A from the blind processing unit 30. Details of these calculations will be described later. The signal indicating the direction of arrival calculated by the direction of arrival calculation unit 50 is sent to the direction of arrival output processing unit 60.
[0036]
The arrival direction output processing unit 60 is configured by a D / A converter. The D / A converter converts the digital signal from the arrival direction calculation unit 50 into an analog signal, and₁~ D₄And output to the outside.
[0037]
Next, the operation of the receiving apparatus according to the first embodiment of the present invention configured as described above will be described.
[0038]
Incident signal S observed at each antenna 1-7_{1, ..., 7}Is sent to the sampling unit 10.
[0039]
The intermediate frequency converter of the sampling unit 10 outputs the incident signal S_{1, ..., 7}, And converts the mixed signal into an intermediate frequency signal, and sends it to the A / D converter 41. The A / D converter converts the received analog signal into a digital signal by sampling using a sampling clock from an oscillator. The converted digital signal is the observation signal x_{1, ..., 7}As (t), it is sent to the signal number estimation processing unit 20.
[0040]
The number-of-signals estimation processing unit 20 calculates the observation signal x from the sampling unit 10._{1, ..., 7}Based on (t), the number k of incident signals is estimated, and observation signals corresponding to the number k of signals are selected from a plurality of observation signals and output to the blind processing unit 30.
[0041]
First, a method of generating a covariance matrix will be described as a first-stage process by the signal number estimation processing unit 20.
[0042]
Now, it is assumed that m (m is an integer of 2 or more) signals are simultaneously incident on an array of n (n is an integer of 2 or more) sensors. This state can be represented by the following equation (1).
[0043]
(Equation 1)
X_{1 ~ n}(T) = AS_{1 ~ n}(T) + N_{1 ~ n}(T) ... (1)
Where X_{1 ~ n}(T) is time-series data composed of complex numbers observed by each sensor, A is an n-row × m-column signal mixing matrix that can be determined by the arrangement and characteristics of the sensors, S_{1 ~ m}(T) is m incident signals including external noise, N_{1 ~ n}(T) is the internal noise in each sensor.
[0044]
When the above equation (1) is expressed in a matrix format, as shown in the following equation (2),
(Equation 2)
X = A · S + N (2)
Becomes Here, X is n rows × time series data sequence, A is n rows × m column, S is m rows × time sequence data sequence, and N is n rows × time sequence data sequence.
[0045]
The covariance matrix Rxx has n rows and n columns, and as shown in the following equation (3),
(Equation 3)
Rxx = XX^TNumber of time series data… (3)
Is represented by Where X^TIs the reverse conjugate transpose of X.
[0046]
The first eigenvalue analysis unit 22 generates a covariance matrix Rxx from the time series data and the number of data of the observation signal input from the sampling unit 10 via the input buffer unit 21 according to the above equation (3). The eigenvalues obtained by eigenvalue analysis of the matrix Rxx are referred to as Ex1, Ex2,..., Exn.
[0047]
FIG. 3A is a diagram showing the magnitudes of the eigenvalues Ex1, Ex2,..., Exn in a bar graph, and shows a state where the boundary between the target signal component and the noise interference component is unclear at the level P1.
[0048]
On the other hand, when an observation signal containing no signal component is represented by Y,
(Equation 4)
Y = A · Snosig + N (4)
Snosig is obtained by removing the target signal component of interest from the state of the incident signal S, and has only a noise interference component including an external noise component and an interference component.
[0049]
The second eigenvalue analysis unit 23 performs the covariance based on the time series data and the number of data of the observation signal Y (noise interference signal) input from the sampling unit 10 via the input buffer unit 21 in the same manner as in the equation (3). A matrix Ryy is generated, and eigenvalues obtained by eigenvalue analysis of the covariance matrix Ryy are referred to as Ey1, Ey2,..., Eyn.
[0050]
FIG. 3B is a diagram showing the magnitudes of the eigenvalues Ey1, Ey2,..., Eyn in a bar graph.
[0051]
Next, a method of removing noise interference components will be described as a second stage process by the signal number estimation processing unit 20.
[0052]
Here, when the difference ΔR between the two eigenvalue groups Rxx and Ryy is obtained,
(Equation 5)
ΔR = eigenvalue of Rxx−eigenvalue of Ryy =
Ex1-Ey1, Ex2-Ey2,..., Exn-Eyn (5)
Becomes
[0053]
The subtracters 24-1 to 24-7 calculate the difference ΔR between the two eigenvalue groups Rxx and Ryy according to the above equation (5).
[0054]
FIG. 3C is a diagram showing a difference ΔR between the two eigenvalue groups Rxx and Ryy, that is, the magnitudes of the eigenvalues Ex1-Ey1, Ex2-Ey2,..., Exn-Eyn in a bar graph. This shows a state in which the boundary between the target signal component and the noise interference component is clear. As shown in FIG. 3C, the target signal component stands out in this distribution.
[0055]
The number-of-signals determination unit 25 counts +1 when the values of the eigenvalues Ex1-Ey1, Ex2-Ey2,..., Exn-Eyn are larger than the reference value ΔRref shown in FIG. By doing so, the number of signals k is determined from each of ΔR from 1 to n, and the selection permission signals EN1 to EN7 that become high level when the eigenvalue is larger than the reference value ΔRref are output to the target signal selection unit 26. . It should be noted that the reference value ΔRref may be, for example, an average value obtained by calculating the sum of the eigenvalues Ex1-Ey1, Ex2-Ey2,..., Exn-Eyn, and dividing the sum by the number n.
[0056]
The target signal selection unit 26 uses the high-level selection permission signals EN1 to EN7 output from the signal number determination unit 25 to generate n observation signals x₁(T)-x_nThe k observation signals selected and permitted by the high-level selection permission signals EN1 to EN7 from (t) are represented by observation signals X_{1, ..., k}It is sent to the blind processing unit 30 as (t).
[0057]
As a result, since the blind processing unit 30 performs the calculation using the k observation signals, the amount of calculation can be reduced as compared with the conventional case where the calculation is performed with n observation signals, and the calculation time can be reduced. Moreover, since only the observation signal whose signal component is larger than the noise component can be used, clear blind processing can be performed.
[0058]
Generally, in azimuth measurement or null steering by super resolution, information corresponding to “A” in the above equation (1) is directly or indirectly estimated.
[0059]
For example, in the MUSIC method, a matrix composed of noise eigenvectors orthogonal to the matrix corresponding to “A” is obtained, and the azimuth is obtained from the correlation between the matrix and the steering vector.
[0060]
On the other hand, in the present embodiment, this “A” is obtained by so-called blind processing. Here, blind processing refers to incident signals S up to the same number as the number of sensors without prior knowledge of the responsiveness of the sensors, the properties of the signals, the direction of arrival, and the like._{1, ..., n}Is separated, and as a result, a matrix A is obtained. This algorithm is called independent component analysis (ICA: Independent Component Analysis) or blind signal separation (BSS: Blind Source Separation), and many known documents have been published.
[0061]
The matrix A indicates an amount corresponding to the amplitude and phase distribution in each sensor, and is an amount related to the arrival direction of the signal. In order to directly determine the direction of arrival of a signal from the matrix A, a steering vector is generally required. In the present invention, the direction of arrival of a signal is determined without using a steering vector.
[0062]
The blind processing unit 30 calculates the observation signal X_{1, ..., 4}Using the four sampling time series data sent as (t), the signal Y_{1, ..., 4}(T) is separated and extracted, and a matrix A is obtained.
[0063]
That is, the blind processing unit 30 outputs the observation signal X_{1, ..., n}The signal y (t) in the selected band is separated and extracted by the blind signal separation algorithm using the n sampling time series transmitted as (t).
[0064]
Hereinafter, a procedure of separation and extraction in the blind processing unit 30 will be described in detail.
[0065]
First, the fundamental problem of blind signal separation is defined here. It is assumed that the signal source is given by the vector of the following equation (6).
[0066]
(Equation 6)
s (t) = (s₁(T), ..., s_n(T))^T{T = 0, 1, 2,...} Equation (6)
Here, it is assumed that s (t) is n incident signals, has an average of “0”, and is independent of each other. T represents transposition.
[0067]
Observation means time-series data observed and band-limited by each antenna 1 to 7,
(Equation 7)
x (t) = (x₁(T), ..., x_n(T))^T{T = 0, 1, 2, ...} Equation (7)
It shall be represented by This can be considered as a signal observed by each of the antennas 1 to 7 of 1,..., N. Generally, the number of antennas does not always match the number of signal sources, but it is assumed here that they match.
[0068]
In a simple ICA problem, between s (t) and x (t),
(Equation 8)
x (t) = As (t) Equation (8)
Assume a simple linear relationship A is a real matrix of a signal mixing matrix (n rows × n columns) determined by the arrangement and characteristics of the antennas 1 to 7. x (t) is separated into independent signal components without knowledge of s (t) or A.
[0069]
That is, by obtaining a certain n × n real matrix,
(Equation 9)
y (t) = Bx (t) Equation (9)
The purpose of ICA is to reconstruct mutually independent y (t) determined by
[0070]
One of the solutions to this problem is a separation method based on the independence of probability distributions. Each s_iUnder the assumption that (t) is (strong) stationary and not Gaussian, y_iVarious methods have been proposed for obtaining B so that (t) is independent of each other, and many of them can be summarized as follows. Assuming that y (t) is a strongly stationary process, the density function of the simultaneous distribution is
(Equation 10)
p (y) = p (y₁, ..., y_n) ... Equation (10)
Then the definition of independence is p (y_i) To p (y) y_iAs a marginal distribution on,
[Equation 11]

And multiply.
[0071]
The distance between the joint distribution and the marginal distribution Kullback-LeiblerｌDivergence (Kullback-Leibler information amount) is:
(Equation 12)

Becomes Here, H (Y; B) is the entropy of the joint distribution p (y), and H (Y_iB) is a marginal distribution p (y_i) Entropy.
[0072]
This is ｛Y_i相互 (i = 1,..., N). Based on the assumption that the signal source is not normally distributed, KL (B) becomes p (y_i) Is “0” only when they are independent of each other. These are determined by p (x) and B.
[0073]
Note that p (y) dy = p (x) dx and p (y) = p (x) / | B | (| B | is the determinant of B).
(Equation 13)

Becomes
[0074]
On the other hand, the entropy of the marginal distribution is
[Equation 14]

It is. Therefore,
[Equation 15]

Becomes
(Equation 16)

Then, the correct B can be obtained as the steepest descent method. The problem in equation (16) is that the inverse matrix (B^T)^-1Is calculated.
[0075]
Regarding the convergence, since this may be multiplied by any positive definite matrix, B^TMultiplying by B (which corresponds to a Rieman metric on a manifold of a regular matrix)
[Equation 17]
ΔB∝ (IE_x[Ψ (y) y^T] B Equation (17)
Becomes a new learning rule. By the assumption of stationarity, p (s), p (x) and p (y) are temporally independent. Under this assumption, the ensemble average can be replaced by a time average.
[0076]
(Equation 18)

Therefore, if η is a positive constant and the parameters are updated according to the following equation (19) every time data is observed, B_tIs obtained.
[0077]
[Equation 19]
B_{t + 1}= B_t+ Η (I-ψ (y) y^T) B_t… Equation (19)
Here, naturally, the problem is that p (y_i) Or ψ (y). Usually, this uses a parametric nonlinear function or a statistical expansion method.
[0078]
A rough idea is that if p (y_i) Is a normal distribution, then ψ (y) is a linear function. On the other hand, when the tail is "heavier" than the normal distribution (sub-Gaussian), it is better to approximate with a polynomial, and when the tail is "lighter" than the normal distribution (super-Gaussian), it is better to approximate with a sigmoid function. It is good. The sigmoid function works well because the tail of the audio signal is "light".
[0079]
The signal Y separated and extracted by the blind processing unit 30_{1, ..., 4}(T) is sent to the separation signal output processing unit 40. In the separation signal output processing unit 40, the signal Y transmitted as a digital signal from the blind processing unit 30_{1, ..., 4}(T) is converted to an analog signal, and the separated signal O₁~ O₄To the outside.
[0080]
Note that the above ICA solution is an example, and the present invention is not limited to this.
[0081]
On the other hand, the matrix A obtained by the blind processing unit 30 is sent to the arrival direction calculation unit 50. The arrival direction calculation unit 50 performs the following inter-sensor phase difference calculation based on the received matrix A.
[0082]
That is, first, a vector for a corresponding signal is extracted from the matrix A. Here, the vector is represented as the following equation (20).
[0083]
(Equation 20)
sig1 = [a11 + jb11, a12 + jb12, a13 + jb13, a14 + jb14] Expression (20)
sig2 = [a21 + jb21, a22 + jb22, a23 + jb23, a24 + jb24]
sig3 = [a31 + jb31, a32 + jb32, a33 + jb33, a34 + jb34]
sig4 = [a41 + jb41, a42 + jb42, a43 + jb43, a44 + jb44]
In the above equation (20), Sig1 is the incident signal S₁Represents the state when only one exists, the complex vector “a11 + jb11” represents the amplitude and phase of the signal observed by the first antenna 11, and “a12 + jb12” represents the amplitude and phase of the signal observed by the second antenna 12. "A13 + jb13" represents the amplitude and phase of the signal observed by the third antenna 13, and "a14 + jb14" represents the amplitude and phase of the signal observed by the fourth antenna 14. The same applies to Sig2, Sig3 and Sig4.
[0084]
From this equation (20), the matrix A reflects the phase of the signal observed by each of the antennas 1 to 7, and from the distribution of this phase, the incident signal S₁~ S₇It can be understood that the arrival direction can be estimated. The interferometer method described later obtains the arrival direction of the signal from the distribution of the phase.
[0085]
From this vector, the phase difference between the sensors (beams) is determined. From the above equation (20), for example, the phase difference between the first antenna 1 and the second antenna 2 can be obtained from “phase angle (a11 + jb11) −phase angle (a12 + jb12)”, and “phase angle (a13 + jb13) −phase angle” (A14 + jb14) ", the phase difference between the third antenna 4 and the fourth antenna 4 can be obtained.
[0086]
Next, the arrival direction calculation unit 50 performs a direction calculation. That is, from any of the phase differences determined above and the physical positional relationship of the antenna, the incident signal S is determined by a normal interferometer method.₁Find the direction of arrival.
[0087]
Hereinafter, the incident signal S₁The method of obtaining the arrival direction of is described in detail. Now, as shown in FIG. 4, the distance between the first antenna 1 and the second antenna 2 is d, and the incident signal S at an angle θ is perpendicular to the line connecting the first antenna 1 and the second antenna 2.₁Shall be incident. In this case, the incident signal S₁Is the distance y from the arrival at the first antenna 1 to the second antenna 2 is “dsin θ”. When this is converted into a phase, the phase difference φ between the first antenna 1 and the second antenna 2 becomes
(Equation 21)
φ = (2π / λ) · dsin θ Expression (21)
Is represented by Here, λ is a wavelength. This phase difference φ is the observed data.
[0088]
Therefore, by calculating θ from the above equation (21), the incident signal S₁Can be obtained. That is,
(Equation 22)
φ = (2π / λ) · dsin θ Expression (22)
Therefore, θ can be obtained from the following equation (23).
[0089]
[Equation 23]
θ = sin^-1(Φ · 1 / ((2π / λ) · d)) Equation (23)
Hereinafter, the incident signal S₂, Incident signal S₃And the incident signal S₄The direction of arrival corresponding to is obtained. The signal Z indicating the direction of arrival determined in this way_{1, ..., 4}(T) is sent to the arrival direction output processing unit 60.
[0090]
In the arrival direction output processing unit 60, the signal Z transmitted from the arrival direction calculation unit 50 as a digital signal is output._{1, ..., 4}(T) is converted to an analog signal, and the azimuth signal D₂~ O₄To the outside.
[0091]
As a result, the four separated signals are output as a time axis waveform together with the signal indicating the direction of arrival. Thus, the observation signal X_{1, ..., 4}When the processing on the sampling data in the predetermined section in (t) is completed, the processing is performed on the sampling data in the next section, and the same processing is repeated thereafter.
[0092]
In the first embodiment, an example in which an antenna is used as a sensor has been described, but the present invention is not limited to this. For example, a microphone that detects voice, a sensor that detects various states of a living body, or the like can be used as the sensor.
[0093]
As described in detail above, according to the receiving apparatus according to the first embodiment of the present invention, the signal number estimation processing unit 20 based on observation signals from a plurality of sensors to which a plurality of incident signals are simultaneously incident. Since the number of observation signals is estimated and the number of observation signals corresponding to the number of observation signals are selected and output from a plurality of observation signals, predetermined signal processing can be performed based on only valid observation signals. The processing speed of the receiving device and the accuracy of the processing result can be improved.
[0094]
Further, in the signal number estimation processing unit 20, the first eigenvalue analysis unit 22 generates a first covariance matrix from the time series data of the observation signal and the number of data, and analyzes the first covariance matrix by eigenvalue analysis. 1 is obtained, a second covariance matrix is generated by the second eigenvalue analysis unit 23 from the time-series data and the number of data of the noise interference signal in the absence of the observation signal, and the second covariance matrix is used as the eigenvalue. The second eigenvalue is obtained by analysis, the difference is obtained from the first and second eigenvalues obtained by the subtracters 24-1 to 24-7, and the difference obtained by the signal number determination unit 25 is smaller than the reference value. A selection permission signal corresponding to the number of signals to be increased is output, and the target signal selection unit 26 selects and outputs the observation signals corresponding to the number of signals from a plurality of observation signals in accordance with the selection permission signal. It is effective to clearly distinguish the boundary from noise Measuring signal can only perform predetermined signal processing based on, it is possible to improve the accuracy of the processing speed and the processing result of the reception apparatus.
[0095]
Further, the blind processing unit 30 (vector calculation unit) separates a plurality of incident signals by blind processing based on the observation signals for the selected number of signals, and obtains the amplitude of each of the plurality of incident signals at the plurality of sensors. And the complex vector amount indicating the relative relationship between the phase distributions, and the arrival direction calculation unit 50 determines the arrival directions of the plurality of incident signals based on the physical positions and the complex vector amounts of the plurality of sensors, and is separated. By outputting the signal and the signal indicating the direction of arrival in association with each other, the direction of arrival of the incident signal can be accurately obtained.
[0096]
Further, the signal separation unit 70 separates eigenvectors obtained simultaneously when the eigenvalue analysis is performed based on the observed signals of the number of signals selected by the signal number estimation processing unit 20 into a signal eigenvector and a noise eigenvector. By mapping only the signal component of the captured signal using the eigenvector and separating a plurality of signals that exist simultaneously by independent component analysis of the mapped data, the incident signal can be accurately separated.
[0097]
Further, in the arrival direction calculation unit 80, when eigenvalue analysis is performed based on the number of observation signals for the number of signals selected by the number-of-signals estimation processing unit 20, the eigenvectors are simultaneously separated into eigenvectors and noise eigenvectors. By detecting the arrival direction of the signal by calculating the eigenvector and the response data of the sensor, the arrival direction of the incident signal can be accurately obtained.
[0098]
(Second embodiment)
FIG. 5 is a block diagram schematically showing a configuration of a receiving apparatus according to the second embodiment of the present invention.
[0099]
The feature of the present embodiment is that a separation method based on independence of probability distribution is adopted as an algorithm of the signal separation unit 70.
[0100]
In this separation method, each component s of the observation signal is_iAssume that (t) follows a probability distribution that is not normal. However, the probability distribution is unknown.
[0101]
Now, let the joint probability density variable of y (t) be
(Equation 24)
y = (y₁, ..., y_n)^T… Equation (24)
If the observation signal x (t) can be correctly separated by the matrix W, each component y of y_iBecomes independent. each component y of y_iSay that independence means y_iEven if you know the value of y_iOther than y_i, ..., y_nThe value of is unknown.
(Equation 25)

Product of simultaneous density function = marginal density function
Becomes
[0102]
Here, as one of the references, the distance between the joint distribution and the peripheral distribution estimated from the restored signal can be considered. There can be various statistical distances between distributions, but in independent component analysis, Kullback-Leibler divergence is often used.
[0103]
The definition of Kullback-Leibler @ divergence is as follows.
[0104]
(Equation 26)

Here, the increased information amount when the distribution of z is replaced by the distribution of x is shown. When the value is large, the distribution is not similar, and when the value is small, the distribution is similar. Represents
[0105]
Therefore, from equation (12),
[Equation 27]

Substituting into {(13) Equation},
[Equation 28]

And y_iAre independent, then equation (9) is zero, that is, it is minimal only when y (t) is independent. To find the matrix W, the gradient of I (y | W) with respect to W is found, and is found by the steepest descent method.
[0106]
(Equation 29)

However,
[Equation 30]

W can be obtained by updating as follows. In calculation, the inverse matrix (W^T)^-1Is a problem. Convergence can be multiplied by a positive definite matrix (W^T)^WMultiply by
(Equation 31)

The calculation amount is smaller and the convergence is faster.
[0107]
This updating rule includes ψ (y) and cannot be calculated unless the form of the density function is known. However, a solution can be correctly obtained even if approximated by an appropriate parametric nonlinear function or statistical expansion method.
[0108]
(Modification 1)
In the first modification, a gradient method is employed as an algorithm of the signal separation unit 70 shown in FIG.
[0109]
The gradient method is a method of obtaining the W that minimizes the mutual information shown in Expression (9) by the steepest descent method, and is a method of obtaining W by Expression (17).
[0110]
Specifically, there is a fixed point method.
[0111]
In the fixed point method, a conditional optimization problem in which a function と す る (w) having a vector w as an argument is maximized or minimized under a condition that a norm is {w} is considered.
[0112]
The problem is solved using Lagrange's undetermined coefficient method,
(Equation 32)
{(w) + λ ({w} −1) → min} Equation (32)
Is rewritten as At the minimum, this derivative is 0, so
[Equation 33]

Solving this for the second term gives
(Equation 34)

Becomes Since w is normalized by {w} = 1, it is sufficient to obtain λ by normalization. The actual calculation is
(Equation 35)

Is repeated until convergence.
[0113]
(Modification 2)
In the second modification, the FastICA method is employed as the algorithm of the signal separation unit 70 shown in FIG.
[0114]
Hyvarinen et al. Discuss the conditions for the convergence speed of the fixed-point method, and in particular, state that their method using kurtosis (FastICA) converges to third order.
[0115]
As φ (w),
[Equation 36]
φ (w) = E ((w · x)⁴) -3E ((wx)²)²… Equation (36)
Using an evaluation function called, find the maximum and minimum values of this using the fixed point method,
This is an algorithm for extracting independent components.
[0116]
Now, the average of the original signal is 0, the variance is 1, and the mixing matrix A is an orthogonal matrix.
(37)
v (t) T = w (t) TA {v} = 1} Equation (37)
Consider the behavior of
[0117]
(The reason why A is limited to an orthogonal matrix is that whitening is performed and the signal is once decorrelated.)
As a result of this rewriting, if w captures one component, v approaches one of the components, and the vector approaches 0.
[0118]
When an evaluation function is written using the original signals s and v,
[Equation 38]
E ((vs)⁴) -3E ((vs)²)²= E ((v · s))⁴-3 {v}⁴Formula (38)
Becomes
[0119]
Given that the components of s are mean 0, variance 1, and independent of each other, consider the ith component,
[Equation 39]

Becomes E (s_i ⁴) = 1, so the last parenthesis is the fourth order cumulant. The normalized vector is the new vector, but considering the ratio,
(Equation 40)

And v_iIs close to 1 and v_jIs close to 0, the above equation gives v_jThis means that the speed of approaching zero to zero is on the order of the third power.
[0120]
(Modification 3)
Modification 3 employs a simultaneous diagonalization method of a correlation function matrix as an algorithm of the signal separation unit 70 shown in FIG.
[0121]
The Jacobi method diagonalizes a real symmetric matrix. That is, it is a classical method for finding an eigenvalue. The basic principle is a method of obtaining orthogonal transform by repeating similarity transform by Givens transform.
[0122]
A proposed method is a simultaneous diagonalization algorithm of a correlation function matrix and a reduced kurtosis.
[0123]
Here, a separation method based on a time structure will be described as a simultaneous diagonalization method of a correlation function matrix.
[0124]
The algorithm using the correlation function matrix considers a correlation signal having a time difference. When the observed signal can be written as x (t) = As (t), assuming the independence of the original signal, its correlation function matrix is
(Equation 41)

Becomes Here, considering a complex-valued time series, * is the complex conjugate transpose, <•> is the averaging operation, and R_si(Τ) is the original signal s_iLet (t) denote the autocorrelation function.
[0125]
Since this is also a symmetric matrix in which the diagonal matrix is similarly transformed by the same matrix, it suffices to search for the matrix W at various time differences at the same time.
[0126]
Specifically, the observation signal x (t) is decorrelated using the correlation matrix at τ = 0, and is calculated by the Jacobi method so as to reduce the off-diagonalization component of the correlation matrix at τ ≠ 0. .
[0127]
(Modification 4)
Modification 4 employs a JADE (Joint \ Approximate \ Diagonization \ of \ Eigen-materials) method, which is an algorithm based on the simultaneous diagonalization of Kurtosis, as the algorithm of the signal separation unit 70 shown in FIG.
[0128]
The simultaneous diagonalization algorithm of the reduced kurtosis is represented as a fourth-order cross-curtosis.
(Equation 42)

Consider the contraction of
[0129]
Assuming that the observed signal is uncorrelated by whitening with an average of 0, the original signal and the observed signal have a certain orthogonal transform U = (u_{1, … , N}), They are linked in a relationship of X = US. From the independence assumption,
[Equation 43]

And the matrix M = (m_ij) As a Kurtosis matrix reduced by
[Equation 44]

Given
[Equation 45]

further,
[Equation 46]

As
[Equation 47]

Becomes Therefore, an appropriate set of ｛N_rTake｝ and ｛C (N_r) It can be seen that it is sufficient to calculate 直交 and obtain an orthogonal matrix for simultaneous diagonalization of｝. This situation
The same applies to matrices.
[0130]
Here, Cardoso {and} A. Souloumic is ｛N_r｝
[Equation 48]

Where e_kProposes a method in which a unit vector in which only the k component is 1 or a kurtosis itself (kurtosis) itself is expanded into an eigenvalue matrix, and an appropriate number of eigenmatrices are used in descending order of eigenvalues.
[0131]
(Modification 5)
Modification 5 employs a model including noise as the algorithm of the signal separation unit 70 shown in FIG.
[0132]
Equation (3) is shown as an ICA model, but external noise may not be negligible in actual measurement. If the external noise is large, it affects the separation and restoration of the signal. Here, the observation noise is
[Equation 49]
ε (t) = [E₁(T), ..., E_n(T)]^T… Equation (49)
Then, the observation signal of the sensor is
[Equation 50]
x (t) = As (t) + ε (t) Equation (50)
Becomes When calculating, the algorithm does not distinguish between signal and noise,
(Equation 51)

If the dimension of the signal is m, A 'is an n × (m + n) matrix, the number of signal sources is larger than the dimension of the observed signal, and there is no linear matrix solution. In addition, decorrelation by whitening becomes meaningless because a signal including noise is decorrelated.
[0133]
Modification 5 employs a factor analysis method as a preprocessing of ICA as a method of reducing noise as an algorithm of the signal separation unit 70 shown in FIG.
[0134]
In this factor analysis method, a factor analysis mixed / restored model is adopted.
[0135]
The factor analysis method divides an observation vector into two parts, a common factor existing between elements and an independent factor between the elements.
[0136]
That is, the observed value x∈R_nTo
(Equation 52)
x = a₁s₁+ A₂s₂+ ... a_ms_m+ Ξ = As + ξ Equation (52)
The common factor s∈R^mAnd the independent factor ξ∈RⁿIt can be decomposed using. Here, A is called a factor loading, s and ξ are independent of each other, and ξ is that each component is independent.
[0137]
Here, the covariance matrix of the original signal s is an m-dimensional unit matrix I_m, Ξ follow a normal distribution, and assume that the covariance is an n-dimensional diagonal matrix 、,
[Equation 53]
Σ = E (xx^T) = AA^T+ Ψ ... Equation (53)
And the covariance determined from the observations is
(Equation 54)

Is estimated based on The simplest method is least squares estimation, which minimizes the error between Σ and r.
[0138]
That is, the loss function L (A, Ψ) is given by r = AA^TBy + Ψ
[Equation 55]
L (A, Ψ) = tr (r-AA^T−Ψ)^T(R-AA^T−Ψ) ... Formula (55)
To minimize. By partially differentiating equation (42) with A and Ψ, the following solution for A and Ψ is obtained for A and に よ る by the least squares method.
[0139]
[Equation 56]
(R−Ψ) A = A (A^TA) ... Equation (56)
[Equation 57]
Ψ = diag (r-AA^T) ... Equation (57)
Next, considering the degrees of freedom, the degree of freedom of Σ is n (n + 1) / 2. On the other hand, the number of parameters of A and Σ is (m + 1) n. However, since the arbitrariness of the orthogonal matrix must be subtracted, subtracting m (m-1) / 2 means that this free parameter is smaller than the degree of freedom of Σ.
[Equation 58]

Here, by adding the condition of m <n and solving this,
[Equation 59]

Becomes This is a necessary and sufficient condition for a solution to exist.
[0140]
From the above, the relation between the observation signal x by the sensor and the signal に よ る by the sensor after the factor analysis for the model including noise is expressed as a generalized inverse matrix,
[Equation 60]
Q = [AΨ^-1A^TΨ^-1]^-1… Equation (60)
Using
[Equation 61]

Is obtained.
[0141]
(Modification 6)
Modification 6 employs a method of maximizing the amount of information (Infomax) as the algorithm of the signal separation unit 70 shown in FIG.
[0142]
Now, suppose that a random variable of an observation signal is X, and a random variable Y is obtained through a communication path G. The relationship between this,
(Equation 62)
Y = G (X) + N Equation (62)
N: Noise
At this time, designing the communication path so that the amount of information on the input X that the output Y has is as large as possible is a basic idea of the maximum information amount (Infomax).
[0143]
More specifically, this is a problem of maximizing the mutual information between X and Y under certain restrictions.
[0144]
In the case of independent component analysis, a combination of linear operation and nonlinear transformation is used.
[Equation 63]

And the amount of mutual information is
[Equation 64]

Becomes
[0145]
In this case, since there is only an ambiguity of the noise component between X and Y, the conditional entropy is equal to the noise entropy. In the concept of the maximum amount of information (Infomax), this part plays an important role, and after this, the discussion will be continued without considering noise.
[0146]
According to the steepest descent method, the derivative of the mutual information is the derivative of the entropy of Y,
[Equation 65]

Becomes In this case, paying attention to the handling of the non-linear conversion, if the relationship between the distribution of X and Y is one-to-one corresponding to y = G (Wx),
[Equation 66]

Can be written. However, J is Jacobian of this correspondence,
[Equation 67]

When specifically calculated, the differential of G is G ′,
[Equation 68]

[Equation 69]

Becomes When the entropy of Y is represented by the distribution of X,
[Equation 70]

Becomes this is,
[Equation 71]

By using
[Equation 72]

Then, as shown in [1], ψ (y_i) Can be updated to obtain W.
[0147]
(Modification 7)
Modification 7 employs the MUSIC (Multiple Signal) Classification (MUSIC) method as the algorithm of the arrival direction calculator 80 shown in FIG.
[0148]
The MUSIC method is the most representative of the Super @ Resolution algorithms. The covariance matrix is eigenvalue-decomposed and divided into signal and noise subspaces, and the direction of arrival is determined using the property that the signal is orthogonal to the noise subspace. This is a type of Noise Subspace Fitting (NSF). Since the orthogonality of space is used, a plurality of arrival directions can be obtained by one direction search. However, correlated incident signals cannot be separated.
[0149]
Hereinafter, an algorithm related to the MUSIC method will be described.
[0150]
(1) Analyze the eigenvalue of R and decompose it into subspaces of observation signals and noise Note that the above-described estimation of the number of incoming waves is required.
[0151]
(2) The MUSIC spectrum in the search direction θ is
[Equation 73]

Becomes Since the observed signal and the noise subspace are complementary spaces orthogonal to each other, P_MUHas a sharp peak in the direction of arrival of the incident signal. A direction search is performed to extract L peaks in descending order.
[0152]
(Modification 8)
Modification 8 employs the cumulant-MUSIC method as the algorithm of the arrival direction calculation unit 80 shown in FIG.
[0153]
The feature is that the algorithm is the same as the MUSIC method, but a fourth-order cumulant is introduced in the calculation of the covariance matrix. The fourth-order cumulant is said to suppress Gaussian noise, and is expected to be more resistant to noise than the MUSIC method.
[0154]
The algorithm is basically the same as that of the MUSIC method. Instead of the second-order covariance matrix R, a fourth-order cumulant matrix C is calculated. Here, the matrix C is
[Equation 74]
C = ｛c_ij{(I, j = 1 ... K)} Equation (74)
c_ij= Σ cum (X_i ^*, X_j, X_k ^*, X_k)
Becomes
[0155]
(Modification 9)
Modification 9 employs the WSF (Weighted Subspace Fitting) method as the algorithm of the arrival direction calculation unit 80 shown in FIG.
[0156]
The feature is that the MUSIC method uses the noise subspace to determine the direction of arrival, whereas the WSF method uses the signal subspace to determine the direction of arrival, and is a type of Signal {Subspace Fitting} (SSF). .
[0157]
Correlated observation signals can also be separated. Since the arrival direction of the L-wave must be multidimensionally searched and the amount of calculation increases, it is preferable to determine the initial value in advance by MUSIC or the like and limit the search range.
[0158]
The algorithm (1) performs eigenvalue analysis of R and decomposes it into a subspace of observation signals and noise. It is necessary to estimate the number of incoming waves in advance. (2) V below_WSFΘ = [θ₁,…, Θ_LSearch for [サ ー チ].
[0159]
[Equation 75]
V_WSF(Θ) = Tr (P_A ^⊥E_sW_sE_s ^H) ... Equation (75)
W_s= (Λ_s− Σ²I)²Λ_s ^-1
(Modification 10)
Modification 10 employs the USF {Unweighted {Subspace} Fitting} method as the algorithm of the arrival direction calculation unit 80 shown in FIG.
[0160]
The feature of the USF method is a kind of the SSF method like the WSF method._sIs replaced by a unit matrix I, and the (unweighted) {MD-MUSIC} (Multidimensional @ MUSIC) method is also called.
[0161]
A multidimensional search of the arrival direction of the L wave is required, and an initial value is preferably determined in advance by the MUSIC method to limit the search range, and an algorithm that converges by iterative calculation may be used.
[0162]
The algorithm (1) performs eigenvalue analysis of R and decomposes it into a subspace of observation signals and noise. At this time, it is necessary to estimate the number of incoming waves. (2) V below_USFΘ = [θ₁,…, Θ_L] Is obtained by a search or a convergence algorithm.
[0163]
[Equation 76]
V_USF(Θ) = Tr (P_A ^⊥E_sE_s ^H) ... Equation (76)
(Modification 11)
Modification 11 employs the ML {(Maximum Likelihood)} method as the algorithm of the arrival direction calculation unit 80 shown in FIG.
[0164]
The feature of the ML method is a general method for estimating an unknown parameter hidden from an observation amount. The parameters are determined so that the likelihood @ function represented by the product of the observed quantity and the probability density function (PDF) is maximized.
[0165]
The ML method includes two methods, {unconditional ML} (UML) (or stoichiometric ML) and {conditional ML} (CML) (or deterministic ML), each of which requires a multidimensional search of the arrival direction of the L wave. It is preferable that the initial value is determined in advance by the MUSIC method or the like to limit the search range, and this is also effective for correlated observation signals.
[0166]
The algorithm estimates the number of arriving waves by (1) the arriving wave number estimation algorithm.
(2) V below_UMLOr V_CMLΘ = [θ₁,…, Θ_L] Is obtained by a search or a convergence algorithm.
[0167]
[Equation 77]
V_UML{(Θ)} = {log} ({det} [P_ARP_A-Tr (P_A ^⊥) P_A ^⊥/ (KL)]) ... Equation (77)
V_CML(Θ) = Tr (P_A ^⊥R)
(Modification 12)
Modification 12 employs the minimum variance method as the algorithm of the arrival direction calculation unit 80 shown in FIG.
[0168]
The feature is also called Capon method, in which the main lobe is directed in one incoming wave direction and null is directed in the other incoming wave direction. The resolution is lower than the MUSIC method, and correlated observation signals cannot be separated. .
[0169]
The algorithm is as follows: (1) Spectrum P in search direction θ_CPIs calculated.
[0170]
[Equation 78]

A direction search is performed to extract L peaks in descending order. Although it is not necessary to directly estimate the number of incoming waves as an algorithm, it is necessary to search for L peaks.
[0171]
(Modification 13)
Modification 13 employs the BF method (beamformer) as the algorithm of the arrival direction calculation unit 80 shown in FIG.
[0172]
The feature is a method of scanning the main lobe to search for a direction in which the output power of the array increases. However, the angular resolution is poor, and when there are a plurality of arriving waves, separation is difficult and the direction becomes inaccurate.
[0173]
The algorithm is as follows: (1) Spectrum P in search direction θ_BFIs calculated.
[0174]
[Expression 79]

A direction search is performed to extract peaks in descending order.
[0175]
(Modification 14)
Modification 14 employs the DOSE {Direction of Source Elimination} method as the algorithm of the arrival direction calculation unit 80 shown in FIG.
[0176]
The feature is to make a null in the direction of arrival and search the spectrum.If there is another peak of the arriving wave, make a null in that direction as well, and gradually make a null in the direction until the residual power becomes sufficiently small. It is a method of specifying each one.
[0177]
The number of arriving waves can be estimated independently of other arriving wave number estimation algorithms. The accuracy equal to or higher than that of the MUSIC method is obtained in a range where the number of snapshots is relatively small. In addition, since the eigenvalue analysis is not performed, the calculation time may be shorter than the MUSIC method when the number of incoming waves is small. In addition, it is considered that since there is no eigenvalue analysis, it is possible to separate even correlated signals. However, the spectrum has poor resolution (same as the BF method when no null is created), and it is considered difficult to separate adjacent arriving waves. It may be suitable for the first-stage search of arrival direction estimation.
[0178]
The algorithm is as follows: (1) {Specify the spectrum P in the search direction θ without first specifying null._DIs calculated, and the maximum direction θ₁Ask for.
[0179]
[Equation 80]

(2) θ₁Create a null in the direction, spectrum P_DIs calculated.
[0180]
(Equation 81)

However,
[Expression 82]
Q = I-A (A^HA)^-1A^H
A: {Steering matrix in the direction to create null} formula (82)
Y = QX R_C= YY^H/ N
(3) Search and P_DFind the maximum value of. If the maximum value is larger than the threshold value, the direction θ₂And θ₁, Θ₂The procedure returns to the step (2) by making a null in the direction of (2).
[0181]
(4) The procedure of (2) and (3) is_DRepeat until the maximum value of falls below the threshold value. If it is repeated L times, θ = [θ₁,…, Θ_L] Is the direction of arrival.
[0182]
(Modification 15)
Modification Example 15 employs a Blind-MAXCOR (Maximum spatial, CORrelation, method, after, Blind, identification, of course, sources, steering, vector) algorithm as the algorithm of the arrival direction calculation unit 80 shown in FIG.
[0183]
The feature of the Blind-MAXCOR method is that the steering vector of the signal calculated by the blind separation algorithm is correlated with the steering vector in the search direction θ, and the direction of arrival is searched by searching the direction in which the correlation coefficient becomes maximum. It is a method of seeking.
[0184]
The algorithm is as follows: (1) Signal separation is performed by a 分離 blind separation algorithm, and a steering vector of the signal is obtained. Let the steering vector of the obtained i-th signal be a '_i(I = 1... L). (2) Correlation coefficient α in the range of search direction θ for arriving waves of 範 囲 1 to L_i(I = 1... L) is calculated, and the azimuth at which the maximum is obtained.
[0185]
[Equation 83]

(Modification 16)
Modification 16 employs the Blind-MUSIC method as the algorithm of the arrival direction calculation unit 80 shown in FIG.
[0186]
A feature of the Blind-MUSIC method is a method of calculating a spectrum by using a steering vector of a signal calculated by a blind separation algorithm and a steering vector in a search direction θ, and searching for a direction in which the signal reaches a maximum to determine an arrival direction. It is.
[0187]
The algorithm is as follows: (1) Signal separation is performed by a 分離 blind separation algorithm, and a steering vector of the signal is obtained. Let the steering vector of the obtained i-th signal be a '_i(I = 1... L) and the steering matrix is A ′ = [a ′₁, ..., 'a'_L].
[0188]
(2) Blind-MUSIC spectrum P in the range of search direction θ_BMUIs calculated, and L peaks and their directions are obtained in descending order.
[0189]
[Equation 84]

(Modification 17)
Modification Example 17 employs the {ESPRIT} (Estimation of Signal Signals via Via Rotational Invariance Technologies) algorithm as the algorithm of the arrival direction calculation unit 80 shown in FIG.
[0190]
The ESPRIT feature applies to arrays containing two identical sub-arrays, similar to the principle of interferometers. Since the two sub-arrays are in a relation of being translated by a known vector Δ, the signals received between the sub-arrays for the respective arriving waves only differ in phase. Arrival direction θ_lIs obtained from the phase difference between the subarrays.
[0191]
Unlike other algorithms, it does not perform a spectrum search. In addition, although there is a restriction of a subarray having two identical patterns, there is a feature that the steering vector may be unknown. However, since it is not suitable for a circular array, ESPRIT is not directly applied.
[0192]
The algorithm is: (1) Eigenvalue analysis of R and the signal subspace E_SAsk for. It is necessary to estimate the number of incoming waves.
[0193]
(2) @E_SFrom the two sub-array parts E_XAnd E_YTake out. (Matrix for extracting two subarrays from the whole array is J₁, J₂E_X= J₁E_S, E_Y= J₂E_S)
(3) E_Y= E_XSolve Ψ for Ψ.
[0194]
(1) LS-ESPRIT (Least Squares ESPRIT)
[Equation 85]
Ψ = (E_X ^HE_X)^-1E_X ^HE_Y… Equation (85)
(2) TLS-ESPRIT (Total Least Squares ESPRIT)
E_HIs calculated.
[0195]
[Equation 86]
E_XY= [E_XE_Y] ... Equation (86)
E_H= E_XY ^HE_XY
E_HIs subjected to eigenvalue decomposition, and eigenvectors are arranged in descending order of eigenvalues. V = [v₁, ..., v_2L]: V₁, ..., v_2L: Eigenvector of R (in order of eigenvalue)
[Equation 87]
V_N= = [V_{L + 1}, ..., v_2L]: G_X: V_NUpper half of G_Y: V_NLower half
Ψ = -G_X G_Y ^-1… Equation (87)
(4) Eigenvalue analysis of Ψ and eigenvalue φ_lFrom ({l} = {1, ..., L}) to θ_l({L} = {1,..., L}).
[0196]
[Equation 88]
θ_l= Sin^-1Arg (φ_l) ... Equation (88)
(Modification 18)
Modification 18 employs a {VESPA} (Virtual {ESPrit} Algorithm) method as an algorithm of the arrival direction calculation unit 80 shown in FIG.
[0197]
Although the ESPRIT method had a limitation of an array including two subarrays having the same pattern, an algorithm extended by applying virtual elements to an array including an antenna having two identical patterns so that ESPRIT can be applied to the array is described. VESPA. Virtual elements are created using cumulant.
[0198]
Similar to ESPRIT, there is a feature that a spectrum search is unnecessary and a steering vector may be unknown if there are two antennas having the same pattern.
[0199]
The algorithm calculates (1) each element of the covariance matrix C by a 4th order cumulant. (It is assumed that the first and Kth antenna patterns are equal.)
[Equation 89]

(2) Apply the ESPRIT algorithm for C instead of the covariance matrix R. The two sub-arrays are antennas 1 to K and antennas K to 2K-1 (K + 1 to 2K-1 are virtual antennas).
[0200]
As described in detail above, the signal separation unit 70 includes a separation method based on independence of the probability distribution, a gradient method, a fixed point method, a FastICA method, a simultaneous diagonalization method of a correlation function matrix, a JADE method, a factor analysis method, By employing any one of the methods for maximizing the amount of information as a processing algorithm, it is possible to separate simultaneously existing incident signals from the observation signals output from the signal number estimation processing unit 20.
[0201]
In addition, the direction-of-arrival detecting unit 80 includes the MUSIC method, the cumulant-MUSIC method, the WSF method, the USF method, the ML method, the minimum dispersion method, the BF method, the DOSE method, the Blind-MAXCOR method, the Blind-MUSIC method, the ESPRIT method, and the VESPA. By adopting any one of the methods as the processing algorithm, it is possible to separate the arrival directions of the incident signals present simultaneously from the observation signals output from the signal number estimation processing unit 20.
[0202]
【The invention's effect】
According to the first aspect of the present invention, the number of incident signals is estimated based on observation signals from a plurality of sensors to which a plurality of incident signals are simultaneously incident, and the number of signals is determined from the plurality of observation signals. Since the minute observation signal is selected and output, predetermined signal processing can be performed based only on the effective observation signal, and the processing speed of the receiving device and the accuracy of the processing result can be improved.
[Brief description of the drawings]
FIG. 1 is a block diagram schematically showing a configuration of a receiving device according to a first embodiment of the present invention.
FIG. 2 is a diagram illustrating an internal configuration of a signal number estimation processing unit 20;
FIG. 3 is a graph (a) showing an eigenvalue of an observation signal, a graph (b) showing an eigenvalue when there is no observation signal, and a diagram showing a difference between the graphs (a) and (b).
FIG. 4 is a diagram for explaining a method of measuring the direction of arrival of an incident signal.
FIG. 5 is a block diagram schematically showing a configuration of a receiving apparatus according to a second embodiment of the present invention.
[Explanation of symbols]
1-7 first to seventh antennas
10 Sampling unit
20 signal number estimation processing unit
21 Input buffer
22 1st eigenvalue analysis unit
23 second eigenvalue analyzer
24-1２４24-7 Subtractor
25 signal number determination unit
26 Target signal selector
30 blind processing unit
40 ° separation signal output processing unit
50 ° arrival direction calculation unit
60 ° arrival direction output processing unit
70 ° signal separation unit
80 ° arrival direction detector

Claims (7)

In a receiving device that performs predetermined signal processing based on the received observation signal,
A plurality of sensors to which a plurality of incident signals are simultaneously incident;
A signal number estimation processing unit for estimating the number of incident signals based on the observation signals from the plurality of sensors and selecting and outputting observation signals corresponding to the number of signals from the plurality of observation signals. A receiving device characterized by the above-mentioned. The number estimation processing unit,
A first eigenvalue analysis unit that generates a first covariance matrix from the time-series data of the observation signal and the number of data, and performs eigenvalue analysis on the first covariance matrix to obtain a first eigenvalue;
A second eigenvalue for generating a second covariance matrix from the time series data and the number of data of the noise interfering signal in the absence of the observation signal and eigenvalue analyzing the second covariance matrix to obtain a second eigenvalue An analysis unit;
A subtractor for calculating a difference between the first and second eigenvalues obtained from the first eigenvalue analysis unit and the second eigenvalue analysis unit;
A signal number determination unit that outputs a selection permission signal for the number of signals in which the difference obtained by the subtracter is larger than a reference value,
And a target signal selection unit that selects and outputs as many observation signals as the number of signals from the plurality of observation signals in accordance with the selection permission signal output from the signal number determination unit. The receiving device according to 1. A separation unit that separates the plurality of incident signals by blind processing based on observation signals for the number of signals selected by the signal number estimation processing unit,
For each of the plurality of incident signals, a vector calculation unit that determines a complex vector amount indicating a relative relationship between amplitude and phase distribution in the plurality of sensors,
Based on the physical position of the plurality of sensors and the complex vector amount determined by the vector calculation unit, an arrival direction calculation unit that determines an arrival direction of the plurality of incident signals;
The receiving apparatus according to claim 1, further comprising an output unit configured to output the signal separated by the separation unit in association with a signal indicating the direction of arrival calculated by the direction-of-arrival calculation unit. When the eigenvalue analysis is performed based on the observed signals for the number of signals selected by the signal number estimation processing unit, the eigenvectors obtained at the same time are separated into a signal eigenvector and a noise eigenvector, and the signal component of the captured signal is obtained by the signal eigenvector. 3. The receiving apparatus according to claim 1, further comprising: a signal separating unit that maps only one of the plurality of signals and separates a plurality of signals that exist simultaneously by independent component analysis of the mapped data. When the eigenvalue analysis is performed based on the observed signals for the number of signals selected by the signal number estimation processing unit, the eigenvectors obtained at the same time are separated into a signal eigenvector and a noise eigenvector, and the eigenvector and the response data of the sensor are separated. The receiving apparatus according to claim 1, further comprising an arrival direction detection unit that detects an arrival direction of a signal by performing the calculation of (1). The signal separation unit,
One of the separation method based on independence of probability distribution, gradient method, fixed point method, FastICA method, simultaneous diagonalization method of correlation function matrix, JADE method, factor analysis method, and method of maximizing information amount 4. The receiving device according to claim 3, wherein the receiving device is adopted as a processing algorithm. The arrival direction detection unit,
Processing algorithm for any one of the MUSIC method, cumulant-MUSIC method, WSF method, USF method, ML method, minimum dispersion method, BF method, DOSE method, Blind-MAXCOR method, Blind-MUSIC method, ESPRIT method, and VESPA method 4. The receiving device according to claim 3, wherein the receiving device is used as a receiving device.

Images