JP4676893B2 - 信号到来方向推定装置、信号到来方向推定方法、信号到来方向推定プログラム及び記録媒体 - Google Patents

Description

本発明は、信号処理の技術分野に属し、特に、それぞれの信号を１つずつ直接観測することはできず、複数の信号が重畳されて観測されるという状況において、全ての信号の到来方向を推定する技術に関する。

ここでは、まず観測信号のモデル化を行い、到来方向推定でよく使われる周波数領域表現を定義した上で、２つの従来技術について簡単に述べる。
［観測信号］
まず、観測信号の定式化を行う。すべての信号はあるサンプリング周波数ｆ_ｓでサンプリングされ、離散的に表現されるものとする。Ｎ個の信号が混合されてＭ個のセンサで観測されたとする。本発明では、信号の発生源からセンサまでの距離により信号が減衰・遅延し、また壁などにより信号が反射して伝送路歪みが発生しうる状況を扱う。このような状況での混合は、源信号ｓ_ｋ（ｔ），ｋ＝１，...，Ｎ、からセンサｘ_q（ｔ）へのインパルス応答ｈ_ｑｋ（ｌ），ｑ＝１，...，Ｍ，ｋ＝１，...，Ｎによる畳込み混合

となる。ここでは、源信号ｓ_ｋ（ｔ）やインパルス応答ｈ_ｑｋ（ｌ）についての情報を事前に得られない状況を考える。この状況で、観測信号ｘ_１（ｔ），...，ｘ_Ｍ（ｔ）のみを用いて源信号ｓ_１（ｔ），...，ｓ_Ｎ（ｔ）の到来方向を推定する。より詳細な到来方向の定義については後述する。

［周波数領域表現］
信号の到来方向の推定は周波数領域で行われることが多い。本発明でも周波数領域において各操作を行う。そのため、センサでの観測信号ｘ_ｑ（ｔ）にＬ点の短時間フーリエ変換を適用して周波数毎の時間系列

を求める。ここで、ｆは周波数であり、ｆ＝０，（１／Ｌ）ｆ_ｓ，...，（（Ｌ-１）／Ｌ）ｆ_ｓと離散化されている（ｆ_ｓはサンプリング周波数）。ｇ（ｌ）は、ハミング窓などの窓関数である。
式（１）で示されている時間領域での畳み込み混合は、周波数領域では、

と各周波数での単純混合に近似表現される。ここで、ｈ_ｑｋ（ｆ）は源信号からｓ_ｋからセンサｘ_ｑまでの周波数応答、ｓ_ｋ（ｆ，τ）は式（２）と同様の式に従って源信号ｓ_ｋ（ｔ）に短時間フーリエ変換を施したものである。式（３）をベクトルを用いて表記すると、

となる。ここで、Ｘ＝［Ｘ_１，...，Ｘ_Ｍ］^Ｔは観測信号ベクトル、ｈ_ｋ＝［ｈ_１ｋ，...，ｈ_Ｍｋ］^Ｔは源信号ｓ_ｋから各センサへの周波数応答をまとめたベクトルである。

［到来方向］
従来技術では、Ｍ個のセンサを直線状に並べるセンサアレイを用いるのが一般的であった。図８（ａ）に従来よく用いられるセンサアレイの構成を示す。センサｑの位置ｄ_ｑは既知であるとし、そこでの観測信号をｘ_ｑ（ｔ）とする。また、センサの並びと垂直な方向を９０°として、源信号ｓ_ｋ（ｔ）の到来方向を０°≦θ_ｋ≦１８０°とする。この到来方向θ_ｋを全ての信号源について推定する。
次に、従来の到来方向推定方法を２つ簡単に述べる。
［従来法１］ＭＵＳＩＣ法
広く用いられる方法としては、ＭＵＳＩＣ（ＭＵｔｉｐｌｅ ＳＩｇｎａｌ Ｃｌａｓｓｉｆｉｃａｔｉｏｎ）法がある（非特許文献１参照）。
これは、周波数ｆ毎に方向推定を行う。まず観測信号Ｘ（ｆ，τ）の相関行列Ｒ（ｆ）＝＜Ｘ（ｆ，τ）・Ｘ（ｆ，τ）^Ｈ＞_τを求め、これを固有値分解Ｒ（ｆ）＝ΣΛΣ^Ｈ（Σ＝［σ_１，...，σ_Ｍ］，Λ＝［λ_１，...，λ_Ｍ］）する。ここで、［＊］^Ｈは共役転置操作を、＜＊＞_τは時間フレームに関する平均操作を、σ_ｑはＲ（ｆ）の固有ベクトル、λ_ｑはσ_ｑに対応する固有値で、大きいものから順にソートされている。Ｎ番目までに大きい固有値に相当する固有ベクトルは信号部分空間を張り、Ｎ＋１〜Ｍ番目の固有値に相当する固有ベクトルは雑音部分空間を張る。

次に、

として（方向θを変数とするベクトル、ｃは信号の速度、ｄ_ｑは直線状に配置されたセンサアレイにおいてある１つのセンサを基準としたときの相対距離を表す。）、雑音部分空間を張る固有ベクトルσ_Ｎ＋１，...，σ_Ｍを用いて、

を計算すると、Ｌ（θ）はＮ個のピークを持つ。これらのＮ個のピークをとるθが、Ｎ個の信号の周波数ｆにおける到来方向θ_ｋ（ｆ）（ｋ＝１，...，Ｎ）となる。すべての周波数の到来方向推定結果の平均値を、信号到来方向とする。
ＭＵＳＩＣ法は雑音部分空間を張るＮ＋１〜Ｍ番目の固有値に相当する固有ベクトルが必要なため、センサ数Ｍ≧信号数Ｎ＋１である必要がある。すなわち、ＭＵＳＩＣ法は、センサ数が信号数より少ない場合（Ｍ＜Ｎ）には適用できない。

［従来法２］ＤＵＥＴ法
センサ数が信号数よりも少ない場合（Ｍ＜Ｎ）にも利用できる方法としてはＤＵＥＴ法がある（非特許文献２参照）。
この手法は信号のスパース性を仮定する。スパースとは、信号がほとんどの時刻τにおいて０であることを示す。信号のスパース性は、例えば音声信号で確認される。信号のスパース性を仮定することで、複数の信号が存在していても、各時間周波数ポイント（ｆ，τ）では互いに重なって観測される確率が低いことを仮定できる。よって、各時間周波数ポイント（ｆ，τ）の各センサにおける観測信号は、その時間周波数ポイント（ｆ，τ）でアクティブな１つの信号ｓ_ｋ（ｆ，τ）のみから成り、その信号ｋの方向情報を含むと仮定できる。このことから、ＤＵＥＴ法では、２つのセンサ（センサ１とセンサ２とする）のみを用いて、２つのセンサにおける観測信号の複数フレームτについて、全ての時間周波数（ｆ，τ）で次式で定義される位相差φ（ｆ，τ）

