JPWO2004079388A1 - 位置情報推定装置、その方法、及びプログラム - Google Patents

Abstract

２次元配置された複数のセンサよりの観測信号ｘ１（ｔ）〜ｘＪ（ｔ）を短時間フーリエ変換し、これらから信号Ｘ１（ω１）〜Ｘ１（ωＮ），〜ＸＪ（ω１）〜ＸＪ（ωＮ）、独立成分分析法により分離行列Ｗ（ω１）〜Ｗ（ωＮ）を生成し、その逆行列Ｈ（ω１）〜Ｈ（ωＮ）を計算し、各ωｎ（ｎ＝１，…，Ｎ）ごとに、Ｈ（ωｎ）の各列の要素対Ｈｊｉ（ωｎ）及びＨｊ’ｉ（ωｎ）について、角度θ＾ｉ，ｊｊ’（ωｎ）＝ｃｏｓ−１（ａｒｇ（Ｈｊｉ（ωｎ）／Ｈｊ’ｉ（ωｎ））／（ωｎｃ−１‖ｄｊ−ｄｊ’‖））を計算し、ａｒｇ（α）はαの偏角、ｃは信号の伝搬速度、‖ｄｊ−ｄｊ’‖はセンサｊとｊ’の間隔である。Ｈ（ω１）〜Ｈ（ωＮ）において各列から得られるθ＾ｉ，ｊｊ’（ωｎ）が昇順になるように、列の入れ替えを行い、列の入れ替えができない列について、‖ｑｉ−ｄｊ’‖／‖ｑｉ−ｄｊ‖＝｜Ｈｊｉ（ωｎ）／Ｈｊ’ｉ（ωｎ）｜＝ＤＲｉ，ｊｊ’（ωｎ）をｑｉについて解いてＲｉ，ｊｊ’（ωｎ）＝‖ＤＲｉ，ｊｊ’（ωｎ）（ｄｊ−ｄｊ’）／（ＤＲ２ｉ，ｊｊ’（ωｎ）−１）‖を計算し、このＲｉ，ｊｊ’（ωｎ）が昇順になるように当該Ｈ（ωｎ）の列を入れ替え、このＨ（ωｎ）を用いてＷ（ω）のパーミュテーションを解決する。

Description

この発明は音源や電波源などの信号源の複数から放射され、互いに混合された信号を複数のセンサで観測して、各信号源の位置情報、正しくは位置を表わすパラメータの少くとも１つを含む情報を推定し、信号の到来方向の検出や、信号源ごとの信号に分離復元に適用される装置、その方法、及びプログラムに関する。

複数の信号源からの各信号が空間内で混合されて複数のセンサに到来し、これらセンサで観測された到来信号から、各源信号の到来方向の推定や各源信号を分離することを、その源信号の混合系の情報を知らずに、独立成分分析（Ｉｎｄｅｐｅｎｄｅｎｔ Ｃｏｍｐｏｎｅｎｔ Ａｎａｌｙｓｉｓ、以後、ＩＣＡと記述する場合もある）により行うことが提案されている。空間内での前記混合は信号源からセンサまでの到達遅延及び減衰度が直接波と障害物などによる複数の反射波で異なるため、ある信号は複数の時間遅れを持って混合された畳み込み混合となる。時間領域で直接分離フィルタを求めるＩＣＡは最終的な解への収束が非常に遅いため、周波数領域で周波数毎にＩＣＡを適用する方法が現実的である。
［到来方向推定］
周波数領域でのＩＣＡを用いて信号源の方向を位置情報として推定する従来技術を図１を参照して簡単に説明する。Ｊ個のセンサ１_１，１_２，…，１_Ｊが直線状に配列されている。センサ１_ｊ（ｊ＝１，２，…，Ｊ）の位置をｄ_ｊ、そのセンサ１_ｊでの観測信号をｘ_ｊ（ｔ）とする。センサ１_１，…，１_ｊの配列方向と垂直な方向を９０°とし、源信号ｓ_ｉ（ｔ）の到来方向を０°≦θ_ｊ≦１８０°とする。Ｉ個の源信号ｓ_１（ｔ），…，ｓ_Ｉ（ｔ）の混合信号をＪ個のセンサ１_１〜１_Ｊにより、観測信号ｘ_１（ｔ），…，ｘ_Ｊ（ｔ）として検出しているものとする。
信号の到来方向の推定は周波数領域で行われることが多い。これには、観測信号ｘ_ｊ（ｔ）を短時間フーリエ変換して周波数領域における時間系列信号Ｘ_ｊ（ω，ｍ）を求める。ωは角周波数（周波数をｆとするとω＝２πｆ）ｍは時刻を表す番号である。源信号ｓ_ｉ（ｔ）（ｉ＝１，…，Ｉ）も、同様に、周波数領域における時間系列Ｓ_ｉ（ω，ｍ）に変換したものとすると、観測信号Ｘ_ｊ（ω，ｍ）は、信号ｓ_ｉの信号源からセンサ１_ｊへの周波数応答Ａ_ｊｉ（ω）を用いＸ_ｊ（ω，ｍ）＝Σ_ｉ＝１ ^ＩＡ_ｊｉ（ω）Ｓ_ｉ（ω，ｍ）
と表現することができる。ベクトルと行列を用いると、

となる。ここで、

はそれぞれ、Ｊ個のセンサによる観測信号およびＩ個の源信号をベクトル表現したものである。Ａ（ω）は周波数応答Ａ_ｊｉ（ω）を要素とするＪ×Ｉ行列であり、信号混合系の周波数応答を表現するので混合行列と呼ばれる。［ａ］^Ｔはベクトル又は行列ａの転置を表わす。
図１において、方向θ_ｉから到来する源信号は、ｄ_１＝０の位置のセンサ１_１に対し、センサ１_ｊにもτ_ｊｉ＝ｃ^−１ｄ_ｊｃｏｓθ_ｉだけ速く到達する。ｃは源信号ｓ_ｉの速度である。従って、直接波のみを考慮すると角周波数ωにおける周波数応答は

とモデル化することができる。方向θを変数とする到来方向ベクトルを
ａ（ω，θ）＝［ｅｘｐ（ｊωｃ^−１ｄ_１ｃｏｓθ）、ｅｘｐ（ｊωｃ^−１ｄ_２ｃｏｓθ），…，ｅｘｐ（ｊωｃ^−１ｄ_Ｊｃｏｓθ）］^Ｔと表現すると、観測信号は

とも近似表現することができる。
独立成分分析を用いて信号源の方向を推定する方法については例えばＳ．Ｋｕｒｉｔａ，Ｈ．Ｓａｒｕｗａｔａｒｉ，Ｓ Ｋａｊｉｔａ，Ｋ．Ｔａｋｅｄａ，ａｎｄ Ｆ．Ｉｔａｋｕｒａ，“Ｅｖａｌｕａｔｉｏｎ ｏｆ ｂｌｉｎｄ ｓｉｇｎａｌ ｓｅｐａｒａｔｉｏｎ ｍｅｔｈｏｄ ｕｓｉｎｇ ｄｉｒｅｃｔｉｖｉｔｙ ｐａｔｔｅｒｎ ｕｎｄｅｒ ｒｅｖｅｒｂｅｒａｎｔ ｃｏｎｄｉｔｉｏｎｓ，”ｉｎ Ｐｒｏｃ．ＩＣＡＳＳＰ２０００，２０００，ｐｐ．３１４０−３１４３（文献１という）に示されている。この方法を以下に簡単に述べる。
観測信号Ｘ（ω，ｍ）はＸ（ω，ｍ）＝Ａ（ω）Ｓ（ω，ｍ）であり、源信号Ｓ（ω，ｍ）が混合されたものであるので、互いに独立ではない。Ｘ（ω，ｍ）に独立成分分析

を適用すると、互に独立な分離信号

が得られる。
Ｗ（ω）は要素がＷ_ｉｊ（ω）であるＩ×Ｊの行列であり、分離行列と呼ばれる。例えば、Ｉ＝Ｊ＝２の場合、独立成分分析は、Ｙ_１（ω，ｍ）とＹ_２（ω，ｍ）が互いに独立になる様に

を満足する分離行列Ｗ（ω）を求める。源信号Ｓ_１（ω，ｍ），…，Ｓ_Ｉ（ω，ｍ）が互いに独立であれば、分離信号Ｙ_１（ω，ｍ），…，Ｙ_Ｉ（ω，ｍ）はそれぞれ、源信号の何れかに対応する。但し、独立成分分析は信号の独立性のみに基づいているところから、分離信号の順序と大きさに関して任意性がある。即ち、分離行列Ｗ（ω）の行を入れ換えても、Ｗ（ω）の行を定数倍しても、独立成分分析の解である。なお後で述べるが順序の任意性はパーミュテーションの問題となり、大きさの任意性はスケーリングの問題となる。
分離行列Ｗ（ω）のｉ行目は、

である。このｗ_ｉ（ω）は分離信号Ｙ_ｉ（ω，ｍ）を作り出していることが分かる。従って、ｗ_ｉ（ω）は、源信号Ｓ_１（ω，ｍ），…，Ｓ_１（ω，ｍ）の何れか１つを強調して取り出し、それ以外を抑圧している。Ｗ_ｉ（ω）が形成する指向特性を解析することにより、何れの方向から到来する信号を取り出し、何れの方向から到来する信号を抑圧しているかを解析することができ、この解析により、源信号ｓ_ｉ（ｔ）の到来方向を推定することができる。これを全てのｗ_ｉ（ω），ｉ＝１，…，Ｉに対して行うと、分離行列のＷ（ω）各行ｗ_ｉ（ω）が取り出している源信号の到来方向Θ（ω）＝［θ_１（ω），…，θ_１（ω）］^Ｔを推定することができる。
Ｗ_ｉ（ω）が形成する指向特性は、到来方向ベクトルａ（ω，θ）を用いてＢ_ｉ（ω，θ）＝ｗ_ｉ（ω）ａ（ω，θ）と表現することができる。このＢ_ｉ（ω，θ）は方向θにある源信号から、分離信号Ｙ_ｉ（ω，ｍ）への周波数応答と考えられる。周波数３１５６Ｈｚにおいて独立成分分析後の指向特性のゲイン｜Ｂ_ｉ（ω，θ）｜特性を図２に示す。図２の横軸はθを、縦軸はゲインを表わす。実線は分離行列の１行目が与える指向特性｜Ｂ_１（ω，θ）｜であり、破線は２行目が与える指向特性｜Ｂ_２（ω，θ）｜である。実線は５５°でゲインが最小となっており、破線は１２１°でゲインが最小となっている。このことから、分離行列の１行目は５５°から到来する信号を抑圧して１２１°から到来する信号を取り出し、逆に分離行列の２行目は１２１°から到来する信号を抑圧して５５°から到来する信号を取り出している。従って、Θ（３１５６Ｈｚ）＝［１２１°，５５°］^Ｔと推定することができる。
複数の信号源方向を複数センサを用いて、センサの観測信号を周波数領域に変換して推定する方法として、ＭＵＳＩＣ（Ｍｕｌｔｉｐｌｅ Ｓｉｇｎａｌ Ｃｌａｓｓｉｆｉｃａｔｉｏｎ）法がある（例えばＳ．Ｕｎｎｉｋｒｉｓｈｎａ Ｐｉｌｌａｉ，“Ａｒｒａｙ Ｓｉｇｎａｌ Ｐｒｏｃｅｓｓｉｎｇ”，Ｓｐｒｉｎｇｅｒ−Ｖｅｒｌａｇ，１９８９，ＩＳＢＮ ０−３８７−９６９５１−９，ＩＳＢＮ ３−５４０−９６９５１−９参照）。この方法はセンサの数Ｊより１個少ない（Ｊ−１）個までしか信号源の方向を推定できない。しかし前記独立成分分析（ＩＣＡ）を用いる方法（以下単にＩＣＡ法ともいう）によれば２つのセンサで２信号の混合に対応することができ、この点で、ＭＵＳＩＣ法より優れている。しかし、上記のＩＣＡを用いる方法では３信号以上の混合に対処するには以下に示すように困難を伴う。また、指向特性のゲインの最小値を求める操作に多くの計算を要することも欠点である。
３信号の混合に対して３つのセンサを用いてＩＣＡ法を適用した状況を説明する。この場合のＩＣＡは分離行列が３×３になるだけで同様に行えるが、指向特性のゲインを解析するのに困難を伴う。この場合のＩＣＡ後の周波数２７３４Ｈｚにおいて指向特性のゲイン｜Ｂ_ｉ（ω，θ）｜を図３に示す。図３の実線は分離行列の１行目が、破線は２行目が、１点鎖線は３行目がそれぞれ与える指向特性である。この場合、各源信号は分離行列のある行によって強調され、他の２つの行によって抑圧される筈である。しかし、抑圧する２つの行が同じ方向で抑圧しているとは限らない。例えば、図３においては４５°付近で｜Ｂ_２｜と｜Ｂ_３｜が共に極小となり、ｗ_１（ω）が４５°付近の源信号を取り出し、ｗ_２（ω）及びｗ_３（ω）はこの源信号を抑圧している。同様に９０°付近で｜Ｂ_１｜及び｜Ｂ_３｜が極小となり、ｗ_２（ω）で９０°付近の源信号を取り出し、ｗ_１（ω）及びｗ_３（ω）でこの信号を抑圧していると読み取れる。しかしｗ_３（ω）が１２０°付近の源信号を取り出し、ｗ_１（ω）およびｗ_２（ω）がこの源信号を抑圧していると考えられるが、１２０°付近での｜Ｂ_１｜及び｜Ｂ_２｜の極小がかなりずれている。この様なずれが大きくなると、何れの方向の抑圧が何れの源信号に対応するか判然としなくなる。従って、ＩＣＡ法による従来技術を３信号以上の状況に対して適用するのは困難であると考えられる。
〔ブラインド信号分離〕
次にＩＣＡを用いたブラインド信号分離の従来技術を説明する。ブラインド信号分離とは、観測された混合信号から源信号を推定する技術である。以下では、Ｉ個の源信号が混在した混合信号がＪ個のセンサで観測される場合を例にとって説明する。
信号源ｉから発生した源信号をｓ_ｉ（ｔ）（ｉ＝１，…，Ｉ、ｔは時刻）、センサｊで観測された混合信号をｘ_ｊ（ｔ）（ｊ＝１，…，Ｊ）とすると、混合信号ｘ_ｊ（ｔ）は次式で表わせる。

