(57)【特許請求の範囲】
【請求項１】
Ｉ箇所の信号源から放射された信号を、Ｊ個のセンサで検出して上記信号源の位置情報を求める信号源位置情報推定装置であって、Ｉは２以上の整数、ＪはＩ以上の整数であり、
上記各センサの観測信号を周波数領域の信号に変換する周波数領域変換手段と、
上記周波数領域の信号から、独立成分分析により周波数ごとに、信号源信号を分離する第１分離行列をそれぞれ算出する分離行列算出手段と、
上記各第１分離行列の逆行列をあるいは擬似逆行列「以下、両者を単に逆行列という」をそれぞれ算出する逆行列算出手段と、
上記周波数ごとの逆行列の少なくとも１つに対し、その各列ごとにその二つの要素比から上記信号源の一つの位置情報を計算する位置情報計算手段とを具備する。
【請求項２】
請求項１記載の装置において、
上記位置情報計算手段は、上記周波数ごとの逆行列の複数について、上記各列ごとの要素比からの位置情報計算により上記各信号源の位置情報を求める手段であり、
上記周波数ごとの信号源ごとの位置情報から、上記複数の逆行列と対応する周波数の分離行列における位置情報と対応する行が予め決めた順序になるように行の入れ替えを行うパーミュテーション行列を生成するパーミュテーション行列生成手段と、
上記パーミュテーション行列と上記第１分離行列とを乗算して行を入れ替えた第２分離行列を求める並替手段とを備える。
【請求項３】
請求項２記載の装置において、
上記Ｊは３以上であり、かつＪ個のセンサは少なくとも２次元に配置され、
上記位置情報は上記センサから上記信号源への方向を含み、その信号源が存在する円錐面であり、
上記位置情報計算手段は、各列ごとの上記要素比から上記円錐面の計算を異なる二つの要素組の複数について行う手段と、周波数ごとの上記複数の円錐面の共通直線の方向を、上記位置情報として推定する到来方向決定手段を含む。
【請求項４】
請求項２記載の装置において、
上記Ｊは３以上であり、上記位置情報は一対のセンサから上記信号源への方向を含み、その信号源が存在する円錐面及びその信号源が存在する曲面であり、
上記位置情報計算手段は、上記二つの要素の比から上記円錐面を計算する手段と、上記円錐面の計算に用いた二つの要素と対応するセンサ間隔よりも大きいセンサ間隔のセンサと対応する二つの要素の比からこれらセンサと上記信号源との距離比を計算する手段と、上記距離比から上記曲面を計算する手段を備え、
上記パーミュテーション行列生成手段は周波数ごとに上記円錐面及び上記曲面に基づいて上記パーミュテーション行列を生成する手段である。
【請求項５】
請求項２、請求項３又は請求項４に記載の装置において、
上記Ｊは３以上であり、上記Ｊ個のセンサは少なくとも２次元に配置され、
上記位置情報は上記センサから上記信号源への方向を含み、その信号源が存在する円錐面であり、
上記円錐面を表わす角度が予め決めた第１角度と第２角度との間であるか否かを判定し、これら第１角度及び第２角度の間の円錐面を有効とする判定手段を備える。
【請求項６】
請求項２記載の装置において、
上記Ｊは３以上であり、上記Ｊ個のセンサは少なくとも２次元に配置され、
上記位置情報は信号源が存在する球面の半径であり、
上記位置情報計算手段は上記二つの要素比から距離比を計算する手段である。
【請求項７】
請求項２に記載の装置において、
上記位置情報は上記センサから上記信号源への方向を含みその信号源が存在する円錐面である。
【請求項８】
請求項２、請求項３、請求項４、請求項6又は請求項７に記載の装置において、
上記観測信号に対し、上記第２分離行列を用いて分離した信号の周波数領域の信号中の、上記パーミュテーション行列生成手段で、パーミュテーション行列が得られない周波数成分と、パーミュテーション行列が得られた周波数成分との相関を計算する相関計算手段と、
上記相関が大きくなるように上記パーミュテーション行列が得られない周波数の分離行列の行を入れ替えるパーミュテーション行列を生成する修正手段とを備える。
【請求項９】
請求項１記載の装置において、
上記位置情報は上記センサから上記信号源への方向情報であり、
上記位置情報計算手段は、上記比の偏角を計算する手段と、単位距離当たりの位相回転量と上記二つの要素とそれぞれ対応する上記センサ間の距離との積を計算する手段と、上記偏角を上記積の結果で除算する手段と、その除算結果の逆コサインを計算して上記方向情報を出力する手段とを備える。
【請求項１０】
請求項９記載の装置において、
上記位置情報計算手段は、一つの信号源の上記方向情報を、周波数ごとに計算する手段であり、上記各一つの信号源ごとに計算された周波数ごとの方向情報を統合して方向情報を確定する統合手段を備える。
【請求項１１】
Ｉ箇所の信号源から放射された信号を、Ｊ個のセンサで検出して上記信号源の位置情報を求める信号源位置情報推定方法であって、Ｉは２以上の整数、ＪはＩ以上の整数であり、
上記各センサの観測信号を周波数領域の信号に変換する周波数領域変換過程と、
上記周波数領域の信号から、独立成分分析により周波数ごとに、信号源信号を分離する第１分離行列をそれぞれ算出する分離行列算出過程と、
上記各第１分離行列の逆行列あるいは擬似逆行列(以下、両者を単に逆行逆列という)をそれぞれ算出する逆行列算出過程と、
上記周波数ごとの逆行列の少なくとも１つにおける列ごとの異なる二つの要素比から上記信号源の一つの位置情報を計算する位置情報計算過程とを有する。
【請求項１２】
請求項１１記載の方法において、
上記位置情報計算過程は、上記周波数ごとの逆行列の複数について、各列ごとの上記要素比からの位置情報計算により上記各信号源の位置情報を求める過程であり、
上記周波数ごとの信号源ごとの位置情報から、上記複数の逆行列と対応する分離行列における位置情報と対応する行が予め決めた順序になるように行の入れ替えを行うパーミュテーション行列を生成するパーミュテーション行列生成過程と、
上記パーミュテーション行列と上記第１分離行列とを乗算して行を入れ替えた第２分離行列を求める分離行列置換過程とを有する。
【請求項１３】
請求項１２記載の方法において、
Ｊ≧３であり、かつＪ個のセンサは少なくとも２次元に配置され、
上記位置情報は上記センサから上記信号源への方向を含み、
その信号源が存在する円錐面であり、
上記位置情報計算過程は、各列ごとの上記要素比から上記円錐面の計算を異なる二つの要素組の複数について行う過程であり、周波数ごとの上記複数の円錐面の共通直線の方向を、上記位置情報として推定する共通線推定過程を含む。
【請求項１４】
請求項１２記載の方法において、
Ｊ≧３であり、上記位置情報は一対のセンサから上記信号源への方向を含み、その信号源が存在する円錐面及びその信号源が存在する曲面であり、
上記位置情報計算過程は、上記２つの要素の比から上記円錐面を計算する過程と、上記円錐面の計算に用いた２つの要素と対応する一対のセンサの間隔より大きい間隔の一対のセンサと対応する２つの要素の比からこれら一対のセンサと一つの信号源との距離比を計算する過程と、上記距離比から上記曲面を計算する過程を有し、
上記パーミュテーション行列生成過程は周波数ごとに上記円錐面及び上記球面に基づいて上記パーミュテーション行列を生成する過程である。
【請求項１５】
請求項１４記載の方法において、
上記パーミュテーション行列生成過程は、円錐面及び上記球面の一方について上記位置情報計算過程を実行させ、得られた上記一方の位置情報に基づいて、各周波数での逆行列の各列と対応する上記一方の位置情報が予め決めた順になるように列の入れ替えを行う第１の入替行列を生成する過程と、
この過程により、列の入れ替えを行うことができない列について上記位置情報の他方について上記位置情報計算過程を実行させ、得られた上記他方の位置情報に基づいて上記逆行列の列の入れ替えを行うように上記第１の入替行列を修正して第２の入替行列を生成する過程と、上記第２の入替行列の逆行列を計算して上記パーミュテーション行列を生成する過程とを有する。
【請求項１６】
請求項１２乃至請求項１５のいずれかに記載の方法において、
Ｊ≧３であり、上記Ｊ個のセンサは少なくとも２次元に配置され、
上記位置情報は上記センサから上記信号源への方向を含み、その信号源が存在する円錐面であり、
上記円錐面を表わす角度が予め決めた第１角度と第２角度との間であるか否かを判定し
、これら角度の間にない円錐面を破棄する判定過程を有する。
【請求項１７】
請求項１２に記載の方法において、
Ｊ≧３であり、上記Ｊ個のセンサは少なくとも２次元に配置され、
上記位置情報は信号源が存在する球面の半径であり、
上記位置情報計算過程は上記二つの要素比から距離比を計算する過程である。
【請求項１８】
請求項１２に記載の方法において、
上記位置情報は上記センサから上記信号源への方向を含み、その信号源が存在する円錐面である。
【請求項１９】
請求項１２〜１５のいずれか又は請求項１７又は１８に記載の方法において、
上記パーミュテーション行列生成過程でパーミュテーション行列を生成できない周波数があれば、
上記観測信号を上記第２分離行列で分離した信号の周波数領域の信号中の上記パーミュテーション行列を生成できた周波数成分と上記パーミュテーション行列を生成できない周波数成分との相関を計算する過程と、
この計算した相関が大きくなるように、上記生成できなかった周波数の分離行列に対するパーミュテーション行列を生成する過程とを有する。
【請求項２０】
請求項１１記載の方法において、
上記位置情報は上記センサから上記信号源への方向情報であり、
上記位置情報計算過程は、上記比の偏角を、単位距離当たりの位相回転量と上記二つの要素とそれぞれ対応する上記センサ間の距離との積で除算し、その除算結果の逆コサインを計算して上記方向情報を出力する過程である。
【請求項２１】
請求項２０記載の方法において、
上記位置情報計算過程は、一つの信号源の上記方向情報を、周波数ごとに計算する過程であり、上記各一つの信号源ごとに計算された周波数ごとの方向情報を統合して方向情報を確定する統合過程を有する。
【請求項２２】
請求項１１乃至請求項１９のいずれかに記載した信号源位置情報推定方法の各過程をコンピュータに実行させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は音源や電波源などの信号源の複数から放射され、互いに混合された信号を複数のセンサで観測して、各信号源の位置情報、正しくは位置を表わすパラメータの少くとも１つを含む情報を推定し、信号の到来方向の検出や、信号源ごとの信号に分離復元に適用される装置、その方法、及びプログラムに関する。
【背景技術】
【０００２】
複数の信号源からの各信号が空間内で混合されて複数のセンサに到来し、これらセンサで観測された到来信号から、各源信号の到来方向の推定や各源信号を分離することを、そ
の源信号の混合系の情報を知らずに、独立成分分析（Independent Component  Analysis、以後、ＩＣＡと記述する場合もある）により行うことが提案されている。空間内での前記混合は信号源からセンサまでの到達遅延及び減衰度が直接波と障害物などによる複数の反射波で異なるため、ある信号は複数の時間遅れを持って混合された畳み込み混合となる。時間領域で直接分離フィルタを求めるＩＣＡは最終的な解への収束が非常に遅いため、周波数領域で周波数毎にＩＣＡを適用する方法が現実的である。
【０００３】
［到来方向推定］
周波数領域でのＩＣＡを用いて信号源の方向を位置情報として推定する従来技術を図１を参照して簡単に説明する。Ｊ個のセンサ１1，１2，…，１Jが直線状に配列されている。センサ１j(ｊ＝１，２，…，Ｊ)の位置をｄj 、そのセンサ１jでの観測信号をｘj(ｔ)とする。センサ１1 ，…，１j の配列方向と垂直な方向を９０゜とし、源信号ｓｉ(ｔ)の到来方向を０゜≦θi ≦１８０゜とする。Ｉ個の源信号ｓ1(ｔ)，…，ｓI(ｔ)の混合信号をＪ個のセンサ１1 〜１J により、観測信号ｘ1(ｔ)，…，ｘJ(ｔ)として検出しているものとする。
信号の到来方向の推定は周波数領域で行われることが多い。これには、観測信号ｘj (ｔ)を短時間フーリエ変換して周波数領域における時間系列信号Ｘj(ω，ｍ)を求める。ωは角周波数（周波数をｆとするとω＝２πｆ）ｍは時刻を表す番号である。源信号ｓi(ｔ)（ｉ＝１，…，Ｉ）も、同様に、周波数領域における時間系列Ｓi (ω，ｍ)に変換したものとすると、観測信号Ｘj (ω，ｍ)は、信号ｓi の信号源からセンサ１j  への周波数応答Ａji(ω)を用いＸj(ω，ｍ)＝Σi=1IＡji（ω）Ｓi(ω，ｍ)
と表現することができる。ベクトルと行列を用いると、
Ｘ（ω，ｍ）＝Ａ（ω）Ｓ（ω，ｍ）                  ・・・（１）
となる。ここで、
Ｘ（ω，ｍ）＝［Ｘ1(ω，ｍ)，…，ＸJ(ω，ｍ)］T       ・・・（２）
Ｓ（ω，ｍ）＝［Ｓ1(ω，ｍ)，…，ＳI(ω，ｍ)］T       ・・・（３）
はそれぞれ、Ｊ個のセンサによる観測信号およびＩ個の源信号をベクトル表現したものである。Ａ（ω）は周波数応答Ａji（ω）を要素とするＪ×Ｉ行列であり、信号混合系の周波数応答を表現するので混合行列と呼ばれる。［ａ］T はベクトル又は行列ａの転置を表わす。
【０００４】
図１において、方向θi から到来する源信号は、ｄ1 ＝０の位置のセンサ１1 に対し、センサ１j にもτji＝ｃ-1ｄj cosθi だけ速く到達する。ｃは源信号ｓi の速度である。従って、直接波のみを考慮すると角周波数ωにおける周波数応答は
Ａji（ω）＝exp（ｊωｃ-1ｄj cosθi）              ・・・（４）
とモデル化することができる。方向θを変数とする到来方向ベクトルを
ａ（ω，θ）=［exp（ｊωｃ-1ｄ1cosθ）、exp（ｊωｃ-1ｄ2cosθ），…，exp（ｊωｃ-1ｄJ cosθ）］T と表現すると、観測信号は
Ｘ(ω，ｍ)＝Σi＝1Ｉａ（ω，θi）Ｓi（ω，ｍ）
とも近似表現することができる。
【０００５】
独立成分分析を用いて信号源の方向を推定する方法については例えばS.Kurita, H.Saruwatari, S Kajita, K.Takeda,and F.Itakura,“Evaluation of blind signal separation
method using directivity pattern under reverberant conditions,”in Proc.ICASSP2000,2000,pp.3140-3143（文献１という）に示されている。この方法を以下に簡単に述べる。
【０００６】
観測信号Ｘ（ω，ｍ）はＸ（ω，ｍ）＝Ａ（ω）Ｓ（ω，ｍ）であり、源信号Ｓ（ω，ｍ）が混合されたものであるので、互いに独立ではない。Ｘ（ω，ｍ）に独立成分分析
Ｙ（ω，ｍ）＝Ｗ（ω）Ｘ（ω，ｍ）                  ・・・（５）
を適用すると、互に独立な分離信号
Ｙ（ω，ｍ）＝［Ｙ1(ω，ｍ)，…，ＹI(ω，ｍ)］T       ・・・（６）
が得られる。
【０００７】
Ｗ（ω）は要素がＷij（ω）であるＩ×Ｊの行列であり、分離行列と呼ばれる。例えば、Ｉ＝Ｊ＝２の場合、独立成分分析は、Ｙ1(ω，ｍ)とＹ2(ω，ｍ)が互いに独立になる様に
【数１】
0003881367.tif 000002
を満足する分離行列Ｗ（ω）を求める。源信号Ｓ1(ω，ｍ)，…，ＳI(ω，ｍ)が互いに独立であれば、分離信号Ｙ1(ω，ｍ)，…，ＹI(ω，ｍ)はそれぞれ、源信号の何れかに対応する。但し、独立成分分析は信号の独立性のみに基づいているところから、分離信号の順序と大きさに関して任意性がある。即ち、分離行列Ｗ（ω）の行を入れ換えても、Ｗ（ω）の行を定数倍しても、独立成分分析の解である。なお後で述べるが順序の任意性はパーミュテーションの問題となり、大きさの任意性はスケーリングの問題となる。
【０００８】
分離行列Ｗ（ω）のｉ行目は、
ｗi(ω)＝［Ｗi1（ω），…，ＷiJ（ω）］
である。このｗi(ω)は分離信号Ｙi(ω，ｍ)を作り出していることが分かる。従って、ｗi(ω)は、源信号Ｓ1(ω，ｍ)，…，ＳI(ω，ｍ)の何れか１つを強調して取り出し、それ以外を抑圧している。ｗi(ω)が形成する指向特性を解析することにより、何れの方向から到来する信号を取り出し、何れの方向から到来する信号を抑圧しているかを解析することができ、この解析により、源信号ｓi(ｔ)の到来方向を推定することができる。これを全てのｗi(ω)，ｉ＝１，…，Ｉに対して行うと、分離行列のＷ（ω）各行ｗi(ω)が取り出している源信号の到来方向Θ（ω）＝［θ1(ω)，…，θI(ω)］T を推定することができる。
【０００９】
ｗi(ω)が形成する指向特性は、到来方向ベクトルａ（ω，θ）を用いてＢi（ω，θ）＝ｗi(ω)ａ（ω，θ）と表現することができる。このＢi（ω，θ）は方向θにある源信号から、分離信号Ｙi(ω，ｍ)への周波数応答と考えられる。周波数３１５６Ｈｚにおいて独立成分分析後の指向特性のゲイン｜Ｂi(ω，θ)｜特性を図２に示す。図２の横軸はθを、縦軸はゲインを表わす。実線は分離行列の１行目が与える指向特性｜Ｂ1(ω，θ)｜であり、破線は２行目が与える指向特性｜Ｂ2(ω，θ)｜である。実線は５５°でゲインが最小となっており、破線は１２１゜でゲインが最小となっている。このことから、分離行列の１行目は５５゜から到来する信号を抑圧して１２１゜から到来する信号を取り出し、逆に分離行列の２行目は１２１゜から到来する信号を抑圧して５５゜から到来する信号を取り出している。従って、Θ（３１５６Ｈｚ）＝［１２１゜，５５゜］T と推定することができる。
【００１０】
複数の信号源方向を複数センサを用いて、センサの観測信号を周波数領域に変換して推定する方法として、ＭＵＳＩＣ（Multiple Signal Classification）法がある（例えばS.Unnikrishna Pillai,“Array Signal Processing”, Springer -Verlag,1989,ISBN 0-387-96951-9,ISBN 3-540-96951-9参照）。この方法はセンサの数Ｊより１個少ない（Ｊ−１）個までしか信号源の方向を推定できない。しかし前記独立成分分析（ＩＣＡ）を用いる方法（以下単にＩＣＡ法ともいう）によれば２つのセンサで２信号の混合に対応すること
ができ、この点で、ＭＵＳＩＣ法より優れている。しかし、上記のＩＣＡを用いる方法では３信号以上の混合に対処するには以下に示すように困難を伴う。また、指向特性のゲインの最小値を求める操作に多くの計算を要することも欠点である。
【００１１】
３信号の混合に対して３つのセンサを用いてＩＣＡ法を適用した状況を説明する。この場合のＩＣＡは分離行列が３×３になるだけで同様に行えるが、指向特性のゲインを解析するのに困難を伴う。この場合のＩＣＡ後の周波数２７３４Ｈｚにおいて指向特性のゲイン｜Ｂi(ω，θ)｜を図３に示す。図３の実線は分離行列の１行目が、破線は２行目が、１点鎖線は３行目がそれぞれ与える指向特性である。この場合、各源信号は分離行列のある行によって強調され、他の２つの行によって抑圧される筈である。しかし、抑圧する２つの行が同じ方向で抑圧しているとは限らない。例えば、図３においては４５°付近で｜Ｂ2｜と｜Ｂ3｜が共に極小となり、ｗ1(ω)が４５°付近の源信号を取り出し、ｗ2(ω)及びｗ3(ω)はこの源信号を抑圧している。同様に９０°付近で｜Ｂ1｜及び｜Ｂ3｜が極小となり、ｗ2(ω)で９０°付近の源信号を取り出し、ｗ1(ω)及びｗ3(ω)でこの信号を抑圧していると読み取れる。しかしｗ3(ω)が１２０°付近の源信号を取り出し、ｗ1(ω)およびｗ2(ω)がこの源信号を抑圧していると考えられるが、１２０°付近での｜Ｂ1 ｜及び｜Ｂ2 ｜の極小がかなりずれている。この様なずれが大きくなると、何れの方向の抑圧が何れの源信号に対応するか判然としなくなる。従って、ＩＣＡ法による従来技術を３信号以上の状況に対して適用するのは困難であると考えられる。
【００１２】
〔ブラインド信号分離〕
次にＩＣＡを用いたブラインド信号分離の従来技術を説明する。ブラインド信号分離とは、観測された混合信号から源信号を推定する技術である。以下では、Ｉ個の源信号が混在した混合信号がＪ個のセンサで観測される場合を例にとって説明する。
【００１３】
信号源ｉから発生した源信号をｓｉ（ｔ）（ｉ＝１,…,Ｉ、ｔは時刻）、センサｊで観測された混合信号をｘj（ｔ）（ｊ＝１,…,Ｊ）とすると、混合信号ｘj（ｔ）は次式で表わせる。
ｘj(ｔ)＝Σi＝1Ｉ（ａji＊ｓi）(t)                    ・・・（７）
ここで、ａjiは信号源ｉからセンサｊへのインパルス応答、＊は畳み込み演算子である。ブラインド信号分離の目的は、観測信号ｘj（ｔ）のみを用いて、分離のためのフィルタｗkj及び分離信号ｙk（ｔ）（ｋ＝１,…,Ｉ）を次式により求めることである。
ｙk(ｔ)＝Σj＝1J（ｗkj＊ｘｊ）(t)                    ・・・（８）
時間領域での畳み込み混合は、周波数領域における複数の瞬時混合に変換することができる。すなわち、前記の式（７）及び式（８）は、それぞれ前記式（１）及び式（５）で表わされ、前述したように、Ｗ（ω）は分離行列であり、ＩＣＡを用いてＹｋ（ω，ｍ）とＹｋ’（ω，ｍ）が互いに独立となるように計算されたもの、つまりＩＣＡの解である。
【００１４】
以上の周波数領域でのブラインド信号分離を行うに際し問題となるのは、パーミュテーションの問題とスケーリングの問題である。
前述したように分離行列Ｗ（ω）は行の入れ替えを行っても独立成分分析の解となる。つまりある角周波数ωにおいて、Ｗ（ω）がＩＣＡの解だとし、複素数を要素とする任意の対角行列をＤ（ω）とし、任意のパーミュテーション行列（この行列を任意の行列の左から掛けた結果は、この任意の行列の行を入れ替えた行列となる）をＰ（ω）とした場合におけるＰ（ω）Ｄ（ω）Ｗ（ω）もまたＩＣＡの解となる。これは、ＩＣＡが源信号間の統計的独立性のみを条件として源信号の分離を行うことに基づくものであり、このＤ（ω）によって与えられる解の自由度をスケーリングの自由度と呼び、Ｐ（ω）によって与えられる解の自由度をパーミュテーションの自由度と呼ぶ。
【００１５】
そのため、適切なブラインド信号分離を行うためには、全てのωについて、ＩＣＡの解から分離行列として適切な解Ｗ（ω）を特定しなければならない。一般に、この適切な解Ｗ（ω）の特定は、任意に求めたＩＣＡの解に、適切なＤ（ω）やＰ（ω）を乗じ、その結果を適切な解Ｗ（ω）とすることによって行われる。そして、全てのωについて、このＤ（ω）を適切に決める問題をスケーリングの問題と呼び、Ｐ（ω）を適切に決める問題をパーミュテーションの問題と呼ぶ。また、パーミュテーションとは、｛1,2,…,I｝から｛1,2,…,I｝への全単射な関数Ｚ：｛1,2,…,I｝→｛1,2,…,I｝であり、パーミュテーション行列と１対１対応する。
【００１６】
スケーリングの自由度は、時間領域において周波数特性を変化させるフィルタの自由度に相当する。それゆえ、時間領域において歪のない分離信号を生成するためには、全てのωについて、Ｄ（ω）を適切に決める必要がある。このスケーリングの問題は、 例えば、Ｄ（ω）＝ｄｉａｇ（Ｗ-1（ω））とすることで容易に解決することができる。ｄｉａｇ（α）は行列αの対角化（対角要素の他は全ての要素を０とする）を表わす。つまり任意に求めたＩＣＡの解Ｗ0（ω）の逆行列を求め、それを対角化した行列をＤ（ω）とし、Ｄ（ω）Ｗ0（ω）を、適切な分離行列Ｗ（ω）として特定する。このことは既に知られている。例えば参考文献：K.Matsuoka and S.Nakashima,“Minimal Distortion Principle for Blind Source Separation”,Proc.ICA 2001,pp.722-727.に示されている。
【００１７】
一方、パーミュテーションの自由度のため、前記の式（５）の演算結果として、例えば、ある角周波数ω1については分離信号Ｙ1（ω1，ｍ）が源信号S1（ω1，ｍ）の推定値として出力され、別の角周波数ω2については分離信号Ｙ1（ω2，ｍ）が源信号Ｓ2（ω2，ｍ）の推定値として出力されることもありうる。このような場合、時間領域の出力信号ｙ1（ｔ）の中に、時間領域の源信号ｓ1（ｔ）の成分とｓ2（ｔ）の成分とが混在してしまい、分離信号を正しく生成することができない。そのため、時間領域における出力信号ｙ1（ｔ）が正しく源信号ｓ1（ｔ）の推定値となるためには、全てのωについて、Ｙ1（ω，ｍ）がＳ1（ω，ｍ）の推定値となるようにＰ（ω）を適切に決める必要がある。
【００１８】
従来のパーミュテーション問題の代表的な解法として、前記文献１に示す信号の到来方向を推定する方法がある。つまり、図２を参照して説明したように、各周波数における分離行列Ｗ（ω）の各行と対応する指向特性を求める（図２ではｆ＝３１５６Ｈｚの指向特性のみを示している）。これら各周波数の指向特性において例えば分離行列Ｗ（ω）の１行目が与える指向特性のゲイン最小値が５５°で、２行目が与える指向特性のゲイン最小値が１２１°になるようにする。つまりある角周波数ωnではその分離行列Ｗ（ωn）の１行目が与える指向特性のゲイン最小値が１２１°で、２行目が与える指向特性のゲイン最小値が５５°であれば、その分離行列Ｗ（ωn）の１行目と２行目とを入れ替える。つまりこの入れ替えを行うパーミュテーション行列Ｐ（ωn）をＷ（ωn）にその左から掛算する。
【００１９】
このパーミュテーション問題の解法は先にも述べたように指向特性のゲイン最小値を求める操作に多くの計算を要し、しかも信号源の数Ｉが３以上では全角周波数のＷ（ω）について適当に並べ替えを行ってみるという試行錯誤が必要となる。更に図３を参照して説明したようにＷ（ω）のある行がある方向の信号Ｓi（ω，ｍ）を取り出している場合は、そのＷ（ω）の他の行は全てその方向の信号Ｓi（ω，ｍ）を抑圧する状態に必ずしもならない。
【００２０】
また指向特性のゲインが低い方向を探索して推定する信号到来方向の推定精度は信号源の位置によって異なり、特に一対のセンサ１jと１j’を結んだ直線（以下、センサペア軸という）に信号到来方向が近い場合に誤差が大きくなる。このことを実験により示す。図４Ａに示すように、２つのセンサ１０１，１０２としてマイクロホンを距離２．８３ｃｍ
離して配置し、センサ１０１と１０２の中点（原点）から一定の距離（約１５０ｃｍ）離れ、かつ角度間隔が２０°離れた２つの信号源１１１，１１２として音源を設け、センサ１０１からセンサ１０２を見た方向を基準（０°）として、原点からみた信号源１１１の角度θが１０°から１５０°となるまで音源１１１と１１２を前記一定距離及び角度間隔を保持した状態で移動させた。
【００２１】
このように音源１１１，１１２を移動させながら行ったブラインド信号分離結果を図４Ｂに例示する。図４Ｂの縦軸は、信号対妨害音比（Signal to Interference Ratio）を示しており、分離音中に含まれる目的信号及び妨害信号を用いて、
ＳＩＲ＝１０log10（目的信号のパワー／妨害信号のパワー）（ｄＢ）
のように計算したものである。また、図４Ｂの横軸は、原点からみた音源１１１の角度θを表し、実線はこの実験結果を示しており、点線はパーミュテーションに正解を与えたときのＳＩＲを示している。
この図より信号源１１１が、センサペア軸に近い方向（０°或いは１８０°付近）に近づくと、実験結果のＳＩＲはパーミュテーションに正解を与えたときのＳＩＲに比べ大きく低下していることがわかる。このことは信号源１１１の方向がセンサペア軸方向と近いとパーミュテーションが誤っていることに基づくと考えられる。
【非特許文献１】S.Unnikrishna Pillai,“Array Signal Processing”,Springer-Verlag,1989,ISBN 0-387-96951-9,ISBN 3-540-96951-9。
【非特許文献２】大賀寿郎，山崎芳男，金田豊，「音響システムとディジタル処理」，電子情報通信学会，1995,ISBN 4-88552-128-9。
【非特許文献３】A.Hyvarinen and J.Karhunen and E.Oja,“Independent Component Analysis,”John Wiley & Sons,2001,ISBN 0-471-40540。
【非特許文献４】S.Kurita, H.Saruwatari, S Kajita, K.Takeda,and F.Itakura,“Evaluation of blind signal separation method using directivity pattern under reverberant conditions,”in Proc.ICASSP2000,2000,pp.3140-3143。
【非特許文献５】H.Sawada, R.Mukai, S.Araki, S.Makino, “A Robust and Precise Method for Solving the Permutation Problem of Frequency-Domain Blind Source Separation,” in Proc. 4th International Symposium on Independent Component Analysis and Blind Signal Separation (ICA2003), 2003, pp.505-510
【特許文献１】特開２００４−１４５１７２
【発明の開示】
【発明が解決しようとする課題】
【００２２】
以上述べたように独立成分分析を用いて分離行列を求め、その各行から指向特性パターンを求めて、その利得が低い方向探索により、信号源の方向（到来信号方向）を求め、更にこれを利用してブラインド信号分離を行う場合は、前記指向特性パターンを求めて、利得が低い方向を探索するために多くの計算時間を必要とした。
この発明の目的は独立成分分析により分離行列を求めて信号源の位置情報を推定するための計算時間が短かい位置情報推定装置、その方法及びプログラムを提供することにある。
【課題を解決するための手段】
【００２３】
この発明によれば周波数領域の分離行列Ｗ（ω1），…，Ｗ（ωN）の逆行列（Ｉ＜Ｊの場合は擬似逆行列）を計算してスケーリングとパーミュテーションの自由度を含む混合行列Ａ（ω1），…，Ａ（ωN）の推定値Ｈ（ω1），…，Ｈ（ωN）を生成し、これらの周波数ごとのＨ（ωn）（ｎ＝１，…，Ｎ）について各列ごとにその２つの要素Ｈji（ωn）とＨj'i（ωn）（ｊ、ｊ’はセンサを表すパラメータ、ｉは信号源を表すパラメータ）の比から、信号源ｉの位置情報の一つのパラメータ、例えば信号源が存在する円錐面あるいは曲面を計算する。
【発明の効果】
【００２４】
このように要素比に基づき表された計算式を演算すればよく、分離された信号の指向性パターンを求め、更にその極小値角度の探索を行う場合より少ない計算量で済む。また要素比をとるため、前述したスケーリングの自由度による影響がなくなる。
【発明を実施するための最良の形態】
【００２５】
まず信号源の位置情報として方向情報を推定する場合にこの発明を適用した実施形態を説明する。なお以下の説明において同一又は対応するものについて図面中に同一の参照番号を付けて、重複説明を省略する。
【００２６】
［第１実施形態］
この第１実施形態は信号源の方向、つまりその信号源から放射された源信号の到来方向を求める。
図５に第１実施形態の機能構成図を、図６にその処理手順の一部の流れ図をそれぞれ示す。
【００２７】
信号源の数Ｉ以上の個数であるＪ個のセンサ１1，１2，…，１Jが、例えば図１に示したように配列されている。隣接センサの間隔は通常は、源信号の最も短かい波長の１／２以下である。これらセンサ１j（ｊ＝１，２，…，Ｊ）で観測された信号ｘj（ｔ）はそれぞれ周波数領域変換部１１jにおいて、例えば短時間フーリエ変換により周波数領域信号Ｘj（ω，ｍ）に変換される（図６、ステップＳ１）。これら周波数領域信号Ｘj（ω，ｍ）に対する各角周波数ωnごとの分離行列Ｗ（ωn）（ｎ＝１，２，…，Ｎ）が分離行列算出部１２において独立成分分析により算出される（図６、ステップＳ２）。
Ｗ（ω）＝（Ｗ（ω1），Ｗ（ω2），…，Ｗ（ωN））
この各周波数ごとの分離行列Ｗ（ωn）の逆行列が逆行列算出部１３で算出され逆行列Ｈ（ωn）が求まる（図６、ステップＳ３）。
Ｈ（ω）＝（Ｈ（ω1），Ｈ（ω2），…，Ｈ（ωN））
なおこの逆行列の計算はＪ＞Ｉの場合は擬似逆行列を計算することになる。擬似逆行列としては例えばＭｏｏｒｅ−Ｐｅｎｒｏｓｅ型一般逆行列が用いられる。
【００２８】
この実施形態においてはＨ（ω）中の少くとも１つの周波数の逆行列Ｈ（ωn）の各列ｉの２つの要素Ｈji（ωn）とＨj’i（ωn）との比の偏角から方向情報、つまり源信号の到来方向が角度計算部１４で計算される（図６、ステップＳ４）。この角度計算部１４の具体的機能構成例と処理手順例を図７及び図８を参照して説明する。選択部１４ａにより角周波数ωnの逆行列Ｈ（ωn）中の未選択の一つの列ｉが選択され（図８、ステップＳ４ａ）、そのｉ列目から２つの要素Ｈji（ωn）とＨj’i（ωn）が選択される（図８、ステップＳ４ｂ）。
【００２９】
偏角計算部１４ｂで、選択された要素の比の偏角arg［Ｈji（ωn）／Ｈj’i（ωn）］が計算され、また間隔計算部１４ｃで、センサ情報格納部１５内のセンサ１jと１j’の位置情報ｄjとｄj’が取り出され、これらの差からセンサ１jと１j’の間隔ｄj−ｄj’が計算される（図８、ステップＳ４ｃ）。位相回転計算部１４ｄで、間隔ｄj−ｄj’と、角周波数ωnでの単位距離当りの信号の位相回転量ωn／ｃ（ｃは源信号の速度）との積が計算されて間隔ｄj−ｄj’での信号の位相回転量が求められ、この位相回転量により偏角arg［Ｈji（ωn）／Ｈj’i（ωn）］が除算部１４ｅで割算される（図８、ステップＳ４ｄ）。なお単位距離当たりの位相回転量ω／ｃは例えばｆ＝６８０Ｈｚの音波（速度ｃ＝３４０ｍ／ｓ）では、１ｍ当たり２π・２の位相が回転することになる。
【００３０】
この除算部１４ｅの割算結果Ｇi（ωn）の絶対値が１以下か否かかが判定部１４ｆで判定され（図８、ステップＳ４ｅ）、１以下であればＧi（ωn）の逆余弦θi（ωn）＝cos-1Ｇi（ωn）が逆余弦計算部１４ｇで計算される（図８、ステップＳ４ｆ）。つまり角度計算部では次の計算が行われる。
【００３１】
Ｇi(ωn)＝arg［Ｈji（ωn）／Ｈj’i（ωn）］／（ω（ｄj−ｄj’）／ｃ））
θi(ωn)＝cos-1Ｇi（ω）  ：  ｜Ｇi(ωn) ｜≦１に対して、  …（９）
ステップＳ４ｅで｜Ｇi（ωn）｜が１以下でなければ角度θi（ωn）が虚数になるので、別の組み合せを選択するため、判定部１４ｆで選択したｉ列目のすべての要素の組み合せを選択したか判定され（図８、ステップＳ４ｇ）、選択していない組み合せがあればステップＳ４ｂに戻り、全ての組み合せを選択していれば、判定部１４ｆですべての列を選択したか判定され（図８、ステップＳ４ｈ）、選択していない列があればステップＳ４ａに戻り、全ての列を選択していれば角度計算処理を終了する。なお図７中の間隔計算部１４ｃは各周波数の角度計算部１４に対して共通に用いられる。
【００３２】
角度計算部１４からＨ（ω）中の選択した角周波数ωnの逆行列Ｈ（ωn）における各列と対応して、その選択順に、もし第１列目から順次選択すればその順に式（７）の計算結果、つまりＩ個の信号源の方向（信号到来方向）Θ（ωn）＝（θ1（ωn），θ2（ωn），…，θＩ（ωn））が出力される。これらθ1（ωn），θ2（ωn），…，θＩ（ωn）は各源信号ｓ1（ｔ），ｓ2（ｔ），…，ｓＩ（ｔ）の到来方向（信号源１，２，…，Ｉの方向）のいずれかの１つと対応する。
【００３３】
この実施形態において信号到来方向が推定できるしくみを以下に説明する。独立成分分析（ＩＣＡ）により分離が達成されていれば、ＩＣＡにより算出した分離行列Ｗ（ω）と真の混合行列Ａ（ω）とはＰ（ω）Ｄ（ω）Ｗ（ω）Ａ（ω）＝Ｉの関係にある。ここで、Ｄ（ω）はスケーリングの自由度を示す対角行列、Ｐ（ω）はパーミュテーションの自由度を示すパーミュテーション行列、Ｉは単位行列である。ＩＣＡを用いても、混合行列Ａ（ω）そのものは一般には算出できない。しかし、Ｗ（ω）の逆行列Ｈ（ω）＝Ｗ-1（ω）＝Ａ（ω）Ｐ（ω）Ｄ（ω）を算出すれば、スケーリングとパーミュテーションの自由度を含む混合行列の推定値が得られる。すなわち、逆行列Ｈ（ω）は、混合行列Ａ（ω）の列をＰ（ω）により並べ替え、さらに各列にそれぞれＤ（ω）の対角要素を掛け合わせたものになる。
【００３４】
この実施形態では、逆行列Ｈ（ω）の同じ列ｉから２つの要素Ｈji（ω）とＨj’i（ω）を取り出してそれらの比Ｈji（ω）／Ｈj’i（ω）を求めることで、算出できないＤ（ω）によるスケーリングの自由度を取り除く。すなわち、
【数２】
0003881367.tif 000003
となる。ここで、Ｚはパーミュテーション行列Ｐ（ω）を右から掛けることに対応する順列である。逆行列Ｈ（ω）の全ての列ｉに対して式（２１）に関わる操作を行うことで、Ｐ（ω）によるパーミュテーションＺにかかわらず、全ての信号の到来方向が推定できる。
【００３５】
背景技術の説明では、混合行列Ａ（ω）の要素をＡji（ω）＝exp（jωｃ-1ｄjcosθi）とモデル化した。しかし、この実施形態では、分離行列Ｗ（ω）の逆行列Ｈ（ω）を用いて、スケーリングとパーミュテーションの自由度を含む混合行列Ａ（ω）の推定値を算出しているので、この単純なモデルでは不十分である。そこで、振幅の減衰度αji（実数）と原点における位相差exp（jφi）を用いて、Ａji（ω）＝αjiexp（ｊφi）exp（ｊωｃ-1ｄjcosθi）とモデル化する。このモデルを用いてＡjZ(i)（ω）／Ａj’Z(i)（ω）を計算すると、式（２１）の関係から、
【数３】
0003881367.tif 000004
が得られる。その結果、上記のＧi（ω）は
Ｇi（ω）＝arg［Ｈji（ω）／Ｈj’i（ω）］／（ωｃ-1（ｄj−ｄj’））＝cosθZ(i)
となり、｜Ｇi（ω）｜＜１であればθZ(i)＝cos-1Ｇi（ω）が実数となり、到来方向を推定することができる。先に述べたように、パーミュテーションＺの自由度により、Ｉ個の方向Θ（ω）＝［θZ(i)（ω），…，θZ(I)（ω）］全体を適当に並べ替えたものが、信号ｓ1，…，ｓIの方向に対応する。
【００３６】
Ｈ（ω）中の複数の周波数又は全ての周波数の各逆行列Ｈ（ωn）（ｎ＝１，…，Ｎ）について先に述べたように各列ごとに角度θi（ωn）を求め、これらの全体により各到来方向を決定してもよい。つまり角度計算部１４で推定された各周波数の到来方向は図５中のソート部３２においてその角度ごとに分けられる（図６、ステップＳ５）。例えばそれぞれ方向（角度）の大きい順に並べられる。角周波数ω1についてのΘ（ω1）中の各成分が大きい順に（θ1（ω1），θ2（ω1），…，θＩ（ω1））と並べ、ω2についてのΘ（ω2）中の各成分が大きい順に（θ1（ω2），θ2（ω2），…，θＩ（ω2））と並べ、…，ωNについてのΘ（ωN）中の各成分が大きい順に（θ1（ωN），θ2（ωN），…，θＩ（ωN））と並べられる。このように大きい角度順に分別された角度は、大きい順が同一であっても周波数間でばらつきがある、つまり角度θi（ω1），θi（ω2），…，θi（ωN）（ｉ＝１，…，Ｉ）はばらついている。
【００３７】
よって統合部３３で、分別角度θi（ω1），θi（ω2），…，θi（ωN）ごとに１つの角度θｉに統合されて１つの角度（到来方向）θiとされる（図６、ステップＳ６）。この統合の方法としては例えば分別角度（θi（ω1），θi（ω2），…，θi（ωN））ごと平均値が統合角度θiとされ、あるいは分別角度（θi（ω1），θi（ω2），…，θi（ωN））ごとの最頻値や中央値などを統合角度θiとしてもよい。このように分別角度ごとに１つの角度に統合することにより、１つのある周波数の逆行列Ｈ（ωn）からのみ到来方向を推定する場合より、正しく到来方向を推定することができる。
【００３８】
なお図８に示した角度（方位）計算推定処理において、その１つでも選択した列について、角度θi（ωn）を求めることができなかった時は、そこでその逆行列Ｈ（ωn）に対する処理を終了し、他の周波数の逆行列Ｈ（ωn）に対する処理を行い、最初に全ての列について角度（方向）を計算することができたら、その時の計算結果を推定（方向）θ1，…，θIとしてもよい。あるいは各周波数の逆行列Ｈ（ωn）中で全ての列について角度計算ができた結果を角度ごとに分別し、統合してもよい。あるいは１つの列について角度θi（ωn）を計算することができないものがあれば、その状態が最初に発生した時に、その後の処理を全て終了し、あらためて観測信号を求める処理からやりなおすことにより、推定結果の信頼性を高めるようにしてもよい。
【００３９】
この第１実施形態の実験例を述べる。残響時間が１９０ｍｓの部屋に３つのマイクロホンを５６．６ｍｍ間隔で１列に配列し、その配列方向を基準として角度４８°、７３°、１１９°の方向に３つの音源を配置した。これら音源からの音響信号の混合時間を６秒、観測信号の標本化周波数を８ｋＨｚ、空間エイリアシングが生じない最大周波数を３ｋＨｚ、短時間フーリエ変換のフレームを１０２４サンプルとした。前述したようにして、各周波数ごとにそれぞれ計算した各角度を図９に、横軸を周波数、縦軸を方向として示す。図９中の◇，＋，□はそれぞれ３つの音源方向の推定計算値を示す。この結果を３つの角度範囲に分別し、各分別された角度平均値は４５°，７４°，１２３°となった。３音源の方向を３つのマイクロホンで推定することはＭＵＳＩＣ法ではできないが、この実施形態によれば可成り正しく各音源方向が推定されていることがわかる。
【００４０】
この実施形態の方法とＭＵＳＩＣ法とを比較するため、音源方向を４８°と１１９°と比較的大きく角差がある状態で２つの音源を配置した場合について同様に実験を行った。図１０Ａにこの実施形態の方法により得られた結果を、図１０ＢにＭＵＳＩＣ法により得られた結果をそれぞれ示す。これらよりいずれの方法も、かなり正しく方向を推定していることがわかる。この実施形態による結果を２つの方向範囲で分別平均した結果は４５°，１２３°であり、ＭＵＳＩＣ法はそれぞれ４５°，１２２°である。次に、音源方向を１０５°と１１９°と方向が比較的接近している状態で２つの音源を配置した場合について同様に実験を行った。この実施形態方法による結果及びＭＵＳＩＣ法による結果を図１１Ａ及び図１１Ｂにそれぞれ示す。ＭＵＳＩＣ法では大部分の周波数で音源方向の推定が
できないが、この実施形態では大部分の周波数で角度計算を行うことができ、しかも、これらを角度範囲で分別して、それぞれ平均した値は１０５°，１２４°であり、ＭＵＳＩＣ法は９４°，１２８°でありこの実施形態による方向推定がかなり正しいことがわかる。
【００４１】
以上述べたようにこの実施形態によれば従来方法において、指向性パターンの利得が低い方向探索する場合と比較して、式（９）に値を代入するだけで推定方向が求まるため計算時間がかなり短かい。なお前述したようにＩＣＡにより求めた分離行列Ｗ（w）はスケーリングの自由度及びパーミュテーションの自由度を含んでいるためこれらを解決した分離行列Ｗ’（ω）の逆行列を計算し、つまり真の混合行列Ａ（w）を計算し、その混合行列Ａ（w）の各列について二つの要素の比を用いて到来方向を推定することが考えられるかもしれない。しかし、真の混合行列Ａ（w）そのものを求めることは、例えば信号源の信号ｓi(ｔ)の平均パワーを１にするなどの制約を与えない限りできない。無線通信分野においては信号源に対しそのような制約を与えることができるかも知れないが、例えば信号源の信号ｓi(ｔ)が人間が直接発声した音声信号である場合は、そのような制約は実質的にはできない。一方この第１実施形態によれば、スケーリング及びパーミュテーションの自由度を含んでいる分離行列Ｗ（w）の逆行列Ｈ（w）の各列について二つの要素の比をとることによりスケーリングの自由度の問題を解決しておりどのような信号源に対しても適用でき、しかも前記両問題を解決した分離行列を計算する必要がなくそれだけ計算時間が短い。更に前述したように、各周波数について求めた推定方向を予め決めた順に分別すればパーミュテーションの問題も簡単に解決できる。センサの数Ｊと同数であっても各信号源の方向を推定することができる。また音源方向が比較的接近していてもかなり正しく推定することができる。
【００４２】
［第２実施形態］
この第２実施形態は信号源の一つの位置情報として方向情報を求めるものであり、かつ少なくとも２次元に配置された少なくとも３個のセンサを用いて、信号源がいずれの方向に位置していても方向を推定でき、従ってブラインド信号分離におけるパーミュテーションの問題も比較的簡単に解決するものである。つまり方向情報に基づく円錐面を推定し、複数の円錐面の共通直線を推定して、方向情報を決定する。
この第２実施形態をブラインド信号分離装置に適用した機能構成を図１２に、その処理手順を図１３に示す。例えば４個のセンサ１1，１2，１3，１4が同一円上で等間隔に配置され、そのいずれの２つのセンサの間隔も、源信号の最短波長の１／２以下とされる。以下の説明はセンサの数はＪ個であり、Ｊ≧３として行う。第１実施形態と同様に各センサｊ（ｊ=１,…,Ｊ）で観測された観測信号ｘj（ｔ）は、それぞれ周波数領域変換部１１で例えば、短時間フーリエ変換によっ
て、周波数領域の信号Ｘj（ω,ｍ）にそれぞれ変換される（ステップＳ１１）。
【００４３】
これら周波数領域の信号Ｘj（ω,ｍ）から、分離行列算出部１２で、独立成分分析により、各周波数ごとの分離行列
【数４】
0003881367.tif 000005
が算出される（ステップＳ１２）。
【００４４】
逆行列算出部１３でこれら周波数ごとの各分離行列Ｗ（ω）の逆行列Ｈ（ω）がそれぞれ算出される（ステップＳ１３）。
【数５】
0003881367.tif 000006
【００４５】
この第２実施形態では、円錐面推定部１４において、周波数ごとの各逆行列Ｈ（ω）における各列ごとに異なる複数の要素対について各要素の比からこれら要素と対応する二つのセンサのセンサペア軸を中心軸とし、何れかの信号源が存在する円錐面、つまり一つの混合行列Ｈ（ω）の各列ごとに複数の円錐面がそれぞれ推定される（ステップＳ１４）。
円錐面推定部１４の機能構成は図１５中の角度計算部１４とほぼ同様であり、その処理手順も図８の手順と似ている。円錐面推定部１４において行われるある一つの周波数の逆行列Ｈ（ω）に対するステップＳ１４の処理の具体例を図１４を参照して説明する。
【００４６】
まず、例えば、円錐面推定部１４内のレジスタ内に格納されている制御パラメータｉとｐを０に初期化する（ステップＳ２０）。ここで、ｉは各信号源の番号と対応し、ｐはｉごとの推定済み円錐面の個数である。
【００４７】
ｉを＋１し（ステップＳ２１）、ｇを＋１し（ステップＳ２２）、制御パラメータｊ,ｊ’はそれぞれ互いに異なるｊ以下の自然数を、例えばランダムに選択する（ステップＳ２３）。この制御パラメータｊ,ｊ’の組は、同一のｉに対しては一度選択されたｊ,ｊ’の組み合わせの選択は行わないようにする。すなわち、ｊ≠ｊ’であり、また、例えば、ｉ＝１について、一度（ｊ,ｊ’）＝（１,２）が選択されると、ｉ＝１については、この処理が終了するまで、再びステップＳ２３において（ｊ,ｊ’）＝（１,２）という選択は行われないものとする。（また、この選択は、選択するｊ,ｊ’によって特定されるセンサペア軸（センサｊとセンサｊ’とを通る直線）が、このルーチン処理においてそれ以前に選択されたｊ,ｊ’によって特定されるセンサペア軸と同一直線上に配置されないように行われることが望ましい。これにより、円錐面推定部１４は、一定の誤差範囲内で中心軸が重ならない複数の円錐面を推定することになる。なお、各センサペア軸が同一直線上となるか否かの判断は、例えば、センサ情報格納部１５に各センサの位置を示すベクトルを格納しておき、ここから各センサの位置を示すベクトルの情報を抽出して行う。）
次に、センサ情報格納部１５から、ステップＳ２３で選択したパラメータｊに対応するｊ番目のセンサｊの位置を示すベクトルｄjと、変数ｊ’に対応するｊ’番目のセンサｊ’の位置を示すベクトルｄj’とを抽出する（ステップＳ２４）。また、混合行列Ｈ（ω）から、ｊ行ｉ列要素Ｈji(ω)と、ｊ’行ｉ列要素Ｈj’i(ω)とを選択して抽出する（ステップＳ２５）。以上の処理は図７中の選択部１４ａにより行われる。従って、ｉ，ｐ，ｊ，ｊ’や選択したｊ，ｊ’などを格納するレジスタも選択部１４ａに設けられることになる。
【００４８】
抽出したこれらの情報を用い、次式（９’）を演算する（ステップ２６）。
【数６】
0003881367.tif 000007
ここで‖ｄj−ｄj’‖はセンサ１jと１j’との間隔(距離)である。この第２実施形態においては複数のセンサが２次元又は３次元に配置されているため各センサの位置情報は例えば図１２中のセンサ１1〜１4が配置されている円の中心を原点とする２又は３要素の座標ベクトルで表わされる。つまり式（９）はセンサが１次元に配置され、信号到来方向の角度が２次元の場合であるが、この式（９’）はセンサが２次元あるいは３次元に配置されていてもよく信号到来方向の角度が３次元空間でもよい場合に拡張されたものである。従って式（９’）は式（９）を包含する。この式（９’）により推定された角度θ＾i,jj’（ω）及びそのｉ，ｊ，ｊ’を円錐面情報として円錐面推定部１４内のレジスタ（メモリ）に一時記録する（ステップＳ２７）。なお図１２中に破線で示すように各周波数に対する間隔計算部１４ｃは共通に用いられる。また式（９）で演算した角度θi（ω）は三次元空間においてセンサペア軸（センサ１jと１j’とを結ぶ直線）に対する角度がθi（ω）となる直線の無数の集合、つまり円錐面に信号源ｉが存在することを推定していることになる。この拡張された式（９’）の演算結果を単なるθi（ω）ではなくθ＾i，jj’（ω）と表わした。ステップＳ２６で行う演算は図７中の偏角計算部１４ｂ、間隔計算部１４ｃ、位相回転計算部１４ｄ、除算部１４ｅ、判定部１４ｆ及び逆余弦計算部１４ｇにより、図８中のステップＳ４ｃ，Ｓ４ｄ，Ｓ４ｅ及びＳ４ｆの処理により行う。
【００４９】
次に、ｐ＝Ｐであるか否かを判断する（ステップＳ２８）。このＰは、ｉごとに推定すべき円錐面の個数であり、このステップでは、そのｉについて円錐面をＰ個推定したか否かを判断する。ｐ＝ＰでなければステップＳ２２に戻り、ｐ＝Ｐであれば、ｉ＝Ｉであるか否かを判断する（ステップＳ２９）。すなわち、すべてのｉについて円錐面の推定が終了したか否かを判断し、ｉ＝ＩでなければステップＳ２１に戻り、ｉ＝Ｉであれば処理を終了する（ステップＳ１４の具体例の説明終了。）。第１実施形態では選択したｉについて角度θi（ω）を１つ推定すれば次のｉについての角度推定に移るがこの第２実施形態では各ｉごとに複数Ｐ個の組の角度（円錐面）θ＾i，jj’（ω）を推定する。ステップＳ２８及びＳ２９の処理は図７中の判定部１４ｆで行われる。
【００５０】
図１２中の到来方向決定部１６で、円錐面推定部１４で推定された複数の円錐面の情報（この例では、ｉ,ｊ,ｊ’、θ＾i，jj’（ω））から源信号の到来方向ｕi（ω）＝（方位角θi（ω），仰角φi（ω））（ｉ=１,…,Ｉ）が決定される（図１２、ステップＳ１５）。具体的には、例えば、ｉについて推定された複数の円錐面のうち、互いに線接触しているとみなされる共通直線の方向を、信号源ｉに対応する信号の到来方向ｕi（ω）として算出する。
【００５１】
この信号の到来方向ｕi（ω）の推定方法を図１５を参照して説明する。
センサ１1と１2の配列方向と直角方向にセンサ１3が配され、センサ１1と１2間、センサ１2と１3間の各間隔は同一とされている。ｉ＝２の信号源２2に対し、センサ１1及び１2によりセンサペア軸３12を中心軸とする円錐面４12（θ＾2，12（ω））が、センサ１2と１3によりセンサペア軸３23を中心軸とする円錐面４23（θ＾2，23（ω））が、センサ１1と１3によりセンサペア軸３13を中心軸とする円錐面４13（θ＾2，13（ω））がそれぞれ推定されている。これらのうち、円錐面４12と４13は、共通直線５2で互いに線接触しているとみなされる。この共通直線５2の方向ｕ2（ω）が、信号源２2から放射された信号の到来方向ｕ2つまり信号源２2の方向と決定する。なお円錐面４23も平行移動すれば円錐面４12，４13とほぼ線接触するが、後述する第４実施形態により円錐面４23は廃棄される。
【００５２】
図１２中の到来方向決定部１６におけるこの共通直線５iの方向θi（ω）を決定するための具体的計算方法の例を説明する。角周波数ωに対して信号源２iが存在すると推定された複数の円錐面を４jj’（１），…，４jj’（Ｐ）とし、これら円錐面４jj’（ｐ）（ｐ＝１，…，Ｐ）の推定に用いられた１対のセンサの位置情報をｄｊ（ｐ），ｄj’（ｐ
）とし、角周波数ωについて推定された円錐面４jj’（ｐ）と対応する角度をθ＾jj’（ω，ｐ）とし、円錐面４jj’（ｐ）を表わすベクトルをｕとする。
【００５３】
正規化軸ベクトル計算部１６ａで前記各一対のセンサの位置間を結ぶ軸ベクトル（ｄj（ｐ）−ｄj’（ｐ））をそれぞれ長さｌに正規化する。つまり
ｖp＝（ｄj（ｐ）−ｄj’（ｐ））／‖ｄj（ｐ）−ｄj’（ｐ）‖
をそれぞれ演算する。このｖpと円錐面ベクトルの内積はこれらベクトルのなす角度の余弦とする。つまり
ｖpT・ｕ／‖ｕ‖＝cosθ＾jj’（ω，ｐ）
関係が成立する。知りたいのは共通直線５ｉの方向のみであるから円錐面ベクトルｕは単位ベクトル、つまり‖ｕ‖＝１とする。全ての円錐面と共通する直線の方向を求めるにはＶ＝(ｖ1…ｖP)T,ｃ＾(ω)＝(cosθ＾jj’(ω，１)…cosθ＾jj’(ω，Ｐ))T
として次式の連立方程式をｕについて解けば良い。
【００５４】
Ｖｕ＝ｃ＾(ω)  （‖ｕ‖＝１）
一般にはこの連立方程式は解が存在しないか、もしくは一意には決まらない。そこで‖Ｖｕ−ｃ＾(ω) ‖を最小化するｕを求めて前記連立方程式の解、つまり共通直線５iの方向ｕi（ω）とする。この誤差を最小化するｕを求める計算が演算部１６ｂで行われる。この方向ｕi（ω）は３次元での方向であるから、極座標表示により方位角θi（ω）及び仰角φi（ω）が求まることになる。
【００５５】
もっと計算を簡便にするには次のようにしてもよい。図１６に示すように正規化軸ベクトル計算部１６ａで円錐面推定に用いた各センサ対についての正規化軸ベクトルｖp（ｐ＝１，…，Ｐ）を求め、逆行列計算部１６ｃでＶ＝（ｖ1，…，ｖP）TのＭｏｏｒｅ−Ｐｅｎｒｏｓｅ型一般逆行列Ｖ＋を計算し、このＶ＋と推定円錐面と対応する余弦ベクトルｃ＾(ω)＝（cosθ＾jj’(ω，１)…cosθ＾jj’(ω，Ｐ)）Tとを用いて最小ノルムな最小２乗誤差解を求め、大きさを正規化して近似解としてもよい。つまり演算部１６ｄで
ｕi（ω）＝Ｖ＋ｃ＾(ω)／‖Ｖ＋ｃ＾(ω)‖
を計算する。
このようにして複数の推定円錐面の共通直線とみなされる直線の方向が各周波数ごとに、かつ各信号源ごとに求められる。
【００５６】
図１２中のパーミュテーション解決部１７は、到来方向決定部１６で決定された到来方向ｕｉ＝（θi，φi）を用いて、分離行列算出部１２で算出された分離行列Ｗ（ω）の行の並び替えを行ってパーミュテーション問題が解決された分離行列を生成する。
【００５７】
パーミュテーション解決部１７で行われる具体例としては、到来方位角θiについて以下のように並べ替えを行い、解決できなかった列について同様に到来仰角φiによる並べ替えを行う。つまり算出決定された到来方位角θi（ω）がどの周波数においても既定の順序、例えば、θ1,θ2,…,θＩが昇順となるように逆行列Ｈ（ω）の列を入れ替え、入れ替えができなかった列について仰角φiがどの周波数でも昇順になるように逆行列Ｈ（ω）の列を入れ替える並替行列が並替行列生成部１７ａで生成され、その並替行列の逆行列Ｐ（ω）が逆行列生成部１７ｂで生成され、並替部１７ｃでその逆行列Ｐ（ω）が、分離行列Ｗ（ω）に左から掛算される。並替行列生成部及び逆行列生成部１７ｂはパーミュテーション行列生成部を構成している。
【００５８】
並替行列生成部１７ａでの処理を具体的に述べる。この例では逆行列Ｈ（ω）の第１列目で（θ1（ω），φ1（ω））が、第２列目で（θ2（ω），φ2（ω））が，…、第Ｉ列目で（θＩ（ω），φＩ（ω））がそれぞれ式（９’）に基づき計算されたとする。いま到来方向決定部１６１から入力された到来方向を、前記昇順に並べたものと区別するため
にθおよびφの各添字にダッシュ「’」を付けて、（θ1’（ω），φ1’（ω）），（θ2’（ω），φ2’（ω）），…，（θＩ’（ω），φＩ’（ω））と表わし、これらを例えばまずθi’（ω）について昇順に並べた時に、例えば（θ3’（ω），φ3’（ω））＞（θＩ’（ω），φＩ’（ω））＞（θ2’（ω），φ2’（ω））＞…であれば、逆行列Ｈ（ω）の３列目を１列目に、Ｉ列目を２列目に、２列目を３列目になるように移動し、残りの列も同様に移動し、θi’（ω）が同一のものについてはφi’（ω）が昇順になるように列を移動する。そのような列の移動並べ替えを行う並替行列を作る。このような並べ替えを行う行列を作ることは従来からよく行われていることである。このようにして全ての周波数について得られた到来方向（θ１（ω），φＩ（ω）），…，（θＩ（ω），φ2（ω））を用いて、並替行列が作られ、更にその逆行列つまりパーミュテーション行列Ｐ（ω）が算出される（図１３、ステップＳ１６）。
【００５９】
この算出したパーミュテーション行列Ｐ（ω）が並替部（１７ｃ）で分離行列Ｗ（ω）の左から掛けられ、その結果Ｗ’（ω）＝Ｐ（ω）Ｗ（ω）がパーミュテーション問題が解決された分離行列として出力される（ステップＳ１７）。つまりこの分離行列Ｗ’（ω）はいずれの周波数のものでも、１行目は信号源２1の信号を分離する要素、２行目は信号源２2の信号を分離する要素、以下同様に同一行の要素は同一信号源からの信号を分離する要素となる。
【００６０】
分離行列Ｗ’（ω）は、時間領域変換部１８で、例えば、逆フーリエ変換によって、時間領域の分離フィルタ係数群
【数７】
0003881367.tif 000008
に変換され、これら分離フィルタ係数群は信号分離部１９に設定される。
【００６１】
信号分離部１９は、各センサからの観測信号ｘ1（ｔ）,…,ｘＪ（ｔ）と分離フィルタ係数群とにより式（８）の演算を行い分離信号ｙ1（ｔ）,…,ｙＩ（ｔ）を出力する。
図１２中に破線で示すように、パーミュテーション解決部１７よりの並替えされた分離行列Ｗ’（ω）と周波数領域変換部１１よりの周波数領域観測信号Ｘ（ω，ｍ）とにより式（５）の演算を周波数領域分離信号生成部１９’で行い、その結果の周波数領域分離信号Ｙ（ω，ｍ）＝Ｗ’（ω）Ｘ（ω，ｍ）を時間領域信号変換部１８’で時間領域信号ｙ1（ｔ），…，ｙＩ（ｔ）を生成してもよい。また到来方向決定部１６で決定された各周波数ごとの到来方向（θ1（ω），φ1（ω）），…，（θＩ（ω），φＩ（ω））を第１実施形態と同様に、つまり図５中のソート部３２で方向範囲ごとに分別し、各分別されたものごとの全周波数における到来方向を統合部３３でそれぞれ統合してもよい。
以上説明したように、この実施形態においても指向特性パターンの低ゲイン方向を探索することなく、式（９’）の演算により円錐面情報（θ^i，jj’（ω），φ^i，jj’（ω））を推定しているため計算量が少ない。しかも同一信号源に対し、円錐面を複数推定してその共通直線から信号の到来方向を算出しているため、０°から３６０°までのどの方向に信号源が存在しても、信号源３iの方向を一意に推定することができる。またその推定方向をパーミュテーション行列Ｐ（ω）の決定に利用しているため、信号源の位置に関わらず、パーミュテーション問題を適切に解決することができる。  
【００６２】
〔第３実施形態〕
第３実施形態は、位置情報として一対のセンサと一つの信号源との各距離の比に基づく
その信号源が存在する曲面を用いる。第１及び第２実施形態においては各信号源はセンサから遠方にあって、信号源からの信号はセンサには平面波として到来することを想定して処理した。しかし信号源とセンサとの距離が比較的短かい場合はセンサには信号が球面波として到来する。このような点から混合行列Ａ（ω）の要素の比Ａji(ω)／Ａj’i(ω)を球面波（近距離場）モデルによって解釈すると、信号源ｉの方向以外の情報を推定することができる。
【００６３】
すなわち、近距離場モデルを用いると、周波数応答Ａji(ω)は以下のように表わせる。
Ａji（ω）＝（１／‖ｑi−ｄj‖）exp（ｊωｃ-1（‖ｑi−ｄj‖））
ここで、ｑiは信号源ｉの位置を示すベクトルである。
このように表現された周波数応答について混合行列の同一列の２つの要素比Ａji（ω）／Ａj’i（ω）をとり、その絶対値を計算すると次式となる。
‖ｑi−ｄj’‖／‖ｑi−ｄj‖＝｜Ａji（ω）／Ａj’i（ω）｜  …（１０）
なお、｜β｜はβの絶対値を表わす。
【００６４】
この式（１０）を満すｑｉの点の無数の集合は信号源ｉが存在する曲面を与え、遠距離場（平面波）モデルを用いて推定された方向（又は円錐面）とあわせて用いることにより、センサから信号源ｉまでの距離を推定することができる。これにより、２つ以上の信号源が同一方向もしくは互いに近い方向にある場合であっても、距離が違っていれば、それらを区別でき、パーミュテーション問題を適切に解決することが可能となる。
【００６５】
この第３実施形態の例では、この信号源が存在する曲面を表わす位置情報と前記方向情報とを用いて分離行列に対するパーミュテーション問題を解決する。
第３実施形態の機能構成を図１７に、その処理手順を図１８にそれぞれ示す。センサは３個以上が２次元又は３次元に配置されるがこの実施形態においては例えば図２０に示すようにセンサ１1と１2の間隔に対しセンサ１2と１3の間隔は１０〜２０倍好ましくは１５倍程度とされる。先の実施形態と同様、各観測信号ｘ1（ｔ），…，ｘJ（ｔ）は、それぞれ周波数領域の信号Ｘ1（ω，ｍ）,…,ＸJ（ω，ｍ）に変換され（ステップＳ１１）、更に独立成分分析により、各周波数ごとに分離行列Ｗ（ω）が算出される（ステップＳ１２）、その各分離行列Ｗ（ω）の逆行列として行列Ｈ（ω）が算出される（ステップＳ１３）。またこの例では第２実施形態と同様に周波数ごとの逆行列Ｈ（ω）の各列から選択した要素対を用い、円錐面が一つ好ましくは複数推定される（ステップＳ１４）。この第３実施形態においては距離比算出部３１において、周波数ごとの逆行列Ｈ（ω）から列ごとに選択した要素対を用いその対応するセンサと１つの信号源ｉとの距離比、つまり式（１０）と式（２１）から次式（１０’）を算出する（ステップＳ３５）。
【００６６】
‖ｑz(i)−ｄj’‖／‖ｑz(i)−ｄj‖＝｜Ａjz(i)（ω）／Ａj’z(i)（ω）｜-＝｜Ｈji（ω）／Ｈj’i（ω）｜＝ＤＲi，jj’（ω）                          …（１０’）
距離比算出部３１において行われるステップＳ３５の処理の具体例を図１９を参照して説明する。この処理は図１４の処理の流れとほぼ同様である。パラメータｉを０に初期化し（ステップＳ２０）、次に、ｉを＋１し（ステップＳ２１）、パラメータｊ,ｊ’としてＪ以下の自然数を、例えばランダムに選択しｊ≠ｊ’、かつ１度選択した組は選択しない（ステップＳ２３）。センサｊの位置ベクトルｄｊと、センサｊ’の位置ベクトルｄj’とを抽出し（ステップＳ２４）、逆行列Ｈ（ω）のｉ列目から、要素Ｈji(ω)と、Ｈj’i(ω)とを、選択する（ステップＳ２５）。
【００６７】
この実施形態ではこれら選択した２つの要素の比ＤＲi，jj’（ω）を計算する（ステップＳ４１）。次に、ｉ＝Ｉであるか否かを判断し（ステップＳ２９）。ｉ＝ＩでなければステップＳ２１に戻り、ｉ＝２であれば処理を終了する。
【００６８】
距離比算出部３１で算出された距離比情報ＤＲi，jj’（ω）は、パーミュテーション解決部１７に送られ、パーミュテーション解決部１７は、到来方向決定部１６で推定された方向情報ｕi（ω）と、距離比算出部３１で算出された距離比情報ＤＲi，jj’（ω）とを用いて、分離行列算出部１２で算出された分離行列を行ってパーミュテーション問題を解決する。
Ｗ（ω）の行替えを行ってパーミュテーションの問題を解決する。例えば方向情報と距離比情報とから信号源２iの距離‖ｑi（ω）‖を距離推定部１７ｄで計算する（ステップＳ３６）。
この距離‖ｑi（ω）‖の計算方法を図２０を参照して説明する。信号源２1及び２2がセンサ１1，１2からみて同一方向Ｂにある。この場合センサ１1，１2及び遠距離場モデルで推定される信号源２1及び２2の方向ｕ1，ｕ2は同一の直線を定める。一方、間隔が大きいセンサ１2，１3及び近距離場モデルを用いて信号源２1が存在すると曲面６1は距離比
ＤＲ1，23（ω）＝｜Ｈ21（ω）／Ｈ31（ω）｜＝‖ｑ1−ｄ3‖／‖ｑ1−ｄ2‖
から推定でき、つまり‖ｑ1（ω）‖が推定できる。また信号源２2が存在する曲面６2は距離比
ＤＲ2，23（ω）＝｜Ｈ22（ω）／Ｈ32（ω）｜＝‖ｑ2−ｄ3‖／‖ｑ2−ｄ2‖
から推定でき、つまり‖ｑ2（ω）‖が推定できる。
【００６９】
信号源２1の位置は、直線ｕ1＝ｕ2と、曲面６1との共通部分上にあると推定でき、信号源２2の位置は、直線ｕ1＝ｕ2と、曲面６2との共通部分上にあると推定できる。例えば直線ｕ1＝ｕ2を表す式と曲面６1及び６2をそれぞれ表わす式との連立方程式を解くことにより‖ｑ1（ω）‖及び‖ｑ2（ω）‖を求めればよい。このようにして信号源の方向が同一や近似する場合であっても、信号源位置を区別することができる。
【００７０】
パーミュテーション解決部１７は、このようにして得られた各周波数ごとの信号源の距離‖ｑi（ω）‖が、所定の順序が例えば昇順になるように、分離行列Ｗ（ω）の行の入れ替えを行う。つまりこの入れ替えを行うためパーミュテーション行列Ｐ（ω）を算出する（ステップＳ３７）。そのパーミュテーション行列Ｐ（ω）の算出は実施形態２中のパーミュテーション解決部で行ったパーミュテーション行列の算出と同様に行えばよい。算出したパーミュテーション行列Ｐ（ω）を分離行列Ｗ（ω）の左から掛け、その結果Ｗ’（ω）＝Ｐ（ω）Ｗ（ω）を分離行列として出力する（ステップＳ３８）。
出力された分離行列Ｗ’（ω）は、時間領域変換部１８に送られ、時間領域変換部１８での信号分離に用いられる。
図２０から理解されるように、センサ１2と１3についてみると、信号源２1及びセンサ１2間の距離と信号源２1及びセンサ１3間の距離との差が大きく、一方信号源２2及びセンサ１2間の距離と信号源２2及びセンサ１3間の距離との差が小さい。従ってＤＲ1，23＝‖ｑ1−ｄ3‖／‖ｑ1−ｄ2‖と数値１との差の絶対値は比較的大きく、ＤＲ2，23＝‖ｑ2−ｄ3‖／‖ｑ2−ｄ2‖と数値１との差の絶対値は小さい。センサ１2と１3の間隔が大きい程、距離比ＤＲ1，23（ω）とＤＲ2，23（ω）との差が大きくなる。このため間隔が大きいセンサ対、この例では１2と１3が設けられる。
【００７１】
パーミュテーション解決部１７において、図１８中のステップＳ３７中に括弧で示すように計算した距離比ＤＲi，jj’（ω）を用いて、すべての周波数において逆行列Ｈ（ω）の１列目，２列目，３列目，…，Ｉ列目の順に、求まる距離比ＤＲi，jj’（ω）が昇順になるように並替行列を生成し、更に前記パーミュテーション行列Ｐ（ω）を生成してもよい。この場合、図１７中の距離推定部１７ｄは省略される。また到来方向決定部１６で方向情報ｕi（ω）を計算できなかったり、逆行列Ｈ（ω）の複数列について同一のｕi（ω）が計算された信号源ｉとセンサｉｊ’についてのみ、距離比ＤＲi，jj’（ω）を計算し、これら自体により、あるいは更に求めた‖ｑi（ω）‖により分離行列Ｗ（ω）の行の入れ替えを行ってもよい。つまり、方向情報ｕi（ω）を用いて分離行列Ｗ（ω）
の行の入れ替えを行い、その後、ｕi（ω）によって入れ替えができなかった行については距離比ＤＲi，jj’（ω）又は‖ｑi（ω）‖を用いてＷ（ω）の行の入れ替えを更に行ってもよい。実際には両入れ替えを同時に行うパーミュテーション行列Ｐ（ω）を生成する。
【００７２】
距離比ＤＲi，jj’（ω）をまず求め、分離行列Ｗ（ω）の行の入れ替えを行い、距離比ＤＲi，jj’（ω）を求めることができなかったものについて方向情報ｕi（ω）を用いて分離行列Ｗ（ω）の入れ替えを更に行ってもよい。この場合も両入れ替えを同時に行うパーミュテーション行列Ｐ（ω）を生成する。一般にはＤＲi，jj’（ω）よりｕi（ω）の方が高い精度で得られるため、ｕi（ω）による入れ替えを主とし、できないものについてＤＲi，jj’（ω）による入れ替えを行う方がよい。
また式（１０'）をその右辺の値を距離比ＤＲi，jj’（ω）と表わして、ｑiについて解くと、中心Ｏi，jj’（ω）及び半径Ｒi，jj’（ω）がそれぞれ下記式（１１）及び（１２）で与えられる球面となる。
【００７３】
Ｏi，jj’（ω）＝ｄj−（ｄj’−ｄj）／（ＤＲ2i，jj’（ω）−１）…（11）
Ｒi，jj’（ω）＝‖ＤＲi，jj’（ω）・（ｄj’−ｄj）／（ＤＲ2i，jj’（ω）−１）‖                                                             …（12）
例えばセンサ１j，１j’の各位置がｄj＝（０，０．１５，０），ｄj’＝（０，−０．１５，０）（単位はメートル）である場合、半径Ｒi，jj’をパラメータとしたとき、式（９）によって決まる球面の様子は図２１に示すようになる。
【００７４】
つまり信号源ｉは位置情報としての式（１１）及び（１２）で与えられる球面上に存在することになる。従って図１７中のパーミュテーション解決部１７中の距離推定部１７ｄを括弧書きで示すように曲面推定部とし、方向情報ｕi（ω）を求めることができなかった各ｉ，ｊｊ’について曲面推定部１７ｄで、式（１１）の半径Ｒi，jj’（ω）及び中心Ｏi，jj’（ω）を、図１８中のステップＳ３６中に括弧書きで示すようにそれぞれ計算し、これらＲi，jj’（ω）及び中心Ｏi，jj’（ω）がどの周波数の逆行列Ｈ（ω）でも同一順になるようにして、パーミュテーション行列Ｐ（ω）で求めてもよい。
【００７５】
この円錐面情報θ＾i，jj’（ω）又は方向情報ｕi（ω）と球面情報ＤＲi，jj’(ω)又はＲi，jj’(ω)を用いてもパーミュテーション行列Ｐ（ω）を求めることができない場合はその周波数について、従来の相関法（例えば、H.Sawada他“A robust and precise
method for solving the permutation problem of frequency-domain blind source separation.”，in Proc.Intl.Symp. on Independent Component Analysis and Blind Signal Separation (ICA 2003),2003,pp.505-510参照）を適用するとよい。図１７中に破線で示すように、信号分離部１９で分離された出力信号ｙ１（ｔ），…，ｙＩ（ｔ）を周波数領域変換部３３で例えば周波数領域変換部１１と同様の方法で周波数領域の信号Ｙ1（ω，ｍ），…，ＹＩ（ω，ｍ）に変換し、相関算出部３４において、これら周波数領域信号Ｙ1（ω，ｍ），…，ＹＩ（ω，ｍ）中のパーミュテーション解決部１７でパーミュテーション行列を作ることができなかった周波数成分ｆank成分と、パーミュテーション行列が得られた周波数中のｆankと隣接した周波数成分との相関を計算し、パーミュテーション解決部１７ではこの相関が大きくなるように周波数ｆankの分離行列の行の入れ替えを行い、この行入れ替えした分離行列に基づく分離出力信号ｙ1（ｔ），…，ｙＩ（ｔ）中のｆank成分について再び相関計算部３４で相関を計算し、同様のことを計算した相関が最大になるまで繰り返す。この相関の最大値が所定値に達しなかった場合は、更にそのｆankの成分と、パーミュテーション行列が得られている周波数成分中のｆankに対する基本波又は高調波(一般には高調波)関係にあるものとの各相関の和を相関計算部３４で計算し、この相関和が大きくなるように、ｆankに対する分離行列の行入れ替えを修正部１７ｅで行うことを、前記相関和が最大になるまで繰り返す。なお第２実施形態で説明したよう
に周波数領域の信号Ｘ（ω，ｍ）に対し分離行列Ｗ（ω）を適用し、周波数領域の分離信号Ｙ（ω，ｍ）を求める場合は、このＹ（ω，ｍ）を相関計算部３４の計算に用いればよい。
【００７６】
以下にこの第３実施形態の実験例を信号源として、室内で実測したインパルス応答を畳み込んで作成した混合音声を用いて分離実験を行った場合について説明する。図２２に示すように残響時間が１３０ミリ秒、３５５ｃｍ×４４５ｃｍ×２５０ｃｍ（高さ）の室内の水平面内でほぼ中心部で高さ１３５ｃｍの所にセンサとしての無指向性マイクロホン１1〜１4を直径４ｃｍの円上に等間隔で配し（これらマイクロホン１1〜１4は図中の左上に拡大して示した）、更に無指向性マイクロホン１5〜１8を直径３０ｃｍの円上に等間隔かつ、マイクロホン１1，１3，１5及び１7が、またマイクロホン１2，１4，１6及び１8がそれぞれ一直線上に配列するように配置した。マイクロホン１1から１5への方向を基準（０°）とし、マイクロホン配置の中心を原点とし、原点から１２０ｃｍ離れ、−３０°方向に音源２1を、＋３０°方向に音源２2を、＋９０°方向に音源２3を、−１５０°方向に音源２6をそれぞれ配置し、１５０°方向で原点からの距離６０ｃｍに音源２4を、距離１８０ｃｍに音源２5をそれぞれ配置した。サンプリング周波数ｋＨｚ、データ長６秒、フレーム長２０４８サンプル（２５６ミリ秒）、フレームシフト５１２サンプル（６４ミリ秒）とした。
【００７７】
周波数領域での分離行列Ｗ（ω）の各行のうち、マイクロホン対１1と１3、１2と１4、１1と１2、１2と１3に相当する行を用いて音源方向（円錐面）推定を行った。推定したヒストグラムを図２２に示す。図２３の横軸は推定方向を、縦軸は信号数である。これより推定方向は５つの群（クラスタ）があり、そのうち１５０°付近のクラスタは他のクラスタより２倍の大きさである。これは６個の音源中の２個が同一方向（１５０°付近）から到来していると推測される。４つの音源については推定した方向をもとにパーミュテーションを解決した。その結果を図２４に示す。図２４の横軸は周波数、縦軸は方向（度）を表わす。
【００７８】
残りの２個の音源については、広い間隔のマイクロホン対１７と１５、１７と１８、１6と１5、１6と１8を用いて音源が存在する球面の半径によって到来信号を区別した。マイクロホン対１7と１5を用いて推定した球面の半径を図２５に示す。図２５の横軸は周波数、縦軸は半径（ｍ）を表わす。
【００７９】
残響や推定誤差の影響により、位置情報だけではパーミュテーションを完全には解決することはできない。従って推定した位置情報に基づいて信号が矛盾なく分類できた周波数はその情報に基づいてパーミュテーション行列を生成し、残りの周波数については相関に基づく方法を用いてパーミュテーションを解決した。最後に時間領域の分離フィルタ係数を求める際にスペクトルの平滑化を行った。スペクトル平滑化については例えばH.Sawada他、“Spectral smoothing for frequency-domain blind source separation”，in Proc.IWAENC 2003, 2003, pp.311-314参照されたい。分離性能（ＳＩＲ）の評価結果を図２６に示す。図２６中の数値の単位はｄＢである。Ｃは相関法のみによってパーミュテーションを解決したもの、Ｄ＋Ｃは方向（円錐面）推定でパーミュテーションを解決し、解決できないものに相関法を適用したものＤ＋Ｓ＋Ｃは方向（円錐面）推定と球面推定とによりパーミュテーション解決を行った後、解決できなかった周波数について、相関法を用いたものである。この後者の方法によれば、６音源すべての分離ができ、ＳＩＲの改善は入力ＳＩＲから平均１７．１ｄＢとなった。
前記した第２実施形態の例ではセンサは２次元配置したが、一対のセンサにより推定する球面はこれらセンサの中心２等分線と対称に現われるため、信号源が３次元に存在している場合は２次元配置センサでは判別できなくなるため、センサも３次元配置とする必要がある。
【００８０】
以上説明したように、この第３実施形態においても式（９’）により円錐面情報を、また式（１０'）により曲面情報をそれぞれ推定しているため計算量が少ない。しかも円錐面と、距離比ＤＲi，jj’（ω）、又は距離‖ｑi（ω）‖あるいは球面の半径Ｒi，jj’（ω）とによりパーミュテーションを解決しているために２つ以上の信号源が同一方向もしくは互いに近い方向にある場合であっても、それらを区別できる。更にこれに相関法を加えれば一層確実に分離が可能である。なお球面情報としては計算が簡単であるからＤＲi，jj’（ω）が好ましい。
【００８１】
〔第４実施形態〕
この第４実施形態では、推定した円錐面の信頼性を検証し、信頼性が高いと判断した円錐面を分離行列のパーミュテーションの解決に用いる。例えば図２７に示すように円錐面推定部１４で推定された円錐面情報θ＾i，jj’（ω）は円錐面検証部４１で許容角情報格納部４２よりの許容角情報に基づき、信頼性があるか否かの検証がなされる。つまり式（９’）により求まる角度θ＾i，jj’（ω）はこれを求めるために用いたセンサ１j及び１j’の配列方向に対する相対角度であり、図４Ｂを参照して説明したように、角度θ＾i，jj’（ω）が０度付近及び１８０度付近の場合、正しいパーミュテーションが行われないと推定される。
【００８２】
従って、パーミュテーションを正しく行うことができると推定される角度の最小値θminと最大値θmaxが許容角情報として許容角情報格納部４２に格納され、円錐面検証部４１で推定円錐面情報θ＾i，jj’（ω）がθminとθmaxとの間にあれば信頼できる円錐面と判定されて出力され、つまりパーミュテーション解決に用いられるが、θminとθmaxとの間になければ、そのθ＾i，jj’（ω）は信頼性がないものとして破棄され、パーミュテーション解決に用いられない。例えば図１５中の円錐面４13は破棄される。
【００８３】
円錐面検証部４１で信頼性があると検証された円錐面情報θ＾i，jj’（ω）は図１２中の到来方向決定部１６へ送られ、又は図１７で説明したように、パーミュテーション解決部１７へ直接送られる。図１３及び図１８の処理手順において、ステップＳ１４の次に破線で示すようにステップＳ１００として推定円錐面θ＾i，jj’（ω）が信頼性があるか否かの検証が行われ、信頼性があるものについてのみ次のステップに移るようにされる。式（９’）において計算された偏角arg（Ｈji／Ｈj’i）に含まれる誤差をΔargＨ＾、推定された角度θ＾iに含まれる誤差をΔθ＾とすると、これらの比｜Δθ＾／ΔargＨ＾｜は式（９’）を偏微分して次式のように近似できる。
｜Δθ＾／ΔargＨ＾｜＝｜１／（ωｃ-1｜ｄj−ｄj’｜sinθ＾i）)｜  （13）
式（１３）をいくつかの周波数について計算した結果を図４Ｃに示す。この図４Ｃから推定された角度θ＾iの方向がセンサ１j，１j’のセンサペア軸に近い場合はarg（Ｈji／Ｈj’i）に含まれる誤差ΔargＨ＾は推定角度θ＾iに対し、大きな誤差を生じさせることがわかる。つまりセンサペア軸に近い方向の推定角度θ＾iを用いて分離行列Ｗ（ω）に対するパーミュテーションの問題を解決した場合、これは正しい解にならない可能性が高い。図４Ｂ及び図４Ｃから例えばθminは２０°程度、θmaxは１６０°程度とされる。図４Ｃからわかるように｜Δθ＾／ΔargＨ＾｜は周波数により可成り異なり、低い周波数ではΔargＨ＾が到来方向の誤差に大きく影響する。従って、低周波数については、信号源が存在する球面の推定に基づく情報ＤＲi，jj’、‖ｑi‖、Ｒi，jj’を用いて、あるいは相関法によりパーミュテーション問題を解決するとよい。
【００８４】
この第４実施形態によれば、推定円錐面情報中の信頼性のないものは破棄されるため、これにより悪影響を受けることなく、それだけ正確に到来方向を推定でき、従って正しいパーミュテーション行列Ｐ（ω）を生成することができ、つまり分離信号のＳＩＲ（性能）が向上する。
【００８５】
〔第５実施形態〕
第５実施形態では位置情報として距離比、又はこれより推定した球面情報を用いる。この機能構成は図１７中に示されている。図２８にこの実施形態の処理手順を簡略に示す。この場合はセンサは比較的広い間隔、例えば図２２に示した信号源が音源の場合、３０ｃｍとして少くとも２次元に配置される。
先の実施形態と同様に時間領域の観測信号は図１７中に示すように周波数領域の信号に変換され、分離行列生成部１２により分離行列Ｗ（ω）が生成され（ステップＳ５１）、これより逆行列生成部１３で逆行列Ｈ（ω）が生成される（ステップＳ５２）。各周波数における逆行列Ｈ（ω）の各列ごとに球面情報が推定される（ステップＳ５３）。この球面情報は第３実施形態で説明した方法と同様に算出される。つまり球面情報は距離算出部３１により算出された距離比ＤＲi，jj’（ω）、あるいは距離推定部又は曲面推定部１７ｄで算出された‖ｑi‖、又は半径Ｒi，jj’（ω）及び中心Ｏi，jj'（ω）である。
【００８６】
これら球面情報を用いて、その配列順が予め決めた順になるよう分離行列Ｗ（ω）に対するパーミュテーション行列が生成され分離行列Ｗ（ω）の行の入れ替えが行われる（ステップＳ５４）。第３実施形態におけるこの処理はパーミュテーション解決部１７で行われるが、ただしこの第５実施形態では球面情報のみが用いられる。この球面情報のみでパーミュテーションを解決できない周波数があれば、その周波数に対しては前述した相関法による解決が行われる（ステップＳ５５）。
【００８７】
図２２に示した室内に配したマイクロホン１6及び１8を用い、１２０°方向において原点から６０ｃｍ及び１５０ｃｍの所に音源２4及び２5を配し、実験室内で測定したインパルス応答に話者４名の音声を畳み込んで作成した１２通りの組合せの混合音声と対する分離実験を行った。パーミュテーション解決は推定された半径Ｒ4，68（ω）とＲ5，68（ω）を比較してＲ4，68（ω）≦Ｒ5，68（ω）となるようにして行った。Ｒ4，68及びＲ5，68中の最大値が、Ｒ4，68及びＲ5，68中の最小値にしきい値Ａthを乗算した値以上、つまりｍａｘ（Ｒ4，68，Ｒ5，68）≧Ａth・ｍｉｎ（Ｒ4，68，Ｒ5，68）を満す周波数は音源位置（Ｒ4，68（ω），Ｒ5，68（ω））によるパーミュテーションの解決が信頼できると判断して、音源位置による方法を適用し、それ以外の周波数では相関による方法を適用した。従ってＡth＝１．０のときは全て音源位置による方法を適用しＡthが無限大のときは全て相関による方法を適用することになる。しきい値Ａthを変化させながら分離性能としてＳＩＲ（信号干渉比）を音声の組合せごとにプロットした結果を図２９に示す。この結果より音源位置のみによる方法では全体的に性能が悪く、相関のみによる方法では分離性能がばらつき、両者を併用すると安定して高い性能が得られることがわかる。またしきい値Ａthとしては比較的大きい値が好ましく、８〜１６で全体の１／５〜１／１０程度の周波数を音源位置によって決定できることがわかる。
【００８８】
また図２２に示した実験条件で音源位置による方法と相関による方法を併用した実験結果を図２６中のＤ＋Ｃの欄に示す。これよりこの併用方法によれば可成り性能よく分離が行われることがわかる。この結果は式（９'）を用いる到来方向（円錐面情報）による方法と相関による方法とを併用した場合と同様の傾向がある。なお、到来方向による方法と相関による方法とを併用した場合の実験結果は、相関による方法の説明で引用した文献に示されている。
【００８９】
この第５実施形態によれば式（１０'）により球面情報を求めているため計算量が少ない。なお球面情報としては距離比ＤＲi，jj’（ω）が好ましい。
【００９０】
[第６実施形態]
第６実施形態は例えば一つの推定円錐面に基づき、パーミュテーションの問題を解決す
る。図１２中に破線で示すように円錐面推定部１４で推定された円錐面θ＾i，jj’（ω１）、・・・    θ＾i，jj’（ωN）、はパーミュテーション解決部１７へ直接入力され、図１３中に破線で示すようにステップｓ１４で設定された円錐面θ＾i，jj’（ωn）（ｎ＝１、・・・、Ｎ）がステップＳ１５の処理をすることなくステップＳ１６でどの周波数においても例えばθ＾i，jj’（ωn）が昇順になるようなパーミュテーション行列Ｐ（ω）が算出される。
この場合、ステップＳ１４で信号源ｉごとに１つの円錐面を推定するのみでもよい。また図に示していないが、行入れ替えができなかったものについては、第３実施形態で説明したように相関法による行入れ替えを得るようにしてもよい。
この第６実施形態によれば逆行列Ｈ（ω）の各列ごとの二つの要素比をとることによりスケーリングの問題が簡単に除去され、式（９'）の演算をすればよく計算時間が短くて済む。
【００９１】
［信号分離のまとめ］
以上述べたブラインド信号分離におけるパーミュテーションを解決した分離行列の生成方法を、逆行列の生成以後の処理を図３０Ａ、図３０Ｂ、図３０Ｃに示す。図３０Ａでは逆行列の各列の要素比から円錐面を推定し（ステップＳ６１）、必要に応じて信頼性のない円錐面を破棄し（ステップＳ６２）、これらの複数の円錐面から共通直線の方向を決定し（ステップＳ６３）、共通直線方向を用いてパーミュテーション行列Ｐ（ω）を生成して分離行列の入れ替えを行い（ステップＳ６４）、パーミュテーションを解決できなかった周波数に対しては相関法による分離行列の行入れ替えを行う（ステップＳ６５）。図３０Ａ中に破線で示すように、ステップＳ６からステップＳ６４に直ちに移り、円錐面を直接用いてパーミュテーション行列Ｐ（ω）を生成してもよい。
【００９２】
図３０Ｂでは前記ステップＳ６１の円錐面推定、必要に応じて前記ステップＳ６２の円錐面破棄を行い、ステップＳ６３の共通直線方向推定を行い、これら共通直線方向を用いて逆行列の列の入れ替えを行う（ステップＳ６６）。この際、図中破線で示すようにステップＳ６１で推定された円錐面又はステップＳ６２で処理後の円錐面を直接、ステップＳ６６の処理に用いてもよい。次に円錐面の推定又は共通直線方向の決定ができなかった、又は円錐面或いは共通直線方向が不確かなものと対応する、或いは同一値となったものについて、その逆行列の列の要素比から球面情報を推定し（ステップＳ６７）、この推定した球面情報を用いて逆行列の入れ替えを行ってパーミュテーション行列Ｐ（ω）を生成して分離行列の行入れ替えを行う（ステップＳ６４）。更にこのパーミュテーション行列Ｐ（ω）を作ることができなかった周波数に対しては相関による方法を適用する（ステップＳ６５）。
【００９３】
図３０Ｃはまず逆行列の各列の要素比から球面情報を推定し（ステップＳ６８）、この球面情報を用いて逆行列の列の入れ替えを行い（ステップＳ６９）、この列の入れ替えを行うことができなかった、又は球面情報が不確かなものについて、その逆行列の列の要素比から円錐面を推定し（ステップＳ７０）、これらの複数の円錐面の共通直線方向を決定し（ステップＳ７１）、この決定した方向又はステップＳ７０で推定した円錐面を直接用いて、逆行列の列の入れ替えを行ってパーミュテーション行列Ｐ（ω）を生成し、分離行列の入れ替えを行う（ステップＳ６４）。パーミュテーション行列を作ることができなかった周波数に対して相関法を適用する（ステップＳ６５）。
【００９４】
これらに対し更に図２８に示した方法がある。
第３実施形態、第５実施形態及び第６実施形態においても、第２実施形態で説明したように、周波数領域の分離行列Ｗ（ω）と観測信号Ｘ（ω）とを用いて信号分離を行い、その後、分離された周波数領域信号Ｙ（ω）を時間領域信号ｙ（ｔ）に変換してもよい。
【００９５】
第１実施形態の説明では２次元における信号の到来方向の推定を行ったが、第２実施形態中でも説明したように、３次元における信号の到来方向の推定にも適用できる。また第２乃至第６実施形態を２次元における信号分離に適用することもできる。この場合は推定円錐面θ＾i，jj’（ω）はその推定に用いたセンサ対のセンサペア軸に対し対称の２方向となり、推定球面Ｒi，jj’（ω）又はＤＲi，jj’（ω）は円の半径又はこれとほぼ対応するものになる。
【００９６】
図５、図１２、図１７にそれぞれ示した装置、また第５実施形態による装置はそれぞれコンピュータに機能させることもできる。この場合、それぞれと対応した処理手順、つまり図６、図１３、図１８などに示した各過程をコンピュータに実行させるためのプログラムを記録した磁気ディスク、半導体メモリ、ＣＤ−ＲＯＭなどから、コンピュータ内のメモリにインストールし、又は通信回線を介してコンピュータ内のメモリに前記プログラムをダウンロードして、コンピュータにそのプログラムを実行させればよい。
【図面の簡単な説明】
【００９７】
【図１】センサアレーと到来信号の各センサへの到来時間差の関係を説明するための図。
【図２】２音源混合信号に対し、ＩＣＡにより計算した分離行列の各行のゲイン指向特性を示す図。
【図３】３音源混合信号に対し、ＩＣＡにより計算した分離行列の各行のゲイン指向特性を示す図。
【図４】図４Ａは予備実験に用いたセンサと信号源の関係を説明するための図、図４Ｂは前記予備実験の結果を示す図、図４Ｃは推定角度とその誤差及び偏角誤差との関係を示す図である。
【図５】この発明を信号到来方向の推定に適用した第１実施形態の機能構成例を示すブロック図。
【図６】第１実施形態の処理手順の例を示す流れ図。
【図７】図５中の角度計算部の具体的機能構成例を示すブロック図。
【図８】図６中のステップＳ４の具体的処理手順の例を示す流れ図。
【図９】この第１実施形態により方向推定した実験結果を示す図。
【図１０】図１０Ａは音源が２個、マイクロホンが３個の場合に第１実施形態により方向推定した実験結果を示す図、図１０Ｂは図１０Ａについての実験と同一条件でＭＵＳＩＣ法により方向推定した実験結果を示す図である。
【図１１】図１１Ａは音源が２個、マイクが３個の場合に第１実施形態により方向推定した実験結果を示す図、図１１Ｂは図１１Ａについての実験と同一条件でＭＵＳＩＣ法により方向推定した実験結果を示す図である。
【図１２】この発明をブラインド信号分離に適用した第２実施形態の機能構成例を示すブロック図。
【図１３】第２実施形態の処理手順の例を示す流れ図。
【図１４】図１３中のステップＳ１４の具体的処理手順の例を示すブロック図。
【図１５】複数円錐面の共通直線方向を説明するための図。
【図１６】図１２中の到来方向決定部１６の具体的機能構成の他の例を示すブロック図。
【図１７】この発明をブラインド信号分離に適用した第３実施形態の機能構成例を示すブロック図。
【図１８】第３実施形態の処理手順の例を示す流れ図。
【図１９】図１８中のステップＳ３５の具体的処理手順の例を示す流れ図。
【図２０】センサ配置と信号源位置と推定曲面との関係を説明するための図。
【図２１】推定球面の例を示す図。
【図２２】実験に用いた室とマイクロホン配置と音源との関係を示す図。
【図２３】小間隔マイクロホン対を用いた推定方向のヒストグラム。
【図２４】その周波数に対する推定方向の分布を示す図。
【図２５】音源２4及び２5に対する推定球面の半径の周波数に対する分布を示す図。
【図２６】各方法での実験結果を示す図。
【図２７】この発明の第４実施形態の要部の機能構成例を示すブロック図。
【図２８】この発明をブラインド信号分離に適用した第５実施形態の処理手順の例を示す流れ図。
【図２９】第５実施形態に対する実験結果を示す図。
【図３０】図３０Ａは複数の円錐面の共通直線方向を求めて分離行列のパーミュテーションを解決する処理手順の要部を示す流れ図、図３０Ｂは円錐面の推定と球面の推定を利用して分離行列のパーミュテーションを解決する処理手順の要部を示す流れ図、図３０Ｃは円錐面の推定と球面の推定を利用して分離行列のパーミュテーションを解決する他の処理手順の要部を示す流れ図である。
【図１】
0003881367.tif 000009
【図２】
0003881367.tif 000010
【図３】
0003881367.tif 000011
0003881367.tif 000012
【図５】
0003881367.tif 000013
【図６】
0003881367.tif 000014
【図７】
0003881367.tif 000015
【図８】
0003881367.tif 000016
【図９】
0003881367.tif 000017
0003881367.tif 000018
0003881367.tif 000019
【図１２】
0003881367.tif 000020
【図１３】
0003881367.tif 000021
【図１４】
0003881367.tif 000022
【図１５】
0003881367.tif 000023
【図１６】
0003881367.tif 000024
【図１７】
0003881367.tif 000025
【図１８】
0003881367.tif 000026
【図１９】
0003881367.tif 000027
【図２０】
0003881367.tif 000028
【図２１】
0003881367.tif 000029
【図２２】
0003881367.tif 000030
【図２３】
0003881367.tif 000031
【図２４】
0003881367.tif 000032
【図２５】
0003881367.tif 000033
【図２６】
0003881367.tif 000034
【図２７】
0003881367.tif 000035
【図２８】
0003881367.tif 000036
【図２９】
0003881367.tif 000037
0003881367.tif 000038