と、同じく次式で定義されるゲイン比α（ｆ，τ）

とを求め、この情報（φ（ｆ，τ），α（ｆ，τ））の全てのｆ、τにおける値のヒストグラムを作成する。そうすると、ヒストグラムはＮ個のピークをもつ。そのｋ番目のピーク（φ_ｋ，α_ｋ）のφ_ｋを用いると、

で求められる。ここでｄは、２つのセンサの間隔である。
この手法は、信号源やセンサの２次元・３次元配置に対応できない。すなわち、例えば、２つの信号源が２センサを結んだ線に対して線対称の位置にある場合に、２つの位置を推定することはできない。また、センサの高さについても推定することはできない。
大賀寿郎、山崎芳男、金田豊、「音響システムとディジタル処理」、電子情報通信学会、１９９５、ＩＳＢＮ ４−８８５５２−１２８−９ S.Rickard and F.Dietrich ‘DOA estimation of many W-disjoint orthogonal sources from two mixtures using DUET,’Proc. SSAP2000,pp.311-314.

従来法１では、直線配置アレイを、従来法２では２センサのみを使うことが一般的であり、０°≦θ_ｋ≦１８０°の範囲の方向推定しかできなかった。すなわち、２次元・３次元的に分散している信号源（０°≦θ_ｋ≦３６０°に分散していたり、高さが異なったりする）の到来方向推定ができなかった。
また、従来法１では、センサ数が信号数より多い場合にしか適用できず、センサ数が信号数より少ない場合には利用できなかった。しかし、実環境においてはセンサ数が信号数に比べて十分である保証はなく、信号数がセンサ数を超えることは容易に起こり得る。例えば、ＴＶ会議システムにおいて、装置のマイク数以上の話者が同時に発言することは起こり得る。
本発明は、以上の点に鑑みてなされたものであり、３個以上のセンサでかつ直線以外の配置（センサの２次元・３次元的な配置）のセンサアレイにも適用可能な到来方向推定方法を提供し、それにより２次元・３次元的に分散している信号源の到来方向推定を可能にする。さらに、センサ数が信号数より少ない場合にも適用できる信号源の到来方向推定を提供する。

本発明では上記課題を解決するために、まず、複数の信号源から発せられた信号を３つ以上のセンサで観測して、周波数領域変換手段において、周波数領域の周波数領域観測信号に変換する。そして、偏角正規化手段において、周波数領域観測信号から構成される観測信号ベクトルの要素の位相を、観測ベクトルの１つの要素の位相を基準として正規化し、第１正規化ベクトルを算出する。そして、ノルム正規化手段において、第１正規化ベクトルを、ノルムが所定の値になる第２正規化ベクトルに正規化する。そして、クラスタリング手段において、第２正規化ベクトルをクラスタリングし、各クラスタのセントロイドベクトルを算出する。その後、角度計算手段において、セントロイドベクトルと前記センサ位置情報とを用いて、各信号の到来方向を算出する。

本発明では、観測信号ベクトルＸ（ｆ，τ）に対して、ある正規化をほどこした特徴量を用いている。その正規化は、センサ数や配置形状によらず、どのような場合にも適用できる。さらに、それをＭ次元複素空間（Ｍはセンサ数）において、クラスタリングをして、そのセントロイドを到来方向推定に用いている。このため、本発明では、従来法とは異なり、２次元・３次元に配置された３個以上のセンサを用いることができるため、２次元・３次元に分布した信号源の到来方向を推定することができる。
また、本発明では、信号のスパース性を用いるので、センサ数が信号数より少ない場合にも利用できる。

以下、本発明の実施の形態を図面を参照して説明する。
まず、本形態の構成について説明する。
＜ハードウェア構成＞
図１は、本形態における信号到来方向推定装置１としてコンピュータに機能させる場合のハードウェア構成を例示したブロック図である。
図１に例示するように、この例の信号到来方向推定装置１は、ＣＰＵ（Central Processing Unit）１０、入力部２０、出力部３０、補助記憶装置４０、ＲＡＭ（Random Access Memory）５０、ＲＯＭ（Read Only Memory）６０及びバス７０を有している。

この例のＣＰＵ１０は、制御部１１、演算部１２及びレジスタ１３を有し、レジスタ１３に読み込まれた各種プログラムに従って様々な演算処理を実行する。また、この例の入力部２０は、データが入力される入力ポート、キーボード、マウス等であり、出力部３０は、データを出力する出力ポート、ディスプレイ等である。補助記憶装置４０は、例えば、ハードディスク、ＭＯ（Magneto-Optical disc）、半導体メモリ等であり、本形態の信号到来方向推定処理の各過程を実行するための信号到来方向推定プログラムを格納した信号到来方向推定プログラム領域４１及びセンサで観測された時間領域の観測信号等の各種データが格納されるデータ領域４２を有している。また、ＲＡＭ５０は、例えば、ＳＲＡＭ (Static Random Access Memory)、ＤＲＡＭ (Dynamic Random Access Memory)等であり、信号到来方向プログラムが書き込まれる信号到来方向推定プログラム領域５１及び各種データが書き込まれるデータ領域５２を有している。また、この例のバス７０は、ＣＰＵ１０、入力部２０、出力部３０、補助記憶装置４０、ＲＡＭ５０及びＲＯＭ６０を通信可能に接続している。

＜ハードウェアとソフトウェアとの協働＞
この例のＣＰＵ１０は、読み込まれたＯＳ（Operating System）プログラムに従い、補助記憶装置４０の信号到来方向推定プログラム領域４１に格納されている信号到来方向推定プログラムを、ＲＡＭ５０の信号到来方向推定プログラム領域５１に書き込む。同様にＣＰＵ１０は、補助記憶装置４０のデータ領域４２に格納されている時間領域の観測信号等の各種データをＲＡＭ５０のデータ領域５２に書き込む。さらに、ＣＰＵ１０は、この信号到来方向推定プログラムや各種データが書き込まれたＲＡＭ５０上のアドレスをレジスタ１３に格納する。そして、ＣＰＵ１０の制御部１１は、レジスタ１３に格納されたこれらのアドレスを順次読み出し、読み出したアドレスが示すＲＡＭ５０上の領域からプログラムやデータを読み出し、そのプログラムが示す演算を演算部１２に順次実行させ、その演算結果をレジスタ１３に格納していく。