ここで、ａ_ｊｉは信号源ｉからセンサｊへのインパルス応答、＊は畳み込み演算子である。ブラインド信号分離の目的は、観測信号ｘ_ｊ（ｔ）のみを用いて、分離のためのフィルタｗ_ｋｊ及び分離信号ｙ_ｋ（ｔ）（ｋ＝１，…，Ｉ）を次式により求めることである。

時間領域での畳み込み混合は、周波数領域における複数の瞬時混合に変換することができる。すなわち、前記の式（７）及び式（８）は、それぞれ前記式（１）及び式（５）で表わされ、前述したように、Ｗ（ω）は分離行列であり、ＩＣＡを用いてＹ_ｋ（ω，ｍ）とＹ_ｋ（ω，ｍ）が互いに独立となるように計算されたもの、つまりＩＣＡの解である。
以上の周波数領域でのブラインド信号分離を行うに際し問題となるのは、パーミュテーションの問題とスケーリングの問題である。
前述したように分離行列Ｗ（ω）は行の入れ替えを行っても独立成分分析の解となる。つまりある角周波数ωにおいて、Ｗ（ω）がＩＣＡの解だとし、複素数を要素とする任意の対角行列をＤ（ω）とし、任意のパーミュテーション行列（この行列を任意の行列の左から掛けた結果は、この任意の行列の行を入れ替えた行列となる）をＰ（ω）とした場合におけるＰ（ω）Ｄ（ω）Ｗ（ω）もまたＩＣＡの解となる。これは、ＩＣＡが源信号間の統計的独立性のみを条件として源信号の分離を行うことに基づくものであり、このＤ（ω）によって与えられる解の自由度をスケーリングの自由度と呼び、Ｐ（ω）によって与えられる解の自由度をパーミュテーションの自由度と呼ぶ。
そのため、適切なブラインド信号分離を行うためには、全てのωについて、ＩＣＡの解から分離行列として適切な解Ｗ（ω）を特定しなければならない。一般に、この適切な解Ｗ（ω）の特定は、任意に求めたＩＣＡの解に、適切なＤ（ω）やＰ（ω）を乗じ、その結果を適切な解Ｗ（ω）とすることによって行われる。そして、全てのωについて、このＤ（ω）を適切に決める問題をスケーリングの問題と呼び、Ｐ（ω）を適切に決める問題をパーミュテーションの問題と呼ぶ。また、パーミュテーションとは、｛１，２，…，Ｉ｝から｛１，２，…，Ｉ｝への全単射な関数Ｚ：｛１，２，…，Ｉ｝→｛１，２，…，Ｉ｝であり、パーミュテーション行列と１対１対応する。
スケーリングの自由度は、時間領域において周波数特性を変化させるフィルタの自由度に相当する。それゆえ、時間領域において歪のない分離信号を生成するためには、全てのωについて、Ｄ（ω）を適切に決める必要がある。このスケーリングの問題は、例えば、Ｄ（ω）＝ｄｉａｇ（Ｗ^−１（ω））とすることで容易に解決することができる。ｄｉａｇ（α）は行列αの対角化（対角要素の他は全ての要素を０とする）を表わす。つまり任意に求めたＩＣＡの解Ｗ_０（ω）の逆行列を求め、それを対角化した行列をＤ（ω）とし、Ｄ（ω）Ｗ_０（ω）を、適切な分離行列Ｗ（ω）として特定する。このことは既に知られている。例えば参考文献：Ｋ．Ｍａｔｓｕｏｋａ ａｎｄ Ｓ．Ｎａｋａｓｈｉｍａ，“Ｍｉｎｉｍａｌ Ｄｉｓｔｏｒｔｉｏｎ Ｐｒｉｎｃｉｐｌｅ ｆｏｒ Ｂｌｉｎｄ Ｓｏｕｒｃｅ Ｓｅｐａｒａｔｉｏｎ”，Ｐｒｏｃ．ＩＣＡ ２００１，ｐｐ．７２２−７２７．に示されている。
一方、パーミュテーションの自由度のため、前記の式（５）の演算結果として、例えば、ある角周波数ω_１については分離信号Ｙ_１（ω_１，ｍ）が源信号Ｓ_１（ω_１，ｍ）の推定値として出力され、別の角周波数ω_２については分離信号Ｙ_１（ω_２，ｍ）が源信号Ｓ_２（ω_２，ｍ）の推定値として出力されることもありうる。このような場合、時間領域の出力信号ｙ_１（ｔ）の中に、時間領域の源信号ｓ_１（ｔ）の成分とｓ_２（ｔ）の成分とが混在してしまい、分離信号を正しく生成することができない。そのため、時間領域における出力信号ｙ_１（ｔ）が正しく源信号ｓ_１（ｔ）の推定値となるためには、全てのωについて、Ｙ_１（ω，ｍ）がＳ_１（ω，ｍ）の推定値となるようにＰ（ω）を適切に決める必要がある。
従来のパーミュテーション問題の代表的な解法として、前記文献１に示す信号の到来方向を推定する方法がある。つまり、図２を参照して説明したように、各周波数における分離行列Ｗ（ω）の各行と対応する指向特性を求める（図２ではｆ＝３１５６Ｈｚの指向特性のみを示している）。これら各周波数の指向特性において例えば分離行列Ｗ（ω）の１行目が与える指向特性のゲイン最小値が５５°で、２行目が与える指向特性のゲイン最小値が１２１°になるようにする。つまりある角周波数ω_ｎではその分離行列Ｗ（ω_ｎ）の１行目が与える指向特性のゲイン最小値が１２１°で、２行目が与える指向特性のゲイン最小値が５５°であれば、その分離行列Ｗ（ω_ｎ）の１行目と２行目とを入れ替える。つまりこの入れ替えを行うパーミュテーション行列Ｐ（ω_ｎ）をＷ（ω_ｎ）にその左から掛算する。
このパーミュテーション問題の解法は先にも述べたように指向特性のゲイン最小値を求める操作に多くの計算を要し、しかも信号源の数Ｉが３以上では全角周波数のＷ（ω）について適当に並べ替えを行ってみるという試行錯誤が必要となる。更に図３を参照して説明したようにＷ（ω）のある行がある方向の信号Ｓ_ｉ（ω，ｍ）を取り出している場合は、そのＷ（ω）の他の行は全てその方向の信号Ｓ_ｉ（ω，ｍ）を抑圧する状態に必ずしもならない。
また指向特性のゲインが低い方向を探索して推定する信号到来方向の推定精度は信号源の位置によって異なり、特に一対のセンサ１_ｊと１_ｊ’を結んだ直線（以下、センサペア軸という）に信号到来方向が近い場合に誤差が大きくなる。このことを実験により示す。図４Ａに示すように、２つのセンサ１０１，１０２としてマイクロホンを距離２．８３ｃｍ離して配置し、センサ１０１と１０２の中点（原点）から一定の距離（約１５０ｃｍ）離れ、かつ角度間隔が２０°離れた２つの信号源１１１，１１２として音源を設け、センサ１０１からセンサ１０２を見た方向を基準（０°）として、原点からみた信号源１１１の角度θが１０°から１５０°となるまで音源１１１と１１２を前記一定距離及び角度間隔を保持した状態で移動させた。
このように音源１１１，１１２を移動させながら行ったブラインド信号分離結果を図４Ｂに例示する。図４Ｂの縦軸は、信号対妨害音比（Ｓｉｇｎａｌ ｔｏ Ｉｎｔｅｒｆｅｒｅｎｃｅ Ｒａｔｉｏ）を示しており、分離音中に含まれる目的信号及び妨害信号を用いて、
ＳＩＲ＝１０ｌｏｇ_１０（目的信号のパワー／妨害信号のパワー）（ｄＢ）
のように計算したものである。また、図４Ｂの横軸は、原点からみた音源１１１の角度θを表し、実線はこの実験結果を示しており、点線はパーミュテーションに正解を与えたときのＳＩＲを示している。
この図より信号源１１１が、センサペア軸に近い方向（０°或いは１８０°付近）に近づくと、実験結果のＳＩＲはパーミュテーションに正解を与えたときのＳＩＲに比べ大きく低下していることがわかる。このことは信号源１１１の方向がセンサペア軸方向と近いとパーミュテーションが誤っていることに基づくと考えられる。
以上述べたように独立成分分析を用いて分離行列を求め、その各行から指向特性パターンを求めて、その利得が低い方向探索により、信号源の方向（到来信号方向）を求め、更にこれを利用してブラインド信号分離を行う場合は、前記指向特性パターンを求めて、利得が低い方向を探索するために多くの計算時間を必要とした。
この発明の目的は独立成分分析により分離行列を求めて信号源の位置情報を推定するための計算時間が短かい位置情報推定装置、その方法及びプログラムを提供することにある。

この発明によれば周波数領域の分離行列Ｗ（ω_１），…，Ｗ（ω_Ｎ）の逆行列（Ｉ＜Ｊの場合は擬似逆行列）を計算してスケーリングとパーミュテーションの自由度を含む混合行列Ａ（ω_１），…，Ａ（ω_Ｎ）の推定値Ｈ（ω_１），…，Ｈ（ω_Ｎ）を生成し、これらの周波数ごとのＨ（ω_ｎ）（ｎ＝１，…，Ｎ）について各列ごとにその２つの要素Ｈ_ｊｉ（ω_ｎ）とＨ_ｊ’ｉ（ω_ｎ）（ｊ、ｊ’はセンサを表すパラメータ、ｉは信号源を表すパラメータ）の比から、信号源ｉの位置情報の一つのパラメータ、例えば信号源が存在する円錐面あるいは曲面を計算する。
このように要素比に基づき表された計算式を演算すればよく、分離された信号の指向性パターンを求め、更にその極小値角度の探索を行う場合より少ない計算量で済む。また要素比をとるため、前述したスケーリングの自由度による影響がなくなる。

図１はセンサアレーと到来信号の各センサへの到来時間差の関係を説明するための図である。
図２は２音源混合信号に対し、ＩＣＡにより計算した分離行列の各行のゲイン指向特性を示す図である。
図３は３音源混合信号に対し、ＩＣＡにより計算した分離行列の各行のゲイン指向特性を示す図である。
図４Ａは予備実験に用いたセンサと信号源の関係を説明するための図である。
図４Ｂは前記予備実験の結果を示す図である。
図４Ｃは推定角度とその誤差及び偏角誤差との関係を示す図である。
図５はこの発明を信号到来方向の推定に適用した第１実施形態の機能構成例を示すブロック図である。
図６は第１実施形態の処理手順の例を示す流れ図である。
図７は図５中の角度計算部の具体的機能構成例を示すブロック図である。
図８は図６中のステップＳ４の具体的処理手順の例を示す流れ図である。
図９はこの第１実施形態により方向推定した実験結果を示す図である。
図１０Ａは音源が２個、マイクロホンが３個の場合に第１実施形態により方向推定した実験結果を示す図である。
図１０Ｂは図１０Ａについての実験と同一条件でＭＵＳＩＣ法により方向推定した実験結果を示す図である。
図１１Ａは音源が２個、マイクが３個の場合に第１実施形態により方向推定した実験結果を示す図である。
図１１Ｂは図１１Ａについての実験と同一条件でＭＵＳＩＣ法により方向推定した実験結果を示す図である。
図１２はこの発明をブラインド信号分離に適用した第２実施形態の機能構成例を示すブロック図である。
図１３は第２実施形態の処理手順の例を示す流れ図である。
図１４は図１３中のステップＳ１４の具体的処理手順の例を示すブロック図である。
図１５は複数円錐面の共通直線方向を説明するための図である。
図１６は図１２中の到来方向決定部１６の具体的機能構成の他の例を示すブロック図である。
図１７はこの発明をブラインド信号分離に適用した第３実施形態の機能構成例を示すブロック図である。
図１８は第３実施形態の処理手順の例を示す流れ図である。
図１９は図１８中のステップＳ３５の具体的処理手順の例を示す流れ図である。
図２０はセンサ配置と信号源位置と推定曲面との関係を説明するための図である。
図２１は推定球面の例を示す図である。
図２２は実験に用いた室とマイクロホン配置と音源との関係を示す図である。
図２３は小間隔マイクロホン対を用いた推定方向のヒストグラムである。
図２４はその周波数に対する推定方向の分布を示す図である。
図２５は音源２_４及び２_５に対する推定球面の半径の周波数に対する分布を示す図である。
図２６は各方法での実験結果を示す図である。
図２７はこの発明の第４実施形態の要部の機能構成例を示すブロック図である。
図２８はこの発明をブラインド信号分離に適用した第５実施形態の処理手順の例を示す流れ図である。
図２９は第５実施形態に対する実験結果を示す図である。
図３０Ａは複数の円錐面の共通直線方向を求めて分離行列のパーミュテーションを解決する処理手順の要部を示す流れ図である。
図３０Ｂは円錐面の推定と球面の推定を利用して分離行列のパーミュテーションを解決する処理手順の要部を示す流れ図である。
図３０Ｃは円錐面の推定と球面の推定を利用して分離行列のパーミュテーションを解決する他の処理手順の要部を示す流れ図である。

