(57)【要約】２次元配置された複数のセンサよりの観測信号ｘ１（ｔ）〜ｘＪ（ｔ）を短時間フーリエ変換し、これらから信号Ｘ１（ω１）〜Ｘ１（ωＮ），〜ＸＪ（ω１）〜ＸＪ（ωＮ）、独立成分分析法により分離行列Ｗ（ω１）〜Ｗ（ωＮ）を生成し、その逆行列Ｈ（ω１）〜Ｈ（ωＮ）を計算し、各ωｎ（ｎ＝１，…，Ｎ）ごとに、Ｈ（ωｎ）の各列の要素対Ｈｊｉ（ωｎ）及びＨｊ’ｉ（ωｎ）について、角度θ＾ｉ，ｊｊ’（ωｎ）＝ｃｏｓ−１（ａｒｇ（Ｈｊｉ（ωｎ）／Ｈｊ’ｉ（ωｎ））／（ωｎｃ−１‖ｄｊ−ｄｊ’‖））を計算し、ａｒｇ（α）はαの偏角、ｃは信号の伝搬速度、‖ｄｊ−ｄｊ’‖はセンサｊとｊ’の間隔である。Ｈ（ω１）〜Ｈ（ωＮ）において各列から得られるθ＾ｉ，ｊｊ’（ωｎ）が昇順になるように、列の入れ替えを行い、列の入れ替えができない列について、‖ｑｉ−ｄｊ’‖／‖ｑｉ−ｄｊ‖＝｜Ｈｊｉ（ωｎ）／Ｈｊ’ｉ（ωｎ）｜＝ＤＲｉ，ｊｊ’（ωｎ）をｑｉについて解いてＲｉ，ｊｊ’（ωｎ）＝‖ＤＲｉ，ｊｊ’（ωｎ）（ｄｊ−ｄｊ’）／（ＤＲ２ｉ，ｊｊ’（ωｎ）−１）‖を計算し、このＲｉ，ｊｊ’（ωｎ）が昇順になるように当該Ｈ（ωｎ）の列を入れ替え、このＨ（ωｎ）を用いてＷ（ω）のパーミュテーションを解決する。
2004079388.tif 000002
【特許請求の範囲】
【請求項１】
Ｉ箇所の信号源から放射された信号を、Ｊ個のセンサで検出して上記信号源の位置情報を求める信号源位置情報推定装置であって、Ｉは２以上の整数、ＪはＩ以上の整数であり、
上記各センサの観測信号を周波数領域の信号に変換する周波数領域変換手段と、
上記周波数領域の信号から、独立成分分析により周波数ごとに、信号源信号を分離する第１分離行列をそれぞれ算出する分離行列算出手段と、
上記各第１分離行列の逆行列をあるいは擬似逆行列「以下、両者を単に逆行列という」をそれぞれ算出する逆行列算出手段と、
上記周波数ごとの逆行列の少なくとも１つに対し、その各列ごとにその二つの要素比から上記信号源の一つの位置情報を計算する位置情報計算手段とを具備する。
【請求項２】
請求の範囲第１項記載の装置において、
上記位置情報計算手段は、上記周波数ごとの逆行列の複数について、上記各列ごとの要素比からの位置情報計算により上記各信号源の位置情報を求める手段であり、
上記周波数ごとの信号源ごとの位置情報から、上記複数の逆行列と対応する周波数の分離行列における位置情報と対応する行が予め決めた順序になるように行の入れ替えを行うパーミュテーション行列を生成するパーミュテーション行列生成手段と、
上記パーミュテーション行列と上記第１の分離行列とを乗算して行を入れ替えた第２の分離行列を求める並替手段とを備える。
【請求項３】
請求の範囲第２項記載の装置において、
上記Ｊは３以上であり、かつＪ個のセンサは少なくとも２次元に配置され、
上記位置情報は上記センサから上記信号源への方向を含み、その信号源が存在する円錐面であり、
上記位置情報計算手段は、各列ごとの上記要素比から上記円錐面の計算を異なる二つの要素組の複数について行う手段と、周波数ごとの上記複数の円錐面の共通直線の方向を、上記位置情報として推定する到来方向決定手段を含む。
【請求項４】
請求の範囲第２項記載の装置において、
上記Ｊは３以上であり、上記位置情報は一対のセンサから上記信号源への方向を含み、その信号源が存在する円錐面及びその信号源が存在する曲面であり、
上記位置情報計算手段は、上記二つの要素の比から上記円錐面を計算する手段と、上記円錐面の計算に用いた二つの要素と対応するセンサ間隔よりも大きいセンサ間隔のセンサと対応する二つの要素の比からこれらセンサと上記信号線との距離比を計算する手段と、上記距離比から上記曲面を計算する手段を備え、
上記パーミュテーション行列生成手段は周波数ごとに上記円錐面及び上記曲面に基づいて上記パーミュテーション行列を生成する手段である。
【請求項５】
請求の範囲第２項、第３項又は第４項に記載の装置において、
上記Ｊは３以上であり、上記Ｊ個のセンサは少なくとも２次元に配置され、
上記位置情報は上記センサから上記信号源への方向を含み、その信号源が存在する円錐面であり、
上記円錐面を表わす角度が予め決めた第１角度と第２角度との間であるか否かを判定し、これら第１角度及び第２角度の間の円錐面を有効とする判定手段を備える。
【請求項６】
請求の範囲第２項記載の装置において、
上記Ｊは３以上であり、上記Ｊ個のセンサは少なくとも２次元に配置され、
上記位置情報は信号源が存在する球面の半径であり、
上記位置情報計算手段は上記二つの要素比から距離比を計算する手段である。
【請求項７】
請求の範囲第２項に記載の装置において、
上記位置情報は上記センサから上記信号源への方向を含みその信号源が存在する円錐面である。
【請求項８】
請求の範囲第２項、第３項、第４項、第６項又は第７項に記載の装置において、
上記観測信号に対し、上記第２分離行列を用いて分離した信号の周波数領域の信号中の、上記パーミュテーション行列生成手段で、パーミュテーション行列が得られない周波数成分と、パーミュテーション行列が得られた周波数成分との相関を計算する相関計算手段と、
上記相関が大きくなるように上記パーミュテーション行列が得られない周波数の分離行列の行を入れ替えるパーミュテーション行列を生成する修正手段とを備える。
【請求項９】
請求の範囲第１項記載の装置において、
上記位置情報は上記センサから上記信号源への方向情報であり、
上記位置情報計算手段は、上記比の偏角を計算する手段と、単位距離当たりの位相回転量と上記二つの要素とそれぞれ対応する上記センサ間の距離との積を計算する手段と、上記偏角を上記積の結果で除算する手段と、その除算結果の逆コサインを計算して上記方向情報を出力する手段とを備える。
【請求項１０】
請求の範囲第９項記載の装置において、
上記位置情報計算手段は、一つの信号源の上記方向情報を、周波数ごとに計算する手段であり、上記各一つの信号源ごとに計算された周波数ごとの方向情報を統合して方向情報を確定する統合手段を備える。
【請求項１１】
Ｉ箇所の信号源から放射された信号を、Ｊ個のセンサで検出して上記信号源の位置情報を求める信号源位置情報推定方法であって、Ｉは２以上の整数、ＪはＩ以上の整数であり、
上記各センサの観測信号を周波数領域の信号に変換する周波数領域変換過程と、
上記周波数領域の信号から、独立成分分析により周波数ごとに、信号源信号を分離する第１分離行列をそれぞれ算出する分離行列算出過程と、
上記各第１分離行列の逆行列あるいは擬似逆行列（以下、両者を単に逆行逆列という）をそれぞれ算出する逆行列算出過程と、
上記周波数ごとの逆行列の少なくとも１つにおける列ごとの異なる二つの要素比から上記信号源の一つの位置情報を計算する位置情報計算過程とを有する。
【請求項１２】
請求の範囲第１１記載の方法において、
上記位置情報計算過程は、上記周波数ごとの逆行列の複数について、各列ごとの上記要素比からの位置情報計算により上記各信号源の位置情報を求める過程であり、
上記周波数ごとの信号源ごとの位置情報から、上記複数の逆行列と対応する分離行列における位置情報と対応する行が予め決めた順序になるように行の入れ替えを行うパーミュテーション行列を生成するパーミュテーション行列生成過程と、
上記パーミュテーション行列と上記第１の分離行列とを乗算して行を入れ替えた第２分離行列を求める分離行列置換過程とを有する。
【請求項１３】
請求の範囲第１２項記載の方法において、
Ｊ≧３であり、かつＪ個のセンサは少なくとも２次元に配置され、
上記位置情報は上記センサから上記信号源への方向を含み、
その信号源が存在する円錐面であり、
上記位置情報計算過程は、各列ごとの上記要素比から上記円錐面の計算を異なる二つの要素組の複数について行う過程であり、周波数ごとの上記複数の円錐面の共通直線の方向を、上記位置情報として推定する共通線推定過程を含む。
【請求項１４】
請求の範囲第１２項記載の方法において、
Ｊ≧３であり、上記位置情報は一対のセンサから上記信号源への方向を含み、その信号源が存在する円錐面及びその信号源が存在する曲面であり、
上記位置情報計算過程は、上記２つの要素の比から上記円錐面を計算する過程と、上記円錐面の計算に用いた２つの要素と対応する一対のセンサの間隔より大きい間隔の一対のセンサと対応する２つの要素の比からこれら一対のセンサと一つの信号源との距離比を計算する過程と、上記距離比から上記曲面を計算する過程を有し、
上記パーミュテーション行列生成過程は周波数ごとに上記円錐面及び上記球面に基づいて上記パーミュテーション行列を生成する過程である。
【請求項１５】
請求の範囲第１４項記載の方法において、
上記パーミュテーション行列生成過程は、円錐面及び上記球面の一方について上記位置情報計算過程を実行させ、得られた上記一方の位置情報に基づいて、各周波数での逆行列の各列と対応する上記一方の位置情報が予め決めた順になるように列の入れ替えを行う第１の入替行列を生成する過程と、
この過程により、列の入れ替えを行うことができない列について上記位置情報の他方について上記位置情報計算過程を実行させ、得られた上記他方の位置情報に基づいて上記逆行列の列の入れ替えを行うように上記第１の入替行列を修正して第２の入替行列を生成する過程と、上記第２の入替行列の逆行列を計算して上記パーミュテーションとを有する。
【請求項１６】
請求の範囲第１２項乃至第１５項のいずれかに記載の方法において、
Ｊ≧３であり、上記Ｊ個のセンサは少なくとも２次元に配置され、
上記位置情報は上記センサから上記信号源への方向を含み、その信号源が存在する円錐面であり、
上記円錐面を表わす角度が予め決めた第１角度と第２角度との間であるか否かを判定し、これら角度の間にない円錐面を破棄する判定過程を有する。
【請求項１７】
請求の範囲第１２項に記載の方法において、
Ｊ≧３であり、上記Ｊ個のセンサは少なくとも２次元に配置され、
上記位置情報は信号源が存在する球面の半径であり、
上記位置情報計算過程は上記二つの要素比から距離比を計算する過程である。
【請求項１８】
請求の範囲第１２項に記載の方法において、
上記位置情報は上記センサから上記信号源への方向を含み、その信号源が存在する円錐面である。
【請求項１９】
請求の範囲第１２項〜第１５項のいずれか又は第１７項又は第１８項に記載の方法において、
上記パーミュテーション行列生成過程でパーミュテーション行列を生成できない周波数があれば、
上記観測信号を上記第２分離行列で分離した信号の周波数領域の信号中の上記パーミュテーション行列を生成できた周波数成分と上記パーミュテーション行列を生成できない周波数成分との相関を計算する過程と、
この計算した相関が大きくなるように、上記生成できなかった周波数の分離行列に対するパーミュテーション行列を生成する過程とを有する。
【請求項２０】
請求の範囲第１１項記載の方法において、
上記位置情報は上記センサから上記信号源への方向情報であり、
上記位置情報計算過程は、上記比の偏角を、単位距離当たりの位相回転量と上記二つの要素とそれぞれ対応する上記センサ間の距離との積で除算し、その除算結果の逆コサインを計算して上記方向情報を出力する過程である。
【請求項２１】
請求の範囲第２０項記載の方法において、
上記位置情報計算過程は、一つの信号源の上記方向情報を、周波数ごとに計算する過程であり、上記各一つの信号源ごとに計算された周波数ごとの方向情報を統合して方向情報を確定する統合過程を有する。
【請求項２２】
請求の範囲第１１項乃至第１９項のいずれかに記載した信号源位置情報推定方法の各過程をコンピュータに実行させるためのプログラム。
【発明の詳細な説明】
【技術分野】  この発明は音源や電波源などの信号源の複数から放射され、互いに混合された信号を複数のセンサで観測して、各信号源の位置情報、正しくは位置を表わすパラメータの少くとも１つを含む情報を推定し、信号の到来方向の検出や、信号源ごとの信号に分離復元に適用される装置、その方法、及びプログラムに関する。
【背景技術】  複数の信号源からの各信号が空間内で混合されて複数のセンサに到来し、これらセンサで観測された到来信号から、各源信号の到来方向の推定や各源信号を分離することを、その源信号の混合系の情報を知らずに、独立成分分析（Ｉｎｄｅｐｅｎｄｅｎｔ  Ｃｏｍｐｏｎｅｎｔ  Ａｎａｌｙｓｉｓ、以後、ＩＣＡと記述する場合もある）により行うことが提案されている。空間内での前記混合は信号源からセンサまでの到達遅延及び減衰度が直接波と障害物などによる複数の反射波で異なるため、ある信号は複数の時間遅れを持って混合された畳み込み混合となる。時間領域で直接分離フィルタを求めるＩＣＡは最終的な解への収束が非常に遅いため、周波数領域で周波数毎にＩＣＡを適用する方法が現実的である。
［到来方向推定］
周波数領域でのＩＣＡを用いて信号源の方向を位置情報として推定する従来技術を図１を参照して簡単に説明する。Ｊ個のセンサ１１，１２，…，１Ｊが直線状に配列されている。センサ１ｊ（ｊ＝１，２，…，Ｊ）の位置をｄｊ、そのセンサ１ｊでの観測信号をｘｊ（ｔ）とする。センサ１１，…，１ｊの配列方向と垂直な方向を９０°とし、源信号ｓｉ（ｔ）の到来方向を０°≦θｊ≦１８０°とする。Ｉ個の源信号ｓ１（ｔ），…，ｓＩ（ｔ）の混合信号をＪ個のセンサ１１〜１Ｊにより、観測信号ｘ１（ｔ），…，ｘＪ（ｔ）として検出しているものとする。
信号の到来方向の推定は周波数領域で行われることが多い。これには、観測信号ｘｊ（ｔ）を短時間フーリエ変換して周波数領域における時間系列信号Ｘｊ（ω，ｍ）を求める。ωは角周波数（周波数をｆとするとω＝２πｆ）ｍは時刻を表す番号である。源信号ｓｉ（ｔ）（ｉ＝１，…，Ｉ）も、同様に、周波数領域における時間系列Ｓｉ（ω，ｍ）に変換したものとすると、観測信号Ｘｊ（ω，ｍ）は、信号ｓｉの信号源からセンサ１ｊへの周波数応答Ａｊｉ（ω）を用いＸｊ（ω，ｍ）＝Σｉ＝１ＩＡｊｉ（ω）Ｓｉ（ω，ｍ）
と表現することができる。ベクトルと行列を用いると、
2004079388.tif 000003
となる。ここで、
2004079388.tif 000004
はそれぞれ、Ｊ個のセンサによる観測信号およびＩ個の源信号をベクトル表現したものである。Ａ（ω）は周波数応答Ａｊｉ（ω）を要素とするＪ×Ｉ行列であり、信号混合系の周波数応答を表現するので混合行列と呼ばれる。［ａ］Ｔはベクトル又は行列ａの転置を表わす。
図１において、方向θｉから到来する源信号は、ｄ１＝０の位置のセンサ１１に対し、センサ１ｊにもτｊｉ＝ｃ−１ｄｊｃｏｓθｉだけ速く到達する。ｃは源信号ｓｉの速度である。従って、直接波のみを考慮すると角周波数ωにおける周波数応答は
2004079388.tif 000005
とモデル化することができる。方向θを変数とする到来方向ベクトルを
ａ（ω，θ）＝［ｅｘｐ（ｊωｃ−１ｄ１ｃｏｓθ）、ｅｘｐ（ｊωｃ−１ｄ２ｃｏｓθ），…，ｅｘｐ（ｊωｃ−１ｄＪｃｏｓθ）］Ｔと表現すると、観測信号は
2004079388.tif 000006
とも近似表現することができる。
独立成分分析を用いて信号源の方向を推定する方法については例えばＳ．Ｋｕｒｉｔａ，Ｈ．Ｓａｒｕｗａｔａｒｉ，Ｓ  Ｋａｊｉｔａ，Ｋ．Ｔａｋｅｄａ，ａｎｄ  Ｆ．Ｉｔａｋｕｒａ，“Ｅｖａｌｕａｔｉｏｎ  ｏｆ  ｂｌｉｎｄ  ｓｉｇｎａｌ  ｓｅｐａｒａｔｉｏｎ  ｍｅｔｈｏｄ  ｕｓｉｎｇ  ｄｉｒｅｃｔｉｖｉｔｙ  ｐａｔｔｅｒｎ  ｕｎｄｅｒ  ｒｅｖｅｒｂｅｒａｎｔ  ｃｏｎｄｉｔｉｏｎｓ，”ｉｎ  Ｐｒｏｃ．ＩＣＡＳＳＰ２０００，２０００，ｐｐ．３１４０−３１４３（文献１という）に示されている。この方法を以下に簡単に述べる。
観測信号Ｘ（ω，ｍ）はＸ（ω，ｍ）＝Ａ（ω）Ｓ（ω，ｍ）であり、源信号Ｓ（ω，ｍ）が混合されたものであるので、互いに独立ではない。Ｘ（ω，ｍ）に独立成分分析
2004079388.tif 000007
を適用すると、互に独立な分離信号
2004079388.tif 000008
が得られる。
Ｗ（ω）は要素がＷｉｊ（ω）であるＩ×Ｊの行列であり、分離行列と呼ばれる。例えば、Ｉ＝Ｊ＝２の場合、独立成分分析は、Ｙ１（ω，ｍ）とＹ２（ω，ｍ）が互いに独立になる様に
2004079388.tif 000009
を満足する分離行列Ｗ（ω）を求める。源信号Ｓ１（ω，ｍ），…，ＳＩ（ω，ｍ）が互いに独立であれば、分離信号Ｙ１（ω，ｍ），…，ＹＩ（ω，ｍ）はそれぞれ、源信号の何れかに対応する。但し、独立成分分析は信号の独立性のみに基づいているところから、分離信号の順序と大きさに関して任意性がある。即ち、分離行列Ｗ（ω）の行を入れ換えても、Ｗ（ω）の行を定数倍しても、独立成分分析の解である。なお後で述べるが順序の任意性はパーミュテーションの問題となり、大きさの任意性はスケーリングの問題となる。
分離行列Ｗ（ω）のｉ行目は、
2004079388.tif 000010
である。このｗｉ（ω）は分離信号Ｙｉ（ω，ｍ）を作り出していることが分かる。従って、ｗｉ（ω）は、源信号Ｓ１（ω，ｍ），…，Ｓ１（ω，ｍ）の何れか１つを強調して取り出し、それ以外を抑圧している。Ｗｉ（ω）が形成する指向特性を解析することにより、何れの方向から到来する信号を取り出し、何れの方向から到来する信号を抑圧しているかを解析することができ、この解析により、源信号ｓｉ（ｔ）の到来方向を推定することができる。これを全てのｗｉ（ω），ｉ＝１，…，Ｉに対して行うと、分離行列のＷ（ω）各行ｗｉ（ω）が取り出している源信号の到来方向Θ（ω）＝［θ１（ω），…，θ１（ω）］Ｔを推定することができる。
Ｗｉ（ω）が形成する指向特性は、到来方向ベクトルａ（ω，θ）を用いてＢｉ（ω，θ）＝ｗｉ（ω）ａ（ω，θ）と表現することができる。このＢｉ（ω，θ）は方向θにある源信号から、分離信号Ｙｉ（ω，ｍ）への周波数応答と考えられる。周波数３１５６Ｈｚにおいて独立成分分析後の指向特性のゲイン｜Ｂｉ（ω，θ）｜特性を図２に示す。図２の横軸はθを、縦軸はゲインを表わす。実線は分離行列の１行目が与える指向特性｜Ｂ１（ω，θ）｜であり、破線は２行目が与える指向特性｜Ｂ２（ω，θ）｜である。実線は５５°でゲインが最小となっており、破線は１２１°でゲインが最小となっている。このことから、分離行列の１行目は５５°から到来する信号を抑圧して１２１°から到来する信号を取り出し、逆に分離行列の２行目は１２１°から到来する信号を抑圧して５５°から到来する信号を取り出している。従って、Θ（３１５６Ｈｚ）＝［１２１°，５５°］Ｔと推定することができる。
複数の信号源方向を複数センサを用いて、センサの観測信号を周波数領域に変換して推定する方法として、ＭＵＳＩＣ（Ｍｕｌｔｉｐｌｅ  Ｓｉｇｎａｌ  Ｃｌａｓｓｉｆｉｃａｔｉｏｎ）法がある（例えばＳ．Ｕｎｎｉｋｒｉｓｈｎａ  Ｐｉｌｌａｉ，“Ａｒｒａｙ  Ｓｉｇｎａｌ  Ｐｒｏｃｅｓｓｉｎｇ”，Ｓｐｒｉｎｇｅｒ−Ｖｅｒｌａｇ，１９８９，ＩＳＢＮ  ０−３８７−９６９５１−９，ＩＳＢＮ  ３−５４０−９６９５１−９参照）。この方法はセンサの数Ｊより１個少ない（Ｊ−１）個までしか信号源の方向を推定できない。しかし前記独立成分分析（ＩＣＡ）を用いる方法（以下単にＩＣＡ法ともいう）によれば２つのセンサで２信号の混合に対応することができ、この点で、ＭＵＳＩＣ法より優れている。しかし、上記のＩＣＡを用いる方法では３信号以上の混合に対処するには以下に示すように困難を伴う。また、指向特性のゲインの最小値を求める操作に多くの計算を要することも欠点である。
３信号の混合に対して３つのセンサを用いてＩＣＡ法を適用した状況を説明する。この場合のＩＣＡは分離行列が３×３になるだけで同様に行えるが、指向特性のゲインを解析するのに困難を伴う。この場合のＩＣＡ後の周波数２７３４Ｈｚにおいて指向特性のゲイン｜Ｂｉ（ω，θ）｜を図３に示す。図３の実線は分離行列の１行目が、破線は２行目が、１点鎖線は３行目がそれぞれ与える指向特性である。この場合、各源信号は分離行列のある行によって強調され、他の２つの行によって抑圧される筈である。しかし、抑圧する２つの行が同じ方向で抑圧しているとは限らない。例えば、図３においては４５°付近で｜Ｂ２｜と｜Ｂ３｜が共に極小となり、ｗ１（ω）が４５°付近の源信号を取り出し、ｗ２（ω）及びｗ３（ω）はこの源信号を抑圧している。同様に９０°付近で｜Ｂ１｜及び｜Ｂ３｜が極小となり、ｗ２（ω）で９０°付近の源信号を取り出し、ｗ１（ω）及びｗ３（ω）でこの信号を抑圧していると読み取れる。しかしｗ３（ω）が１２０°付近の源信号を取り出し、ｗ１（ω）およびｗ２（ω）がこの源信号を抑圧していると考えられるが、１２０°付近での｜Ｂ１｜及び｜Ｂ２｜の極小がかなりずれている。この様なずれが大きくなると、何れの方向の抑圧が何れの源信号に対応するか判然としなくなる。従って、ＩＣＡ法による従来技術を３信号以上の状況に対して適用するのは困難であると考えられる。
〔ブラインド信号分離〕
次にＩＣＡを用いたブラインド信号分離の従来技術を説明する。ブラインド信号分離とは、観測された混合信号から源信号を推定する技術である。以下では、Ｉ個の源信号が混在した混合信号がＪ個のセンサで観測される場合を例にとって説明する。
信号源ｉから発生した源信号をｓｉ（ｔ）（ｉ＝１，…，Ｉ、ｔは時刻）、センサｊで観測された混合信号をｘｊ（ｔ）（ｊ＝１，…，Ｊ）とすると、混合信号ｘｊ（ｔ）は次式で表わせる。
2004079388.tif 000011
ここで、ａｊｉは信号源ｉからセンサｊへのインパルス応答、＊は畳み込み演算子である。ブラインド信号分離の目的は、観測信号ｘｊ（ｔ）のみを用いて、分離のためのフィルタｗｋｊ及び分離信号ｙｋ（ｔ）（ｋ＝１，…，Ｉ）を次式により求めることである。
2004079388.tif 000012
時間領域での畳み込み混合は、周波数領域における複数の瞬時混合に変換することができる。すなわち、前記の式（７）及び式（８）は、それぞれ前記式（１）及び式（５）で表わされ、前述したように、Ｗ（ω）は分離行列であり、ＩＣＡを用いてＹｋ（ω，ｍ）とＹｋ（ω，ｍ）が互いに独立となるように計算されたもの、つまりＩＣＡの解である。
以上の周波数領域でのブラインド信号分離を行うに際し問題となるのは、パーミュテーションの問題とスケーリングの問題である。
前述したように分離行列Ｗ（ω）は行の入れ替えを行っても独立成分分析の解となる。つまりある角周波数ωにおいて、Ｗ（ω）がＩＣＡの解だとし、複素数を要素とする任意の対角行列をＤ（ω）とし、任意のパーミュテーション行列（この行列を任意の行列の左から掛けた結果は、この任意の行列の行を入れ替えた行列となる）をＰ（ω）とした場合におけるＰ（ω）Ｄ（ω）Ｗ（ω）もまたＩＣＡの解となる。これは、ＩＣＡが源信号間の統計的独立性のみを条件として源信号の分離を行うことに基づくものであり、このＤ（ω）によって与えられる解の自由度をスケーリングの自由度と呼び、Ｐ（ω）によって与えられる解の自由度をパーミュテーションの自由度と呼ぶ。
そのため、適切なブラインド信号分離を行うためには、全てのωについて、ＩＣＡの解から分離行列として適切な解Ｗ（ω）を特定しなければならない。一般に、この適切な解Ｗ（ω）の特定は、任意に求めたＩＣＡの解に、適切なＤ（ω）やＰ（ω）を乗じ、その結果を適切な解Ｗ（ω）とすることによって行われる。そして、全てのωについて、このＤ（ω）を適切に決める問題をスケーリングの問題と呼び、Ｐ（ω）を適切に決める問題をパーミュテーションの問題と呼ぶ。また、パーミュテーションとは、｛１，２，…，Ｉ｝から｛１，２，…，Ｉ｝への全単射な関数Ｚ：｛１，２，…，Ｉ｝→｛１，２，…，Ｉ｝であり、パーミュテーション行列と１対１対応する。
スケーリングの自由度は、時間領域において周波数特性を変化させるフィルタの自由度に相当する。それゆえ、時間領域において歪のない分離信号を生成するためには、全てのωについて、Ｄ（ω）を適切に決める必要がある。このスケーリングの問題は、例えば、Ｄ（ω）＝ｄｉａｇ（Ｗ−１（ω））とすることで容易に解決することができる。ｄｉａｇ（α）は行列αの対角化（対角要素の他は全ての要素を０とする）を表わす。つまり任意に求めたＩＣＡの解Ｗ０（ω）の逆行列を求め、それを対角化した行列をＤ（ω）とし、Ｄ（ω）Ｗ０（ω）を、適切な分離行列Ｗ（ω）として特定する。このことは既に知られている。例えば参考文献：Ｋ．Ｍａｔｓｕｏｋａ  ａｎｄ  Ｓ．Ｎａｋａｓｈｉｍａ，“Ｍｉｎｉｍａｌ  Ｄｉｓｔｏｒｔｉｏｎ  Ｐｒｉｎｃｉｐｌｅ  ｆｏｒ  Ｂｌｉｎｄ  Ｓｏｕｒｃｅ  Ｓｅｐａｒａｔｉｏｎ”，Ｐｒｏｃ．ＩＣＡ  ２００１，ｐｐ．７２２−７２７．に示されている。
一方、パーミュテーションの自由度のため、前記の式（５）の演算結果として、例えば、ある角周波数ω１については分離信号Ｙ１（ω１，ｍ）が源信号Ｓ１（ω１，ｍ）の推定値として出力され、別の角周波数ω２については分離信号Ｙ１（ω２，ｍ）が源信号Ｓ２（ω２，ｍ）の推定値として出力されることもありうる。このような場合、時間領域の出力信号ｙ１（ｔ）の中に、時間領域の源信号ｓ１（ｔ）の成分とｓ２（ｔ）の成分とが混在してしまい、分離信号を正しく生成することができない。そのため、時間領域における出力信号ｙ１（ｔ）が正しく源信号ｓ１（ｔ）の推定値となるためには、全てのωについて、Ｙ１（ω，ｍ）がＳ１（ω，ｍ）の推定値となるようにＰ（ω）を適切に決める必要がある。
従来のパーミュテーション問題の代表的な解法として、前記文献１に示す信号の到来方向を推定する方法がある。つまり、図２を参照して説明したように、各周波数における分離行列Ｗ（ω）の各行と対応する指向特性を求める（図２ではｆ＝３１５６Ｈｚの指向特性のみを示している）。これら各周波数の指向特性において例えば分離行列Ｗ（ω）の１行目が与える指向特性のゲイン最小値が５５°で、２行目が与える指向特性のゲイン最小値が１２１°になるようにする。つまりある角周波数ωｎではその分離行列Ｗ（ωｎ）の１行目が与える指向特性のゲイン最小値が１２１°で、２行目が与える指向特性のゲイン最小値が５５°であれば、その分離行列Ｗ（ωｎ）の１行目と２行目とを入れ替える。つまりこの入れ替えを行うパーミュテーション行列Ｐ（ωｎ）をＷ（ωｎ）にその左から掛算する。
このパーミュテーション問題の解法は先にも述べたように指向特性のゲイン最小値を求める操作に多くの計算を要し、しかも信号源の数Ｉが３以上では全角周波数のＷ（ω）について適当に並べ替えを行ってみるという試行錯誤が必要となる。更に図３を参照して説明したようにＷ（ω）のある行がある方向の信号Ｓｉ（ω，ｍ）を取り出している場合は、そのＷ（ω）の他の行は全てその方向の信号Ｓｉ（ω，ｍ）を抑圧する状態に必ずしもならない。
また指向特性のゲインが低い方向を探索して推定する信号到来方向の推定精度は信号源の位置によって異なり、特に一対のセンサ１ｊと１ｊ’を結んだ直線（以下、センサペア軸という）に信号到来方向が近い場合に誤差が大きくなる。このことを実験により示す。図４Ａに示すように、２つのセンサ１０１，１０２としてマイクロホンを距離２．８３ｃｍ離して配置し、センサ１０１と１０２の中点（原点）から一定の距離（約１５０ｃｍ）離れ、かつ角度間隔が２０°離れた２つの信号源１１１，１１２として音源を設け、センサ１０１からセンサ１０２を見た方向を基準（０°）として、原点からみた信号源１１１の角度θが１０°から１５０°となるまで音源１１１と１１２を前記一定距離及び角度間隔を保持した状態で移動させた。
このように音源１１１，１１２を移動させながら行ったブラインド信号分離結果を図４Ｂに例示する。図４Ｂの縦軸は、信号対妨害音比（Ｓｉｇｎａｌ  ｔｏ  Ｉｎｔｅｒｆｅｒｅｎｃｅ  Ｒａｔｉｏ）を示しており、分離音中に含まれる目的信号及び妨害信号を用いて、
ＳＩＲ＝１０ｌｏｇ１０（目的信号のパワー／妨害信号のパワー）（ｄＢ）
のように計算したものである。また、図４Ｂの横軸は、原点からみた音源１１１の角度θを表し、実線はこの実験結果を示しており、点線はパーミュテーションに正解を与えたときのＳＩＲを示している。
この図より信号源１１１が、センサペア軸に近い方向（０°或いは１８０°付近）に近づくと、実験結果のＳＩＲはパーミュテーションに正解を与えたときのＳＩＲに比べ大きく低下していることがわかる。このことは信号源１１１の方向がセンサペア軸方向と近いとパーミュテーションが誤っていることに基づくと考えられる。
以上述べたように独立成分分析を用いて分離行列を求め、その各行から指向特性パターンを求めて、その利得が低い方向探索により、信号源の方向（到来信号方向）を求め、更にこれを利用してブラインド信号分離を行う場合は、前記指向特性パターンを求めて、利得が低い方向を探索するために多くの計算時間を必要とした。
この発明の目的は独立成分分析により分離行列を求めて信号源の位置情報を推定するための計算時間が短かい位置情報推定装置、その方法及びプログラムを提供することにある。
【発明の開示】
この発明によれば周波数領域の分離行列Ｗ（ω１），…，Ｗ（ωＮ）の逆行列（Ｉ＜Ｊの場合は擬似逆行列）を計算してスケーリングとパーミュテーションの自由度を含む混合行列Ａ（ω１），…，Ａ（ωＮ）の推定値Ｈ（ω１），…，Ｈ（ωＮ）を生成し、これらの周波数ごとのＨ（ωｎ）（ｎ＝１，…，Ｎ）について各列ごとにその２つの要素Ｈｊｉ（ωｎ）とＨｊ’ｉ（ωｎ）（ｊ、ｊ’はセンサを表すパラメータ、ｉは信号源を表すパラメータ）の比から、信号源ｉの位置情報の一つのパラメータ、例えば信号源が存在する円錐面あるいは曲面を計算する。
このように要素比に基づき表された計算式を演算すればよく、分離された信号の指向性パターンを求め、更にその極小値角度の探索を行う場合より少ない計算量で済む。また要素比をとるため、前述したスケーリングの自由度による影響がなくなる。
【図面の簡単な説明】  図１はセンサアレーと到来信号の各センサへの到来時間差の関係を説明するための図である。
図２は２音源混合信号に対し、ＩＣＡにより計算した分離行列の各行のゲイン指向特性を示す図である。
図３は３音源混合信号に対し、ＩＣＡにより計算した分離行列の各行のゲイン指向特性を示す図である。
図４Ａは予備実験に用いたセンサと信号源の関係を説明するための図である。
図４Ｂは前記予備実験の結果を示す図である。
図４Ｃは推定角度とその誤差及び偏角誤差との関係を示す図である。
図５はこの発明を信号到来方向の推定に適用した第１実施形態の機能構成例を示すブロック図である。
図６は第１実施形態の処理手順の例を示す流れ図である。
図７は図５中の角度計算部の具体的機能構成例を示すブロック図である。
図８は図６中のステップＳ４の具体的処理手順の例を示す流れ図である。
図９はこの第１実施形態により方向推定した実験結果を示す図である。
図１０Ａは音源が２個、マイクロホンが３個の場合に第１実施形態により方向推定した実験結果を示す図である。
図１０Ｂは図１０Ａについての実験と同一条件でＭＵＳＩＣ法により方向推定した実験結果を示す図である。
図１１Ａは音源が２個、マイクが３個の場合に第１実施形態により方向推定した実験結果を示す図である。
図１１Ｂは図１１Ａについての実験と同一条件でＭＵＳＩＣ法により方向推定した実験結果を示す図である。
図１２はこの発明をブラインド信号分離に適用した第２実施形態の機能構成例を示すブロック図である。
図１３は第２実施形態の処理手順の例を示す流れ図である。
図１４は図１３中のステップＳ１４の具体的処理手順の例を示すブロック図である。
図１５は複数円錐面の共通直線方向を説明するための図である。
図１６は図１２中の到来方向決定部１６の具体的機能構成の他の例を示すブロック図である。
図１７はこの発明をブラインド信号分離に適用した第３実施形態の機能構成例を示すブロック図である。
図１８は第３実施形態の処理手順の例を示す流れ図である。
図１９は図１８中のステップＳ３５の具体的処理手順の例を示す流れ図である。
図２０はセンサ配置と信号源位置と推定曲面との関係を説明するための図である。
図２１は推定球面の例を示す図である。
図２２は実験に用いた室とマイクロホン配置と音源との関係を示す図である。
図２３は小間隔マイクロホン対を用いた推定方向のヒストグラムである。
図２４はその周波数に対する推定方向の分布を示す図である。
図２５は音源２４及び２５に対する推定球面の半径の周波数に対する分布を示す図である。
図２６は各方法での実験結果を示す図である。
図２７はこの発明の第４実施形態の要部の機能構成例を示すブロック図である。
図２８はこの発明をブラインド信号分離に適用した第５実施形態の処理手順の例を示す流れ図である。
図２９は第５実施形態に対する実験結果を示す図である。
図３０Ａは複数の円錐面の共通直線方向を求めて分離行列のパーミュテーションを解決する処理手順の要部を示す流れ図である。
図３０Ｂは円錐面の推定と球面の推定を利用して分離行列のパーミュテーションを解決する処理手順の要部を示す流れ図である。
図３０Ｃは円錐面の推定と球面の推定を利用して分離行列のパーミュテーションを解決する他の処理手順の要部を示す流れ図である。
【発明を実施するための最良の形態】  まず信号源の位置情報として方向情報を推定する場合にこの発明を適用した実施形態を説明する。なお以下の説明において同一又は対応するものについて図面中に同一の参照番号を付けて、重複説明を省略する。
［第１実施形態］
この第１実施形態は信号源の方向、つまりその信号源から放射された源信号の到来方向を求める。
図５に第１実施形態の機能構成図を、図６にその処理手順の一部の流れ図をそれぞれ示す。
信号源の数Ｉ以上の個数であるＪ個のセンサ１１，１２，…，１Ｊが、例えば図１に示したように配列されている。隣接センサの間隔は通常は、源信号の最も短かい波長の１／２以下である。これらセンサ１ｊ（ｊ＝１，２，…，Ｊ）で観測された信号ｘｊ（ｔ）はそれぞれ周波数領域変換部１１ｊにおいて、例えば短時間フーリエ変換により周波数領域信号Ｘｊ（ω，ｍ）に変換される（図６、ステップＳ１）。これら周波数領域信号Ｘｊ（ω，ｍ）に対する各角周波数ωｎごとの分離行列Ｗ（ωｎ）（ｎ＝１，２，…，Ｎ）が分離行列算出部１２において独立成分分析により算出される（図６、ステップＳ２）。
2004079388.tif 000013
この各周波数ごとの分離行列Ｗ（ωｎ）の逆行列が逆行列算出部１３で算出され逆行列Ｈ（ωｎ）が求まる（図６、ステップＳ３）。
2004079388.tif 000014
なおこの逆行列の計算はＪ＞Ｉの場合は疑似逆行列を計算することになる。疑似逆行列としては例えばＭｏｏｒｅ−Ｐｅｎｒｏｓｅ型一般逆行列が用いられる。
この実施形態においてはＨ（ω）中の少くとも１つの周波数の逆行列Ｈ（ωｎ）の各列ｉの２つの要素Ｈｊｉ（ωｎ）とＨｊ’ｉ（ωｎ）との比の偏角から方向情報、つまり源信号の到来方向が角度計算部１４で計算される（図６、ステップＳ４）。この角度計算部１４の具体的機能構成例と処理手順例を図７及び図８を参照して説明する。選択部１４ａにより角周波数ωｎの逆行列Ｈ（ωｎ）中の未選択の一つの列ｉが選択され（図８、ステップＳ４ａ）、そのｉ列目から２つの要素Ｈｊｉ（ωｎ）とＨｊ’ｉ（ωｎ）が選択される（図８、ステップＳ４ｂ）。
偏角計算部１４ｂで、選択された要素の比の偏角ａｒｇ［Ｈｊｉ（ωｎ）／Ｈｊ’ｉ（ωｎ）］が計算され、また間隔計算部１４ｃで、センサ情報格納部１５内のセンサ１ｊと１ｊ’の位置情報ｄｊとｄｊ’が取り出され、これらの差からセンサ１ｊと１ｊ’の間隔ｄｊ−ｄｊ’が計算される（図８、ステップＳ４ｃ）。位相回転計算部１４ｄで、間隔ｄｊ−ｄｊ’と、角周波数ωｎでの単位距離当りの信号の位相回転量ωｎ／ｃ（ｃは源信号の速度）との積が計算されて間隔ｄｊ−ｄｊ’での信号の位相回転量が求められ、この位相回転量により偏角ａｒｇ［Ｈｊｉ（ωｎ）／Ｈｊ’ｉ（ωｎ）］が除算部１４ｅで割算される（図８、ステップＳ４ｄ）。なお単位距離当たりの位相回転量ω／ｃは例えばｆ＝６８０Ｈｚの音波（速度ｃ＝３４０ｍ／ｓ）では、１ｍ当たり２π・２の位相が回転することになる。
この除算部１４ｅの割算結果Ｇｉ（ωｎ）の絶対値が１以下か否かかが判定部１４ｆで判定され（図８、ステップＳ４ｅ）、１以下であればＧｉ（ωｎ）の逆余弦θｉ（ωｎ）＝ｃｏｓ−１Ｇｉ（ωｎ）が逆余弦計算部１４ｇで計算される（図８、ステップＳ４ｆ）。つまり角度計算部では次の計算が行われる。
2004079388.tif 000015
ステップＳ４ｅで｜Ｇｉ（ωｎ）｜が１以下でなければ角度θｉ（ωｎ）が虚数になるので、別の組み合せを選択するため、判定部１４ｆで選択したｉ列目のすべての要素の組み合せを選択したか判定され（図８、ステップＳ４ｇ）、選択していない組み合せがあればステップＳ４ｂに戻り、全ての組み合せを選択していれば、判定部１４ｆですべての列を選択したか判定され（図８、ステップＳ４ｈ）、選択していない列があればステップＳ４ａに戻り、全ての列を選択していれば角度計算処理を終了する。なお図７中の間隔計算部１４ｃは各周波数の角度計算部１４に対して共通に用いられる。
角度計算部１４からＨ（ω）中の選択した角周波数ωｎの逆行列Ｈ（ωｎ）における各列と対応して、その選択順に、もし第１列目から順次選択すればその順に式（７）の計算結果、つまりＩ個の信号源の方向（信号到来方向）Θ（ωｎ）＝（θ１（ωｎ），θ２（ωｎ），…，θＩ（ωｎ））が出力される。これらθ１（ωｎ），θ２（ωｎ），…，θＩ（ωｎ）は各源信号ｓ１（ｔ），ｓ２（ｔ），…，ｓＩ（ｔ）の到来方向（信号源１，２，…，Ｉの方向）のいずれかの１つと対応する。
この実施形態において信号到来方向が推定できるしくみを以下に説明する。独立成分分析（ＩＣＡ）により分離が達成されていれば、ＩＣＡにより算出した分離行列Ｗ（ω）と真の混合行列Ａ（ω）とはＰ（ω）Ｄ（ω）Ｗ（ω）Ａ（ω）＝Ｉの開係にある。ここで、Ｄ（ω）はスケーリングの自由度を示す対角行列、Ｐ（ω）はパーミュテーションの自由度を示すパーミュテーション行列、Ｉは単位行列である。ＩＣＡを用いても、混合行列Ａ（ω）そのものは一般には算出できない。しかし、Ｗ（ω）の逆行列Ｈ（ω）＝Ｗ−１（ω）＝Ａ（ω）Ｐ（ω）Ｄ（ω）を算出すれば、スケーリングとパーミュテーションの自由度を含む混合行列の推定値が得られる。すなわち、逆行列Ｈ（ω）は、混合行列Ａ（ω）の列をＰ（ω）により並べ替え、さらに各列にそれぞれＤ（ω）の対角要素を掛け合わせたものになる。
この実施形態では、逆行列Ｈ（ω）の同じ列ｉから２つの要素Ｈｊｉ（ω）とＨｊ’ｉ（ω）を取り出してそれらの比Ｈｊｉ（ω）／Ｈｊ’ｉ（ω）を求めることで、算出できないＤ（ω）によるスケーリングの自由度を取り除く。すなわち、
2004079388.tif 000016
となる。ここで、Ｚはパーミュテーション行列Ｐ（ω）を右から掛けることに対応する順列である。逆行列Ｈ（ω）の全ての列ｉに対して式（２１）に関わる操作を行うことで、Ｐ（ω）によるパーミュテーションＺにかかわらず、全ての信号の到来方向が推定できる。
背景技術の説明では、混合行列Ａ（ω）の要素をＡｊｉ（ω）＝ｅｘｐ（ｊωｃ−１ｄｊｃｏｓθｉ）とモデル化した。しかし、この実施形態では、分離行列Ｗ（ω）の逆行列Ｈ（ω）を用いて、スケーリングとパーミュテーションの自由度を含む混合行列Ａ（ω）の推定値を算出しているので、この単純なモデルでは不十分である。そこで、振幅の減衰度αｊｉ（実数）と原点における位相差ｅｘｐ（ｊφｉ）を用いて、Ａｊｉ（ω）＝αｊｉｅｘｐ（ｊφｉ）ｅｘｐ（ｊωｃ−１ｄｊｃｏｓθｉ）とモデル化する。このモデルを用いてＡｊＺ（ｉ）（ω）／Ａｊ’Ｚ（ｉ）（ω）を計算すると、式（２１）の関係から、
2004079388.tif 000017
が得られる。その結果、上記のＧｊ（ω）は
2004079388.tif 000018
となり、｜Ｇｉ（ω）｜＜１であればθＺ（ｉ）＝ｃｏｓ−１Ｇｉ（ω）が実数となり、到来方向を推定することができる。先に述べたように、パーミュテーションＺの自由度により、Ｉ個の方向Θ（ω）＝［θＺ（ｉ）（ω），…，θＺ（Ｉ）（ω）］全体を適当に並べ替えたものが、信号ｓ１，…，ｓＩの方向に対応する。
Ｈ（ω）中の複数の周波数又は全ての周波数の各逆行列Ｈ（ωｎ）（ｎ＝１，…，Ｎ）について先に述べたように各列ごとに角度θｉ（ωｎ）を求め、これらの全体により各到来方向を決定してもよい。つまり角度計算部１４で推定された各周波数の到来方向は図５中のソート部３２においてその角度ごとに分けられる（図６、ステップＳ５）。例えばそれぞれ方向（角度）の大きい順に並べられる。角周波数ω１についてのΘ（ω１）中の各成分が大きい順に（θ１（ω１），θ２（ω１），…，θＩ（ω１））と並べ、ω２についてのΘ（ω２）中の各成分が大きい順に（θ１（ω２），θ２（ω２），…，θＩ（ω２））と並べ、…，ωＮについてのΘ（ωＮ）中の各成分が大きい順に（θ１（ωＮ），θ２（ωＮ），…，θＩ（ωＮ））と並べられる。このように大きい角度順に分別された角度は、大きい順が同一であっても周波数間でばらつきがある、つまり角度θｉ（ω１），θｉ（ω２），…，θｉ（ωＮ）（ｉ＝１，…，Ｉ）はばらついている。
よって統合部３３で、分別角度θｉ（ω１），θｉ（ω２），…，θｉ（ωＮ）ごとに１つの角度θｉに統合されて１つの角度（到来方向）θｉとされる（図６、ステップＳ６）。この統合の方法としては例えば分別角度（θｉ（ω１），θｉ（ω２），…，θｉ（ωＮ））ごと平均値が統合角度θｉとされ、あるいは分別角度（θｉ（ω１），θｉ（ω２），…，θｉ（ωＮ））ごとの最頻値や中央値などを統合角度θｉとしてもよい。このように分別角度ごとに１つの角度に統合することにより、１つのある周波数の逆行列Ｈ（ωｎ）からのみ到来方向を推定する場合より、正しく到来方向を推定することができる。
なお図８に示した角度（方位）計算推定処理において、その１つでも選択した列について、角度θｉ（ωｎ）を求めることができなかった時は、そこでその逆行列Ｈ（ωｎ）に対する処理を終了し、他の周波数の逆行列Ｈ（ωｎ）に対する処理を行い、最初に全ての列について角度（方向）を計算することができたら、その時の計算結果を推定（方向）θ１，…，θ１としてもよい。あるいは各周波数の逆行列Ｈ（ωｎ）中で全ての列について角度計算ができた結果を角度ごとに分別し、統合してもよい。あるいは１つの列について角度θｉ（ωｎ）を計算することができないものがあれば、その状態が最初に発生した時に、その後の処理を全て終了し、あらためて観測信号を求める処理からやりなおすことにより、推定結果の信頼性を高めるようにしてもよい。
この第１実施形態の実験例を述べる。残響時間が１９０ｍｓの部屋に３つのマイクロホンを５６．６ｍｍ間隔で１列に配列し、その配列方向を基準として角度４８°、７３°、１１９°の方向に３つの音源を配置した。これら音源からの音響信号の混合時間を６秒、観測信号の標本化周波数を８ｋＨｚ、空間エイリアシングが生じない最大周波数を３ｋＨｚ、短時間フーリエ変換のフレームを１０２４サンプルとした。前述したようにして、各周波数ごとにそれぞれ計算した各角度を図９に、横軸を周波数、縦軸を方向として示す。図９中の◇，＋，□はそれぞれ３つの音源方向の推定計算値を示す。この結果を３つの角度範囲に分別し、各分別された角度平均値は４５°，７４°，１２３°となった。３音源の方向を３つのマイクロホンで推定することはＭＵＳＩＣ法ではできないが、この実施形態によれば可成り正しく各音源方向が推定されていることがわかる。
この実施形態の方法とＭＵＳＩＣ法とを比較するため、音源方向を４８°と１１９°と比較的大きく角差がある状態で２つの音源を配置した場合について同様に実験を行った。図１０Ａにこの実施形態の方法により得られた結果を、図１０ＢにＭＵＳＩＣ法により得られた結果をそれぞれ示す。これらよりいずれの方法も、かなり正しく方向を推定していることがわかる。この実施形態による結果を２つの方向範囲で分別平均した結果は４５°，１２３°であり、ＭＵＳＩＣ法はそれぞれ４５°，１２２°である。次に、音源方向を１０５°と１１９°と方向が比較的接近している状態で２つの音源を配置した場合について同様に実験を行った。この実施形態方法による結果及びＭＵＳＩＣ法による結果を図１１Ａ及び図１１Ｂにそれぞれ示す。ＭＵＳＩＣ法では大部分の周波数で音源方向の推定ができないが、この実施形態では大部分の周波数で角度計算を行うことができ、しかも、これらを角度範囲で分別して、それぞれ平均した値は１０５°，１２４°であり、ＭＵＳＩＣ法は９４°，１２８°でありこの実施形態による方向推定がかなり正しいことがわかる。
以上述べたようにこの実施形態によれば従来方法において、指向性パターンの利得が低い方向探索する場合と比較して、式（９）に値を代入するだけで推定方向が求まるため計算時間がかなり短かい。なお前述したようにＩＣＡにより求めた分離行列Ｗ（ω）はスケーリングの自由度及びパーミュテーションの自由度を含んでいるためこれらを解決した分離行列Ｗ’（ω）の逆行列を計算し、つまり真の混合行列Ａ（ω）を計算し、その混合行列Ａ（ω）の各列について二つの要素の比を用いて到来方向を推定することが考えられるかもしれない。しかし、真の混合行列Ａ（ω）そのものを求めることは、例えば信号源の信号ｓｉ（ｔ）の平均パワーを１にするなどの制約を与えない限りできない。無線通信分野においては信号源に対しそのような制約を与えることができるかも知れないが、例えば信号源の信号ｓｉ（ｔ）が人間が直接発声した音声信号である場合は、そのような制約は実質的にはできない。一方この第１実施形態によれば、スケーリング及びパーミュテーションの自由度を含んでいる分離行列Ｗ（ω）の逆行列Ｈ（ω）の各列について二つの要素の比をとることによりスケーリングの自由度の問題を解決しておりどのような信号源に対しても適用でき、しかも前記両問題を解決した分離行列を計算する必要がなくそれだけ計算時間が短い。更に前述したように、各周波数について求めた推定方向を予め決めた順に分別すればパーミュテーションの問題も簡単に解決できる。センサの数Ｊと同数であっても各信号源の方向を推定することができる。また音源方向が比較的接近していてもかなり正しく推定することができる。
［第２実施形態］
この第２実施形態は信号源の一つの位置情報として方向情報を求めるものであり、かつ少なくとも２次元に配置された少なくとも３個のセンサを用いて、信号源がいずれの方向に位置していても方向を推定でき、従ってブラインド信号分離におけるパーミュテーションの問題も比較的簡単に解決するものである。つまり方向情報に基づく円錐面を推定し、複数の円錐面の共通直線を推定して、方向情報を決定する。
この第２実施形態をブラインド信号分離装置に適用した機能構成を図１２に、その処理手順を図１３に示す。例えば４個のセンサ１１，１２，１３，１４が同一円上で等間隔に配置され、そのいずれの２つのセンサの間隔も、源信号の最短波長の１／２以下とされる。以下の説明はセンサの数はＪ個であり、Ｊ≧３として行う。第１実施形態と同様に各センサｊ（ｊ＝１，…，Ｊ）で観測された観測信号ｘｊ（ｔ）は、それぞれ周波数領域変換部１１で例えば、短時間フーリエ変換によって、周波数領域の信号Ｘｊ（ω，ｍ）にそれぞれ変換される（ステップＳ１１）。
これら周波数領域の信号Ｘｊ（ω，ｍ）から、分離行列算出部１２で、独立成分分析により、各周波数ごとの分離行列
2004079388.tif 000019
が算出される（ステップＳ１２）。
逆行列算出部１３でこれら周波数ごとの各分離行列Ｗ（ω）の逆行列Ｈ（ω）がそれぞれ算出される（ステップＳ１３）。
2004079388.tif 000020
この第２実施形態では、円錐面推定部１４において、周波数ごとの各逆行列Ｈ（ω）における各列ごとに異なる複数の要素対について各要素の比からこれら要素と対応する二つのセンサのセンサペア軸を中心軸とし、何れかの信号源が存在する円錐面、つまり一つの混合行列Ｈ（ω）の各列ごとに複数の円錐面がそれぞれ推定される（ステップＳ１４）。
円錐面推定部１４の機能構成は図１５中の角度計算部１４とほぼ同様であり、その処理手順も図８の手順と似ている。円錐面推定部１４において行われるある一つの周波数の逆行列Ｈ（ω）に対するステップＳ１４の処理の具体例を図１４を参照して説明する。
まず、例えば、円錐面推定部１４内のレジスタ内に格納されている制御パラメータｉとｐを０に初期化する（ステップＳ２０）。ここで、ｉは各信号源の番号と対応し、ｐはｉごとの推定済み円錐面の個数である。
ｉを＋１し（ステップＳ２１）、ｇを＋１し（ステップＳ２２）、制御パラメータｊ，ｊ’はそれぞれ互いに異なるｊ以下の自然数を、例えばランダムに選択する（ステップＳ２３）。この制御パラメータｊ，ｊ’の組は、同一のｉに対しては一度選択されたｊ，ｊ’の組み合わせの選択は行わないようにする。すなわち、ｊ≠ｊ’であり、また、例えば、ｉ＝１について、一度（ｊ，ｊ’）＝（１，２）が選択されると、ｉ＝１については、この処理が終了するまで、再びステップＳ２３において（ｊ，ｊ’）＝（１，２）という選択は行われないものとする。（また、この選択は、選択するｊ，ｊ’によって特定されるセンサペア軸（センサｊとセンサｊ’とを通る直線）が、このルーチン処理においてそれ以前に選択されたｊ，ｊ’によって特定されるセンサペア軸と同一直線上に配置されないように行われることが望ましい。これにより、円錐面推定部１４は、一定の誤差範囲内で中心軸が重ならない複数の円錐面を推定することになる。なお、各センサペア軸が同一直線上となるか否かの判断は、例えば、センサ情報格納部１５に各センサの位置を示すベクトルを格納しておき、ここから各センサの位置を示すベクトルの情報を抽出して行う。）
次に、センサ情報格納部１５から、ステップＳ２３で選択したパラメータｊに対応するｊ番目のセンサｊの位置を示すベクトルｄｊと、変数ｊ’に対応するｊ’番目のセンサｊ’の位置を示すベクトルｄｊ’とを抽出する（ステップＳ２４）。また、混合行列Ｈ（ω）から、ｊ行ｉ列要素Ｈｊｉ（ω）と、ｊ’行ｉ列要素Ｈｊ’ｉ（ω）とを選択して抽出する（ステップＳ２５）。以上の処理は図７中の選択部１４ａにより行われる。従って、ｉ，ｐ，ｊ，ｊ’や選択したｊ，ｊ’などを格納するレジスタも選択部１４ａに設けられることになる。
抽出したこれらの情報を用い、次式（９’）を演算する（ステップ２６）。
2004079388.tif 000021
ここで‖ｄｊ−ｄｊ’‖はセンサ１ｊと１ｊ’との間隔（距離）である。この第２実施形態においては複数のセンサが２次元又は３次元に配置されているため各センサの位置情報は例えば図１２中のセンサ１１〜１４が配置されている円の中心を原点とする２又は３要素の座標ベクトルで表わされる。つまり式（９）はセンサが１次元に配置され、信号到来方向の角度が２次元の場合であるが、この式（９’）はセンサが２次元あるいは３次元に配置されていてもよく信号到来方向の角度が３次元空間でもよい場合に拡張されたものである。従って式（９’）は式（９）を包含する。この式（９’）により推定された角度θ＾ｉ，ｊｊ’（ω）及びそのｉ，ｊ，ｊ’を円錐面情報として円錐面推定部１４内のレジスタ（メモリ）に一時記録する（ステップＳ２７）。なお図１２中に破線で示すように各周波数に対する間隔計算部１４ｃは共通に用いられる。また式（９）で演算した角度θｉ（ω）は三次元空間においてセンサペア軸（センサ１ｊと１ｊ’とを結ぶ直線）に対する角度がθｉ（ω）となる直線の無数の集合、つまり円錐面に信号源ｉが存在することを推定していることになる。この拡張された式（９’）の演算結果を単なるθｉ（ω）ではなくθ＾ｉ，ｊｊ’（ω）と表わした。ステップＳ２６で行う演算は図７中の偏角計算部１４ｂ、間隔計算部１４ｃ、位相回転計算部１４ｄ、除算部１４ｅ、判定部１４ｆ及び逆余弦計算部１４ｇにより、図８中のステップＳ４ｃ，Ｓ４ｄ，Ｓ４ｅ及びＳ４ｆの処理により行う。
次に、ｐ＝Ｐであるか否かを判断する（ステップＳ２８）。このＰは、ｉごとに推定すべき円錐面の個数であり、このステップでは、そのｉについて円錐面をＰ個推定したか否かを判断する。ｐ＝ＰでなければステップＳ２２に戻り、ｐ＝Ｐであれば、ｉ＝Ｉであるか否かを判断する（ステップＳ２９）。すなわち、すべてのｉについて円錐面の推定が終了したか否かを判断し、ｉ＝ＩでなければステップＳ２１に戻り、ｉ＝Ｉであれば処理を終了する（ステップＳ１４の具体例の説明終了。）。第１実施形態では選択したｉについて角度θｉ（ω）を１つ推定すれば次のｉについての角度推定に移るがこの第２実施形態では各ｉごとに複数Ｐ個の組の角度（円錐面）θ＾ｉ，ｊｊ’（ω）を推定する。ステップＳ２８及びＳ２９の処理は図７中の判定部１４ｆで行われる。
図１２中の到来方向決定部１６で、円錐面推定部１４で推定された複数の円錐面の情報（この例では、ｉ，ｊ，ｊ’、θ＾ｉ，ｊｊ’（ω））から源信号の到来方向ｕｉ（ω）＝（方位角θｉ（ω），仰角φｉ（ω））（ｉ＝１，…，Ｉ）が決定される（図１２、ステップＳ１５）。具体的には、例えば、ｉについて推定された複数の円錐面のうち、互いに線接触しているとみなされる共通直線の方向を、信号源ｉに対応する信号の到来方向ｕｉ（ω）として算出する。
この信号の到来方向ｕｉ（ω）の推定方法を図１５を参照して説明する。
センサ１１と１２の配列方向と直角方向にセンサ１３が配され、センサ１１と１２間、センサ１２と１３間の各間隔は同一とされている。ｉ＝２の信号源２２に対し、センサ１１及び１２によりセンサペア軸３１２を中心軸とする円錐面４１２（θ＾２，１２（ω））が、センサ１２と１３によりセンサペア軸３２３を中心軸とする円錐面４２３（θ＾２，２３（ω））が、センサ１１と１３によりセンサペア軸３１３を中心軸とする円錐面４１３（θ＾２，１３（ω））がそれぞれ推定されている。これらのうち、円錐面４１２と４１３は、共通直線５２で互いに線接触しているとみなされる。この共通直線５２の方向ｕ２（ω）が、信号源２２から放射された信号の到来方向ｕ２つまり信号源２２の方向と決定する。なお円錐面４２３も平行移動すれば円錐面４１２，４１３とほぼ線接触するが、後述する第４実施形態により円錐面４２３は廃棄される。
図１２中の到来方向決定部１６におけるこの共通直線５ｉの方向θｉ（ω）を決定するための具体的計算方法の例を説明する。角周波数ωに対して信号源２ｉが存在すると推定された複数の円錐面を４ｊｊ’（１），…，４ｊｊ’（Ｐ）とし、これら円錐面４ｊｊ’（ｐ）（ｐ＝１，…，Ｐ）の推定に用いられた１対のセンサの位置情報をｄｊ（ｐ），ｄｊ’（ｐ）とし、角周波数ωについて推定された円錐面４ｊｊ’（ｐ）と対応する角度をθ＾ｊｊ’（ω，ｐ）とし、円錐面４ｊｊ’（ｐ）を表わすベクトルをｕとする。
正規化軸ベクトル計算部１６ａで前記各一対のセンサの位置間を結ぶ軸ベクトル（ｄｊ（ｐ）−ｄｊ’（ｐ））をそれぞれ長さｌに正規化する。つまり
2004079388.tif 000022
をそれぞれ演算する。このｖｐと円錐面ベクトルの内積はこれらベクトルのなす角度の余弦とする。つまり
2004079388.tif 000023
関係が成立する。知りたいのは共通直線５ｉの方向のみであるから円錐面ベクトルｕは単位ベクトル、つまり‖ｕ‖＝１とする。全ての円錐面と共通する直線の方向を求めるには
2004079388.tif 000024
として次式の連立方程式をｕについて解けば良い。
2004079388.tif 000025
一般にはこの連立方程式は解が存在しないか、もしくは一意には決まらない。そこで‖Ｖｕ−ｃ＾（ω）‖を最小化するｕを求めて前記連立方程式の解、つまり共通直線５ｉの方向ｕｉ（ω）とする。この誤差を最小化するｕを求める計算が演算部１６ｂで行われる。この方向ｕｉ（ω）は３次元での方向であるから、極座標表示により方位角θｉ（ω）及び仰角φｉ（ω）が求まることになる。
もっと計算を簡便にするには次のようにしてもよい。図１６に示すように正規化軸ベクトル計算部１６ａで円錐面推定に用いた各センサ対についての正規化軸ベクトルｖｐ（ｐ＝１，…，Ｐ）を求め、逆行列計算部１６ｃでＶ＝（ｖ１，…，ｖＰ）ＴのＭｏｏｒｅ−Ｐｅｎｒｏｓｅ型一般逆行列Ｖ＋を計算し、このＶ＋と推定円錐面と対応する余弦ベクトルｃ＾（ω）＝（ｃｏｓθ＾ｊｊ’（ω，１）…ｃｏｓθ＾ｊｊ’（ω，Ｐ））Ｔとを用いて最小ノルムな最小２乗誤差解を求め、大きさを正規化して近似解としてもよい。つまり演算部１６ｄで
2004079388.tif 000026
を計算する。
このようにして複数の推定円錐面の共通直線とみなされる直線の方向が各周波数ごとに、かつ各信号源ごとに求められる。
図１２中のパーミュテーション解決部１７は、到来方向決定部１６で決定された到来方向ｕｉ＝（θｉ，φｉ）を用いて、分離行列算出部１２で算出された分離行列Ｗ（ω）の行の並び替えを行ってパーミュテーション問題が解決された分離行列を生成する。
パーミュテーション解決部１７で行われる具体例としては、到来方位角θｉについて以下のように並べ替えを行い、解決できなかった列について同様に到来仰角φｉによる並べ替えを行う。つまり算出決定された到来方位角θｉ（ω）がどの周波数においても既定の順序、例えば、θ１，θ２，…，θＩが昇順となるように逆行列Ｈ（ω）の列を入れ替え、入れ替えができなかった列について仰角φｉがどの周波数でも昇順になるように逆行列Ｈ（ω）の列を入れ替える並替行列が並替行列生成部１７ａで生成され、その並替行列の逆行列Ｐ（ω）が逆行列生成部１７ｂで生成され、並替部１７ｃでその逆行列Ｐ（ω）が、分離行列Ｗ（ω）に左から掛算される。並替行列生成部及び逆行列生成部１７ｂはパーミュテーション行列生成部を構成している。
並替行列生成部１７ａでの処理を具体的に述べる。この例では逆行列Ｈ（ω）の第１列目で（θ１（ω），φ１（ω））が、第２列目で（θ２（ω），φ２（ω））が，…、第Ｉ列目で（θＩ（ω），φＩ（ω））がそれぞれ式（９’）に基づき計算されたとする。いま到来方向決定部１６１から入力された到来方向を、前記昇順に並べたものと区別するためにθおよびφの各添字にダッシュ「’」を付けて、（θ１’（ω），φ１’（ω）），（θ２’（ω），φ２’（ω）），…，（θＩ’（ω），φＩ’（ω））と表わし、これらを例えばまずθｉ’（ω）について昇順に並べた時に、例えば（θ３’（ω），φ３’（ω））＞（θＩ’（ω），φＩ’（ω））＞（θ２’（ω），φ２’（ω））＞…であれば、逆行列Ｈ（ω）の３列目を１列目に、Ｉ列目を２列目に、２列目を３列目になるように移動し、残りの列も同様に移動し、θｉ’（ω）が同一のものについてはφｉ’（ω）が昇順になるように列を移動する。そのような列の移動並べ替えを行う並替行列を作る。このような並べ替えを行う行列を作ることは従来からよく行われていることである。このようにして全ての周波数について得られた到来方向（θ１（ω），φＩ（ω）），…，（θＩ（ω），φ２（ω））を用いて、並替行列が作られ、更にその逆行列つまりパーミュテーション行列Ｐ（ω）が算出される（図１３、ステップＳ１６）。
この算出したパーミュテーション行列Ｐ（ω）が並替部（１７ｃ）で分離行列Ｗ（ω）の左から掛けられ、その結果Ｗ’（ω）＝Ｐ（ω）Ｗ（ω）がパーミュテーション問題が解決された分離行列として出力される（ステップＳ１７）。つまりこの分離行列Ｗ’（ω）はいずれの周波数のものでも、１行目は信号源２１の信号を分離する要素、２行目は信号源２２の信号を分離する要素、以下同様に同一行の要素は同一信号源からの信号を分離する要素となる。
分離行列Ｗ’（ω）は、時間領域変換部１８で、例えば、逆フーリエ変換によって、時間領域の分離フィルタ係数群
2004079388.tif 000027
に変換され、これら分離フィルタ係数群は信号分離部１９に設定される。
信号分離部１９は、各センサからの観測信号ｘ１（ｔ），…，ｘＪ（ｔ）と分離フィルタ係数群とにより式（８）の演算を行い分離信号ｙ１（ｔ），…，ｙＩ（ｔ）を出力する。
図１２中に破線で示すように、パーミュテーション解決部１７よりの並替えされた分離行列Ｗ’（ω）と周波数領域変換部１１よりの周波数領域観測信号Ｘ（ω，ｍ）とにより式（５）の演算を周波数領域分離信号生成部１９’で行い、その結果の周波数領域分離信号Ｙ（ω，ｍ）＝Ｗ’（ω）Ｘ（ω，ｍ）を時間領域信号変換部１８’で時間領域信号ｙ１（ｔ），…，ｙＩ（ｔ）を生成してもよい。また到来方向決定部１６で決定された各周波数ごとの到来方向（θ１（ω），φ１（ω）），…，（θＩ（ω），φＩ（ω））を第１実施形態と同様に、つまり図５中のソート部３２で方向範囲ごとに分別し、各分別されたものごとの全周波数における到来方向を統合部３３でそれぞれ統合してもよい。
以上説明したように、この実施形態においても指向特性パターンの低ゲイン方向を探索することなく、式（９’）の演算により円錐面情報（θ＾ｉ，ｊｊ’（ω），φ＾ｉ，ｊｊ’（ω））を推定しているため計算量が少ない。しかも同一信号源に対し、円錐面を複数推定してその共通直線から信号の到来方向を算出しているため、０°から３６０°までのどの方向に信号源が存在しても、信号源３ｉの方向を一意に推定することができる。またその推定方向をパーミュテーション行列Ｐ（ω）の決定に利用しているため、信号源の位置に関わらず、パーミュテーション問題を適切に解決することができる。
〔第３実施形態〕
第３実施形態は、位置情報として一対のセンサと一つの信号源との各距離の比に基づくその信号源が存在する曲面を用いる。第１及び第２実施形態においては各信号源はセンサから遠方にあって、信号源からの信号はセンサには平面波として到来することを想定して処理した。しかし信号源とセンサとの距離が比較的短かい場合はセンサには信号が球面波として到来する。このような点から混合行列Ａ（ω）の要素の比Ａｊｉ（ω）／Ａｊ’ｉ（ω）を球面波（近距離場）モデルによって解釈すると、信号源ｉの方向以外の情報を推定することができる。
すなわち、近距離場モデルを用いると、周波数応答Ａｊｉ（ω）は以下のように表わせる。
2004079388.tif 000028
ここで、ｑｉは信号源ｉの位置を示すベクトルである。
このように表現された周波数応答について混合行列の同一列の２つの要素比Ａｊｉ（ω）／Ａｊ’ｉ（ω）をとり、その絶対値を計算すると次式となる。
2004079388.tif 000029
なお、｜β｜はβの絶対値を表わす。
この式（１０）を満すｑｉの点の無数の集合は信号源ｉが存在する曲面を与え、遠距離場（平面波）モデルを用いて推定された方向（又は円錐面）とあわせて用いることにより、センサから信号源ｉまでの距離を推定することができる。これにより、２つ以上の信号源が同一方向もしくは互いに近い方向にある場合であっても、距離が違っていれば、それらを区別でき、パーミュテーション問題を適切に解決することが可能となる。
この第３実施形態の例では、この信号源が存在する曲面を表わす位置情報と前記方向情報とを用いて分離行列に対するパーミュテーション問題を解決する。
第３実施形態の機能構成を図１７に、その処理手順を図１８にそれぞれ示す。センサは３個以上が２次元又は３次元に配置されるがこの実施形態においては例えば図２０に示すようにセンサ１１と１２の間隔に対しセンサ１２と１３の間隔は１０〜２０倍好ましくは１５倍程度とされる。先の実施形態と同様、各観測信号ｘ１（ｔ），…，ＸＪ（ｔ）は、それぞれ周波数領域の信号Ｘ１（ω，ｍ），…，ＸＪ（ω，ｍ）に変換され（ステップＳ１１）、更に独立成分分析により、各周波数ごとに分離行列Ｗ（ω）が算出される（ステップＳ１２）、その各分離行列Ｗ（ω）の逆行列として行列Ｈ（ω）が算出される（ステップＳ１３）。またこの例では第２実施形態と同様に周波数ごとの逆行列Ｈ（ω）の各列から選択した要素対を用い、円錐面が一つ好ましくは複数推定される（ステップＳ１４）。この第３実施形態においては距離比算出部３１において、周波数ごとの逆行列Ｈ（ω）から列ごとに選択した要素対を用いその対応するセンサと１つの信号源ｉとの距離比、つまり式（１０）と式（２１）から次式（１０’）を算出する（ステップＳ３５）。
2004079388.tif 000030
距離比算出部３１において行われるステップＳ３５の処理の具体例を図１９を参照して説明する。この処理は図１４の処理の流れとほぼ同様である。パラメータｉを０に初期化し（ステップＳ２０）、次に、ｉを＋１し（ステップＳ２１）、パラメータｊ，ｊ’としてＪ以下の自然数を、例えばランダムに選択しｊ≠ｊ’、かつ１度選択した組は選択しない（ステップＳ２３）。センサｊの位置ベクトルｄｊと、センサｊ’の位置ベクトルｄｊ’とを抽出し（ステップＳ２４）、逆行列Ｈ（ω）のｉ列目から、要素Ｈｊｉ（ω）と、Ｈｊ’ｉ（ω）とを、選択する（ステップＳ２５）。
この実施形態ではこれら選択した２つの要素の比ＤＲｉ，ｊｊ’（ω）を計算する（ステップＳ４１）。次に、ｉ＝Ｉであるか否かを判断し（ステップＳ２９）。ｉ＝ＩでなければステップＳ２１に戻り、ｉ＝２であれば処理を終了する。
距離比算出部３１で算出された距離比情報ＤＲｉ，ｊｊ’（ω）は、パーミュテーション解決部１７に送られ、パーミュテーション解決部１７は、到来方向決定部１６で推定された方向情報ｕｉ（ω）と、距離比算出部３１で算出された距離比情報ＤＲｉ，ｊｊ’（ω）とを用いて、分離行列算出部１２で算出された分離行列を行ってパーミュテーション問題を解決する。
Ｗ（ω）の行替えを行ってパーミュテーションの問題を解決する。例えば方向情報と距離比情報とから信号源２ｉの距離‖ｑｉ（ω）‖を距離推定部１７ｄで計算する（ステップＳ３６）。
この距離‖ｑｉ（ω）‖の計算方法を図２０を参照して説明する。信号源２１及び２２がセンサ１１，１２からみて同一方向Ｂにある。この場合センサ１１，１２及び遠距離場モデルで推定される信号源２１及び２２の方向ｕ１，ｕ２は同一の直線を定める。一方、間隔が大きいセンサ１２，１３及び近距離場モデルを用いて信号源２１が存在すると曲面６１は距離比
2004079388.tif 000031
から推定でき、つまり‖ｑ１（ω）‖が推定できる。また信号源２２が存在する曲面６２は距離比
2004079388.tif 000032
から推定でき、つまり‖ｑ２（ω）‖が推定できる。
信号源２１の位置は、直線ｕ１＝ｕ２と、曲面６１との共通部分上にあると推定でき、信号源２２の位置は、直線ｕ１＝ｕ２と、曲面６２との共通部分上にあると推定できる。例えば直線ｕ１＝ｕ２を表す式と曲面６１及び６２をそれぞれ表わす式との連立方程式を解くことにより‖ｑ１（ω）‖及び‖ｑ２（ω）‖を求めればよい。このようにして信号源の方向が同一や近似する場合であっても、信号源位置を区別することができる。
パーミュテーション解決部１７は、このようにして得られた各周波数ごとの信号源の距離‖ｑｉ（ω）‖が、所定の順序が例えば昇順になるように、分離行列Ｗ（ω）の行の入れ替えを行う。つまりこの入れ替えを行うためパーミュテーション行列Ｐ（ω）を算出する（ステップＳ３７）。そのパーミュテーション行列Ｐ（ω）の算出は実施形態２中のパーミュテーション解決部で行ったパーミュテーション行列の算出と同様に行えばよい。算出したパーミュテーション行列Ｐ（ω）を分離行列Ｗ（ω）の左から掛け、その結果Ｗ’（ω）＝Ｐ（ω）Ｗ（ω）を分離行列として出力する（ステップＳ３８）。
出力された分離行列Ｗ’（ω）は、時間領域変換部１８に送られ、時間領域変換部１８での信号分離に用いられる。
図２０から理解されるように、センサ１２と１３についてみると、信号源２１及びセンサ１２間の距離と信号源２１及びセンサ１３間の距離との差が大きく、一方信号源２２及びセンサ１２間の距離と信号源２２及びセンサ１３間の距離との差が小さい。従ってＤＲ１，２３＝‖ｑ１−ｄ３‖／‖ｑ１−ｄ２‖と数値１との差の絶対値は比較的大きく、ＤＲ２，２３＝‖ｑ２−ｄ３‖／‖ｑ２−ｄ２‖と数値１との差の絶対値は小さい。センサ１２と１３の間隔が大きい程、距離比ＤＲ１，２３（ω）とＤＲ２，２３（ω）との差が大きくなる。このため間隔が大きいセンサ対、この例では１２と１３が設けられる。
パーミュテーション解決部１７において、図１８中のステップＳ３７中に括弧で示すように計算した距離比ＤＲｉ，ｊｊ’（ω）を用いて、すべての周波数において逆行列Ｈ（ω）の１列目，２列目，３列目，…，Ｉ列目の順に、求まる距離比ＤＲｉ，ｊｊ’（ω）が昇順になるように並替行列を生成し、更に前記パーミュテーション行列Ｐ（ω）を生成してもよい。この場合、図１７中の距離推定部１７ｄは省略される。また到来方向決定部１６で方向情報ｕｉ（ω）を計算できなかったり、逆行列Ｈ（ω）の複数列について同一のｕｉ（ω）が計算された信号源ｉとセンサｉｊ’についてのみ、距離比ＤＲｉ，ｊｊ’（ω）を計算し、これら自体により、あるいは更に求めた‖ｑｉ（ω）‖により分離行列Ｗ（ω）の行の入れ替えを行ってもよい。つまり、方向情報ｕｉ（ω）を用いて分離行列Ｗ（ω）の行の入れ替えを行い、その後、ｕｉ（ω）によって入れ替えができなかった行については距離比ＤＲｉ，ｊｊ’（ω）又は‖ｑｉ（ω）‖を用いてＷ（ω）の行の入れ替えを更に行ってもよい。実際には両入れ替えを同時に行うパーミュテーション行列Ｐ（ω）を生成する。
距離比ＤＲｉ，ｊｊ’（ω）をまず求め、分離行列Ｗ（ω）の行の入れ替えを行い、距離比ＤＲｉ，ｊｊ’（ω）を求めることができなかったものについて方向情報ｕｉ（ω）を用いて分離行列Ｗ（ω）の入れ替えを更に行ってもよい。この場合も両入れ替えを同時に行うパーミュテーション行列Ｐ（ω）を生成する。一般にはＤＲｉ，ｊｊ’（ω）よりｕｉ（ω）の方が高い精度で得られるため、ｕｉ（ω）による入れ替えを主とし、できないものについてＤＲｉ，ｊｊ’（ω）による入れ替えを行う方がよい。
また式（１０’）をその右辺の値を距離比ＤＲｉ，ｊｊ’（ω）と表わして、ｑｉについて解くと、中心Ｏｉ，ｊｊ’（ω）及び半径Ｒｉ，ｊｊ’（ω）がそれぞれ下記式（１１）及び（１２）で与えられる球面となる。
2004079388.tif 000033
例えばセンサ１ｊ，１ｊ’の各位置がｄｊ＝（０，０．１５，０），ｄｊ’＝（０，−０．１５，０）（単位はメートル）である場合、半径Ｒｉ，ｊｊ’をパラメータとしたとき、式（９）によって決まる球面の様子は図２１に示すようになる。
つまり信号源ｉは位置情報としての式（１１）及び（１２）で与えられる球面上に存在することになる。従って図１７中のパーミュテーション解決部１７中の距離推定部１７ｄを括弧書きで示すように曲面推定部とし、方向情報ｕｉ（ω）を求めることができなかった各ｉ，ｊｊ’について曲面推定部１７ｄで、式（１１）の半径Ｒｉ，ｊｊ’（ω）及び中心Ｏｉ，ｊｊ’（ω）を、図１８中のステップＳ３６中に括弧書きで示すようにそれぞれ計算し、これらＲｉ，ｊｊ’（ω）及び中心Ｏｉ，ｊｊ’（ω）がどの周波数の逆行列Ｈ（ω）でも同一順になるようにして、パーミュテーション行列Ｐ（ω）で求めてもよい。
この円錐面情報θ＾ｉ，ｊｊ’（ω）又は方向情報ｕｉ（ω）と球面情報ＤＲｉ，ｊｊ’（ω）又はＲｉ，ｊｊ’（ω）を用いてもパーミュテーション行列Ｐ（ω）を求めることができない場合はその周波数について、従来の相関法（例えば、Ｈ．Ｓａｗａｄａ他“Ａ  ｒｏｂｕｓｔ  ａｎｄ  ｐｒｅｃｉｓｅ  ｍｅｔｈｏｄ  ｆｏｒ  ｓｏｌｖｉｎｇ  ｔｈｅ  ｐｅｒｍｕｔａｔｉｏｎ  ｐｒｏｂｌｅｍ  ｏｆ  ｆｒｅｑｕｅｎｃｙ−ｄｏｍａｉｎ  ｂｌｉｎｄ  ｓｏｕｒｃｅ  ｓｅｐａｒａｔｉｏｎ．”，ｉｎ  Ｐｒｏｃ．Ｉｎｔｌ．Ｓｙｍｐ．ｏｎ  Ｉｎｄｅｐｅｎｄｅｎｔ  Ｃｏｍｐｏｎｅｎｔ  Ａｎａｌｙｓｉｓ  ａｎｄ  Ｂｌｉｎｄ  Ｓｉｇｎａｌ  Ｓｅｐａｒａｔｉｏｎ（ＩＣＡ  ２００３），２００３，ｐｐ．５０５−５１０参照）を適用するとよい。図１７中に破線で示すように、信号分離部１９で分離された出力信号ｙ１（ｔ），…，ｙＩ（ｔ）を周波数領域変換部３３で例えば周波数領域変換部１１と同様の方法で周波数領域の信号Ｙ１（ω，ｍ），…，ＹＩ（ω，ｍ）に変換し、相関算出部３４において、これら周波数領域信号Ｙ１（ω，ｍ），…，ＹＩ（ω，ｍ）中のパーミュテーション解決部１７でパーミュテーション行列を作ることができなかった周波数成分ｆａｎｋ成分と、パーミュテーション行列が得られた周波数中のｆａｎｋと隣接した周波数成分との相関を計算し、パーミュテーション解決部１７ではこの相関が大きくなるように周波数ｆａｎｋの分離行列の行の入れ替えを行い、この行入れ替えした分離行列に基づく分離出力信号ｙ１（ｔ），…，ｙＩ（ｔ）中のｆａｎｋ成分について再び相関計算部３４で相関を計算し、同様のことを計算した相関が最大になるまで繰り返す。この相関の最大値が所定値に達しなかった場合は、更にそのｆａｎｋの成分と、パーミュテーション行列が得られている周波数成分中のｆａｎｋに対する基本波又は高調波（一般には高調波）関係にあるものとの各相関の和を相関計算部３４で計算し、この相関和が大きくなるように、ｆａｎｋに対する分離行列の行入れ替えを修正部１７ｅで行うことを、前記相関和が最大になるまで繰り返す。なお第２実施形態で説明したように周波数領域の信号Ｘ（ω，ｍ）に対し分離行列Ｗ（ω）を適用し、周波数領域の分離信号Ｙ（ω，ｍ）を求める場合は、このＹ（ω，ｍ）を相関計算部３４の計算に用いればよい。
以下にこの第３実施形態の実験例を信号源として、室内で実測したインパルス応答を畳み込んで作成した混合音声を用いて分離実験を行った場合について説明する。図２２に示すように残響時間が１３０ミリ秒、３５５ｃｍ×４４５ｃｍ×２５０ｃｍ（高さ）の室内の水平面内でほぼ中心部で高さ１３５ｃｍの所にセンサとしての無指向性マイクロホン１１〜１４を直径４ｃｍの円上に等間隔で配し（これらマイクロホン１１〜１４は図中の左上に拡大して示した）、更に無指向性マイクロホン１５〜１８を直径３０ｃｍの円上に等間隔かつ、マイクロホン１１，１３，１５及び１７が、またマイクロホン１２，１４，１６及び１８がそれぞれ一直線上に配列するように配置した。マイクロホン１１から１５への方向を基準（０°）とし、マイクロホン配置の中心を原点とし、原点から１２０ｃｍ離れ、−３０°方向に音源２１を、＋３０°方向に音源２２を、＋９０°方向に音源２３を、−１５０°方向に音源２６をそれぞれ配置し、１５０°方向で原点からの距離６０ｃｍに音源２４を、距離１８０ｃｍに音源２５をそれぞれ配置した。サンプリング周波数ｋＨｚ、データ長６秒、フレーム長２０４８サンプル（２５６ミリ秒）、フレームシフト５１２サンプル（６４ミリ秒）とした。
周波数領域での分離行列Ｗ（ω）の各行のうち、マイクロホン対１１と１３、１２と１４、１１と１２、１２と１３に相当する行を用いて音源方向（円錐面）推定を行った。推定したヒストグラムを図２２に示す。図２３の横軸は推定方向を、縦軸は信号数である。これより推定方向は５つの群（クラスタ）があり、そのうち１５０°付近のクラスタは他のクラスタより２倍の大きさである。これは６個の音源中の２個が同一方向（１５０°付近）から到来していると推測される。４つの音源については推定した方向をもとにパーミュテーションを解決した。その結果を図２４に示す。図２４の横軸は周波数、縦軸は方向（度）を表わす。
残りの２個の音源については、広い間隔のマイクロホン対１７と１５、１７と１８、１６と１５、１６と１８を用いて音源が存在する球面の半径によって到来信号を区別した。マイクロホン対１７と１５を用いて推定した球面の半径を図２５に示す。図２５の横軸は周波数、縦軸は半径（ｍ）を表わす。
残響や推定誤差の影響により、位置情報だけではパーミュテーションを完全には解決することはできない。従って推定した位置情報に基づいて信号が矛盾なく分類できた周波数はその情報に基づいてパーミュテーション行列を生成し、残りの周波数については相関に基づく方法を用いてパーミュテーションを解決した。最後に時間領域の分離フィルタ係数を求める際にスペクトルの平滑化を行った。スペクトル平滑化については例えばＨ．Ｓａｗａｄａ他、“Ｓｐｅｃｔｒａｌ  ｓｍｏｏｔｈｉｎｇ  ｆｏｒ  ｆｒｅｑｕｅｎｃｙ−ｄｏｍａｉｎ  ｂｌｉｎｄ  ｓｏｕｒｃｅ  ｓｅｐａｒａｔｉｏｎ”，ｉｎ  Ｐｒｏｃ．ＩＷＡＥＮＣ  ２００３，２００３，ｐｐ．３１１−３１４参照されたい。分離性能（ＳＩＲ）の評価結果を図２６に示す。図２６中の数値の単位はｄＢである。Ｃは相関法のみによってパーミュテーションを解決したもの、Ｄ＋Ｃは方向（円錐面）推定でパーミュテーションを解決し、解決できないものに相関法を適用したものＤ＋Ｓ＋Ｃは方向（円錐面）推定と球面推定とによりパーミュテーション解決を行った後、解決できなかった周波数について、相関法を用いたものである。この後者の方法によれば、６音源すべての分離ができ、ＳＩＲの改善は入力ＳＩＲから平均１７．１ｄＢとなった。
前記した第２実施形態の例ではセンサは２次元配置したが、一対のセンサにより推定する球面はこれらセンサの中心２等分線と対称に現われるため、信号源が３次元に存在している場合は２次元配置センサでは判別できなくなるため、センサも３次元配置とする必要がある。
以上説明したように、この第３実施形態においても式（９’）により円錐面情報を、また式（１０’）により曲面情報をそれぞれ推定しているため計算量が少ない。しかも円錐面と、距離比ＤＲｉ，ｊｊ’（ω）、又は距離‖ｑｉ（ω）‖あるいは球面の半径Ｒｉ，ｊｊ’（ω）とによりパーミュテーションを解決しているために２つ以上の信号源が同一方向もしくは互いに近い方向にある場合であっても、それらを区別できる。更にこれに相関法を加えれば一層確実に分離が可能である。なお球面情報としては計算が簡単であるからＤＲｉ，ｊｊ’（ω）が好ましい。
〔第４実施形態〕
この第４実施形態では、推定した円錐面の信頼性を検証し、信頼性が高いと判断した円錐面を分離行列のパーミュテーションの解決に用いる。例えば図２７に示すように円錐面推定部１４で推定された円錐面情報θ＾ｉ，ｊｊ’（ω）は円錐面検証部４１で許容角情報格納部４２よりの許容角情報に基づき、信頼性があるか否かの検証がなされる。つまり式（９’）により求まる角度θ＾ｉ，ｊｊ’（ω）はこれを求めるために用いたセンサ１ｊ及び１ｊ’の配列方向に対する相対角度であり、図４Ｂを参照して説明したように、角度θ＾ｉ，ｊｊ’（ω）が０度付近及び１８０度付近の場合、正しいパーミュテーションが行われないと推定される。
従って、パーミュテーションを正しく行うことができると推定される角度の最小値θｍｉｎと最大値θｍａｘが許容角情報として許容角情報格納部４２に格納され、円錐面検証部４１で推定円錐面情報θ＾ｉ，ｊｊ’（ω）がθｍｉｎとθｍａｘとの間にあれば信頼できる円錐面と判定されて出力され、つまりパーミュテーション解決に用いられるが、θｍｉｎとθｍａｘとの間になければ、そのθ＾ｉ，ｊｊ’（ω）は信頼性がないものとして破棄され、パーミュテーション解決に用いられない。例えば図１５中の円錐面４１３は破棄される。
円錐面検証部４１で信頼性があると検証された円錐面情報θ＾ｉ，ｊｊ’（ω）は図１２中の到来方向決定部１６へ送られ、又は図１７で説明したように、パーミュテーション解決部１７へ直接送られる。図１３及び図１８の処理手順において、ステップＳ１４の次に破線で示すようにステップＳ１００として推定円錐面θ＾ｉ，ｊｊ’（ω）が信頼性があるか否かの検証が行われ、信頼性があるものについてのみ次のステップに移るようにされる。式（９’）において計算された偏角ａｒｇ（Ｈｊｉ／Ｈｊ’ｉ）に含まれる誤差をΔａｒｇ＾、推定された角度θ＾ｉに含まれる誤差をΔθ＾とすると、これらの比｜Δθ＾／ΔａｒｇＨ＾｜は式（９’）を偏微分して次式のように近似できる。
2004079388.tif 000034
式（１３）をいくつかの周波数について計算した結果を図４Ｃに示す。この図４Ｃから推定された角度θ＾ｉの方向がセンサ１ｊ，１ｊ’のセンサペア軸に近い場合はａｒｇ（Ｈｊｉ／Ｈｊ’ｉ）に含まれる誤差ΔａｒｇＨ＾は推定角度θ＾ｉに対し、大きな誤差を生じさせることがわかる。つまりセンサペア軸に近い方向の推定角度θ＾ｉを用いて分離行列Ｗ（ω）に対するパーミュテーションの問題を解決した場合、これは正しい解にならない可能性が高い。図４Ｂ及び図４Ｃから例えばθｍｉｎは２０°程度、θｍａｘは１６０°程度とされる。図４Ｃからわかるように｜Δθ＾／ΔａｒｇＨ＾｜は周波数により可成り異なり、低い周波数ではΔａｒｇＨ＾が到来方向の誤差に大きく影響する。従って、低周波数については、信号源が存在する球面の推定に基づく情報ＤＲｉ，ｊｊ’、‖ｑｉ‖、Ｒｉ，ｊｊ’を用いて、あるいは相関法によりパーミュテーション問題を解決するとよい。
この第４実施形態によれば、推定円錐面情報中の信頼性のないものは破棄されるため、これにより悪影響を受けることなく、それだけ正確に到来方向を推定でき、従って正しいパーミュテーション行列Ｐ（ω）を生成することができ、つまり分離信号のＳＩＲ（性能）が向上する。
〔第５実施形態〕
第５実施形態では位置情報として距離比、又はこれより推定した球面情報を用いる。この機能構成は図１７中に示されている。図２８にこの実施形態の処理手順を簡略に示す。この場合はセンサは比較的広い間隔、例えば図２２に示した信号源が音源の場合、３０ｃｍとして少くとも２次元に配置される。
先の実施形態と同様に時間領域の観測信号は図１７中に示すように周波数領域の信号に変換され、分離行列生成部１２により分離行列Ｗ（ω）が生成され（ステップＳ５１）、これより逆行列生成部１３で逆行列Ｈ（ω）が生成される（ステップＳ５２）。各周波数における逆行列Ｈ（ω）の各列ごとに球面情報が推定される（ステップＳ５３）。この球面情報は第３実施形態で説明した方法と同様に算出される。つまり球面情報は距離算出部３１により算出された距離比ＤＲｉ，ｊｊ’（ω）、あるいは距離推定部又は曲面推定部１７ｄで算出された‖ｑｉ‖、又は半径Ｒｉ，ｊｊ’（ω）及び中心Ｏｉ，ｊｊ’（ω）である。
これら球面情報を用いて、その配列順が予め決めた順になるよう分離行列Ｗ（ω）に対するパーミュテーション行列が生成され分離行列Ｗ（ω）の行の入れ替えが行われる（ステップＳ５４）。第３実施形態におけるこの処理はパーミュテーション解決部１７で行われるが、ただしこの第５実施形態では球面情報のみが用いられる。この球面情報のみでパーミュテーションを解決できない周波数があれば、その周波数に対しては前述した相関法による解決が行われる（ステップＳ５５）。
図２２に示した室内に配したマイクロホン１６及び１８を用い、１２０°方向において原点から６０ｃｍ及び１５０ｃｍの所に音源２４及び２５を配し、実験室内で測定したインパルス応答に話者４名の音声を畳み込んで作成した１２通りの組合せの混合音声と対する分離実験を行った。パーミュテーション解決は推定された半径Ｒ４，６８（ω）とＲ５，６８（ω）を比較してＲ４，６８（ω）≦Ｒ５，６８（ω）となるようにして行った。Ｒ４，６８及びＲ５，６８中の最大値が、Ｒ４，６８及びＲ５，６８中の最小値にしきい値Ａｔｈを乗算した値以上、つまりｍａｘ（Ｒ４，６８，Ｒ５，６８）≧Ａｔｈ・ｍｉｎ（Ｒ４，６８，Ｒ５，６８）を満す周波数は音源位置（Ｒ４，６８（ω），Ｒ５，６８（ω））によるパーミュテーションの解決が信頼できると判断して、音源位置による方法を適用し、それ以外の周波数では相関による方法を適用した。従ってＡｔｈ＝１．０のときは全て音源位置による方法を適用しＡｔｈが無限大のときは全て相関による方法を適用することになる。しきい値Ａｔｈを変化させながら分離性能としてＳＩＲ（信号干渉比）を音声の組合せごとにプロットした結果を図２９に示す。この結果より音源位置のみによる方法では全体的に性能が悪く、相関のみによる方法では分離性能がばらつき、両者を併用すると安定して高い性能が得られることがわかる。またしきい値Ａｔｈとしては比較的大きい値が好ましく、８〜１６で全体の１／５〜１／１０程度の周波数を音源位置によって決定できることがわかる。
また図２２に示した実験条件で音源位置による方法と相関による方法を併用した実験結果を図２６中のＤ＋Ｃの欄に示す。これよりこの併用方法によれば可成り性能よく分離が行われることがわかる。この結果は式（９’）を用いる到来方向（円錐面情報）による方法と相関による方法とを併用した場合と同様の傾向がある。なお、到来方向による方法と相関による方法とを併用した場合の実験結果は、相関による方法の説明で引用した文献に示されている。
この第５実施形態によれば式（１０’）により球面情報を求めているため計算量が少ない。なお球面情報としては距離比ＤＲｉ，ｊｊ’（ω）が好ましい。
［第６実施形態］
第６実施形態は例えば一つの推定円錐面に基づき、パーミュテーションの問題を解決する。図１２中に破線で示すように円錐面推定部１４で推定された円錐面θ＾ｉ，ｊｊ’（ω１）、・・・    θ＾ｉ，ｊｊ’（ωＮ）、はパーミュテーション解決部１７へ直接入力され、図１３中に破線で示すようにステップｓ１４で設定された円錐面θ＾ｉ，ｊｊ’（ωｎ）（ｎ＝１、・・・、Ｎ）がステップＳ１５の処理をすることなくステップＳ１６でどの周波数においても例えばθ＾ｉ，ｊｊ’（ωｎ）が昇順になるようなパーミュテーション行列Ｐ（ω）が算出される。
この場合、ステップＳ１４で信号源ｉごとに１つの円錐面を推定するのみでもよい。また図に示していないが、行入れ替えができなかったものについては、第３実施形態で説明したように相関法による行入れ替えを得るようにしてもよい。
この第６実施形態によれば逆行列Ｈ（ω）の各列ごとの二つの要素比をとることによりスケーリングの問題が簡単に除去され、式（９’）の演算をすればよく計算時間が短くて済む。
［信号分離のまとめ］
以上述べたブラインド信号分離におけるパーミュテーションを解決した分離行列の生成方法を、逆行列の生成以後の処理を図３０Ａ、図３０Ｂ、図３０Ｃに示す。図３０Ａでは逆行列の各列の要素比から円錐面を推定し（ステップＳ６１）、必要に応じて信頼性のない円錐面を破棄し（ステップＳ６２）、これらの複数の円錐面から共通直線の方向を決定し（ステップＳ６３）、共通直線方向を用いてパーミュテーション行列Ｐ（ω）を生成して分離行列の入れ替えを行い（ステップＳ６４）、パーミュテーションを解決できなかった周波数に対しては相関法による分離行列の行入れ替えを行う（ステップＳ６５）。図３０Ａ中に破線で示すように、ステップＳ６からステップＳ６４に直ちに移り、円錐面を直接用いてパーミュテーション行列Ｐ（ω）を生成してもよい。
図３０Ｂでは前記ステップＳ６１の円錐面推定、必要に応じて前記ステップＳ６２の円錐面破棄を行い、ステップＳ６３の共通直線方向推定を行い、これら共通直線方向を用いて逆行列の列の入れ替えを行う（ステップＳ６６）。この際、図中破線で示すようにステップＳ６１で推定された円錐面又はステップＳ６２で処理後の円錐面を直接、ステップＳ６６の処理に用いてもよい。次に円錐面の推定又は共通直線方向の決定ができなかった、又は円錐面或いは共通直線方向が不確かなものと対応する、或いは同一値となったものについて、その逆行列の列の要素比から球面情報を推定し（ステップＳ６７）、この推定した球面情報を用いて逆行列の入れ替えを行ってパーミュテーション行列Ｐ（ω）を生成して分離行列の行入れ替えを行う（ステップＳ６４）。更にこのパーミュテーション行列Ｐ（ω）を作ることができなかった周波数に対しては相関による方法を適用する（ステップＳ６５）。
図３０Ｃはまず逆行列の各列の要素比から球面情報を推定し（ステップＳ６８）、この球面情報を用いて逆行列の列の入れ替えを行い（ステップＳ６９）、この列の入れ替えを行うことができなかった、又は球面情報が不確かなものについて、その逆行列の列の要素比から円錐面を推定し（ステップＳ７０）、これらの複数の円錐面の共通直線方向を決定し（ステップＳ７１）、この決定した方向又はステップＳ７０で推定した円錐面を直接用いて、逆行列の列の入れ替えを行ってパーミュテーション行列Ｐ（ω）を生成し、分離行列の入れ替えを行う（ステップＳ６４）。パーミュテーション行列を作ることができなかった周波数に対して相関法を適用する（ステップＳ６５）。
これらに対し更に図２８に示した方法がある。
第３実施形態、第５実施形態及び第６実施形態においても、第２実施形態で説明したように、周波数領域の分離行列Ｗ（ω）と観測信号Ｘ（ω）とを用いて信号分離を行い、その後、分離された周波数領域信号Ｙ（ω）を時間領域信号ｙ（ｔ）に変換してもよい。
第１実施形態の説明では２次元における信号の到来方向の推定を行ったが、第２実施形態中でも説明したように、３次元における信号の到来方向の推定にも適用できる。また第２乃至第６実施形態を２次元における信号分離に適用することもできる。この場合は推定円錐面θ＾ｉ，ｊｊ’（ω）はその推定に用いたセンサ対のセンサペア軸に対し対称の２方向となり、推定球面Ｒｉ，ｊｊ’（ω）又はＤＲｉ，ｊｊ’（ω）は円の半径又はこれとほぼ対応するものになる。
図５、図１２、図１７にそれぞれ示した装置、また第５実施形態による装置はそれぞれコンピュータに機能させることもできる。この場合、それぞれと対応した処理手順、つまり図６、図１３、図１８などに示した各過程をコンピュータに実行させるためのプログラムを記録した磁気ディスク、半導体メモリ、ＣＤ−ＲＯＭなどから、コンピュータ内のメモリにインストールし、又は通信回線を介してコンピュータ内のメモリに前記プログラムをダウンロードして、コンピュータにそのプログラムを実行させればよい。
【図１】
2004079388.tif 000035
【図２】
2004079388.tif 000036
【図３】
2004079388.tif 000037
2004079388.tif 000038
【図５】
2004079388.tif 000039
【図６】
2004079388.tif 000040
【図７】
2004079388.tif 000041
【図８】
2004079388.tif 000042
【図９】
2004079388.tif 000043
2004079388.tif 000044
2004079388.tif 000045
【図１２】
2004079388.tif 000046
【図１３】
2004079388.tif 000047
【図１４】
2004079388.tif 000048
【図１５】
2004079388.tif 000049
【図１６】
2004079388.tif 000050
【図１７】
2004079388.tif 000051
【図１８】
2004079388.tif 000052
【図１９】
2004079388.tif 000053
【図２０】
2004079388.tif 000054
【図２１】
2004079388.tif 000055
【図２２】
2004079388.tif 000056
【図２３】
2004079388.tif 000057
【図２４】
2004079388.tif 000058
【図２５】
2004079388.tif 000059
【図２６】
2004079388.tif 000060
【図２７】
2004079388.tif 000061
【図２８】
2004079388.tif 000062
【図２９】
2004079388.tif 000063
2004079388.tif 000064
【国際調査報告】
2004079388.tif 000065
2004079388.tif 000066
2004079388.tif 000067
2004079388.tif 000068