図２は、このようにＣＰＵ１０に信号到来方向推定プログラムが読み込まれることにより構成される信号到来方向推定装置１の機能構成の例示である。なお、この図における実線の矢印は実際のデータの流れを示し、破線の矢印は理論的な情報の流れを示す。また、制御部１４０に出入りするデータの流れに対応する矢印は省略してある。
図２に例示するように、信号到来方向推定装置１は、メモリ１００、周波数領域変換部１１０、到来方向推定部１２０及び制御部１４０を有している。この到来方向推定部１２０は、観測信号ベクトル正規化部１２１、クラスタリング部１２６及び角度計算部１２７を有している。制御部１４０は一時メモリ１４１を有している。この観測信号ベクトル正規化部１２１は、偏角正規化部１２２及びノルム正規化部１２３を有している。角度計算部１２７は、逆行列計算部１２８及び積計算部１２９、ノルム正規化部１３０を有している。
ここでメモリ１００は補助記憶装置４０及びＲＡＭ５０に、一時メモリ１４１はレジスタ１３に相当する。また、時間周波数領域変換部１１０、到来方向推定部１２０及び制御部１４０は、ＣＰＵ１０にＯＳプログラムや信号到来方向推定プログラムが読み込まれ、ＣＰＵ１０がこれらを実行することにより機能として構成されるものである。

＜処理＞
次に、本形態の信号到来方向推定装置１の処理について説明する。なお、以下では、Ｎ個の源信号が混合され、Ｍ個のセンサで観測された状況を取り扱う。ここで、Ｎは２以上の整数、Ｍは３以上の整数であり、Ｍ個のセンサは２次元又は３次元に配置されているものとする。また、前処理において、各センサで観測された時間領域の観測信号ｘ_ｑ(t)(q=1,...,M)がメモリ１００の記憶領域１０１に格納され、実数α、信号の伝達速度ｃ、Ｍ以下の自然数から選択された基準値Ｑ（Ｍ個のセンサから選択された基準センサの添字番号）の各パラメータ、センサ位置ｄ_ｑ(ｑ=1,...,M)及びセントロイドの初期値ｃ^０ _ｋ（ｋ＝1,...,N）の各パラメータが記憶領域１０５に格納されているものとする。
ここで、センサ間の距離は、測定対象となる音源の波長の長さ以上あるように、センサｑが配置され、また、センサ位置ｄ_ｑ(ｑ=1,...,M)が記憶領域１０５に格納されているものとする。なお、センサ位置とは、センサの位置情報のことである。
図３は、本形態における信号到来方向推定装置１の処理の全体を説明するためのフローチャートである。以下、この図に沿って、本形態における信号到来方向推定装置１の処理を説明していく。

［処理の全体］
まず、周波数領域変換部１１０において、メモリ１００の記憶領域１０１から時間領域の観測信号x_q(t)を読み出し、これらを短時間離散フーリエ変換等によって、短時間ごとに周波数領域信号に変換された周波数ごとの時系列信号（「周波数領域観測信号」と呼ぶ。）X_q(f,τ)に変換し、メモリ１００の記憶領域１０２に格納する（ステップＳ１）。ここで、上記、X_q(f,τ)において、ｑ=1,...,M、f=0,f_s/L,…,f_s(L-1)/L、ｆ_ｓはサンプリング周波数、τは上記短時間変換をする時刻である。
次に、観測信号ベクトル正規化部１２１の偏角正規化部１２２において、メモリ１００の記憶領域１０２から周波数領域観測信号X_q(f,τ)を読み出す。周波数領域観測信号X_q(f,τ)を読み出した偏角正規化部１２２は、これらからなる観測信号ベクトルX(f,τ)=[X₁(f,τ),...,X_M(f,τ)]^Tを、観測信号ベクトルの１つの要素の位相を基準にして第１正規化ベクトルX’(f,τ)に正規化する（ステップＳ２）。生成された第１正規化ベクトルX’(f,τ)はメモリ１００の記憶領域１０３に格納される。なお、ステップＳ２の処理の詳細については後述する。
次に、観測信号ベクトル正規化部１２１のノルム正規化部１２３において、メモリ１００の記憶領域１０３から第１正規化ベクトルX’(f,τ)を読み込み、これらをノルムが所定の値（例えば１）になる第２正規化ベクトルX''(f,τ)に正規化する。そして、ノルム正規化部１２３は、生成した第２正規化ベクトルX''(f,τ)をメモリ１００の記憶領域１０４に格納する（ステップＳ３）。なお、この詳細については後述する。

次に、到来方向推定部１２０のクラスタリング部１２６において、メモリ１００の記憶領域１０４から第２正規化ベクトルX''(f,τ)を読み込み、これらをクラスタリングして、各クラスタのセントロイドベクトルｃ_ｋ(k=1,…,N)を求め、これをメモリ１００の記憶領域１０６に格納する（ステップＳ４）。なお、この詳細については後述する。
次に、角度計算部１２７の逆行列計算部１２８において、メモリ１００の記憶領域１０５からセンサ位置情報ｄ_ｑ及び基準値Ｑを読み込む。そして、逆行列計算部１２８は、センサ位置情報ｄ_ｑ及び基準値Ｑを用い、Ｑ番目のセンサ位置ｄ_Ｑを原点として他のセンサ位置を表した行列Ｄを計算し、その行列Ｄの一般化逆行列Ｄ^＋を求める。そして、逆行列計算部１２８は、計算した一般化逆行列Ｄ^＋をメモリ１００の記憶領域１０７に格納する。
次に、角度計算部１２７の積計算部１２９において、メモリ１００の記憶領域１０６からセントロイドベクトルｃ_ｋ及びメモリ１００の記憶領域１０７から一般化逆行列Ｄ^＋を読み込み、このセントロイドベクトルｃ_ｋ及び一般化逆行列Ｄ^＋から、信号到来方向ベクトルｖ_ｋを計算する。角度計算部１２７は、計算した信号到来方向ベクトルｖ_ｋをメモリ１００の記憶領域１０９に格納する（ステップＳ５）。なお、この詳細については後述する。
次に、各処理の例を詳細に説明する。

［観測信号ベクトル正規化部１２１の処理の詳細］
観測信号ベクトル正規化部１２１は、偏角正規化部１２２及びノルム正規化部１２３から構成され、全ての観測信号ベクトルX(f,τ)=[X₁(f,τ),...,X_M(f,τ)]^T（f=0,f_s/L,…,f_s(L-1)/L)を全ての時間周波数の要素について正規化する。この正規化は、クラスタリング部１２６でクラスタを形成させるために必要な手順であり、短時間フーリエ変換で用いる窓の位置や、信号原ｓ_ｋの位相や振幅、周波数特性ｈ_ｑｋ（ｆ）の周波数依存性を除去するためのものであり、クラスタリングを行うために必要である。正規化を適切に行わないと、クラスタが形成されない。
前述した通り、本形態における正規化は二段階からなる。第１の正規化は、偏角正規化部１２２において、観測信号ベクトルX(f,τ)の偏角を、短時間フーリエ変換で用いる窓の位置、信号原ｓ_ｋの位相及び周波数特性ｈ_ｑｋ（ｆ）の周波数依存性に依存しない第１正規化ベクトルX’(f,τ)に正規化するものである。第２の正規化は、ノルム正規化部１２３において、第１正規化ベクトルX’(f,τ)をノルムが所定の値（この例では１）をとる第２正規化ベクトルX''(f,τ)に正規化するものである。以下にこれらの正規化の詳細を説明する。