まず信号源の位置情報として方向情報を推定する場合にこの発明を適用した実施形態を説明する。なお以下の説明において同一又は対応するものについて図面中に同一の参照番号を付けて、重複説明を省略する。
［第１実施形態］
この第１実施形態は信号源の方向、つまりその信号源から放射された源信号の到来方向を求める。
図５に第１実施形態の機能構成図を、図６にその処理手順の一部の流れ図をそれぞれ示す。
信号源の数Ｉ以上の個数であるＪ個のセンサ１_１，１_２，…，１_Ｊが、例えば図１に示したように配列されている。隣接センサの間隔は通常は、源信号の最も短かい波長の１／２以下である。これらセンサ１_ｊ（ｊ＝１，２，…，Ｊ）で観測された信号ｘ_ｊ（ｔ）はそれぞれ周波数領域変換部１１_ｊにおいて、例えば短時間フーリエ変換により周波数領域信号Ｘ_ｊ（ω，ｍ）に変換される（図６、ステップＳ１）。これら周波数領域信号Ｘ_ｊ（ω，ｍ）に対する各角周波数ω_ｎごとの分離行列Ｗ（ω_ｎ）（ｎ＝１，２，…，Ｎ）が分離行列算出部１２において独立成分分析により算出される（図６、ステップＳ２）。

この各周波数ごとの分離行列Ｗ（ω_ｎ）の逆行列が逆行列算出部１３で算出され逆行列Ｈ（ω_ｎ）が求まる（図６、ステップＳ３）。

なおこの逆行列の計算はＪ＞Ｉの場合は疑似逆行列を計算することになる。疑似逆行列としては例えばＭｏｏｒｅ−Ｐｅｎｒｏｓｅ型一般逆行列が用いられる。
この実施形態においてはＨ（ω）中の少くとも１つの周波数の逆行列Ｈ（ω_ｎ）の各列ｉの２つの要素Ｈ_ｊｉ（ω_ｎ）とＨ_ｊ’ｉ（ω_ｎ）との比の偏角から方向情報、つまり源信号の到来方向が角度計算部１４で計算される（図６、ステップＳ４）。この角度計算部１４の具体的機能構成例と処理手順例を図７及び図８を参照して説明する。選択部１４ａにより角周波数ω_ｎの逆行列Ｈ（ω_ｎ）中の未選択の一つの列ｉが選択され（図８、ステップＳ４ａ）、そのｉ列目から２つの要素Ｈ_ｊｉ（ω_ｎ）とＨ_ｊ’ｉ（ω_ｎ）が選択される（図８、ステップＳ４ｂ）。
偏角計算部１４ｂで、選択された要素の比の偏角ａｒｇ［Ｈ_ｊｉ（ω_ｎ）／Ｈ_ｊ’ｉ（ω_ｎ）］が計算され、また間隔計算部１４ｃで、センサ情報格納部１５内のセンサ１_ｊと１_ｊ’の位置情報ｄ_ｊとｄ_ｊ’が取り出され、これらの差からセンサ１_ｊと１_ｊ’の間隔ｄ_ｊ−ｄ_ｊ’が計算される（図８、ステップＳ４ｃ）。位相回転計算部１４ｄで、間隔ｄ_ｊ−ｄ_ｊ’と、角周波数ω_ｎでの単位距離当りの信号の位相回転量ω_ｎ／ｃ（ｃは源信号の速度）との積が計算されて間隔ｄ_ｊ−ｄ_ｊ’での信号の位相回転量が求められ、この位相回転量により偏角ａｒｇ［Ｈ_ｊｉ（ω_ｎ）／Ｈ_ｊ’ｉ（ω_ｎ）］が除算部１４ｅで割算される（図８、ステップＳ４ｄ）。なお単位距離当たりの位相回転量ω／ｃは例えばｆ＝６８０Ｈｚの音波（速度ｃ＝３４０ｍ／ｓ）では、１ｍ当たり２π・２の位相が回転することになる。
この除算部１４ｅの割算結果Ｇ_ｉ（ω_ｎ）の絶対値が１以下か否かかが判定部１４ｆで判定され（図８、ステップＳ４ｅ）、１以下であればＧ_ｉ（ω_ｎ）の逆余弦θ_ｉ（ω_ｎ）＝ｃｏｓ^−１Ｇ_ｉ（ω_ｎ）が逆余弦計算部１４ｇで計算される（図８、ステップＳ４ｆ）。つまり角度計算部では次の計算が行われる。

ステップＳ４ｅで｜Ｇ_ｉ（ω_ｎ）｜が１以下でなければ角度θ_ｉ（ω_ｎ）が虚数になるので、別の組み合せを選択するため、判定部１４ｆで選択したｉ列目のすべての要素の組み合せを選択したか判定され（図８、ステップＳ４ｇ）、選択していない組み合せがあればステップＳ４ｂに戻り、全ての組み合せを選択していれば、判定部１４ｆですべての列を選択したか判定され（図８、ステップＳ４ｈ）、選択していない列があればステップＳ４ａに戻り、全ての列を選択していれば角度計算処理を終了する。なお図７中の間隔計算部１４ｃは各周波数の角度計算部１４に対して共通に用いられる。
角度計算部１４からＨ（ω）中の選択した角周波数ω_ｎの逆行列Ｈ（ω_ｎ）における各列と対応して、その選択順に、もし第１列目から順次選択すればその順に式（７）の計算結果、つまりＩ個の信号源の方向（信号到来方向）Θ（ω_ｎ）＝（θ_１（ω_ｎ），θ_２（ω_ｎ），…，θ_Ｉ（ω_ｎ））が出力される。これらθ_１（ω_ｎ），θ_２（ω_ｎ），…，θ_Ｉ（ω_ｎ）は各源信号ｓ_１（ｔ），ｓ_２（ｔ），…，ｓ_Ｉ（ｔ）の到来方向（信号源１，２，…，Ｉの方向）のいずれかの１つと対応する。
この実施形態において信号到来方向が推定できるしくみを以下に説明する。独立成分分析（ＩＣＡ）により分離が達成されていれば、ＩＣＡにより算出した分離行列Ｗ（ω）と真の混合行列Ａ（ω）とはＰ（ω）Ｄ（ω）Ｗ（ω）Ａ（ω）＝Ｉの開係にある。ここで、Ｄ（ω）はスケーリングの自由度を示す対角行列、Ｐ（ω）はパーミュテーションの自由度を示すパーミュテーション行列、Ｉは単位行列である。ＩＣＡを用いても、混合行列Ａ（ω）そのものは一般には算出できない。しかし、Ｗ（ω）の逆行列Ｈ（ω）＝Ｗ^−１（ω）＝Ａ（ω）Ｐ（ω）Ｄ（ω）を算出すれば、スケーリングとパーミュテーションの自由度を含む混合行列の推定値が得られる。すなわち、逆行列Ｈ（ω）は、混合行列Ａ（ω）の列をＰ（ω）により並べ替え、さらに各列にそれぞれＤ（ω）の対角要素を掛け合わせたものになる。
この実施形態では、逆行列Ｈ（ω）の同じ列ｉから２つの要素Ｈ_ｊｉ（ω）とＨ_{ｊ’ ｉ}（ω）を取り出してそれらの比Ｈ_ｊｉ（ω）／Ｈ_ｊ’ｉ（ω）を求めることで、算出できないＤ（ω）によるスケーリングの自由度を取り除く。すなわち、

となる。ここで、Ｚはパーミュテーション行列Ｐ（ω）を右から掛けることに対応する順列である。逆行列Ｈ（ω）の全ての列ｉに対して式（２１）に関わる操作を行うことで、Ｐ（ω）によるパーミュテーションＺにかかわらず、全ての信号の到来方向が推定できる。
背景技術の説明では、混合行列Ａ（ω）の要素をＡ_ｊｉ（ω）＝ｅｘｐ（ｊωｃ^−１ｄ_ｊｃｏｓθ_ｉ）とモデル化した。しかし、この実施形態では、分離行列Ｗ（ω）の逆行列Ｈ（ω）を用いて、スケーリングとパーミュテーションの自由度を含む混合行列Ａ（ω）の推定値を算出しているので、この単純なモデルでは不十分である。そこで、振幅の減衰度α_ｊｉ（実数）と原点における位相差ｅｘｐ（ｊφ_ｉ）を用いて、Ａ_ｊｉ（ω）＝α_ｊｉｅｘｐ（ｊφ_ｉ）ｅｘｐ（ｊωｃ^−１ｄ_ｊｃｏｓθ_ｉ）とモデル化する。このモデルを用いてＡ_{ｊＺ（ｉ）}（ω）／Ａ_{ｊ’Ｚ（ｉ）}（ω）を計算すると、式（２１）の関係から、

が得られる。その結果、上記のＧ_ｊ（ω）は