［偏角正規化部１２２における処理の詳細（ステップＳ２の処理の詳細）］
図４（ａ）は、図３に示したステップＳ２の処理の詳細を説明するためのフローチャートである。以下、この図に沿ってステップＳ２の処理の例を詳細に説明する。
まず、制御部１４０（図２）においてパラメータｑを１に設定し、これを一時メモリ１４１に格納する（ステップＳ１１）。次に、偏角正規化部１２２（図２）において、メモリ１００の記憶領域１０５から前述のパラメータα，ｃ，Ｑを読み込み、記憶領域１０２から各(f,τ)に対応する観測信号ベクトルX(f,τ)の要素X_q(f,τ)を読み込み、一時メモリ１４１からパラメータｑを読み込む。そして、偏角正規化部１２２は、

の演算を行い、当該演算結果を第１正規化ベクトルX'(f,τ)=[X₁'(f,τ),...,X_M'(f,τ)]^Tの各要素として、メモリ１００の記憶領域１０３に格納する（ステップＳ１２）。ここで、Ｘ_Ｑ（ｆ，τ）はある選択された基準センサＱの観測信号、αはある正数、ｃは信号の伝達速度、ａｒｇ［＊］は＊に対する偏角、ｊは虚数単位を意味する。αの選び方としては、４ｄ_max「ただし、ｄ_maxは、選択された基準センサＱと他のセンサとの距離の最大値」が最も良いが、他の数、例えば、距離単位として、ｍ（メートル）を用いた場合には１や２πなどでも良い。
次に、制御部１４０において、一時メモリ１４１に格納されたパラメータｑがｑ＝Ｍを満たすか否かを判断する（ステップＳ１３）。ここでｑ＝Ｍでなければ、制御部１４０は、ｑ＋１の演算結果を新たなパラメータｑの値とし、これを一時メモリ１４１に格納し（ステップＳ１４）、処理をステップＳ１２へ戻す。一方、ｑ＝Ｍであれば、制御部１４０はステップＳ２の処理を終了させ、以下のステップＳ３の処理を実行させる。

［ノルム正規化部１２３における処理の詳細（ステップＳ３の詳細）］
図４（ｂ）は、図３に示したステップＳ３の処理の詳細を説明するためのフローチャートである。以下、この図に沿ってステップＳ３の処理の詳細を説明する。
ノルム正規化部１２３（図２）において、メモリ１００の記憶領域１０３から、各(f,τ)に対応する第１正規化ベクトルX'(f,τ)=[X₁'(f,τ),...,X_M'(f,τ)]^Tを読み込む。そして、ノルム正規化部１２３は、

の演算を行ってそれらのノルム‖X'(f,τ)‖を求め、ノルム‖X'(f,τ)‖を一時メモリ１４１に格納する（ステップＳ２１）。
次に、ノルム正規化部１２３は、メモリ１００の記憶領域１０３から各(f,τ)に対応する第１正規化ベクトルX'(f,τ)及び、一時メモリ１４１からノルム‖X'(f,τ)‖を読み出し、

の演算を行って第２正規化ベクトルX''(f,τ)を求める（ステップＳ２２）。求められた第２正規化ベクトルX''(f,τ)は、メモリ１００の記憶領域１０４に格納される。これによりステップＳ３の処理が終了する。
このように生成された第２正規化ベクトルX''(f,τ)は、短時間フーリエ変換で用いる窓の位置や信号源Ｓ_ｋの位相や振幅、周波数特性ｈ_ｑｋ（ｆ）の周波数特性には依存せず、信号源の位置のみに依存するベクトルとなる。その結果、この第２正規化ベクトルX''(f,τ)はクラスタを形成することになる。以下にこの理由を説明する。

［第２正規化ベクトルX''(f,τ)がクラスタを形成する理由］
まず、信号源ｓ_ｋからセンサｑまでの周波数領域でのインパルス応答、すなわち、周波数応答ｈ_ｑｋ（ｆ）を以下のように取り扱う。前述した式（３）に示される畳み込み混合モデルを平面波（遠距離場）混合モデル

を用いて近似する。すなわち、ここではセンサ間における減衰差は考えず、音源ｋから原点までの遅延と、音源ｋからセンサｑまでの遅延の差ｃ^−１ｄ_ｑ ^Ｔｑ_ｋのみを考える。ここで、ｄ_ｑ ^Ｔはセンサｑの位置を表す３次元ベクトルである。また、ｖ_ｋは、原点から見た信号源ｓ_ｋの方向を表す大きさ１の３次元ベクトル（到来方向ベクトル）であり、図８（ｂ）に示される方位θ_ｋと仰角φ_ｋを用いて、

と書くことができる。
また、信号のスパース性を仮定する。すなわち、各時間周波数ポイント（ｆ，τ）の各センサにおける観測信号は、その時間周波数ポイント（ｆ，τ）でアクティブな１つの信号ｓ_ｋ（ｆ，τ）のみからなり、その信号ｋの方向情報を含むと仮定する。このとき、各時間周波数(f,τ)において以下の関係が成り立つ。

式（１３）と（１５）を式（６）に代入して偏角を取ると、

となる。この式から分かるように、第１正規化ベクトルX'(f,τ)の偏角は、信号源の振幅や位相に依らなくなることが分かる。これが偏角正規化部１２２の役割である。
さらに、式（１０）〜式（１３）より、第２正規化ベクトルX'' (f,τ)は、

又は

となる。ここで、ｃｏｓΘ_ｋ ^ｑＱは、センサｑとＱを結ぶ軸に対する信号ｓ_ｋ（ｔ）の到来方向ベクトルであり、図９（ａ）において太いベクトルで示されているものである。上記式（１７）及び式（１８）より、第２正規化ベクトルX'' (f,τ)は、周波数応答ｈ_ｑｋ（ｆ）が持っていた周波数依存性が無くなっていることがわかる。さらに、第２正規化ベクトルX''(f,τ)は、信号源kとセンサqの位置のみに依存することも分かる。このため、第２正規化ベクトルX''(f,τ)をクラスタリングすると、同じ信号源に対応するものが、その位置によってクラスタを形成するのである。（［観測信号ベクトル正規化部１２１の処理の詳細・偏角正規化部１２２・ノルム正規化部１２３の処理の詳細］の説明終わり）。

［クラスタリング部１２６の処理の詳細（ステップＳ４の詳細）］
前述のようにクラスタリング部１２６は、メモリ１００の記憶領域１０４から第２正規化ベクトルX''(f,τ)を読み込み、これらをクラスタリングしてＮ個のクラスタＣ_１，...Ｃ_Ｎを生成する。すなわち、Ｍ次元複素ベクトルである正規化された観測信号ベクトルX''(f,τ)をＭ次元複素空間で直接クラスタリングする。このクラスタリングは、例えば、k-means法を用いて効果的に行うことができる（［参考文献］R.O.Duda,P.E.Hart,and D.G.Stock,Pattern Classification, Wiley Interscience, ２nd edition,2000.）。以下、本発明によるクラスタリング手順を説明する。
１．クラスタのセントロイドの初期値ｃ^ｉ _ｋ（ｉ＝０、ｋ＝１，...，Ｎ）をメモリ１００の記憶領域１０５から読み込む。セントロイドの初期値ｃ^ｉ _ｋは、第２正規化ベクトルX''(f,τ)と同じ次元のベクトル（Ｍ次元複素ベクトル）である。なお、セントロイドの初期値ｃ^０ _ｋの選び方については後述する。
２．ｉ＋１を新たなｉとする。
３．すべての時間周波数(f,τ)における第２正規化ベクトルX''(f,τ)を、最も近いセントロイドｃ^ｉ−１ _ｋで代表されるクラスタＣ_ｋに割り当てる。すなわち、各第２正規化ベクトルX''(f,τ)に対して、‖X''(f,τ)−ｃ^ｉ−１ _ｋ‖が最も小さくなるようにｋを選ぶ。
４．各クラスタＣ_ｋに割りあてられた第２正規化ベクトルX''(f,τ)の平均値を計算し、そのノルムを１にすることでセントロイドを更新する。すなわち、各クラスタＣ_ｋに割りあてられた第２正規化ベクトルX''(f,τ)に対して、

の演算を行うことにより、セントロイドを更新する。ここで、Ｅ［・］_ｋは、クラスタＣ_ｋのメンバに対する平均操作を表す。
５．セントロイドｃ^ｉ _ｋが収束するまで、手順２−５を繰り返す。最後に収束したセントロイドを、ｃ_ｋ（ｋ＝１，...，Ｎ）として、メモリ１００の記憶領域１０６に格納する。
以上が、本発明によるクラスタリング手順である。

次に、セントロイドの初期値の選び方の例を説明する。
《初期値設定方法１》
第２正規化ベクトルX''(f,τ)の中からランダムにＮ個のベクトルを選び、それをセントロイドの初期値ｃ^０ _ｋ（ｋ＝１，...，Ｎ）とする。
《初期値設定方法２》
セントロイドは、後述するように、式（２３）のように書けるので、これを用いる。具体的には、式（２３）において、センサ位置ｄ_ｑ（ｑ＝１，...，Ｍ）は、メモリ１００の記憶領域１０５において保持されている値を、ｖ_ｋとしては、式（１４）のように考え、方位θ_ｋと仰角φ_ｋ（ｋ＝１，...，Ｎ）を適当に与える。ここで、方位θ_ｋと仰角φ_ｋは初期値であるので、適当な値で良い。例えば、

とすると、空間的に散らばった初期値が得られる（［クラスタリング部１２６の処理の詳細（ステップＳ４の詳細）］の説明終わり）。

［角度計算部１２７・逆行列計算部１２８・積計算部１２９・ノルム正規化部１３０における処理の詳細（ステップＳ５の処理の詳細）］
図５は、図３に示したステップＳ５の処理の例の詳細を説明するためのフローチャートである。以下、この図に沿ってステップＳ５の処理の例の詳細を説明する。
角度計算部１２７において、到来方向ベクトルｖ_ｋを計算する。この例では、３個以上のセンサが２次元又は３次元に任意に配置されている。まず、角度計算部１２７内の逆行列計算部１２８は、メモリ１００の記憶領域１０５から、センサ位置情報ｄ_ｑ（ｑ＝１，...，Ｍ）及び基準値Ｑを読み込み、Ｑ番目のセンサ位置ｄ_Ｑを原点にして正規化した行列Ｄを計算する。そして、この行列Ｄの一般化逆行列Ｄ^＋を計算する。ここで、一般化逆行列Ｄ^＋として、例えばMoore-Penroseの一般化逆行列を用いることができる（ステップＳ３１）。
なお、後述のように、正立方体のある頂点と、それと隣り合う３つの頂点にセンサを配置した場合には、逆行列計算部１２８の処理を省略することができる。

次に、角度計算部１２７内の積計算部１２９において、メモリ１００の記憶領域１０６からセントロイドベクトルｃ_ｋ及び記憶領域１０７から一般化逆行列Ｄ^＋を読み出し、

の演算をする。そして、積計算部１２９は、この演算により求められたｖ_ｋ’をメモリ１００の記憶領域１０８に格納する（ステップＳ３２）。理想的には、この演算により、到来方向ベクトルが得られる。なお、その理由については、後述する。また、後述のように、正立方体のある頂点と、それと隣り合う３つの頂点にセンサを配置した場合には、積計算部１２９の処理を省略することができる。

しかし、上記式（２１）の右辺は計算誤差を含むことが多いので、角度計算部１２７内のノルム正規化部１３０において、第１到来方向ベクトルｖ_ｋ’のノルムを１にする正規化を行う。すなわち、ノルム正規化部１３０は、メモリ１００の記憶領域１０８から第１到来方向ベクトルｖ_ｋ’を読み出し、

の演算を行い、第２到来方向ベクトルｖ_ｋを求め、これを信号到来方向ベクトルとする（ステップＳ３３）。

《到来方向ベクトルが求まる理由》
上述した式（１７）と式（１８）により、第２正規化ベクトルの各要素が、角度情報である到来方向ベクトルｖ_ｋもしくは、２つのセンサｑ、Ｑに対する到来方向Θ_ｋ ^ｑＱを含むことが分かる。また、式（１９）から分かるように、クラスタのセントロイドｃ_ｋは、第２正規化ベクトルの各クラスタの平均値を与えることが分かる。このため、式（１７）と式（１８）により、クラスタのセントロイドｃ_ｋのｑ番目の要素｛ｃ_ｋ｝_ｑは、次のように信号の角度情報を持つことが分かる。

式（２３）の偏角を考えると、

となり、これを全てのセンサペア（ｑ＝１，...，Ｍ）について考えると、