となり、｜Ｇ_ｉ（ω）｜＜１であればθ_Ｚ（ｉ）＝ｃｏｓ^−１Ｇ_ｉ（ω）が実数となり、到来方向を推定することができる。先に述べたように、パーミュテーションＺの自由度により、Ｉ個の方向Θ（ω）＝［θ_Ｚ（ｉ）（ω），…，θ_Ｚ（Ｉ）（ω）］全体を適当に並べ替えたものが、信号ｓ_１，…，ｓ_Ｉの方向に対応する。
Ｈ（ω）中の複数の周波数又は全ての周波数の各逆行列Ｈ（ω_ｎ）（ｎ＝１，…，Ｎ）について先に述べたように各列ごとに角度θ_ｉ（ω_ｎ）を求め、これらの全体により各到来方向を決定してもよい。つまり角度計算部１４で推定された各周波数の到来方向は図５中のソート部３２においてその角度ごとに分けられる（図６、ステップＳ５）。例えばそれぞれ方向（角度）の大きい順に並べられる。角周波数ω_１についてのΘ（ω_１）中の各成分が大きい順に（θ_１（ω_１），θ_２（ω_１），…，θ_Ｉ（ω_１））と並べ、ω_２についてのΘ（ω_２）中の各成分が大きい順に（θ_１（ω_２），θ_２（ω_２），…，θ_Ｉ（ω_２））と並べ、…，ω_ＮについてのΘ（ω_Ｎ）中の各成分が大きい順に（θ_１（ω_Ｎ），θ_２（ω_Ｎ），…，θ_Ｉ（ω_Ｎ））と並べられる。このように大きい角度順に分別された角度は、大きい順が同一であっても周波数間でばらつきがある、つまり角度θ_ｉ（ω_１），θ_ｉ（ω_２），…，θ_ｉ（ω_Ｎ）（ｉ＝１，…，Ｉ）はばらついている。
よって統合部３３で、分別角度θ_ｉ（ω_１），θ_ｉ（ω_２），…，θ_ｉ（ω_Ｎ）ごとに１つの角度θ_ｉに統合されて１つの角度（到来方向）θ_ｉとされる（図６、ステップＳ６）。この統合の方法としては例えば分別角度（θ_ｉ（ω_１），θ_ｉ（ω_２），…，θ_ｉ（ω_Ｎ））ごと平均値が統合角度θ_ｉとされ、あるいは分別角度（θ_ｉ（ω_１），θ_ｉ（ω_２），…，θ_ｉ（ω_Ｎ））ごとの最頻値や中央値などを統合角度θ_ｉとしてもよい。このように分別角度ごとに１つの角度に統合することにより、１つのある周波数の逆行列Ｈ（ω_ｎ）からのみ到来方向を推定する場合より、正しく到来方向を推定することができる。
なお図８に示した角度（方位）計算推定処理において、その１つでも選択した列について、角度θ_ｉ（ω_ｎ）を求めることができなかった時は、そこでその逆行列Ｈ（ω_ｎ）に対する処理を終了し、他の周波数の逆行列Ｈ（ω_ｎ）に対する処理を行い、最初に全ての列について角度（方向）を計算することができたら、その時の計算結果を推定（方向）θ_１，…，θ_１としてもよい。あるいは各周波数の逆行列Ｈ（ω_ｎ）中で全ての列について角度計算ができた結果を角度ごとに分別し、統合してもよい。あるいは１つの列について角度θ_ｉ（ω_ｎ）を計算することができないものがあれば、その状態が最初に発生した時に、その後の処理を全て終了し、あらためて観測信号を求める処理からやりなおすことにより、推定結果の信頼性を高めるようにしてもよい。
この第１実施形態の実験例を述べる。残響時間が１９０ｍｓの部屋に３つのマイクロホンを５６．６ｍｍ間隔で１列に配列し、その配列方向を基準として角度４８°、７３°、１１９°の方向に３つの音源を配置した。これら音源からの音響信号の混合時間を６秒、観測信号の標本化周波数を８ｋＨｚ、空間エイリアシングが生じない最大周波数を３ｋＨｚ、短時間フーリエ変換のフレームを１０２４サンプルとした。前述したようにして、各周波数ごとにそれぞれ計算した各角度を図９に、横軸を周波数、縦軸を方向として示す。図９中の◇，＋，□はそれぞれ３つの音源方向の推定計算値を示す。この結果を３つの角度範囲に分別し、各分別された角度平均値は４５°，７４°，１２３°となった。３音源の方向を３つのマイクロホンで推定することはＭＵＳＩＣ法ではできないが、この実施形態によれば可成り正しく各音源方向が推定されていることがわかる。
この実施形態の方法とＭＵＳＩＣ法とを比較するため、音源方向を４８°と１１９°と比較的大きく角差がある状態で２つの音源を配置した場合について同様に実験を行った。図１０Ａにこの実施形態の方法により得られた結果を、図１０ＢにＭＵＳＩＣ法により得られた結果をそれぞれ示す。これらよりいずれの方法も、かなり正しく方向を推定していることがわかる。この実施形態による結果を２つの方向範囲で分別平均した結果は４５°，１２３°であり、ＭＵＳＩＣ法はそれぞれ４５°，１２２°である。次に、音源方向を１０５°と１１９°と方向が比較的接近している状態で２つの音源を配置した場合について同様に実験を行った。この実施形態方法による結果及びＭＵＳＩＣ法による結果を図１１Ａ及び図１１Ｂにそれぞれ示す。ＭＵＳＩＣ法では大部分の周波数で音源方向の推定ができないが、この実施形態では大部分の周波数で角度計算を行うことができ、しかも、これらを角度範囲で分別して、それぞれ平均した値は１０５°，１２４°であり、ＭＵＳＩＣ法は９４°，１２８°でありこの実施形態による方向推定がかなり正しいことがわかる。
以上述べたようにこの実施形態によれば従来方法において、指向性パターンの利得が低い方向探索する場合と比較して、式（９）に値を代入するだけで推定方向が求まるため計算時間がかなり短かい。なお前述したようにＩＣＡにより求めた分離行列Ｗ（ω）はスケーリングの自由度及びパーミュテーションの自由度を含んでいるためこれらを解決した分離行列Ｗ’（ω）の逆行列を計算し、つまり真の混合行列Ａ（ω）を計算し、その混合行列Ａ（ω）の各列について二つの要素の比を用いて到来方向を推定することが考えられるかもしれない。しかし、真の混合行列Ａ（ω）そのものを求めることは、例えば信号源の信号ｓ_ｉ（ｔ）の平均パワーを１にするなどの制約を与えない限りできない。無線通信分野においては信号源に対しそのような制約を与えることができるかも知れないが、例えば信号源の信号ｓ_ｉ（ｔ）が人間が直接発声した音声信号である場合は、そのような制約は実質的にはできない。一方この第１実施形態によれば、スケーリング及びパーミュテーションの自由度を含んでいる分離行列Ｗ（ω）の逆行列Ｈ（ω）の各列について二つの要素の比をとることによりスケーリングの自由度の問題を解決しておりどのような信号源に対しても適用でき、しかも前記両問題を解決した分離行列を計算する必要がなくそれだけ計算時間が短い。更に前述したように、各周波数について求めた推定方向を予め決めた順に分別すればパーミュテーションの問題も簡単に解決できる。センサの数Ｊと同数であっても各信号源の方向を推定することができる。また音源方向が比較的接近していてもかなり正しく推定することができる。
［第２実施形態］
この第２実施形態は信号源の一つの位置情報として方向情報を求めるものであり、かつ少なくとも２次元に配置された少なくとも３個のセンサを用いて、信号源がいずれの方向に位置していても方向を推定でき、従ってブラインド信号分離におけるパーミュテーションの問題も比較的簡単に解決するものである。つまり方向情報に基づく円錐面を推定し、複数の円錐面の共通直線を推定して、方向情報を決定する。
この第２実施形態をブラインド信号分離装置に適用した機能構成を図１２に、その処理手順を図１３に示す。例えば４個のセンサ１_１，１_２，１_３，１_４が同一円上で等間隔に配置され、そのいずれの２つのセンサの間隔も、源信号の最短波長の１／２以下とされる。以下の説明はセンサの数はＪ個であり、Ｊ≧３として行う。第１実施形態と同様に各センサｊ（ｊ＝１，…，Ｊ）で観測された観測信号ｘ_ｊ（ｔ）は、それぞれ周波数領域変換部１１で例えば、短時間フーリエ変換によって、周波数領域の信号Ｘ_ｊ（ω，ｍ）にそれぞれ変換される（ステップＳ１１）。
これら周波数領域の信号Ｘ_ｊ（ω，ｍ）から、分離行列算出部１２で、独立成分分析により、各周波数ごとの分離行列

が算出される（ステップＳ１２）。
逆行列算出部１３でこれら周波数ごとの各分離行列Ｗ（ω）の逆行列Ｈ（ω）がそれぞれ算出される（ステップＳ１３）。

この第２実施形態では、円錐面推定部１４において、周波数ごとの各逆行列Ｈ（ω）における各列ごとに異なる複数の要素対について各要素の比からこれら要素と対応する二つのセンサのセンサペア軸を中心軸とし、何れかの信号源が存在する円錐面、つまり一つの混合行列Ｈ（ω）の各列ごとに複数の円錐面がそれぞれ推定される（ステップＳ１４）。
円錐面推定部１４の機能構成は図１５中の角度計算部１４とほぼ同様であり、その処理手順も図８の手順と似ている。円錐面推定部１４において行われるある一つの周波数の逆行列Ｈ（ω）に対するステップＳ１４の処理の具体例を図１４を参照して説明する。
まず、例えば、円錐面推定部１４内のレジスタ内に格納されている制御パラメータｉとｐを０に初期化する（ステップＳ２０）。ここで、ｉは各信号源の番号と対応し、ｐはｉごとの推定済み円錐面の個数である。
ｉを＋１し（ステップＳ２１）、ｇを＋１し（ステップＳ２２）、制御パラメータｊ，ｊ’はそれぞれ互いに異なるｊ以下の自然数を、例えばランダムに選択する（ステップＳ２３）。この制御パラメータｊ，ｊ’の組は、同一のｉに対しては一度選択されたｊ，ｊ’の組み合わせの選択は行わないようにする。すなわち、ｊ≠ｊ’であり、また、例えば、ｉ＝１について、一度（ｊ，ｊ’）＝（１，２）が選択されると、ｉ＝１については、この処理が終了するまで、再びステップＳ２３において（ｊ，ｊ’）＝（１，２）という選択は行われないものとする。（また、この選択は、選択するｊ，ｊ’によって特定されるセンサペア軸（センサｊとセンサｊ’とを通る直線）が、このルーチン処理においてそれ以前に選択されたｊ，ｊ’によって特定されるセンサペア軸と同一直線上に配置されないように行われることが望ましい。これにより、円錐面推定部１４は、一定の誤差範囲内で中心軸が重ならない複数の円錐面を推定することになる。なお、各センサペア軸が同一直線上となるか否かの判断は、例えば、センサ情報格納部１５に各センサの位置を示すベクトルを格納しておき、ここから各センサの位置を示すベクトルの情報を抽出して行う。）
次に、センサ情報格納部１５から、ステップＳ２３で選択したパラメータｊに対応するｊ番目のセンサｊの位置を示すベクトルｄ_ｊと、変数ｊ’に対応するｊ’番目のセンサｊ’の位置を示すベクトルｄ_ｊ’とを抽出する（ステップＳ２４）。また、混合行列Ｈ（ω）から、ｊ行ｉ列要素Ｈ_ｊｉ（ω）と、ｊ’行ｉ列要素Ｈ_ｊ’ｉ（ω）とを選択して抽出する（ステップＳ２５）。以上の処理は図７中の選択部１４ａにより行われる。従って、ｉ，ｐ，ｊ，ｊ’や選択したｊ，ｊ’などを格納するレジスタも選択部１４ａに設けられることになる。
抽出したこれらの情報を用い、次式（９’）を演算する（ステップ２６）。

ここで‖ｄ_ｊ−ｄ_ｊ’‖はセンサ１_ｊと１_ｊ’との間隔（距離）である。この第２実施形態においては複数のセンサが２次元又は３次元に配置されているため各センサの位置情報は例えば図１２中のセンサ１_１〜１_４が配置されている円の中心を原点とする２又は３要素の座標ベクトルで表わされる。つまり式（９）はセンサが１次元に配置され、信号到来方向の角度が２次元の場合であるが、この式（９’）はセンサが２次元あるいは３次元に配置されていてもよく信号到来方向の角度が３次元空間でもよい場合に拡張されたものである。従って式（９’）は式（９）を包含する。この式（９’）により推定された角度θ＾_{ｉ，ｊｊ’}（ω）及びそのｉ，ｊ，ｊ’を円錐面情報として円錐面推定部１４内のレジスタ（メモリ）に一時記録する（ステップＳ２７）。なお図１２中に破線で示すように各周波数に対する間隔計算部１４ｃは共通に用いられる。また式（９）で演算した角度θ_ｉ（ω）は三次元空間においてセンサペア軸（センサ１_ｊと１_ｊ’とを結ぶ直線）に対する角度がθ_ｉ（ω）となる直線の無数の集合、つまり円錐面に信号源ｉが存在することを推定していることになる。この拡張された式（９’）の演算結果を単なるθ_ｉ（ω）ではなくθ＾_{ｉ，ｊｊ’}（ω）と表わした。ステップＳ２６で行う演算は図７中の偏角計算部１４ｂ、間隔計算部１４ｃ、位相回転計算部１４ｄ、除算部１４ｅ、判定部１４ｆ及び逆余弦計算部１４ｇにより、図８中のステップＳ４ｃ，Ｓ４ｄ，Ｓ４ｅ及びＳ４ｆの処理により行う。
次に、ｐ＝Ｐであるか否かを判断する（ステップＳ２８）。このＰは、ｉごとに推定すべき円錐面の個数であり、このステップでは、そのｉについて円錐面をＰ個推定したか否かを判断する。ｐ＝ＰでなければステップＳ２２に戻り、ｐ＝Ｐであれば、ｉ＝Ｉであるか否かを判断する（ステップＳ２９）。すなわち、すべてのｉについて円錐面の推定が終了したか否かを判断し、ｉ＝ＩでなければステップＳ２１に戻り、ｉ＝Ｉであれば処理を終了する（ステップＳ１４の具体例の説明終了。）。第１実施形態では選択したｉについて角度θ_ｉ（ω）を１つ推定すれば次のｉについての角度推定に移るがこの第２実施形態では各ｉごとに複数Ｐ個の組の角度（円錐面）θ＾_{ｉ，ｊｊ’}（ω）を推定する。ステップＳ２８及びＳ２９の処理は図７中の判定部１４ｆで行われる。
図１２中の到来方向決定部１６で、円錐面推定部１４で推定された複数の円錐面の情報（この例では、ｉ，ｊ，ｊ’、θ＾_{ｉ，ｊｊ’}（ω））から源信号の到来方向ｕ_ｉ（ω）＝（方位角θ_ｉ（ω），仰角φ_ｉ（ω））（ｉ＝１，…，Ｉ）が決定される（図１２、ステップＳ１５）。具体的には、例えば、ｉについて推定された複数の円錐面のうち、互いに線接触しているとみなされる共通直線の方向を、信号源ｉに対応する信号の到来方向ｕ_ｉ（ω）として算出する。
この信号の到来方向ｕ_ｉ（ω）の推定方法を図１５を参照して説明する。
センサ１_１と１_２の配列方向と直角方向にセンサ１_３が配され、センサ１_１と１_２間、センサ１_２と１_３間の各間隔は同一とされている。ｉ＝２の信号源２_２に対し、センサ１_１及び１_２によりセンサペア軸３_１２を中心軸とする円錐面４_１２（θ＾_２，１２（ω））が、センサ１_２と１_３によりセンサペア軸３_２３を中心軸とする円錐面４_２３（θ＾_２，２３（ω））が、センサ１_１と１_３によりセンサペア軸３_１３を中心軸とする円錐面４_１３（θ＾_２，１３（ω））がそれぞれ推定されている。これらのうち、円錐面４_１２と４_１３は、共通直線５_２で互いに線接触しているとみなされる。この共通直線５_２の方向ｕ_２（ω）が、信号源２_２から放射された信号の到来方向ｕ_２つまり信号源２_２の方向と決定する。なお円錐面４_２３も平行移動すれば円錐面４_１２，４_１３とほぼ線接触するが、後述する第４実施形態により円錐面４_２３は廃棄される。
図１２中の到来方向決定部１６におけるこの共通直線５_ｉの方向θ_ｉ（ω）を決定するための具体的計算方法の例を説明する。角周波数ωに対して信号源２_ｉが存在すると推定された複数の円錐面を４_ｊｊ’（１），…，４_ｊｊ’（Ｐ）とし、これら円錐面４_ｊｊ’（ｐ）（ｐ＝１，…，Ｐ）の推定に用いられた１対のセンサの位置情報をｄ_ｊ（ｐ），ｄ_ｊ’（ｐ）とし、角周波数ωについて推定された円錐面４_ｊｊ’（ｐ）と対応する角度をθ＾_ｊｊ’（ω，ｐ）とし、円錐面４_ｊｊ’（ｐ）を表わすベクトルをｕとする。
正規化軸ベクトル計算部１６ａで前記各一対のセンサの位置間を結ぶ軸ベクトル（ｄ_ｊ（ｐ）−ｄ_ｊ’（ｐ））をそれぞれ長さｌに正規化する。つまり

をそれぞれ演算する。このｖ_ｐと円錐面ベクトルの内積はこれらベクトルのなす角度の余弦とする。つまり

関係が成立する。知りたいのは共通直線５_ｉの方向のみであるから円錐面ベクトルｕは単位ベクトル、つまり‖ｕ‖＝１とする。全ての円錐面と共通する直線の方向を求めるには

として次式の連立方程式をｕについて解けば良い。

一般にはこの連立方程式は解が存在しないか、もしくは一意には決まらない。そこで‖Ｖｕ−ｃ＾（ω）‖を最小化するｕを求めて前記連立方程式の解、つまり共通直線５_ｉの方向ｕ_ｉ（ω）とする。この誤差を最小化するｕを求める計算が演算部１６ｂで行われる。この方向ｕ_ｉ（ω）は３次元での方向であるから、極座標表示により方位角θ_ｉ（ω）及び仰角φ_ｉ（ω）が求まることになる。
もっと計算を簡便にするには次のようにしてもよい。図１６に示すように正規化軸ベクトル計算部１６ａで円錐面推定に用いた各センサ対についての正規化軸ベクトルｖ_ｐ（ｐ＝１，…，Ｐ）を求め、逆行列計算部１６ｃでＶ＝（ｖ_１，…，ｖ_Ｐ）^ＴのＭｏｏｒｅ−Ｐｅｎｒｏｓｅ型一般逆行列Ｖ^＋を計算し、このＶ^＋と推定円錐面と対応する余弦ベクトルｃ＾（ω）＝（ｃｏｓθ＾_ｊｊ’（ω，１）…ｃｏｓθ＾_ｊｊ’（ω，Ｐ））^Ｔとを用いて最小ノルムな最小２乗誤差解を求め、大きさを正規化して近似解としてもよい。つまり演算部１６ｄで

を計算する。
このようにして複数の推定円錐面の共通直線とみなされる直線の方向が各周波数ごとに、かつ各信号源ごとに求められる。
図１２中のパーミュテーション解決部１７は、到来方向決定部１６で決定された到来方向ｕ_ｉ＝（θ_ｉ，φ_ｉ）を用いて、分離行列算出部１２で算出された分離行列Ｗ（ω）の行の並び替えを行ってパーミュテーション問題が解決された分離行列を生成する。
パーミュテーション解決部１７で行われる具体例としては、到来方位角θ_ｉについて以下のように並べ替えを行い、解決できなかった列について同様に到来仰角φ_ｉによる並べ替えを行う。つまり算出決定された到来方位角θ_ｉ（ω）がどの周波数においても既定の順序、例えば、θ_１，θ_２，…，θ_Ｉが昇順となるように逆行列Ｈ（ω）の列を入れ替え、入れ替えができなかった列について仰角φ_ｉがどの周波数でも昇順になるように逆行列Ｈ（ω）の列を入れ替える並替行列が並替行列生成部１７ａで生成され、その並替行列の逆行列Ｐ（ω）が逆行列生成部１７ｂで生成され、並替部１７ｃでその逆行列Ｐ（ω）が、分離行列Ｗ（ω）に左から掛算される。並替行列生成部及び逆行列生成部１７ｂはパーミュテーション行列生成部を構成している。
並替行列生成部１７ａでの処理を具体的に述べる。この例では逆行列Ｈ（ω）の第１列目で（θ_１（ω），φ_１（ω））が、第２列目で（θ_２（ω），φ_２（ω））が，…、第Ｉ列目で（θ_Ｉ（ω），φ_Ｉ（ω））がそれぞれ式（９’）に基づき計算されたとする。いま到来方向決定部１６１から入力された到来方向を、前記昇順に並べたものと区別するためにθおよびφの各添字にダッシュ「’」を付けて、（θ_１’（ω），φ_１’（ω）），（θ_２’（ω），φ_２’（ω）），…，（θ_Ｉ’（ω），φ_Ｉ’（ω））と表わし、これらを例えばまずθ_ｉ’（ω）について昇順に並べた時に、例えば（θ_３’（ω），φ_３’（ω））＞（θ_Ｉ’（ω），φ_Ｉ’（ω））＞（θ_２’（ω），φ_２’（ω））＞…であれば、逆行列Ｈ（ω）の３列目を１列目に、Ｉ列目を２列目に、２列目を３列目になるように移動し、残りの列も同様に移動し、θ_ｉ’（ω）が同一のものについてはφ_ｉ’（ω）が昇順になるように列を移動する。そのような列の移動並べ替えを行う並替行列を作る。このような並べ替えを行う行列を作ることは従来からよく行われていることである。このようにして全ての周波数について得られた到来方向（θ_１（ω），φ_Ｉ（ω）），…，（θ_Ｉ（ω），φ_２（ω））を用いて、並替行列が作られ、更にその逆行列つまりパーミュテーション行列Ｐ（ω）が算出される（図１３、ステップＳ１６）。
この算出したパーミュテーション行列Ｐ（ω）が並替部（１７ｃ）で分離行列Ｗ（ω）の左から掛けられ、その結果Ｗ’（ω）＝Ｐ（ω）Ｗ（ω）がパーミュテーション問題が解決された分離行列として出力される（ステップＳ１７）。つまりこの分離行列Ｗ’（ω）はいずれの周波数のものでも、１行目は信号源２_１の信号を分離する要素、２行目は信号源２_２の信号を分離する要素、以下同様に同一行の要素は同一信号源からの信号を分離する要素となる。
分離行列Ｗ’（ω）は、時間領域変換部１８で、例えば、逆フーリエ変換によって、時間領域の分離フィルタ係数群

に変換され、これら分離フィルタ係数群は信号分離部１９に設定される。
信号分離部１９は、各センサからの観測信号ｘ_１（ｔ），…，ｘ_Ｊ（ｔ）と分離フィルタ係数群とにより式（８）の演算を行い分離信号ｙ_１（ｔ），…，ｙ_Ｉ（ｔ）を出力する。
図１２中に破線で示すように、パーミュテーション解決部１７よりの並替えされた分離行列Ｗ’（ω）と周波数領域変換部１１よりの周波数領域観測信号Ｘ（ω，ｍ）とにより式（５）の演算を周波数領域分離信号生成部１９’で行い、その結果の周波数領域分離信号Ｙ（ω，ｍ）＝Ｗ’（ω）Ｘ（ω，ｍ）を時間領域信号変換部１８’で時間領域信号ｙ_１（ｔ），…，ｙ_Ｉ（ｔ）を生成してもよい。また到来方向決定部１６で決定された各周波数ごとの到来方向（θ_１（ω），φ_１（ω）），…，（θ_Ｉ（ω），φ_Ｉ（ω））を第１実施形態と同様に、つまり図５中のソート部３２で方向範囲ごとに分別し、各分別されたものごとの全周波数における到来方向を統合部３３でそれぞれ統合してもよい。
以上説明したように、この実施形態においても指向特性パターンの低ゲイン方向を探索することなく、式（９’）の演算により円錐面情報（θ＾_{ｉ，ｊｊ’}（ω），φ＾_{ｉ， ｊｊ’}（ω））を推定しているため計算量が少ない。しかも同一信号源に対し、円錐面を複数推定してその共通直線から信号の到来方向を算出しているため、０°から３６０°までのどの方向に信号源が存在しても、信号源３_ｉの方向を一意に推定することができる。またその推定方向をパーミュテーション行列Ｐ（ω）の決定に利用しているため、信号源の位置に関わらず、パーミュテーション問題を適切に解決することができる。
〔第３実施形態〕
第３実施形態は、位置情報として一対のセンサと一つの信号源との各距離の比に基づくその信号源が存在する曲面を用いる。第１及び第２実施形態においては各信号源はセンサから遠方にあって、信号源からの信号はセンサには平面波として到来することを想定して処理した。しかし信号源とセンサとの距離が比較的短かい場合はセンサには信号が球面波として到来する。このような点から混合行列Ａ（ω）の要素の比Ａ_ｊｉ（ω）／Ａ_ｊ’ｉ（ω）を球面波（近距離場）モデルによって解釈すると、信号源ｉの方向以外の情報を推定することができる。
すなわち、近距離場モデルを用いると、周波数応答Ａ_ｊｉ（ω）は以下のように表わせる。

ここで、ｑ_ｉは信号源ｉの位置を示すベクトルである。
このように表現された周波数応答について混合行列の同一列の２つの要素比Ａ_ｊｉ（ω）／Ａ_ｊ’ｉ（ω）をとり、その絶対値を計算すると次式となる。

なお、｜β｜はβの絶対値を表わす。
この式（１０）を満すｑ_ｉの点の無数の集合は信号源ｉが存在する曲面を与え、遠距離場（平面波）モデルを用いて推定された方向（又は円錐面）とあわせて用いることにより、センサから信号源ｉまでの距離を推定することができる。これにより、２つ以上の信号源が同一方向もしくは互いに近い方向にある場合であっても、距離が違っていれば、それらを区別でき、パーミュテーション問題を適切に解決することが可能となる。
この第３実施形態の例では、この信号源が存在する曲面を表わす位置情報と前記方向情報とを用いて分離行列に対するパーミュテーション問題を解決する。
第３実施形態の機能構成を図１７に、その処理手順を図１８にそれぞれ示す。センサは３個以上が２次元又は３次元に配置されるがこの実施形態においては例えば図２０に示すようにセンサ１_１と１_２の間隔に対しセンサ１_２と１_３の間隔は１０〜２０倍好ましくは１５倍程度とされる。先の実施形態と同様、各観測信号ｘ_１（ｔ），…，Ｘ_Ｊ（ｔ）は、それぞれ周波数領域の信号Ｘ_１（ω，ｍ），…，Ｘ_Ｊ（ω，ｍ）に変換され（ステップＳ１１）、更に独立成分分析により、各周波数ごとに分離行列Ｗ（ω）が算出される（ステップＳ１２）、その各分離行列Ｗ（ω）の逆行列として行列Ｈ（ω）が算出される（ステップＳ１３）。またこの例では第２実施形態と同様に周波数ごとの逆行列Ｈ（ω）の各列から選択した要素対を用い、円錐面が一つ好ましくは複数推定される（ステップＳ１４）。この第３実施形態においては距離比算出部３１において、周波数ごとの逆行列Ｈ（ω）から列ごとに選択した要素対を用いその対応するセンサと１つの信号源ｉとの距離比、つまり式（１０）と式（２１）から次式（１０’）を算出する（ステップＳ３５）。

距離比算出部３１において行われるステップＳ３５の処理の具体例を図１９を参照して説明する。この処理は図１４の処理の流れとほぼ同様である。パラメータｉを０に初期化し（ステップＳ２０）、次に、ｉを＋１し（ステップＳ２１）、パラメータｊ，ｊ’としてＪ以下の自然数を、例えばランダムに選択しｊ≠ｊ’、かつ１度選択した組は選択しない（ステップＳ２３）。センサｊの位置ベクトルｄ_ｊと、センサｊ’の位置ベクトルｄ_ｊ’とを抽出し（ステップＳ２４）、逆行列Ｈ（ω）のｉ列目から、要素Ｈ_ｊｉ（ω）と、Ｈ_ｊ’ｉ（ω）とを、選択する（ステップＳ２５）。
この実施形態ではこれら選択した２つの要素の比ＤＲ_{ｉ，ｊｊ’}（ω）を計算する（ステップＳ４１）。次に、ｉ＝Ｉであるか否かを判断し（ステップＳ２９）。ｉ＝ＩでなければステップＳ２１に戻り、ｉ＝２であれば処理を終了する。
距離比算出部３１で算出された距離比情報ＤＲ_{ｉ，ｊｊ’}（ω）は、パーミュテーション解決部１７に送られ、パーミュテーション解決部１７は、到来方向決定部１６で推定された方向情報ｕ_ｉ（ω）と、距離比算出部３１で算出された距離比情報ＤＲ_{ｉ，ｊｊ’}（ω）とを用いて、分離行列算出部１２で算出された分離行列を行ってパーミュテーション問題を解決する。
Ｗ（ω）の行替えを行ってパーミュテーションの問題を解決する。例えば方向情報と距離比情報とから信号源２_ｉの距離‖ｑ_ｉ（ω）‖を距離推定部１７ｄで計算する（ステップＳ３６）。
この距離‖ｑ_ｉ（ω）‖の計算方法を図２０を参照して説明する。信号源２_１及び２_２がセンサ１_１，１_２からみて同一方向Ｂにある。この場合センサ１_１，１_２及び遠距離場モデルで推定される信号源２_１及び２_２の方向ｕ_１，ｕ_２は同一の直線を定める。一方、間隔が大きいセンサ１_２，１_３及び近距離場モデルを用いて信号源２_１が存在すると曲面６_１は距離比