が得られる。この式を、ｖ_ｋについて整理すると、

となる。このため、一般化逆行列Ｄ^＋とセントロイドベクトルｃ_ｋの偏角の積を計算することにより、到来方向ベクトルが得られるのである（［角度計算部１２７・逆行列計算部１２８・積計算部１２９・ノルム正規化部１３０における処理の詳細（ステップＳ５の処理の詳細）］の説明終わり）。

［実験結果］
上記実施例の効果を示すために、到来方向推定実験を行った。実験条件は図９（ｂ）に示す通りである。長辺が８８０ｃｍ、短辺が３７５ｃｍ、高さが２４０ｃｍ、残響は１２０ｍｓの室内において、底面の長辺から２００ｃｍ、短辺から２８２ｃｍの位置にセンサを配置した。センサ位置を中心に半径５０ｃｍの円上に３つの音源を、半径８０ｃｍの円上に２つの音源を配した。長辺と平行軸をｘ、短辺と平行軸をｙとし、図中［Setup１］として示すように、３つのセンサをｙ軸に２個、ｘ軸に１個、辺の長さ４ｃｍの正三角形の頂点につまり２次元に配した場合と、［Setup２］として示すように、４つのセンサを原点に１個、原点からそれぞれ４ｃｍ離してｘ軸、ｙ軸、ｚ軸上にそれぞれ１個配した３次元配置とした場合の２種類の実験を行う。音源としては５秒間の英語の音声を用い、またセンサとしてはマイクロホンを用いた。［Setup1］および［ Setup2］の結果をそれぞれ図１０（１）の表１及び表２に示す。ここでは、MUSIC法にて、予めそれぞれ1つずつ音源を鳴らして音源到来方向を推定しそれを正解とした。表１ではMUSIC(each)として示している。(注意：５音源が同時に存在する場合の推定ではない。MUSIC法はその場合適用不可。)。また、［Setup1］では、全ての仰角φ_ｋ＝０であるため、方位θ_kの値のみを示している。
図１０（１）の表１及び表２より、この実施例の方法によれば、いずれの音源に対してもＭＵＳＩＣ（ｅａｃｈ）の測定結果とよく一致している。この発明法を用いることで信号源数＞センサ数の場合でも、精度よく到来方向推定ができることが分かる、従来のMUSIC法ではこのような信号源数＞マイク数の場合の到来方向推定は不可能であった。
［Setup１］(図１０（１）の表1)では、センサのある平面上の全ての方向の方向推定ができている。また、［Setup２］(図１０（１）の表2)では、仰角φ_ｋまで含めた全ての方向の方向推定ができている。これは従来のDUET法では不可能であった。

さらに、上記実施例では、従来最もよく用いられる方法であるMUSIC法よりも優れた角度分解能を持つことが確かめられた。図１０（２）において(a)(c)はMUSICスペクトルを、(b)(d)は上記実施法における到来方向推定を全ての時間周波数（ｆ，τ）について行った場合に得られる到来方向のヒストグラムを示している。双方とも、ピーク位置が音源方向を示す。
まず、図１０（２）(a)(b)は、２つの信号の角度差が大きい(θ_１＝117°,θ₂=217°)場合のMUSICスペクトルと上記実施例における方法のヒストグラムを示している。この場合、両者ともピークは２つであり２方向を推定できる。
一方、図１０（２）(c)(d)は、２信号の角度差が小さい(θ₁=132°,θ₂=154°)場合のMUSICスペクトルと上記実施法のヒストグラムを示している。この場合、MUSIC法では、２信号が存在するにもかかわらず、スペクトルのピークが1つしか確認できない。一方、上記実施法では、ピークが２つ存在し、２信号の角度差が小さい場合にも正しく方向を推定できることが分かる。なお、この場合の到来方向推定値は、θ₁=128°,θ₂=156°であった。

上記実施例の効果をまとめると以下のようになる。
・ 信号のスパース性に基づいた到来方向推定であるため、信号源数＞センサ数の場合にも信号の到来方向推定が可能である.これは従来のMUSIC法ではできなかった。
（２）全センサの配置情報を知る必要のない正規化による特徴量を用いたクラスタリングに基づく方法であるため、不規則なセンサ配置の採用が容易である。
（３）（２）に述べたように,不規則なセンサ配置を採用できるため、センサの２次元・３次元配置が可能となり、信号が２次元・３次元に配置している場合でも全方位の到来方向推定が可能である。これは従来の２センサによるDUET法や線形アレイによる方法ではできなかった。
（４）信号のスパース性に基づいた到来方向推定であるため、従来のMUSIC法に比べて、優れた角度分解能を持つことが確認された。(なお、角度分解能は、信号のスパース性の程度に依存する。)

[変形例等]
《変形例１》
変形例１では、偏角正規化部１２２において、

とする代わりに、観測信号ベクトルの大きさ|X_q(f,τ)|を用いずに、位相を表す項のみを用いる。すなわち、偏角正規化部１２２において、

の演算を行うことにより、第１正規化ベクトルX'(f,τ)を求めることもできる。

《変形例２》
変形例２では、偏角正規化部１２２において、式（１０）の代わりに、式（１０）の指数の肩の項のみを用いて、

の演算を行うことにより、第１正規化ベクトルX'(f,τ)を求めることもできる。

また、変形例２は、図６に示す角度計算部１２７ａを有する。すなわち、変形例２においては、図６に示すように、角度計算部１２７ａの積計算部１２９ａにおいて、式（２１）の演算ではなく、

の演算を行い、第１到来方向ベクトルｖ_ｋ’を求める。これに続けてノルム正規化部１３０において、

の演算を行い、第２到来方向ベクトルｖ_ｋを推定する。

《変形例３》
変形例３では、偏角正規化部１２２において、式（１０）の代わりに、

の演算を行うことにより、第１正規化ベクトルX'(f,τ)を求める。
観測信号ベクトル正規化部１２１内の偏角正規化部１２２において、式（２９）ではなく式（３２）の演算を行う点のみが、変形例２と異なり、その他の機能構成・処理は変形例２の場合と同一である。

《変形例４》
変形例４においては、正立方体のある頂点と、それに隣り合う３つの頂点に、すなわち、例えば、図９（ｃ）に示すように、［ｄ_１，ｄ_２，ｄ_３，ｄ_４］^Ｔ＝Ａ［（０，０，０）^Ｔ，（１，０，０）^Ｔ，（０，１，０）^Ｔ，（０，０，１）^Ｔ］^Ｔ（Ａ：ある定数）のようにセンサを配置する。
この場合、原点（０，０，０）^Ｔを基準Ｑに他のセンサ位置を正規化したベクトルＤの一般化逆行列Ｄ^＋は、

となる。この場合、任意のＭ次元ベクトルａに対して、Ｄ^＋ａ＝ａ／Ａとなるため、角度計算部１２７における逆行列演算部１２８と積計算部１２９を省略することができる。
変形例４では、観測信号ベクトル正規化部１２１において式（１０）の後に式（１２）の正規化を行った後、クラスタリング部１２６でクラスタリングされた結果得られたセントロイドベクトルｃ_ｋについて、図７（ａ）に示す角度計算部１２７ｂの偏角計算部１３１ａにおいて、

の演算をして、偏角ａｒｇ［ｃ_ｋ］を求めて、その結果を第１到来方向ベクトルｖ_ｋ’とする。そして、ノルム正規化部１３０において、

の演算を行い、第２到来方向ベクトルｖ_ｋを求める。
なお、上記の変形例４では、式（１０）を用いて観測信号ベクトルを正規化したが、式（２８）を用いて正規化しても良い。

《変形例５》
変形例５においては、変形例４の場合と同様に、正立方体のある頂点と、それに隣り合う３つの頂点に、すなわち、［ｄ_１，ｄ_２，ｄ_３，ｄ_４］^Ｔ＝Ａ［（０，０，０）^Ｔ，（１，０，０）^Ｔ，（０，１，０）^Ｔ，（０，０，１）^Ｔ］^Ｔ（Ａ：ある定数）のようにセンサを配置する。
その後、観測信号ベクトル正規化部１２１において式（２９）又は式（３２）の正規化を行った後、式（１２）のノルム正規化を行い、クラスタリング部１２６でクラスタリングされた結果得られたセントロイドベクトルｃ_ｋについて、図７（ｂ）に示す角度計算部１２７ｃの偏角計算部１３１ｂにおいて、

の演算を行い、その結果を第１到来方向ベクトルｖ_ｋ’とする。そして、ノルム正規化部１３０において、式（３５）の演算を行い、第２到来方向ベクトルｖ_ｋを求める。

《変形例６》
ここまでは、全ての時間周波数（ｆ，τ）における正規化された観測信号ベクトルをクラスタリング部１２６におけるクラスタリングの対象としていたが、このクラスタリングを周波数ｆ毎に行い、その後全ての周波数での到来方向推定の結果を統合しても良い。これには例えば、全ての周波数での到来方向推定結果を再度クラスタリングし、そのセントロイドを最終的な到来方向推定結果として用いれば良い。

《その他の変形例》
上述では、センサ位置情報を表す行列Ｄの一般化逆行列Ｄ^＋を用いる場合と、これを用いない場合としたが、例えば角度計算部１２７にセンサ配置法判定部１３２を設け、センサの配置状態が正立方体の頂点に配置されたものか否かを判定し、その判定結果に応じて、その後の処理を、前記一般化逆行列を求める場合と、これを用いない場合とに区別するようにしても良い。センサ配置法判定部１３２としては、例えば、予めセンサ配置法が知られ、正立方体頂点位置に配されているか否かを示す情報が入力部２０（図１）から入力され、メモリ１００のいずれかの領域に格納され、この情報からセンサ配置法を判断する。あるいは、各センサ位置情報ｄ_ｑ（ｑ＝１，...，Ｍ）から、これらが正立方体頂点位置に配された関係にあるかを演算により判断しても良い。
また、本発明は上述の実施の形態に限定されるものではない。例えば、上述の実施の形態では、ノルム正規化部１２２においてノルムが１になるような正規化を行ったが、ノルムが１以外の規定値となるように正規化を行ってもよい。

また、上述の構成をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよいが、具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ
−Ｒ（Recordable）／ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto-Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。

また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
また、上述した実施形態とは別の実行形態として、コンピュータが可搬型記録媒体から直接このプログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。

また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。

本技術により、例えば、様々な妨害信号・干渉音が存在する環境下においても、それぞれの信号の到来方向を精度よく推定することが可能となる。
例えば、オーディオ分野に応用として、ロボット聴覚やテレビ会議システムなどがある。これらにおいては、複数の音声や音楽などが混在する状況が起こりやすく、また同時に入力マイクロホンと話者が離れた位置になりやすいため、マイクロホンが複数の音を同時に収音してしまう。このような状況でも、それぞれの音の方向を推定できることで、ロボットにおいてはロボットカメラを向けることが、テレビ会議システムにおいては話者の方向情報をも含めた情報の通信が可能となり、ステレオ受聴や会議システムのカメラを話者へ向けることが可能となる。

本発明の実施形態における信号到来方向推定装置のハードウェア構成を例示したブロック図。 本発明の実施形態における信号到来方向推定装置のブロック図。 本発明の実施形態における信号到来方向推定装置の処理の全体を説明するためのフローチャート。 （ａ）は、図３中のステップＳ２の詳細例を説明するためのフローチャート、（ｂ）は、図３中のステップＳ３の詳細例を説明するためのフローチャート。 図３中のステップＳ５の詳細例を説明するためのフローチャート。 変形例２及び変形例３における角度計算部の詳細例を示したブロック図。 （ａ）は、変形例４における角度計算部の詳細例を示したブロック図、（ｂ）は、変形例５における角度計算部の詳細例を示したブロック図。 （ａ）は、従来の到来方向推定問題について説明する図、（ｂ）は、本発明における信号の到来方向を定義する図。 （ａ）は、本発明法の信号の到来方向推定問題を定義する図、（ｂ）は、実験条件を示した図、（ｃ）は、変形例４及び変形例５におけるセンサ位置を表した図。 （１）及び（２）は、実験結果を示した図。

符号の説明

１ 信号到来方向推定装置

Claims (14)

1. 複数の信号源から発せられた信号を３つ以上のセンサで観測し、観測された観測信号を用いて個々の信号の到来方向を推定する装置であって、
   前記観測信号を周波数領域の周波数領域観測信号に変換する周波数領域変換手段と、
   前記周波数領域観測信号から構成される観測信号ベクトルの要素の位相を、前記観測信号ベクトルの１つの要素の位相を基準として正規化し、第１正規化ベクトルを算出する偏角正規化手段と、
   前記第１正規化ベクトルを、ノルムが所定の値になる第２正規化ベクトルに正規化するノルム正規化手段と、
   前記第２正規化ベクトルをクラスタリングし、各クラスタのセントロイドベクトルを算出するクラスタリング手段と、
   前記各セントロイドベクトルと前記センサの位置情報とを用いて、各信号の第１到来方向ベクトルを算出する角度計算手段と、
   前記第１到来方向ベクトルを、ノルムが１になる第２到来方向ベクトルに正規化し、当該第２到来方向ベクトルを到来方向ベクトルとするノルム正規化手段と、
   を有することを特徴とする信号到来方向推定装置。