から推定でき、つまり‖ｑ_１（ω）‖が推定できる。また信号源２_２が存在する曲面６_２は距離比

から推定でき、つまり‖ｑ_２（ω）‖が推定できる。
信号源２_１の位置は、直線ｕ_１＝ｕ_２と、曲面６_１との共通部分上にあると推定でき、信号源２_２の位置は、直線ｕ_１＝ｕ_２と、曲面６_２との共通部分上にあると推定できる。例えば直線ｕ_１＝ｕ_２を表す式と曲面６_１及び６_２をそれぞれ表わす式との連立方程式を解くことにより‖ｑ_１（ω）‖及び‖ｑ_２（ω）‖を求めればよい。このようにして信号源の方向が同一や近似する場合であっても、信号源位置を区別することができる。
パーミュテーション解決部１７は、このようにして得られた各周波数ごとの信号源の距離‖ｑ_ｉ（ω）‖が、所定の順序が例えば昇順になるように、分離行列Ｗ（ω）の行の入れ替えを行う。つまりこの入れ替えを行うためパーミュテーション行列Ｐ（ω）を算出する（ステップＳ３７）。そのパーミュテーション行列Ｐ（ω）の算出は実施形態２中のパーミュテーション解決部で行ったパーミュテーション行列の算出と同様に行えばよい。算出したパーミュテーション行列Ｐ（ω）を分離行列Ｗ（ω）の左から掛け、その結果Ｗ’（ω）＝Ｐ（ω）Ｗ（ω）を分離行列として出力する（ステップＳ３８）。
出力された分離行列Ｗ’（ω）は、時間領域変換部１８に送られ、時間領域変換部１８での信号分離に用いられる。
図２０から理解されるように、センサ１_２と１_３についてみると、信号源２_１及びセンサ１_２間の距離と信号源２_１及びセンサ１_３間の距離との差が大きく、一方信号源２_２及びセンサ１_２間の距離と信号源２_２及びセンサ１_３間の距離との差が小さい。従ってＤＲ_１，２３＝‖ｑ_１−ｄ_３‖／‖ｑ_１−ｄ_２‖と数値１との差の絶対値は比較的大きく、ＤＲ_２，２３＝‖ｑ_２−ｄ_３‖／‖ｑ_２−ｄ_２‖と数値１との差の絶対値は小さい。センサ１_２と１_３の間隔が大きい程、距離比ＤＲ_１，２３（ω）とＤＲ_２，２３（ω）との差が大きくなる。このため間隔が大きいセンサ対、この例では１_２と１_３が設けられる。
パーミュテーション解決部１７において、図１８中のステップＳ３７中に括弧で示すように計算した距離比ＤＲ_{ｉ，ｊｊ’}（ω）を用いて、すべての周波数において逆行列Ｈ（ω）の１列目，２列目，３列目，…，Ｉ列目の順に、求まる距離比ＤＲ_{ｉ，ｊｊ’}（ω）が昇順になるように並替行列を生成し、更に前記パーミュテーション行列Ｐ（ω）を生成してもよい。この場合、図１７中の距離推定部１７ｄは省略される。また到来方向決定部１６で方向情報ｕ_ｉ（ω）を計算できなかったり、逆行列Ｈ（ω）の複数列について同一のｕ_ｉ（ω）が計算された信号源ｉとセンサｉｊ’についてのみ、距離比ＤＲ_{ｉ，ｊｊ’}（ω）を計算し、これら自体により、あるいは更に求めた‖ｑ_ｉ（ω）‖により分離行列Ｗ（ω）の行の入れ替えを行ってもよい。つまり、方向情報ｕ_ｉ（ω）を用いて分離行列Ｗ（ω）の行の入れ替えを行い、その後、ｕ_ｉ（ω）によって入れ替えができなかった行については距離比ＤＲ_{ｉ，ｊｊ’}（ω）又は‖ｑ_ｉ（ω）‖を用いてＷ（ω）の行の入れ替えを更に行ってもよい。実際には両入れ替えを同時に行うパーミュテーション行列Ｐ（ω）を生成する。
距離比ＤＲ_{ｉ，ｊｊ’}（ω）をまず求め、分離行列Ｗ（ω）の行の入れ替えを行い、距離比ＤＲ_{ｉ，ｊｊ’}（ω）を求めることができなかったものについて方向情報ｕ_ｉ（ω）を用いて分離行列Ｗ（ω）の入れ替えを更に行ってもよい。この場合も両入れ替えを同時に行うパーミュテーション行列Ｐ（ω）を生成する。一般にはＤＲ_{ｉ，ｊｊ’}（ω）よりｕ_ｉ（ω）の方が高い精度で得られるため、ｕ_ｉ（ω）による入れ替えを主とし、できないものについてＤＲ_{ｉ，ｊｊ’}（ω）による入れ替えを行う方がよい。
また式（１０’）をその右辺の値を距離比ＤＲ_{ｉ，ｊｊ’}（ω）と表わして、ｑ_ｉについて解くと、中心Ｏ_{ｉ，ｊｊ’}（ω）及び半径Ｒ_{ｉ，ｊｊ’}（ω）がそれぞれ下記式（１１）及び（１２）で与えられる球面となる。

例えばセンサ１_ｊ，１_ｊ’の各位置がｄ_ｊ＝（０，０．１５，０），ｄ_ｊ’＝（０，−０．１５，０）（単位はメートル）である場合、半径Ｒ_{ｉ，ｊｊ’}をパラメータとしたとき、式（９）によって決まる球面の様子は図２１に示すようになる。
つまり信号源ｉは位置情報としての式（１１）及び（１２）で与えられる球面上に存在することになる。従って図１７中のパーミュテーション解決部１７中の距離推定部１７ｄを括弧書きで示すように曲面推定部とし、方向情報ｕ_ｉ（ω）を求めることができなかった各ｉ，ｊｊ’について曲面推定部１７ｄで、式（１１）の半径Ｒ_{ｉ，ｊｊ’}（ω）及び中心Ｏ_{ｉ，ｊｊ’}（ω）を、図１８中のステップＳ３６中に括弧書きで示すようにそれぞれ計算し、これらＲ_{ｉ，ｊｊ’}（ω）及び中心Ｏ_{ｉ，ｊｊ’}（ω）がどの周波数の逆行列Ｈ（ω）でも同一順になるようにして、パーミュテーション行列Ｐ（ω）で求めてもよい。
この円錐面情報θ＾_{ｉ，ｊｊ’}（ω）又は方向情報ｕ_ｉ（ω）と球面情報ＤＲ_{ｉ，ｊｊ’}（ω）又はＲ_{ｉ，ｊｊ’}（ω）を用いてもパーミュテーション行列Ｐ（ω）を求めることができない場合はその周波数について、従来の相関法（例えば、Ｈ．Ｓａｗａｄａ他“Ａ ｒｏｂｕｓｔ ａｎｄ ｐｒｅｃｉｓｅ ｍｅｔｈｏｄ ｆｏｒ ｓｏｌｖｉｎｇ ｔｈｅ ｐｅｒｍｕｔａｔｉｏｎ ｐｒｏｂｌｅｍ ｏｆ ｆｒｅｑｕｅｎｃｙ−ｄｏｍａｉｎ ｂｌｉｎｄ ｓｏｕｒｃｅ ｓｅｐａｒａｔｉｏｎ．”，ｉｎ Ｐｒｏｃ．Ｉｎｔｌ．Ｓｙｍｐ．ｏｎ Ｉｎｄｅｐｅｎｄｅｎｔ Ｃｏｍｐｏｎｅｎｔ Ａｎａｌｙｓｉｓ ａｎｄ Ｂｌｉｎｄ Ｓｉｇｎａｌ Ｓｅｐａｒａｔｉｏｎ（ＩＣＡ ２００３），２００３，ｐｐ．５０５−５１０参照）を適用するとよい。図１７中に破線で示すように、信号分離部１９で分離された出力信号ｙ_１（ｔ），…，ｙ_Ｉ（ｔ）を周波数領域変換部３３で例えば周波数領域変換部１１と同様の方法で周波数領域の信号Ｙ_１（ω，ｍ），…，Ｙ_Ｉ（ω，ｍ）に変換し、相関算出部３４において、これら周波数領域信号Ｙ_１（ω，ｍ），…，Ｙ_Ｉ（ω，ｍ）中のパーミュテーション解決部１７でパーミュテーション行列を作ることができなかった周波数成分ｆ_ａｎｋ成分と、パーミュテーション行列が得られた周波数中のｆ_ａｎｋと隣接した周波数成分との相関を計算し、パーミュテーション解決部１７ではこの相関が大きくなるように周波数ｆ_ａｎｋの分離行列の行の入れ替えを行い、この行入れ替えした分離行列に基づく分離出力信号ｙ_１（ｔ），…，ｙ_Ｉ（ｔ）中のｆ_ａｎｋ成分について再び相関計算部３４で相関を計算し、同様のことを計算した相関が最大になるまで繰り返す。この相関の最大値が所定値に達しなかった場合は、更にそのｆ_ａｎｋの成分と、パーミュテーション行列が得られている周波数成分中のｆ_ａｎｋに対する基本波又は高調波（一般には高調波）関係にあるものとの各相関の和を相関計算部３４で計算し、この相関和が大きくなるように、ｆ_ａｎｋに対する分離行列の行入れ替えを修正部１７ｅで行うことを、前記相関和が最大になるまで繰り返す。なお第２実施形態で説明したように周波数領域の信号Ｘ（ω，ｍ）に対し分離行列Ｗ（ω）を適用し、周波数領域の分離信号Ｙ（ω，ｍ）を求める場合は、このＹ（ω，ｍ）を相関計算部３４の計算に用いればよい。
以下にこの第３実施形態の実験例を信号源として、室内で実測したインパルス応答を畳み込んで作成した混合音声を用いて分離実験を行った場合について説明する。図２２に示すように残響時間が１３０ミリ秒、３５５ｃｍ×４４５ｃｍ×２５０ｃｍ（高さ）の室内の水平面内でほぼ中心部で高さ１３５ｃｍの所にセンサとしての無指向性マイクロホン１_１〜１_４を直径４ｃｍの円上に等間隔で配し（これらマイクロホン１_１〜１_４は図中の左上に拡大して示した）、更に無指向性マイクロホン１_５〜１_８を直径３０ｃｍの円上に等間隔かつ、マイクロホン１_１，１_３，１_５及び１_７が、またマイクロホン１_２，１_４，１_６及び１_８がそれぞれ一直線上に配列するように配置した。マイクロホン１_１から１_５への方向を基準（０°）とし、マイクロホン配置の中心を原点とし、原点から１２０ｃｍ離れ、−３０°方向に音源２_１を、＋３０°方向に音源２_２を、＋９０°方向に音源２_３を、−１５０°方向に音源２_６をそれぞれ配置し、１５０°方向で原点からの距離６０ｃｍに音源２_４を、距離１８０ｃｍに音源２_５をそれぞれ配置した。サンプリング周波数ｋＨｚ、データ長６秒、フレーム長２０４８サンプル（２５６ミリ秒）、フレームシフト５１２サンプル（６４ミリ秒）とした。
周波数領域での分離行列Ｗ（ω）の各行のうち、マイクロホン対１_１と１_３、１_２と１_４、１_１と１_２、１_２と１_３に相当する行を用いて音源方向（円錐面）推定を行った。推定したヒストグラムを図２２に示す。図２３の横軸は推定方向を、縦軸は信号数である。これより推定方向は５つの群（クラスタ）があり、そのうち１５０°付近のクラスタは他のクラスタより２倍の大きさである。これは６個の音源中の２個が同一方向（１５０°付近）から到来していると推測される。４つの音源については推定した方向をもとにパーミュテーションを解決した。その結果を図２４に示す。図２４の横軸は周波数、縦軸は方向（度）を表わす。
残りの２個の音源については、広い間隔のマイクロホン対１_７と１_５、１_７と１_８、１_６と１_５、１_６と１_８を用いて音源が存在する球面の半径によって到来信号を区別した。マイクロホン対１_７と１_５を用いて推定した球面の半径を図２５に示す。図２５の横軸は周波数、縦軸は半径（ｍ）を表わす。
残響や推定誤差の影響により、位置情報だけではパーミュテーションを完全には解決することはできない。従って推定した位置情報に基づいて信号が矛盾なく分類できた周波数はその情報に基づいてパーミュテーション行列を生成し、残りの周波数については相関に基づく方法を用いてパーミュテーションを解決した。最後に時間領域の分離フィルタ係数を求める際にスペクトルの平滑化を行った。スペクトル平滑化については例えばＨ．Ｓａｗａｄａ他、“Ｓｐｅｃｔｒａｌ ｓｍｏｏｔｈｉｎｇ ｆｏｒ ｆｒｅｑｕｅｎｃｙ−ｄｏｍａｉｎ ｂｌｉｎｄ ｓｏｕｒｃｅ ｓｅｐａｒａｔｉｏｎ”，ｉｎ Ｐｒｏｃ．ＩＷＡＥＮＣ ２００３，２００３，ｐｐ．３１１−３１４参照されたい。分離性能（ＳＩＲ）の評価結果を図２６に示す。図２６中の数値の単位はｄＢである。Ｃは相関法のみによってパーミュテーションを解決したもの、Ｄ＋Ｃは方向（円錐面）推定でパーミュテーションを解決し、解決できないものに相関法を適用したものＤ＋Ｓ＋Ｃは方向（円錐面）推定と球面推定とによりパーミュテーション解決を行った後、解決できなかった周波数について、相関法を用いたものである。この後者の方法によれば、６音源すべての分離ができ、ＳＩＲの改善は入力ＳＩＲから平均１７．１ｄＢとなった。
前記した第２実施形態の例ではセンサは２次元配置したが、一対のセンサにより推定する球面はこれらセンサの中心２等分線と対称に現われるため、信号源が３次元に存在している場合は２次元配置センサでは判別できなくなるため、センサも３次元配置とする必要がある。
以上説明したように、この第３実施形態においても式（９’）により円錐面情報を、また式（１０’）により曲面情報をそれぞれ推定しているため計算量が少ない。しかも円錐面と、距離比ＤＲ_{ｉ，ｊｊ’}（ω）、又は距離‖ｑ_ｉ（ω）‖あるいは球面の半径Ｒ_{ｉ，ｊｊ’}（ω）とによりパーミュテーションを解決しているために２つ以上の信号源が同一方向もしくは互いに近い方向にある場合であっても、それらを区別できる。更にこれに相関法を加えれば一層確実に分離が可能である。なお球面情報としては計算が簡単であるからＤＲ_{ｉ，ｊｊ’}（ω）が好ましい。
〔第４実施形態〕
この第４実施形態では、推定した円錐面の信頼性を検証し、信頼性が高いと判断した円錐面を分離行列のパーミュテーションの解決に用いる。例えば図２７に示すように円錐面推定部１４で推定された円錐面情報θ＾_{ｉ，ｊｊ’}（ω）は円錐面検証部４１で許容角情報格納部４２よりの許容角情報に基づき、信頼性があるか否かの検証がなされる。つまり式（９’）により求まる角度θ＾_{ｉ，ｊｊ’}（ω）はこれを求めるために用いたセンサ１_ｊ及び１_ｊ’の配列方向に対する相対角度であり、図４Ｂを参照して説明したように、角度θ＾_{ｉ，ｊｊ’}（ω）が０度付近及び１８０度付近の場合、正しいパーミュテーションが行われないと推定される。
従って、パーミュテーションを正しく行うことができると推定される角度の最小値θ_ｍｉｎと最大値θ_ｍａｘが許容角情報として許容角情報格納部４２に格納され、円錐面検証部４１で推定円錐面情報θ＾_{ｉ，ｊｊ’}（ω）がθ_ｍｉｎとθ_ｍａｘとの間にあれば信頼できる円錐面と判定されて出力され、つまりパーミュテーション解決に用いられるが、θ_ｍｉｎとθ_ｍａｘとの間になければ、そのθ＾_{ｉ，ｊｊ’}（ω）は信頼性がないものとして破棄され、パーミュテーション解決に用いられない。例えば図１５中の円錐面４_１３は破棄される。
円錐面検証部４１で信頼性があると検証された円錐面情報θ＾_{ｉ，ｊｊ’}（ω）は図１２中の到来方向決定部１６へ送られ、又は図１７で説明したように、パーミュテーション解決部１７へ直接送られる。図１３及び図１８の処理手順において、ステップＳ１４の次に破線で示すようにステップＳ１００として推定円錐面θ＾_{ｉ，ｊｊ’}（ω）が信頼性があるか否かの検証が行われ、信頼性があるものについてのみ次のステップに移るようにされる。式（９’）において計算された偏角ａｒｇ（Ｈ_ｊｉ／Ｈ_ｊ’ｉ）に含まれる誤差をΔａｒｇ＾、推定された角度θ＾_ｉに含まれる誤差をΔθ＾とすると、これらの比｜Δθ＾／ΔａｒｇＨ＾｜は式（９’）を偏微分して次式のように近似できる。

式（１３）をいくつかの周波数について計算した結果を図４Ｃに示す。この図４Ｃから推定された角度θ＾_ｉの方向がセンサ１_ｊ，１_ｊ’のセンサペア軸に近い場合はａｒｇ（Ｈ_ｊｉ／Ｈ_ｊ’ｉ）に含まれる誤差ΔａｒｇＨ＾は推定角度θ＾_ｉに対し、大きな誤差を生じさせることがわかる。つまりセンサペア軸に近い方向の推定角度θ＾_ｉを用いて分離行列Ｗ（ω）に対するパーミュテーションの問題を解決した場合、これは正しい解にならない可能性が高い。図４Ｂ及び図４Ｃから例えばθ_ｍｉｎは２０°程度、θ_ｍａｘは１６０°程度とされる。図４Ｃからわかるように｜Δθ＾／ΔａｒｇＨ＾｜は周波数により可成り異なり、低い周波数ではΔａｒｇＨ＾が到来方向の誤差に大きく影響する。従って、低周波数については、信号源が存在する球面の推定に基づく情報ＤＲ_{ｉ，ｊｊ’}、‖ｑ_ｉ‖、Ｒ_{ｉ，ｊｊ’}を用いて、あるいは相関法によりパーミュテーション問題を解決するとよい。
この第４実施形態によれば、推定円錐面情報中の信頼性のないものは破棄されるため、これにより悪影響を受けることなく、それだけ正確に到来方向を推定でき、従って正しいパーミュテーション行列Ｐ（ω）を生成することができ、つまり分離信号のＳＩＲ（性能）が向上する。
〔第５実施形態〕
第５実施形態では位置情報として距離比、又はこれより推定した球面情報を用いる。この機能構成は図１７中に示されている。図２８にこの実施形態の処理手順を簡略に示す。この場合はセンサは比較的広い間隔、例えば図２２に示した信号源が音源の場合、３０ｃｍとして少くとも２次元に配置される。
先の実施形態と同様に時間領域の観測信号は図１７中に示すように周波数領域の信号に変換され、分離行列生成部１２により分離行列Ｗ（ω）が生成され（ステップＳ５１）、これより逆行列生成部１３で逆行列Ｈ（ω）が生成される（ステップＳ５２）。各周波数における逆行列Ｈ（ω）の各列ごとに球面情報が推定される（ステップＳ５３）。この球面情報は第３実施形態で説明した方法と同様に算出される。つまり球面情報は距離算出部３１により算出された距離比ＤＲ_{ｉ，ｊｊ’}（ω）、あるいは距離推定部又は曲面推定部１７ｄで算出された‖ｑ_ｉ‖、又は半径Ｒ_{ｉ，ｊｊ’}（ω）及び中心Ｏ_{ｉ，ｊｊ’}（ω）である。
これら球面情報を用いて、その配列順が予め決めた順になるよう分離行列Ｗ（ω）に対するパーミュテーション行列が生成され分離行列Ｗ（ω）の行の入れ替えが行われる（ステップＳ５４）。第３実施形態におけるこの処理はパーミュテーション解決部１７で行われるが、ただしこの第５実施形態では球面情報のみが用いられる。この球面情報のみでパーミュテーションを解決できない周波数があれば、その周波数に対しては前述した相関法による解決が行われる（ステップＳ５５）。
図２２に示した室内に配したマイクロホン１_６及び１_８を用い、１２０°方向において原点から６０ｃｍ及び１５０ｃｍの所に音源２_４及び２_５を配し、実験室内で測定したインパルス応答に話者４名の音声を畳み込んで作成した１２通りの組合せの混合音声と対する分離実験を行った。パーミュテーション解決は推定された半径Ｒ_４，６８（ω）とＲ_５，６８（ω）を比較してＲ_４，６８（ω）≦Ｒ_５，６８（ω）となるようにして行った。Ｒ_４，６８及びＲ_５，６８中の最大値が、Ｒ_４，６８及びＲ_５，６８中の最小値にしきい値Ａｔｈを乗算した値以上、つまりｍａｘ（Ｒ_４，６８，Ｒ_５，６８）≧Ａｔｈ・ｍｉｎ（Ｒ_４，６８，Ｒ_５，６８）を満す周波数は音源位置（Ｒ_４，６８（ω），Ｒ_５，６８（ω））によるパーミュテーションの解決が信頼できると判断して、音源位置による方法を適用し、それ以外の周波数では相関による方法を適用した。従ってＡｔｈ＝１．０のときは全て音源位置による方法を適用しＡｔｈが無限大のときは全て相関による方法を適用することになる。しきい値Ａｔｈを変化させながら分離性能としてＳＩＲ（信号干渉比）を音声の組合せごとにプロットした結果を図２９に示す。この結果より音源位置のみによる方法では全体的に性能が悪く、相関のみによる方法では分離性能がばらつき、両者を併用すると安定して高い性能が得られることがわかる。またしきい値Ａｔｈとしては比較的大きい値が好ましく、８〜１６で全体の１／５〜１／１０程度の周波数を音源位置によって決定できることがわかる。
また図２２に示した実験条件で音源位置による方法と相関による方法を併用した実験結果を図２６中のＤ＋Ｃの欄に示す。これよりこの併用方法によれば可成り性能よく分離が行われることがわかる。この結果は式（９’）を用いる到来方向（円錐面情報）による方法と相関による方法とを併用した場合と同様の傾向がある。なお、到来方向による方法と相関による方法とを併用した場合の実験結果は、相関による方法の説明で引用した文献に示されている。
この第５実施形態によれば式（１０’）により球面情報を求めているため計算量が少ない。なお球面情報としては距離比ＤＲ_{ｉ，ｊｊ’}（ω）が好ましい。
［第６実施形態］
第６実施形態は例えば一つの推定円錐面に基づき、パーミュテーションの問題を解決する。図１２中に破線で示すように円錐面推定部１４で推定された円錐面θ＾_{ｉ，ｊｊ’}（ω_１）、・・・ θ＾_{ｉ，ｊｊ’}（ω_Ｎ）、はパーミュテーション解決部１７へ直接入力され、図１３中に破線で示すようにステップｓ１４で設定された円錐面θ＾_{ｉ，ｊｊ’}（ω_ｎ）（ｎ＝１、・・・、Ｎ）がステップＳ１５の処理をすることなくステップＳ１６でどの周波数においても例えばθ＾_{ｉ，ｊｊ’}（ω_ｎ）が昇順になるようなパーミュテーション行列Ｐ（ω）が算出される。
この場合、ステップＳ１４で信号源ｉごとに１つの円錐面を推定するのみでもよい。また図に示していないが、行入れ替えができなかったものについては、第３実施形態で説明したように相関法による行入れ替えを得るようにしてもよい。
この第６実施形態によれば逆行列Ｈ（ω）の各列ごとの二つの要素比をとることによりスケーリングの問題が簡単に除去され、式（９’）の演算をすればよく計算時間が短くて済む。
［信号分離のまとめ］
以上述べたブラインド信号分離におけるパーミュテーションを解決した分離行列の生成方法を、逆行列の生成以後の処理を図３０Ａ、図３０Ｂ、図３０Ｃに示す。図３０Ａでは逆行列の各列の要素比から円錐面を推定し（ステップＳ６１）、必要に応じて信頼性のない円錐面を破棄し（ステップＳ６２）、これらの複数の円錐面から共通直線の方向を決定し（ステップＳ６３）、共通直線方向を用いてパーミュテーション行列Ｐ（ω）を生成して分離行列の入れ替えを行い（ステップＳ６４）、パーミュテーションを解決できなかった周波数に対しては相関法による分離行列の行入れ替えを行う（ステップＳ６５）。図３０Ａ中に破線で示すように、ステップＳ６からステップＳ６４に直ちに移り、円錐面を直接用いてパーミュテーション行列Ｐ（ω）を生成してもよい。
図３０Ｂでは前記ステップＳ６１の円錐面推定、必要に応じて前記ステップＳ６２の円錐面破棄を行い、ステップＳ６３の共通直線方向推定を行い、これら共通直線方向を用いて逆行列の列の入れ替えを行う（ステップＳ６６）。この際、図中破線で示すようにステップＳ６１で推定された円錐面又はステップＳ６２で処理後の円錐面を直接、ステップＳ６６の処理に用いてもよい。次に円錐面の推定又は共通直線方向の決定ができなかった、又は円錐面或いは共通直線方向が不確かなものと対応する、或いは同一値となったものについて、その逆行列の列の要素比から球面情報を推定し（ステップＳ６７）、この推定した球面情報を用いて逆行列の入れ替えを行ってパーミュテーション行列Ｐ（ω）を生成して分離行列の行入れ替えを行う（ステップＳ６４）。更にこのパーミュテーション行列Ｐ（ω）を作ることができなかった周波数に対しては相関による方法を適用する（ステップＳ６５）。
図３０Ｃはまず逆行列の各列の要素比から球面情報を推定し（ステップＳ６８）、この球面情報を用いて逆行列の列の入れ替えを行い（ステップＳ６９）、この列の入れ替えを行うことができなかった、又は球面情報が不確かなものについて、その逆行列の列の要素比から円錐面を推定し（ステップＳ７０）、これらの複数の円錐面の共通直線方向を決定し（ステップＳ７１）、この決定した方向又はステップＳ７０で推定した円錐面を直接用いて、逆行列の列の入れ替えを行ってパーミュテーション行列Ｐ（ω）を生成し、分離行列の入れ替えを行う（ステップＳ６４）。パーミュテーション行列を作ることができなかった周波数に対して相関法を適用する（ステップＳ６５）。
これらに対し更に図２８に示した方法がある。
第３実施形態、第５実施形態及び第６実施形態においても、第２実施形態で説明したように、周波数領域の分離行列Ｗ（ω）と観測信号Ｘ（ω）とを用いて信号分離を行い、その後、分離された周波数領域信号Ｙ（ω）を時間領域信号ｙ（ｔ）に変換してもよい。
第１実施形態の説明では２次元における信号の到来方向の推定を行ったが、第２実施形態中でも説明したように、３次元における信号の到来方向の推定にも適用できる。また第２乃至第６実施形態を２次元における信号分離に適用することもできる。この場合は推定円錐面θ＾_{ｉ，ｊｊ’}（ω）はその推定に用いたセンサ対のセンサペア軸に対し対称の２方向となり、推定球面Ｒ_{ｉ，ｊｊ’}（ω）又はＤＲ_{ｉ，ｊｊ’}（ω）は円の半径又はこれとほぼ対応するものになる。
図５、図１２、図１７にそれぞれ示した装置、また第５実施形態による装置はそれぞれコンピュータに機能させることもできる。この場合、それぞれと対応した処理手順、つまり図６、図１３、図１８などに示した各過程をコンピュータに実行させるためのプログラムを記録した磁気ディスク、半導体メモリ、ＣＤ−ＲＯＭなどから、コンピュータ内のメモリにインストールし、又は通信回線を介してコンピュータ内のメモリに前記プログラムをダウンロードして、コンピュータにそのプログラムを実行させればよい。