2. 請求項１記載の信号到来方向推定装置であって、
   前記偏角正規化手段が、
   Ｍを前記センサの数とし、ｑ＝１，...，Ｍとし、前記観測信号ベクトルの各要素をX_q(f,τ)とし、arg（*）を*に対する偏角とし、ｊを虚数単位とし、ｃを信号の伝達速度とし、ＱをＭ以下の自然数から選択された値とし、αを正数とし、ｆを周波数とし、τを離散時間とした場合における、
   

   

   の演算を行い、
   当該演算結果を前記第１正規化ベクトルの各要素とする手段である、
   ことを特徴とする信号到来方向推定装置。
3. 請求項１記載の信号到来方向推定装置であって、
   前記偏角正規化手段が、
   Ｍを前記センサの数とし、ｑ＝１，...，Ｍとし、前記観測信号ベクトルの各要素をX_q(f,τ)とし、arg（*）を*に対する偏角とし、ｊを虚数単位とし、ｃを信号の伝達速度とし、ＱをＭ以下の自然数から選択された値とし、αを正数とし、ｆを周波数とし、τを離散時間とした場合における、
   

   

   の演算を行い、
   当該演算結果を前記第１正規化ベクトルの各要素とする手段である、
   ことを特徴とする信号到来方向推定装置。
4. 請求項１記載の信号到来方向推定装置であって、
   前記偏角正規化手段が、
   Ｍを前記センサの数とし、ｑ＝１，...，Ｍとし、前記観測信号ベクトルの各要素をX_q(f,τ)とし、arg（*）を*に対する偏角とし、ｊを虚数単位とし、ｃを信号の伝達速度とし、ＱをＭ以下の自然数から選択された値とし、αを正数とし、ｆを周波数とし、τを離散時間とした場合における、
   

   

   の演算を行い、
   当該演算結果を前記第１正規化ベクトルの各要素とする手段である、
   ことを特徴とする信号到来方向推定装置。
5. 請求項１記載の信号到来方向推定装置であって、
   前記偏角正規化手段が、
   Ｍを前記センサの数とし、ｑ＝１，...，Ｍとし、前記観測信号ベクトルの各要素をX_q(f,τ)とし、arg（*）を*に対する偏角とし、ＱをＭ以下の自然数から選択された値とし、ｆを周波数とし、τを離散時間とした場合における、
   

   

   の演算を行い、
   当該演算結果を前記第１正規化ベクトルの各要素とする手段である、
   ことを特徴とする信号到来方向推定装置。
6. 請求項２又は３記載の信号到来方向推定装置であって、
   前記角度計算手段が、
   前記各セントロイドベクトルの偏角と前記センサの位置情報を表す行列の一般化逆行列との積を計算して各信号の第１到来方向ベクトルを求める手段である、
   ことを特徴とする信号到来方向推定装置。
7. 請求項４又は５記載の信号到来方向推定装置であって、
   前記角度計算手段が、
   前記各セントロイドベクトルと前記センサの位置情報を表す行列の一般化逆行列との積を計算して各信号の第１到来方向ベクトルを求める手段である、
   ことを特徴とする信号到来方向推定装置。
8. 請求項２又は３記載の信号到来方向推定装置であって、
   （１）前記周波数領域変換手段が、
   正立方体のある頂点とそれにとなりあう３つの頂点に配置された複数のセンサにおいて観測された観測信号を周波数領域の周波数領域観測信号に変換する周波数領域変換手段であり、
   （２）前記角度計算手段が、
   前記各セントロイドベクトルの偏角を計算して各信号の第１到来方向ベクトルを求める手段である、
   ことを特徴とする信号到来方向推定装置。
9. 請求項４又は５記載の信号到来方向推定装置であって、
   （１）前記周波数領域変換手段が、
   正立方体のある頂点とそれにとなりあう３つの頂点に配置された複数のセンサにおいて観測された観測信号を周波数領域の周波数領域観測信号に変換する周波数領域変換手段であり、
   （２）前記角度計算手段が、
   前記各セントロイドベクトルを第１到来方向ベクトルとする手段である、
   ことを特徴とする信号到来方向推定装置。
10. 請求項２，３，６，８のいずれかに記載の信号到来方向推定装置であって、
    前記クラスタリング手段が、
    Ｍを前記センサの数とし、ｄを正数とし、ｄ_ｑをセンサｑの位置情報とし、ｖ_ｋ＝［ｃｏｓθ_ｋｃｏｓφ_ｋ，ｓｉｎθ_ｋｃｏｓφ_ｋ，ｓｉｎφ_ｋ］とし、θ_ｋを０≦θ_ｋ≦３６０°、φ_ｋを０≦φ_ｋ≦１８０°の範囲から選んだ値とした場合における、
    

    

    をｋ＝１，...，Ｎ，ｑ＝１，...，Ｍについて計算し、
    当該計算結果であるｃ_ｋをセントロイドベクトルの初期値としてクラスタリングを行う手段である、
    ことを特徴とする信号到来方向推定装置。
11. 複数の信号源から発せられた信号を３つ以上のセンサで観測し、観測された観測信号を用いて個々の信号の到来方向を推定する方法であって、
    周波数領域変換手段において、前記観測信号を周波数領域の周波数領域観測信号に変換する過程と、
    偏角正規化手段において、前記周波数領域観測信号から構成される観測信号ベクトルの要素の位相を、前記観測ベクトルの１つの要素の位相を基準として正規化し、第１正規化ベクトルを算出する過程と、
    ノルム正規化手段において、前記第１正規化ベクトルを、ノルムが１になる第２正規化ベクトルに正規化する過程と、
    クラスタリング手段において、前記第２正規化ベクトルをクラスタリングし、各クラスタのセントロイドベクトルを算出する過程と、
    角度計算手段において、前記各セントロイドベクトルと前記センサの位置情報とを用いて、各信号の第１到来方向ベクトルを算出する過程と、
    ノルム正規化手段において、前記第１到来方向ベクトルを、ノルムが１になる第２到来方向ベクトルに正規化し、当該第２到来方向ベクトルを到来方向ベクトルとする過程と、
    を有することを特徴とする信号到来方向推定方法。
12. 請求項１１記載の信号到来方向推定方法であって、
    前記周波数領域変換手段が、
    正立方体のある頂点とそれにとなりあう３つの頂点に配置された複数のセンサにおいて観測された観測信号を周波数領域の周波数領域観測信号に変換する過程である、
    ことを特徴とする信号到来方向推定方法。
13. 請求項１から１０の何れかに記載の信号到来方向推定装置としてコンピュータを機能させるための信号到来方向推定プログラム。
14. 請求項１３記載の信号到来方向推定プログラムを格納したコンピュータ読み取り可能な記録媒体。

Images