Claims (22)

1. Ｉ箇所の信号源から放射された信号を、Ｊ個のセンサで検出して上記信号源の位置情報を求める信号源位置情報推定装置であって、Ｉは２以上の整数、ＪはＩ以上の整数であり、
   上記各センサの観測信号を周波数領域の信号に変換する周波数領域変換手段と、
   上記周波数領域の信号から、独立成分分析により周波数ごとに、信号源信号を分離する第１分離行列をそれぞれ算出する分離行列算出手段と、
   上記各第１分離行列の逆行列をあるいは擬似逆行列「以下、両者を単に逆行列という」をそれぞれ算出する逆行列算出手段と、
   上記周波数ごとの逆行列の少なくとも１つに対し、その各列ごとにその二つの要素比から上記信号源の一つの位置情報を計算する位置情報計算手段とを具備する。
2. 請求の範囲第１項記載の装置において、
   上記位置情報計算手段は、上記周波数ごとの逆行列の複数について、上記各列ごとの要素比からの位置情報計算により上記各信号源の位置情報を求める手段であり、
   上記周波数ごとの信号源ごとの位置情報から、上記複数の逆行列と対応する周波数の分離行列における位置情報と対応する行が予め決めた順序になるように行の入れ替えを行うパーミュテーション行列を生成するパーミュテーション行列生成手段と、
   上記パーミュテーション行列と上記第１の分離行列とを乗算して行を入れ替えた第２の分離行列を求める並替手段とを備える。
3. 請求の範囲第２項記載の装置において、
   上記Ｊは３以上であり、かつＪ個のセンサは少なくとも２次元に配置され、
   上記位置情報は上記センサから上記信号源への方向を含み、その信号源が存在する円錐面であり、
   上記位置情報計算手段は、各列ごとの上記要素比から上記円錐面の計算を異なる二つの要素組の複数について行う手段と、周波数ごとの上記複数の円錐面の共通直線の方向を、上記位置情報として推定する到来方向決定手段を含む。
4. 請求の範囲第２項記載の装置において、
   上記Ｊは３以上であり、上記位置情報は一対のセンサから上記信号源への方向を含み、その信号源が存在する円錐面及びその信号源が存在する曲面であり、
   上記位置情報計算手段は、上記二つの要素の比から上記円錐面を計算する手段と、上記円錐面の計算に用いた二つの要素と対応するセンサ間隔よりも大きいセンサ間隔のセンサと対応する二つの要素の比からこれらセンサと上記信号線との距離比を計算する手段と、上記距離比から上記曲面を計算する手段を備え、
   上記パーミュテーション行列生成手段は周波数ごとに上記円錐面及び上記曲面に基づいて上記パーミュテーション行列を生成する手段である。
5. 請求の範囲第２項、第３項又は第４項に記載の装置において、
   上記Ｊは３以上であり、上記Ｊ個のセンサは少なくとも２次元に配置され、
   上記位置情報は上記センサから上記信号源への方向を含み、その信号源が存在する円錐面であり、
   上記円錐面を表わす角度が予め決めた第１角度と第２角度との間であるか否かを判定し、これら第１角度及び第２角度の間の円錐面を有効とする判定手段を備える。
6. 請求の範囲第２項記載の装置において、
   上記Ｊは３以上であり、上記Ｊ個のセンサは少なくとも２次元に配置され、
   上記位置情報は信号源が存在する球面の半径であり、
   上記位置情報計算手段は上記二つの要素比から距離比を計算する手段である。
7. 請求の範囲第２項に記載の装置において、
   上記位置情報は上記センサから上記信号源への方向を含みその信号源が存在する円錐面である。
8. 請求の範囲第２項、第３項、第４項、第６項又は第７項に記載の装置において、
   上記観測信号に対し、上記第２分離行列を用いて分離した信号の周波数領域の信号中の、上記パーミュテーション行列生成手段で、パーミュテーション行列が得られない周波数成分と、パーミュテーション行列が得られた周波数成分との相関を計算する相関計算手段と、
   上記相関が大きくなるように上記パーミュテーション行列が得られない周波数の分離行列の行を入れ替えるパーミュテーション行列を生成する修正手段とを備える。
9. 請求の範囲第１項記載の装置において、
   上記位置情報は上記センサから上記信号源への方向情報であり、
   上記位置情報計算手段は、上記比の偏角を計算する手段と、単位距離当たりの位相回転量と上記二つの要素とそれぞれ対応する上記センサ間の距離との積を計算する手段と、上記偏角を上記積の結果で除算する手段と、その除算結果の逆コサインを計算して上記方向情報を出力する手段とを備える。
10. 請求の範囲第９項記載の装置において、
    上記位置情報計算手段は、一つの信号源の上記方向情報を、周波数ごとに計算する手段であり、上記各一つの信号源ごとに計算された周波数ごとの方向情報を統合して方向情報を確定する統合手段を備える。
11. Ｉ箇所の信号源から放射された信号を、Ｊ個のセンサで検出して上記信号源の位置情報を求める信号源位置情報推定方法であって、Ｉは２以上の整数、ＪはＩ以上の整数であり、
    上記各センサの観測信号を周波数領域の信号に変換する周波数領域変換過程と、
    上記周波数領域の信号から、独立成分分析により周波数ごとに、信号源信号を分離する第１分離行列をそれぞれ算出する分離行列算出過程と、
    上記各第１分離行列の逆行列あるいは擬似逆行列（以下、両者を単に逆行逆列という）をそれぞれ算出する逆行列算出過程と、
    上記周波数ごとの逆行列の少なくとも１つにおける列ごとの異なる二つの要素比から上記信号源の一つの位置情報を計算する位置情報計算過程とを有する。
12. 請求の範囲第１１記載の方法において、
    上記位置情報計算過程は、上記周波数ごとの逆行列の複数について、各列ごとの上記要素比からの位置情報計算により上記各信号源の位置情報を求める過程であり、
    上記周波数ごとの信号源ごとの位置情報から、上記複数の逆行列と対応する分離行列における位置情報と対応する行が予め決めた順序になるように行の入れ替えを行うパーミュテーション行列を生成するパーミュテーション行列生成過程と、
    上記パーミュテーション行列と上記第１の分離行列とを乗算して行を入れ替えた第２分離行列を求める分離行列置換過程とを有する。
13. 請求の範囲第１２項記載の方法において、
    Ｊ≧３であり、かつＪ個のセンサは少なくとも２次元に配置され、
    上記位置情報は上記センサから上記信号源への方向を含み、
    その信号源が存在する円錐面であり、
    上記位置情報計算過程は、各列ごとの上記要素比から上記円錐面の計算を異なる二つの要素組の複数について行う過程であり、周波数ごとの上記複数の円錐面の共通直線の方向を、上記位置情報として推定する共通線推定過程を含む。
14. 請求の範囲第１２項記載の方法において、
    Ｊ≧３であり、上記位置情報は一対のセンサから上記信号源への方向を含み、その信号源が存在する円錐面及びその信号源が存在する曲面であり、
    上記位置情報計算過程は、上記２つの要素の比から上記円錐面を計算する過程と、上記円錐面の計算に用いた２つの要素と対応する一対のセンサの間隔より大きい間隔の一対のセンサと対応する２つの要素の比からこれら一対のセンサと一つの信号源との距離比を計算する過程と、上記距離比から上記曲面を計算する過程を有し、
    上記パーミュテーション行列生成過程は周波数ごとに上記円錐面及び上記球面に基づいて上記パーミュテーション行列を生成する過程である。
15. 請求の範囲第１４項記載の方法において、
    上記パーミュテーション行列生成過程は、円錐面及び上記球面の一方について上記位置情報計算過程を実行させ、得られた上記一方の位置情報に基づいて、各周波数での逆行列の各列と対応する上記一方の位置情報が予め決めた順になるように列の入れ替えを行う第１の入替行列を生成する過程と、
    この過程により、列の入れ替えを行うことができない列について上記位置情報の他方について上記位置情報計算過程を実行させ、得られた上記他方の位置情報に基づいて上記逆行列の列の入れ替えを行うように上記第１の入替行列を修正して第２の入替行列を生成する過程と、上記第２の入替行列の逆行列を計算して上記パーミュテーションとを有する。
16. 請求の範囲第１２項乃至第１５項のいずれかに記載の方法において、
    Ｊ≧３であり、上記Ｊ個のセンサは少なくとも２次元に配置され、
    上記位置情報は上記センサから上記信号源への方向を含み、その信号源が存在する円錐面であり、
    上記円錐面を表わす角度が予め決めた第１角度と第２角度との間であるか否かを判定し、これら角度の間にない円錐面を破棄する判定過程を有する。
17. 請求の範囲第１２項に記載の方法において、
    Ｊ≧３であり、上記Ｊ個のセンサは少なくとも２次元に配置され、
    上記位置情報は信号源が存在する球面の半径であり、
    上記位置情報計算過程は上記二つの要素比から距離比を計算する過程である。
18. 請求の範囲第１２項に記載の方法において、
    上記位置情報は上記センサから上記信号源への方向を含み、その信号源が存在する円錐面である。
19. 請求の範囲第１２項〜第１５項のいずれか又は第１７項又は第１８項に記載の方法において、
    上記パーミュテーション行列生成過程でパーミュテーション行列を生成できない周波数があれば、
    上記観測信号を上記第２分離行列で分離した信号の周波数領域の信号中の上記パーミュテーション行列を生成できた周波数成分と上記パーミュテーション行列を生成できない周波数成分との相関を計算する過程と、
    この計算した相関が大きくなるように、上記生成できなかった周波数の分離行列に対するパーミュテーション行列を生成する過程とを有する。
20. 請求の範囲第１１項記載の方法において、
    上記位置情報は上記センサから上記信号源への方向情報であり、
    上記位置情報計算過程は、上記比の偏角を、単位距離当たりの位相回転量と上記二つの要素とそれぞれ対応する上記センサ間の距離との積で除算し、その除算結果の逆コサインを計算して上記方向情報を出力する過程である。
21. 請求の範囲第２０項記載の方法において、
    上記位置情報計算過程は、一つの信号源の上記方向情報を、周波数ごとに計算する過程であり、上記各一つの信号源ごとに計算された周波数ごとの方向情報を統合して方向情報を確定する統合過程を有する。
22. 請求の範囲第１１項乃至第１９項のいずれかに記載した信号源位置情報推定方法の各過程をコンピュータに実行させるためのプログラム。

Images