(57)【特許請求の範囲】
【請求項１】
複数の信号源から発せられた源信号の混合からなる混合信号を、当該源信号の推定信号である分離信号に分離する信号分離装置であって、
複数のセンサで観測された前記混合信号を、周波数領域の混合信号に変換する周波数領域変換部と、
前記周波数領域の混合信号を用いて、周波数毎に分離行列を算出する分離行列算出部と、
前記分離行列の一般化逆行列を算出する逆行列算出部と、
前記一般化逆行列を構成する基底ベクトルの正規化を行い、正規化基底ベクトルを算出する基底ベクトル正規化部と、
前記正規化基底ベクトルをクラスタリングし、クラスタを生成するクラスタリング部と、
前記クラスタの中心ベクトルと前記正規化基底ベクトルとを用い、前記分離行列の要素を並び替えるための順列を算出する順列算出部と、
前記順列に従って前記分離行列の行を並び替えた第２の分離行列を生成する並び替え部と、
前記第２の分離行列に、各センサで観測された前記混合信号から変換された前記周波数領域の混合信号を要素とする混合信号ベクトルを乗じ、周波数領域の各分離信号を要素とする分離信号ベクトルを生成する分離信号生成部と、
を有する信号分離装置。
【請求項２】
請求項１記載の信号分離装置であって、
前記基底ベクトル正規化部は、
前記基底ベクトルの周波数依存性を排除する正規化を行うことを特徴とする信号分離装置。
【請求項３】
請求項２記載の信号分離装置であって、
前記基底ベクトルの周波数依存性を排除する正規化は、
前記基底ベクトルが具備する特定の１つの要素を基準として当該基底ベクトルの各要素の偏角を正規化し、さらに、各要素の偏角を周波数に比例した値で除算する正規化であることを特徴とする信号分離装置。
【請求項４】
請求項２記載の信号分離装置であって、
前記基底ベクトルの周波数依存性を排除する正規化は、
前記基底ベクトルAp（f）（ｐ＝１，...，Ｎ、Ｎは前記信号源の数）の各要素Aqp（f）（ｑ＝１，...，Ｍ、Ｍは前記混合信号を観測するセンサの数）に対し、ｅｘｐをネピア数とし、ａｒｇ［・］を偏角とし、ｆを周波数とし、ｊを虚数単位とし、ｃを信号の伝達速度とし、ＱをＭ以下の自然数から選択された基準値とし、ｄを実数とした場合における、
【数６４】
0004406428.tif 000067
の演算によって行われる正規化であることを特徴とする信号分離装置。
【請求項５】
請求項４記載の信号分離装置であって、
前記ｄは、
前記要素AQp（f）に対応する基準センサと他のセンサとの最大距離ｄmaxであることを特徴とする信号分離装置。
【請求項６】
請求項１記載の信号分離装置であって、
前記基底ベクトル正規化部は、
前記基底ベクトルの周波数依存性を排除する正規化と、ノルムを規定値にする正規化とを行うことを特徴とする信号分離装置。
【請求項７】
複数の信号源から発せられた源信号の混合からなる混合信号を、当該源信号の推定信号である分離信号に分離する信号分離装置であって、
複数のセンサで観測された前記混合信号を周波数領域の混合信号に変換する周波数領域変換部と、
前記周波数領域の混合信号を用いて、周波数毎に分離行列を算出する分離行列算出部と、
前記分離行列の一般化逆行列を算出する逆行列算出部と、
前記一般化逆行列を構成する基底ベクトルの正規化を行い、正規化基底ベクトルを算出する基底ベクトル正規化部と、
前記正規化基底ベクトルをクラスタリングし、クラスタを生成するクラスタリング部と、
前記周波数領域の混合信号と前記分離行列或いはその行を並び替えた分離行列とから求めた分離信号のエンベロープと、前記クラスタの中心ベクトルと、前記正規化基底ベクトルとを用い、前記分離行列の要素を並び替えるための順列を算出する順列算出部と、
前記順列に従って前記分離行列の行を並び替えた第２の分離行列を生成する並び替え部と、
周波数毎に前記順列の信頼度を評価し、前記順列の信頼度が低いと評価された周波数に対し、前記分離信号のエンベロープを用いて第２の順列を算出し、この周波数のみについて第２の分離行列の行を前記第２の順列に従って並び替えた第３の分離行列を生成する順列評価部と、
前記第３の分離行列に、各センサで観測された前記混合信号から変換された前記周波数領域の混合信号を要素とする混合信号ベクトルを乗じ、各信号源に対応する周波数領域の分離信号を要素とする分離信号ベクトルを生成する分離信号生成部と、
を有する信号分離装置。
【請求項８】
複数の信号源から発せられた源信号の混合からなる混合信号を、当該源信号の推定信号である分離信号に分離する信号分離装置であって、
複数のセンサで観測された前記混合信号を周波数領域の混合信号に変換する周波数領域変換部と、
前記周波数領域の混合信号を用い、周波数毎に、分離行列と、周波数領域の各分離信号を要素とする分離信号ベクトルと、を算出する信号分離部と、
前記各分離信号に対応する前記分離行列の一般化逆行列の列である基底ベクトルを生成する逆行列算出部と、
前記基底ベクトルを正規化し、正規化された基底ベクトルを生成する正規化部と、
前記正規化された基底ベクトルをクラスタリングしてクラスタを生成するクラスタリング部と、
前記各クラスタの分散を指標として１個以上のクラスタを選択する分散判定部と、
選択されたクラスタの中心ベクトルとの距離が最も近い正規化された基底ベクトルに対応する分離信号を、目的信号を含む選択信号として選択する選択部と、を有する信号分離装置。
【請求項９】
請求項８に記載の信号分離装置であって、
更に、前記周波数領域の混合信号と前記基底ベクトルとを用いて時間周波数マスクを生成するマスク生成部と、
前記時間周波数マスクを前記選択部で選択された前記選択信号に適用し、マスク処理選択信号を生成するマスキング部と、を有することを特徴とする信号分離装置。
【請求項１０】
請求項９に記載の信号分離装置であって、
前記マスク生成部は、
前記周波数領域の混合信号を用いて白色化行列を生成する白色化行列生成部と、
前記白色化行列を用い、前記周波数領域の混合信号を要素とする混合信号ベクトルを変換した白色化混合信号ベクトルと、前記基底ベクトルを変換した白色化基底ベクトルとを算出する白色化部と、
前記白色化混合信号ベクトルと前記白色化基底ベクトルとがなす角度を時間周波数毎に算出する角度計算部と、
前記角度を要素とする関数である前記時間周波数マスクを生成する関数演算部と、を有することを特徴とする信号分離装置。
【請求項１１】
請求項１０に記載の信号分離装置であって、
前記白色化行列は、
周波数をfとし、離散時間をτとし、前記混合信号ベクトルをX（f,τ）とし、ベクトル＊を時間平均したベクトルを<＊>τとし、ベクトル＊の複素共役転置ベクトルを＊Hとし、R（f）=<X（f,τ）・X（f,τ）H>τとした場合におけるV（f）=R（f）−1/2であり、
前記白色化部は、
Z（f,τ）=V（f）・X（f,τ）の演算により、前記白色化混合信号ベクトルZ（f,τ）を算出し、前記基底ベクトルをA（f）とした場合におけるB（f）=V（f）・A（f）の演算により、前記白色化基底ベクトルB（f）を算出し、
前記角度計算部は、
ベクトル＊の絶対値を｜＊｜とし、ベクトル＊のノルムを‖＊‖とした場合におけるθ（f,τ）=cos−1（|B H（f）・Z（f,τ）|／‖B（f）‖・‖Z（f,τ）‖）の演算により、前記角度θ（f,τ）を算出し、
前記関数演算部は、
α，g，θTを実数とした場合におけるロジスティック関数M（θ（f,τ））=α／（1+eg・（θ（f,τ）−θT））を前記時間周波数マスクとして算出することを特徴とする信号分離装置。
【請求項１２】
請求項９に記載の信号分離装置であって、
前記マスク生成部は、
前記周波数領域の混合信号を用いて生成された混合信号ベクトルX（f,τ）を、周波数に依存しない周波数正規化ベクトルX'（f,τ）に正規化する周波数正規化部と、
前記周波数正規化ベクトルX'（f,τ）を、ノルムが所定の値になるノルム正規化ベクトルX''（f,τ）に正規化する第１ノルム正規化部と、
前記選択信号に対応するセントロイドηιを抽出するセントロイド選択部と、
前記選択信号に対応するセントロイドηιを、そのノルムが前記所定の値になるノルム正規化セントロイドηι’に正規化する第２ノルム正規化部と、
前記周波数正規化ベクトルX'（f,τ）と前記ノルム正規化セントロイドηι’との距離の自乗DS（f,τ）を算出する自乗距離算出部と、
前記距離の自乗DS（f,τ）を要素とする関数を用いた前記時間周波数マスクを生成する関数生成部と、を有することを特徴とする信号分離装置。
【請求項１３】
請求項８に記載の信号分離装置であって、
前記正規化部は、
前記基底ベクトルの周波数依存性を排除する正規化を行うことを特徴とする信号分離装置。
【請求項１４】
請求項１３記載の信号分離装置であって、
前記基底ベクトルの周波数依存性を排除する正規化は、
前記基底ベクトルが具備する特定の１つの要素を基準として当該基底ベクトルの各要素の偏角を正規化し、さらに、各要素の偏角を周波数に比例した値で除算する正規化であることを特徴とする信号分離装置。
【請求項１５】
請求項１４記載の信号分離装置であって、
前記基底ベクトルの周波数依存性を排除する正規化は、
前記基底ベクトルAp（f）（ｐは自然数）の各要素Aqp（f）（ｑ＝１，...，Ｍ、Ｍは前記混合信号を観測するセンサの数）に対し、ｅｘｐをネピア数とし、ａｒｇ［・］を偏角とし、ｆを周波数とし、ｊを虚数単位とし、ｃを信号の伝達速度とし、ＱをＭ以下の自然数から選択された基準値とし、ｄを実数とした場合における、
【数６５】
0004406428.tif 000068
の演算によって行われる正規化であることを特徴とする信号分離装置。
【請求項１６】
請求項１５記載の信号分離装置であって、
前記実数ｄは、
前記基準値Ｑに対応する基準センサと他のセンサとの最大距離ｄmaxであることを特徴とする信号分離装置。
【請求項１７】
請求項８記載の信号分離装置であって、
前記正規化部は、
前記基底ベクトルの周波数依存性を排除する正規化と、ノルムを規定値にする正規化とを行うことを特徴とする信号分離装置。
【請求項１８】
請求項８に記載の信号分離装置であって、
前記分散判定部は、
分散が最小となるクラスタを選択し、選択されたクラスタに対応する前記分離信号を前記選択信号として選択することを特徴とする信号分離装置。
【請求項１９】
複数の信号源から発せられた源信号の混合からなる混合信号を、当該源信号の推定信号である分離信号に分離する信号分離装置であって、
複数のセンサで観測された前記混合信号を周波数領域の混合信号に変換する周波数領域変換部と、
前記周波数領域の混合信号から構成される混合信号ベクトルの正規化を行い、正規化ベクトルを算出するベクトル正規化部と、
前記正規化ベクトルをクラスタリングし、クラスタを生成するクラスタリング部と、
ｋ番目の前記クラスタに属する前記正規化ベクトルの時間周波数に対応する前記混合信号ベクトルから所定番目の要素を抽出し、これをｋ番目の要素とした分離信号ベクトルを生成する分離信号生成部と、を有し、
前記分離信号生成部で生成された前記分離信号ベクトルの各要素が、周波数領域の各分離信号であることを特徴とする信号分離装置。
【請求項２０】
請求項１９に記載の信号分離装置であって、
前記ベクトル正規化部は、
前記周波数領域の混合信号から構成される混合信号ベクトルの周波数依存性を排除する正規化を行うことを特徴とする信号分離装置。
【請求項２１】
請求項２０記載の信号分離装置であって、
前記混合信号ベクトルの周波数依存性を排除する正規化は、
前記混合信号ベクトルが具備する特定の１つの要素を基準として当該基底ベクトルの各要素の偏角を正規化し、さらに、各要素の偏角を周波数に比例した値で除算する正規化であることを特徴とする信号分離装置。
【請求項２２】
請求項２１記載の信号分離装置であって、
前記混合信号ベクトルの周波数依存性を排除する正規化は、
Ｍを前記センサの数とし、ｑ＝１，...，Ｍとし、前記混合信号ベクトルの各要素をXq（f,τ）とし、ｅｘｐをネピア数とし、ａｒｇ［・］を偏角とし、ｊを虚数単位とし、ｃを信号の伝達速度とし、ＱをＭ以下の自然数から選択された値とし、ｄを実数とし、ｆを周波数とし、τを離散時間とした場合における、
【数６６】
0004406428.tif 000069
の演算によって行われることを特徴とする信号分離装置。
【請求項２３】
請求項２２記載の信号分離装置であって、
前記ｄは、
前記要素XQ（f,τ）に対応するセンサと他のセンサとの最大距離ｄmaxであることを特徴とする信号分離装置。
【請求項２４】
請求項１９に記載の信号分離装置であって、
前記ベクトル正規化部は、
前記混合信号ベクトルの周波数依存性を排除する正規化と、ノルムを規定値にする正規化とを行うことを特徴とする信号分離装置。
【請求項２５】
複数の信号源から発せられた源信号の混合からなる混合信号を、当該源信号の推定信号である分離信号に分離する信号分離方法であって、
複数のセンサで観測された前記混合信号を周波数領域の混合信号に変換して出力する過程と、
入力された前記周波数領域の混合信号を用い、周波数毎に分離行列を算出して出力する過程と、
入力された前記分離行列の一般化逆行列を算出して出力する過程と、
入力された前記一般化逆行列を構成する基底ベクトルの正規化を行い、正規化基底ベクトルを算出して出力する過程と、
入力された前記正規化基底ベクトルをクラスタリングし、クラスタを生成して出力する過程と、
入力された前記クラスタの中心ベクトルと前記正規化基底ベクトルとを用い、前記分離行列の要素を並び替えるための順列を算出して出力する過程と、
前記順列に従って前記分離行列の行を並び替えた第２の分離行列を生成する過程と、
前記第２の分離行列に、各センサで観測された前記混合信号から変換された前記周波数領域の混合信号を要素とする混合信号ベクトルを乗じ、周波数領域の各分離信号を要素とする分離信号ベクトルを生成する過程と、
を有する信号分離方法。
【請求項２６】
複数の信号源から発せられた源信号の混合からなる混合信号を、当該源信号の推定信号である分離信号に分離する信号分離方法であって、
複数のセンサで観測された前記混合信号を周波数領域の混合信号に変換して出力する過程と、
入力された前記周波数領域の混合信号を用い、周波数毎に、分離行列と、周波数領域の各分離信号を要素とする分離信号ベクトルと、を算出して出力する過程と、
入力された前記分離行列の一般化逆行列の列である基底ベクトルを生成する過程と、
前記基底ベクトルを正規化し、正規化された基底ベクトルを生成する過程と、
前記基底ベクトルを正規化し、正規化された基底ベクトルを生成する過程と、
前記正規化された基底ベクトルをクラスタリングしてクラスタを生成する過程と、
前記各クラスタの分散を指標として１個以上のクラスタを選択する過程と、
選択されたクラスタの中心ベクトルとの距離が最も近い正規化された基底ベクトルに対応する分離信号を、目的信号を含む選択信号として選択する過程と、
を有する信号分離方法。
【請求項２７】
複数の信号源から発せられた源信号の混合からなる混合信号を、当該源信号の推定信号である分離信号に分離する信号分離方法であって、
複数のセンサで観測された前記混合信号を周波数領域の混合信号に変換して出力する過程と、
入力された前記周波数領域の混合信号から構成される混合信号ベクトルの正規化を行い、正規化ベクトルを算出して出力する過程と、
入力された前記正規化ベクトルをクラスタリングし、クラスタを生成して出力する過程と、
ｋ番目の前記クラスタに属する前記正規化ベクトルの時間周波数に対応する前記混合信号ベクトルから所定番目の要素を抽出し、これをｋ番目の要素とした分離信号ベクトルを生成して出力する過程と、を有し、
生成された前記分離信号ベクトルの各要素が、周波数領域の各分離信号であることを特徴とする信号分離方法。
【請求項２８】
請求項１から２４の何れかに記載の信号分離装置としてコンピュータを機能させるための信号分離プログラム。
【請求項２９】
請求項２８記載の信号分離プログラムを格納したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、信号処理の技術分野に属し、特に複数の源信号が空間内で混合された混合信号から源信号を抽出する技術に関する。
【背景技術】
【０００２】
複数のセンサを用いることによりある信号を抽出してその他の信号を抑圧する従来の技術として、ビームフォーマ（beamformer）（ビームフォーミング（beamforming）とも呼ぶ）が広く知られている（例えば、非特許文献１参照）。しかし、ビームフォーマでは、目的信号の方向に関する情報を必要とするため、そのような情報が得られない（または推定できない）状況では利用しにくいという欠点がある。
また、より新しい別の技術としては、ブラインド信号分離（BSS：Blind Signal Separation）が挙げられる（例えば、非特許文献２参照）。これは、上記ビームフォーマが必要とした情報を必要としないという点で優れており、様々な状況での利用が期待されている。以下、このブラインド信号分離による信号分離について説明する。【０００３】
［ブラインド信号分離］
まず、ブラインド信号分離の定式化を行う。すべての信号はあるサンプリング周波数ｆｓでサンプリングされ、離散的に表現されるものとする。また、Ｎ個の信号が混合されてＭ個のセンサで観測されたとする。以下では、信号の発生源からセンサまでの距離により信号が減衰・遅延し、また壁などにより信号が反射して伝送路の歪みが発生しうる状況を扱う。このような状況で混合される信号は、信号源ｋからセンサｑ（ｑはセンサの番号を示す〔ｑ＝１，…，Ｍ〕。また、ｋは信号源の番号を示す〔ｋ＝１，…，Ｎ〕。）へのインパルス応答hqk（r）による畳み込み混合【０００４】
【数１】
0004406428.tif 000002
で表現できる。ここでｔはサンプリング時刻を示している。また、sk（t）は、サンプリング時刻ｔにおいて信号源から発せられる源信号を示しており、xq（t）は、サンプリング時刻ｔにおいてセンサｑで観測される信号を示している。また、ｒは掃引のための変数を示している。
一般的なインパルス応答hqk（r）は、適当な時間経過後にパルス的な強い応答を持ち、時間と共に減衰していく。ブラインド信号分離の目的は、源信号s1（t）,…,sN（t）やインパルス応答h11（r）,…,h1N（r）,…,hM1（r）,…,hMN（r）を知らずに、観測信号（以下「混合信号」と呼ぶ）x1（t）,…,xM（t）のみから、源信号s1（t）,…,sN（t）にそれぞれ対応する分離信号y1（t）,…,yN（t）を求めることにある。【０００５】
［周波数領域］
次に、従来のブラインド信号分離の手順について説明する。
ここでは周波数領域において分離の操作を行う。そのためにセンサｑでの混合信号ｘｑ（ｔ）にＬ点の短時間離散フーリエ変換（STFT：Short-Time Fourier Transform）を適用し、周波数ごとの時間系列【０００６】
【数２】
0004406428.tif 000003
を求める。ここでｆは周波数でありf=0,fs/L,…,fs（L-1）/Lと離散化されている（ｆｓはサンプリング周波数）。また、τは離散時間であり、ｊは虚数単位である。さらにg（r）は窓関数である。窓関数としては、例えば、ハニング窓【０００７】
【数３】
0004406428.tif 000004
などのg（0）にパワーの中心を持つ窓関数を用いる。この場合、Xq（f,τ）は時刻t=τを中心とする混合信号xq（t）の周波数特性を表現する。なお、Xq（f,τ）はＬサンプルにわたる情報を含んでいるため、すべてのτに対してXq（f,τ）を求める必要はなく、適当な間隔のτごとにXq（f,τ）を求める。
周波数領域で処理を行うと、式（１）で示される時間領域での畳み込み混合が、【０００８】
【数４】
0004406428.tif 000005
と各周波数での単純混合に近似表現でき、分離の操作が単純になる。ここで、Hqk（f）は源信号ｋからセンサｑまでの周波数応答であり、Sk（f,τ）は式（２）と同様な式に従って源信号sk（t）に短時間離散フーリエ変換を施したものである。ベクトルを用いて式（３）を表記すると、【０００９】
【数５】
0004406428.tif 000006
となる。ここで、X（f,τ）=[X1（f,τ）,…,XM（f,τ）]Tは混合信号ベクトル、Hk（f）=[H1k（f）,…,HMk（f）]Tは信号源ｋから各センサへの周波数応答をまとめたベクトルである。なお、[*]Tは[*]の転置ベクトルを示す。
［独立成分分析による信号分離］
ブラインド信号分離手法の１つとして、独立成分分析（ICA:Independent Component Analysis） による信号分離がある。このＩＣＡによる手法では、混合信号ベクトルX（f,τ）のみから、Ｎ行Ｍ列の分離行列W（f）及び分離信号ベクトル
Y（f,τ）=W（f）X（f,τ）  …（5）
を算出する。ここで、分離行列W（f）は、分離信号ベクトルY（f,τ）=[Y1（f,τ）,...,YN（f,τ）]Tの各要素（分離信号）Y1（f,τ）,...,YN（f,τ）が互いに独立になるように算出される。そのためのアルゴリズムには、非特許文献４に記載されているものなどがある。【００１０】
ＩＣＡでは信号の独立性に着目して分離を行うため、得られる分離信号Y1（f,τ）,...,YN（f,τ）には、順序の任意性がある。信号の順序が入れ替わってもそれらの独立性は保たれるからである。この順序の任意性の問題は、パーミュテーション（permutation）の問題と呼ばれ、周波数領域での信号分離において非常に重要な問題である。そして、このパーミュテーションの問題は、同じ源信号Sk（f,τ）に対応する分離信号Yp（f,τ）の添字ｐが、すべての周波数ｆで同じになるように解決されなければならない。
このパーミュテーション問題を解決する従来手法として、例えば非特許文献５に示されるものがある。この手法では、選択された２つのセンサ（センサペア）の位置を基準として信号源の位置に関する情報（方向と距離比）を推定する。そして、複数のセンサペアでの推定値を統合することで、より詳細な位置情報が得られる。そして、これらの位置情報としての推定値をクラスタリングし、同じクラスタに属する推定値は同じ信号源に対応するものとみなしてパーミュテーション問題を解決する。【００１１】
［時間周波数マスクによる信号分離］
また、他のブラインド信号分離手法として、時間周波数マスクによる方法がある。この手法は、信号源の数Ｎとセンサ数ＭがＭ＜Ｎの関係にある場合でも有効な信号分離抽出手法である。
この手法では信号のスパース性を仮定する。スパースとは、信号が殆どの離散時間τにおいて０であることを指す。信号のスパース性は、例えば、周波数領域での音声信号で確認される。信号のスパース性と相互独立性を仮定することで、複数の信号が同時に存在していても、各時間周波数ポイント（f,τ）では互いに重なって観測される確率が低いことを仮定できる。よって、各時間周波数ポイント（f,τ）の各センサにおける混合信号は、その時間周波数ポイント（f,τ）でアクティブな１つの信号sp（f,τ）のみから成ると仮定できる。よって、混合信号ベクトルを適当な特徴量によりクラスタリングし、各クラスタCｋのメンバの時間周波数（f,τ）に対応する混合信号Ｘ（f,τ）を抽出する時間周波数マスクMk（f,τ）を推定し、
Yk（f,τ）=Mk（f,τ）XQ'（f,τ）
により、各信号を分離抽出する。ここで、XQ'（f,τ）は混合信号のうちの１つであり、Ｑ’∈{1,…,M}である。【００１２】
クラスタリングに用いる特徴量としては、例えば、２つのセンサ（センサｑと基準センサＱ〔なお、Ｑを基準値と呼び、基準値Ｑに対応するセンサを基準センサＱと表記する。〕）における混合信号の位相差【００１３】
【数６】
0004406428.tif 000007
から計算される信号の推定到来方向（Direction of Arrival : DOA）
【００１４】
【数７】
0004406428.tif 000008
を例示できる（例えば、非特許文献３参照）。なお、ｄはセンサｑと基準センサＱとの距離であり、ｃは信号速度である。また、クラスタリングには、k-means法（例えば、非特許文献６参照）等を用いることができる。また、時間周波数マスクMk（f,τ）としては、例えばそれぞれのクラスタCｋに属するメンバの平均値θ１〜，θ2〜，…，θN〜を求め、
【００１５】
【数８】
0004406428.tif 000009
のようにして生成したものを用いることができる。ここでΔは信号を抽出する範囲を与える。この方法では、Δを小さくすると、よい分離抽出性能が得られるが非線形型歪みは大きくなる。また、Δを大きくすると、非線形型歪みは減少するが分離性能が劣化する。
その他、クラスタリングの特徴量として、２つのセンサ（センサｑと基準センサＱ）における混合信号の位相差（式（８））や両者のゲイン比【００１６】
【数９】
0004406428.tif 000010
を用いてもよい。
【非特許文献１】B. D. Van Veen and K. M. Buckley, "Beamforming: a versatile approach to spacial filtering," IEEE ASSP Magazine, pp. 4-24, April 1988
【非特許文献２】S. Haykin, eds, "Unsupervised Adaptive Filtering," John-Wiley & Sons, 2000, ISBN 0-471-29412-8
【非特許文献３】S. Araki, S. Makino, A. Blin, R. Mukai, and H. Sawada, "Underdetermined blind separation for speech in real environments with sparseness and ICA," in Proc. ICASSP 2004, vol. III, May 2004, pp. 881-884
【非特許文献４】A. Hyvarinen and J. Karhunen and E. Oja, "Independent Component Analysis," John Wiley & Sons, 2001, ISBN 0-471-40540
【非特許文献５】R. Mukai, H. Sawada, S. Araki and S. Makino, "Frequency Domain Blind Source Separation using Small and Large Spacing Sensor Pairs," in Proc. of ISCAS 2004, vol. V, pp.1-4, May 2004.
【非特許文献６】R. O. Duda, P. E. Hart, and D. G. Stork, Pattern Classification, Wiley Interscience, 2nd edition, 2000
【発明の開示】
【発明が解決しようとする課題】
【００１７】
しかし、上述した従来の技術では、複数のセンサによる観測信号から得られた情報を容易かつ効率的に利用して信号分離処理を行うことができないという問題点があった。
例えば、独立成分分析による信号分離の場合、正確にパーミュテーション問題を解決するためには繁雑な操作が必要となるという問題がある。すなわち、従来のパーミュテーション問題の解決手法では、センサペア毎に方向と距離比とが推定される。従って、正確にパーミュテーション問題を解決するためには複数のセンサペアでの推定値を統合する必要があった。また、この推定値には誤差が含まれているが、より誤差が少ないと思われるセンサペアを優先的に用いたり、誤差をうまく吸収できるように統合方法を工夫したりする必要があった。さらに、この従来手法では、信号源の位置に関する情報を推定する必要性から、センサの位置情報をあらかじめ取得しておかなければならないという問題点もあった。これはセンサを不規則に配置する場合に不利である。また、たとえセンサを規則的に配置するとしても、配置情報を厳密に知ることは困難であり、より正確なパーミュテーション問題の解決のためにはキャリブレーションなどの操作を必要とする。【００１８】
また、従来の時間周波数マスクによる信号分離手法の場合、センサが２個の場合のための手法しか提案されておらず、センサが３個以上ある場合であっても、そのうち特定の２つのセンサｑ，Ｑのみの情報しか用いずに特徴量を計算していた。これは次元の低下を意味し、全てのセンサを用いる場合に比べて情報量が欠落する。そのため、全てのセンサの情報を効率的に使えず性能に限度があった。また、全てのセンサの情報を有効利用するには、例えば、非特許文献５にならって、複数のセンサペアで求めた特徴量を統合することも可能であろうが、この統合のためには、特徴量抽出のための更なる処理が必要であり、より統合誤差が少ないと思われるセンサペアを優先的に用いるなど、統合の際に工夫する必要があると思われる。さらに、この手法でも、センサの厳密な配置情報をあらかじめ知っておかなければならないという問題点がある。これは、センサを自由に配置する場合に不利である。また、たとえ規則的に配置するとしても、配置情報を厳密に知ることは困難であり、より正確な信号抽出のためにはキャリブレーションなどの操作を必要とする。【００１９】
さらに、ブラインド信号分離の処理の基本は、センサで観測された混合信号を分離して、複数の分離信号を出力することである。しかし、すべての分離信号が重要なわけではなく、一部の分離信号のみに目的とする目的信号が含まれていることもある。そのような場合、目的信号が含まれている分離信号を選択する必要があるが、従来のブラインド信号分離では、どの分離信号に目的信号が含まれているかという情報までは提供されない。従って、何らかの別の手段によって、どの分離信号に目的信号が含まれているかを判定しなければならない。【００２０】
本発明はこのような点に鑑みてなされたものであり、複数のセンサによる観測信号から得られた情報を容易かつ効率的に利用して信号分離処理を行うことが可能な技術を提供することを目的とする。
【課題を解決するための手段】
【００２１】
本発明では、上記課題を解決するために、まず、周波数領域変換部が、複数のセンサで観測された混合信号を、周波数領域の混合信号に変換する。次に、正規化部が、当該周波数領域の混合信号を用いて生成された複素ベクトルを正規化し、当該複素ベクトルの周波数依存性を排除した正規化ベクトルを生成する。次に、クラスタリング部が、当該正規化ベクトルをクラスタリングし、クラスタを生成する。そして、このクラスタを用いて信号分離処理を行う。
ここで、当該クラスタの生成には、混合信号を観測する各センサの厳密な位置情報を、入力情報として直接用いることは必要としない。また、各クラスタは、各信号源の位置に依存する情報に基づき形成される。これらより、本発明では、各センサの厳密な位置情報を用いることなく、信号分離処理が可能となる。【００２２】
また、本発明において好ましくは、正規化部は、複素ベクトルが具備する特定の１つの要素を基準として当該複素ベクトルの各要素の偏角を正規化する第１正規化部と、第１正規化部で正規化された各要素の偏角を周波数に比例した値で除算する第２正規化部とを有する。
これらの正規化後の複素ベクトルは、信号源の位置に依存したクラスタを形成する。これにより、各センサの厳密な位置情報を用いることなく、信号分離処理を行うことが可能となる。【００２３】
また、本発明において好ましくは、正規化部は、第２正規化部で正規化された各要素からなるベクトルのノルムを所定の値に正規化する第３正規化部をさらに有する。
これらの正規化後の複素ベクトルは、信号源の位置に依存したクラスタを形成する。また、第２正規化された各要素からなるベクトルのノルムを正規化することにより、クラスタリングの処理が簡略化される。
また、好ましい第１の本発明では、まず、周波数領域変換部が、複数のセンサで観測された混合信号を、周波数領域の混合信号に変換する。次に、分離行列算出部が、周波数領域の混合信号を用いて、周波数毎に分離行列を算出し、逆行列算出部が、分離行列の一般化逆行列を算出する。次に、基底ベクトル正規化部が、一般化逆行列を構成する基底ベクトルの正規化を行い、正規化基底ベクトルを算出する。次に、クラスタリング部が、正規化基底ベクトルをクラスタリングし、クラスタを生成する。その後、順列算出部が、クラスタの中心ベクトルと正規化基底ベクトルとを用い、分離行列の要素を並び替えるための順列を算出する。なお、基底ベクトルは、複素ベクトルに含まれる概念である。【００２４】
この第１の本発明では、基底ベクトルを正規化してクラスタリングし、パーミュテーション問題解決のための順列を算出する。そのため、このクラスタリングに際し、センサの位置情報をあらかじめ取得しておく必要がない。さらに、この好ましい本発明では、正規化基底ベクトルのすべての要素を対象としてクラスタリングを行い、その結果からパーミュテーション問題解決のための順列を算出する。そのため、従来のように推定結果を統合するための操作を必要としない。
また、当該第１の本発明において、より好ましくは、基底ベクトル正規化部は、基底ベクトルの周波数依存性を排除する正規化を行う。そして、より好ましくは、当該基底ベクトルの周波数依存性を排除する正規化は、基底ベクトルが具備する特定の１つの要素を基準として当該基底ベクトルの各要素の偏角を正規化し、さらに、各要素の偏角を周波数に比例した値で除算する正規化である。このような正規化により、信号源の位置に依存したクラスタの形成が可能となる。【００２５】
また、当該第１の本発明において、より好ましくは、基底ベクトルの周波数依存性を排除する正規化は、基底ベクトルAp（f）（ｐ＝１，...，Ｎ、Ｎは信号源の数）の各要素Aqp（f）（ｑ＝１，...，Ｍ、Ｍは混合信号を観測するセンサの数）に対し、【００２６】
【数１０】
0004406428.tif 000011
の演算を行う正規化である。ただし、ｅｘｐをネピア数とし、ａｒｇ［・］を偏角とし、ｆを周波数とし、ｊを虚数単位とし、ｃを信号の伝達速度とし、ＱをＭ以下の自然数から選択された基準値とし、ｄを実数とする。すなわち、式（１０）の演算によって行われる正規化は、基底ベクトルが具備する特定の１つの要素を基準として当該基底ベクトルの各要素の偏角を正規化し、さらに、各要素の偏角を周波数に比例した値で除算する正規化である。この正規化により、周波数への依存性を排除できる。また、当該正規化には、センサの厳密な配置情報を必要としない。【００２７】
また、上述の式（１０）の実数ｄは、要素AQp（f）に対応する基準センサＱと他のセンサとの最大距離ｄmaxであることが望ましい。クラスタリングの精度が一般的に向上するからである。この詳細は後述する。
また、第１の本発明において好ましくは、基底ベクトルを周波数に依存しない周波数正規化ベクトルに正規化し、この周波数正規化ベクトルをノルムが規定値をとる正規化基底ベクトルに正規化する。以上の二段階の正規化によって生成された正規化基底ベクトルは、周波数に依存せず、信号源の位置のみに依存する。なお、ノルムの正規化により、クラスタリングの処理が簡略化できる。【００２８】
さらにまた、第１の本発明において好ましくは、分離信号のエンベロープ（分離信号の絶対値の包絡線）とクラスタの中心ベクトルと正規化基底ベクトルとを用いて順列を算出する。これにより、より精度よくパーミュテーション問題を解決することができる。
また、好ましい第２の本発明では、周波数領域変換部が、複数のセンサで観測された混合信号を周波数領域の混合信号に変換し、信号分離部が、周波数領域の混合信号を用い、周波数毎に分離行列と分離信号とを算出する。そして、目的信号選択部が、複数の分離信号から、目的信号を含むものを選択信号として選択する。その手順は、分離行列の一般化逆行列の列である基底ベクトルを正規化し、その正規化された基底ベクトルをクラスタリングし、そのクラスタの分散を指標として選択信号を決定するものである。なお、分離行列が正方行列の場合、「一般化逆行列」は逆行列に相当する。すなわち、ここでの一般化逆行列は、通常の逆行列をも含む概念である。【００２９】
クラスタの分散を指標とすることで、センサに近い信号を目的信号として位置付け、そのような目的信号を含む分離信号を選択信号として選び出すことができる。その理由を以下に述べる。ここで、基底ベクトルの正規化は、複数の信号源から発せられた信号の畳み込み混合を所定のモデル（例えば、近距離場モデル）に近似した場合に、正規化された基底ベクトルが信号源の位置のみに依存したクラスタを形成するように行われる。しかし、実環境では、このようなモデルに反映されない様々な因子が存在する。例えば、近距離場モデルでは、信号が壁などに反射することによって生じる伝送歪みが反映されない。このような実環境とモデルとの乖離は、信号源からセンサまでの距離が大きいほど大きくなり、センサに近い信号ほど小さくなる。従って、センサに近い信号ほど実環境に近い条件で正規化ができ、実環境とモデルとの乖離に起因するクラスタの分散を小さくできる。好ましい第２の本発明では、この関係に着目し、クラスタの分散を指標として、センサに近い信号である目的信号を含む選択信号を抽出する。以上の操作により、目的信号を抽出し、他の妨害信号をある程度抑圧することができる。
【００３０】
しかし、分離行列及び分離信号の算出を独立成分分析（ICA：Independent Component Analysis）で行う場合、以上の処理によって完全に抑制できる妨害信号の数はセンサの数から１を引いた数までである。すなわち、妨害信号の数がそれ以上の場合は抑圧しきれない残留成分が残る。そのため、本発明において好ましくは、更に、マスク生成部が、周波数領域の混合信号と基底ベクトルとを用いて時間周波数マスクを生成し、マスキング部が、この時間周波数マスクを、選択された選択信号に適用する。これにより、信号源の数がセンサ数より多い場合でも、選択信号に残留した妨害信号をより良く抑圧することができる。【００３１】
更に、第２の本発明において好ましくは、マスク生成部は、周波数領域の混合信号を用いて白色化行列を生成し、白色化行列を用い、周波数領域の混合信号を要素とする混合信号ベクトルを変換した白色化混合信号ベクトルと、基底ベクトルを変換した白色化基底ベクトルとを算出し、白色化混合信号ベクトルと白色化基底ベクトルとがなす角度を時間周波数毎に算出し、この角度を要素とする関数を用いて時間周波数マスクを生成する。この時間周波数マスクを選択信号に適用することにより、選択信号に残留した妨害信号を更に抑圧できる。【００３２】
また更に、第２の本発明において好ましくは、白色化行列が、周波数をfとし、離散時間をτとし、混合信号ベクトルをX（f,τ）とし、ベクトル＊を時間平均したベクトルを<＊>τとし、ベクトル＊の複素共役転置ベクトル（各要素の複素共役を採り転置したベクトル）を＊Hとし、R（f）=<X（f,τ）・X（f,τ）H>τとした場合におけるV（f）=R（f）−1/2である。そして、Z（f,τ）=V（f）・X（f,τ）の演算により、白色化混合信号ベクトルZ（f,τ）を算出し、基底ベクトルをA（f）とした場合におけるB（f）=V（f）・A（f）の演算により、白色化基底ベクトルB（f）を算出する。また、ベクトル＊の絶対値を｜＊｜とし、ベクトル＊のノルムを‖＊‖とした場合におけるθ（f,τ）=cos−1（|BH（f）・Z（f,τ）|／‖B（f）‖・‖Z（f,τ）‖）の演算により、角度θ（f,τ）を算出し、α，g，θTを実数とした場合におけるロジスティック関数M（θ（f,τ））=α／（1+eg・（θ（f,τ）−θT））を時間周波数マスクとして算出する。この時間周波数マスクを、抽出した選択信号に適用することにより、選択信号に残留した妨害信号を更に抑圧できる。【００３３】
また、第２の本発明において好ましくは、目的信号選択部は、基底ベクトルの周波数依存性を排除する正規化を行う。第２の本発明において、より好ましくは、基底ベクトルの周波数依存性を排除する正規化は、基底ベクトルが具備する特定の１つの要素を基準として当該基底ベクトルの各要素の偏角を正規化し、さらに、各要素の偏角を周波数に比例した値で除算する正規化である。またさらに、第２の本発明において、好ましくは、基底ベクトルの周波数依存性を排除する正規化は、基底ベクトルAp（f）（ｐは自然数）の各要素Aqp（f）（ｑ＝１，...，Ｍ、Ｍは混合信号を観測するセンサの数）に対し、
【００３４】
【数１１】
0004406428.tif 000012
の演算を行う正規化である。ただし、ｅｘｐをネピア数とし、ａｒｇ［・］を偏角とし、ｆを周波数とし、ｊを虚数単位とし、ｃを信号の伝達速度とし、ＱをＭ以下の自然数から選択された基準値とし、ｄを実数とする。このような正規化により、複数の信号源から発せられた信号の畳み込み混合を所定のモデルに近似した場合に、正規化された基底ベクトルが信号源の位置に依存したクラスタを形成することとなる。その結果、上述のようにクラスタの分散の大きさを指標として、目的信号が含まれる分離信号を選択することが可能となる。なお、当該正規化には、センサの厳密な配置情報を必要としない。【００３５】
またさらに、上述の式（１１）の実数ｄは、基準センサＱと他のセンサとの最大距離ｄmaxであることが望ましい。クラスタリングの精度が一般的に向上するからである。この詳細は後述する。
また、第２の本発明において好ましくは、目的信号選択部は、分散が最小となるクラスタを選択し、選択されたクラスタに対応する分離信号を目的信号が含まれるものとして選択する。これにより、モデルからの乖離が最も少ない信号（例えば、センサに最も近い信号）を目的信号として抽出することができる。【００３６】
また、好ましい第３の本発明では、まず、周波数領域変換部が、複数のセンサで観測された混合信号を周波数領域の混合信号に変換する。次に、ベクトル正規化部が、周波数領域の混合信号から構成される混合信号ベクトルの正規化を行い、正規化ベクトルを算出する。そして、クラスタリング部が、正規化ベクトルをクラスタリングし、クラスタを生成する。その後、分離信号生成部が、ｋ番目のクラスタに属する正規化ベクトルの時間周波数に対応する混合信号ベクトルから所定番目の要素を抽出し、これをｋ番目の要素とした分離信号ベクトルを生成する。【００３７】
ここで、第３の本発明では、全てのセンサで観測された混合信号を正規化してクラスタリングを行い、各クラスタの情報を用いて分離信号ベクトルを生成している。これは、全てのセンサの情報を同時に用いて分離信号を抽出していることを意味する。また、この処理では、センサの厳密な配置情報を必要としない。以上より、第３の本発明では、センサの厳密な配置情報を必要とすることなく、全ての観測信号から得られた情報を容易かつ効率的に利用して信号分離を行うことができる。
また、第３の本発明において、好ましくは、ベクトル正規化部は、周波数領域の混合信号から構成される混合信号ベクトルの周波数依存性を排除する正規化を行う。また、さらに好ましくは、混合信号ベクトルの周波数依存性を排除する正規化は、混合信号ベクトルが具備する特定の１つの要素を基準として当該基底ベクトルの各要素の偏角を正規化し、さらに、各要素の偏角を周波数に比例した値で除算する正規化である。またさらに好ましくは、混合信号ベクトルの周波数依存性を排除する正規化は、混合信号ベクトルの各要素Xq（f,τ） （ｑ＝１，...，Ｍ、Ｍは混合信号を観測するセンサの数）に対し、【００３８】
【数１２】
0004406428.tif 000013
の演算を行う正規化である。ただし、ｅｘｐをネピア数とし、ａｒｇ［・］を偏角とし、ｊを虚数単位とし、ｃを信号の伝達速度とし、ＱをＭ以下の自然数から選択された値とし、ｄを実数とし、ｆを周波数とし、τを離散時間とする。これにより、周波数の依存性を排除できる。これにより、信号源の位置に依存したクラスタを形成できる。なお、当該正規化には、センサの厳密な配置情報を必要としない。【００３９】
また、上述の式（１２）の実数ｄは、要素ＸＱ（ｆ，τ）に対応するセンサと他のセンサとの最大距離ｄmaxであることが望ましい。クラスタリングの精度が一般的に向上するからである。この詳細は後述する。
また、第３の本発明において、好ましくは、ベクトル正規化部は、混合信号ベクトルの周波数依存性を排除する正規化と、ノルムを規定値にする正規化とを行う。これにより、クラスタリングの処理が簡略化される。
【発明の効果】
【００４０】
以上のように、本発明では、複数のセンサによる観測信号から得られた情報を容易かつ効率的に利用して信号分離処理を行うことができる。
例えば、第１の発明では、厳密なセンサ配置に係る情報の事前取得や煩雑な操作を必要とせずに正確にパーミュテーション問題を解決できる。また、第２の発明では、目的信号の方向に関する情報を知ることなく、複数の信号源から発せられた信号が混合された混合信号から目的信号を抽出することができる（たとえＮ＞Ｍであったとしても）。また、第３の発明では、センサの厳密な配置情報を必要とすることなく、全ての観測信号から得られた情報を容易かつ効率的に利用して信号分離を行うことができる（たとえＮ＞Ｍであったとしても）。
【図面の簡単な説明】
【００４１】
【図１】図１は、本発明の原理を備えた信号分離装置の機能構成を例示したブロック図である。
【図２】図２は、第１の実施の形態における信号分離装置のハードウェア構成を例示したブロック図である。
【図３】図３は、第１の実施の形態における信号分離装置のブロック図の例示である。
【図４】図４Ａは、図３におけるパーミュテーション問題解決部の詳細を例示したブロック図である。図４Ｂは、図４Ａの基底ベクトル正規化部の詳細を例示したブロック図である。
【図５】図５は、第１の実施の形態における信号分離装置の処理の全体を説明するためのフローチャートである。
【図６】図６は、第１の実施の形態の正規化処理の詳細を説明するためのフローチャートである。
【図７】図７Ａは、ｄmax／２≧dである場合における、パラメータｄ毎の正規化基底ベクトルの要素Aqp''（f）とその偏角arg[Aqp''（f）]との関係を説明するための複素平面図である。図７Ｂは、ｄmax／２＜ｄ＜ｄmaxである場合における、パラメータｄ毎の正規化基底ベクトルの要素Aqp''（f）とその偏角arg[Aqp''（f）]との関係を説明するための複素平面図である。
【図８】図８Ａは、ｄ＝ｄmaxである場合における、パラメータｄ毎の正規化基底ベクトルの要素Aqp''（f）とその偏角arg[Aqp''（f）]との関係を説明するための複素平面図である。図８Ｂは、ｄ＞ｄmaxである場合における、パラメータｄ毎の正規化基底ベクトルの要素Aqp''（f）とその偏角arg[Aqp''（f）]との関係を説明するための複素平面図である。
【図９】図９は、第２の実施の形態における信号分離装置のブロック図の例示である。
【図１０】図１０Ａは、図９におけるパーミュテーション問題解決部の詳細を例示したブロック図である。図１０Ｂは、図１０Ａの順列修正部の詳細を例示したブロック図である。
【図１１】図１１は、第２の実施の形態における信号分離装置の処理の全体を説明するためのフローチャートである。
【図１２】図１２は、図１１におけるステップＳ５８の例を説明するためのフローチャートである。
【図１３】図１３は、図１１におけるステップＳ５８の例を説明するためのフローチャートである。
【図１４】図１４Ａは、第１の実施の形態及び第２の実施の形態による音源分離実験条件を示した図である。図１４Ｂは、第１の実施の形態及び第２の実施の形態による音源分離実験結果を示した図である。
【図１５】図１５Ａは、第１の実施の形態及び第２の実施の形態による音源分離実験条件を示した図である。図１５Ｂは、第１の実施の形態及び第２の実施の形態による音源分離実験結果を示した図である。
【図１６】図１６は、第３の実施の形態における信号分離装置の構成を例示したブロック図である。
【図１７】図１７Ａは、図１６における目的信号選択部の詳細構成を例示したブロック図である。図１７Ｂは、図１７Ａにおける基底ベクトルクラスタリング部の詳細構成を例示したブロック図である。
【図１８】図１８Ａは、図１６における時間周波数マスキング部の詳細構成を例示したブロック図である。図１８Ｂは、図１８Ａにおけるマスク生成部の詳細を例示したブロック図である。
【図１９】図１９は、第３の実施の形態における信号分離処理の全体を説明するためのフローチャートである。
【図２０】図２０は、第３の実施の形態における目的信号選択部の処理の詳細を説明するためのフローチャートである。
【図２１】図２１Ａは、ステップＳ１１２の周波数正規化の詳細を説明するためのフローチャートである。図２１Ｂは、ステップＳ１１３のノルム正規化の詳細を説明するためのフローチャートである。
【図２２】図２２は、選択信号の選択手順（ステップＳ１１５）の詳細を例示したフローチャートである。
【図２３】図２３は、図１９におけるステップＳ１０４の詳細を説明するためのフローチャートである。
【図２４】図２４Ａは、二種類の実数パラメータθＴ，ｇに対し、式（４６）に従って算出された時間周波数マスクM（f,τ）を例示した図である。図２４Ｂは、ある時間周波数位置（f,τ）において目的信号（I（f）=1とする）に対応するベクトルV（f）・H1（f）の他に、妨害信号に対応ベクトルするV（f）・H2（f），V（f）・H3（f）が共存することを示した図である。
【図２５】図２５は、第４の実施の形態における信号分離装置のブロック図の例示である。
【図２６】図２６は、第４の実施の形態における信号分離装置の処理を説明するためのフローチャートである。
【図２７】図２７は、第５の実施の形態における信号分離装置のブロック図の例示である。
【図２８】図２８Ａは、図２７における時間周波数マスキング部の詳細構成を示すブロック図である。図２８Ｂは、図２８Ａのマスク生成部の詳細構成を示すブロック図である。
【図２９】図２９は、第５の実施の形態における時間周波数マスクの生成処理を説明するためのフローチャートである。
【図３０】図３０Ａは、図２９におけるステップＳ１７１の詳細を説明するためのフローチャートである。図３０Ｂは、図２９におけるステップＳ１７２の詳細を説明するためのフローチャートである。
【図３１】図３１Ａは、第３の実施の形態及び第４の実施の形態による効果を示すための実験条件を示した図である。図３１Ｂは、ICAのみの場合（第４の実施の形態）と、ICAと時間周波数マスキングの両方を組み合わせた場合（第３の実施の形態）のSIRの平均改善量を示す表である。
【図３２】図３２は、第６の実施の形態における信号分離装置のブロック図の例示である。
【図３３】図３３は、図３２における信号分離部の詳細を例示したブロック図である。
【図３４】図３４は、第６の実施の形態における信号分離装置の処理の全体を説明するためのフローチャートである。
【図３５】図３５Ａは、図３４に示したステップＳ２０２の処理の詳細を説明するためのフローチャートである。図３５Ｂは、図３４に示したステップＳ２０３の処理の詳細を説明するためのフローチャートである。
【図３６】図３６は、図３４に示したステップＳ２０５の処理の詳細を説明するためのフローチャートである。
【図３７】図３７Ａは、ｄmax／２≧dの場合における、パラメータｄ毎のノルム正規化ベクトルX''（f,τ）の要素Xq''（f,τ）と、その偏角arg[Xq'' （f,τ）]との関係を説明するための複素平面図である。図３７Ｂは、ｄmax／２＜ｄ＜ｄmaxの場合における、パラメータｄ毎のノルム正規化ベクトルX''（f,τ）の要素Xq'' （f,τ）と、その偏角arg[Xq''（f,τ）]との関係を説明するための複素平面図である。
【図３８】図３８Ａは、ｄ＝ｄmaxの場合における、パラメータｄ毎のノルム正規化ベクトルX''（f,τ）の要素Xq''（f,τ）と、その偏角arg[Xq''（f,τ）]との関係を説明するための複素平面図である。図３８Ｂは、ｄ＞ｄmaxの場合における、パラメータｄ毎のノルム正規化ベクトルX''（f,τ）の要素Xq''（f,τ）と、その偏角arg[Xq''（f,τ）]との関係を説明するための複素平面図である。
【図３９】図３９Ａは、第６の実施の形態による音源分離実験条件を示した図である。図３９Ｂは、第６の実施の形態による音源分離実験結果を示した図である。
【図４０】図４０Ａは、第６の実施の形態による音源分離実験条件を示した図である。図４０Ｂは、第６の実施の形態による音源分離実験結果を示した図である。
【図４１】図４１Ａは、第６の実施の形態による音源分離実験条件を示した図である。図４１Ｂは、第６の実施の形態による音源分離実験結果を示した図である。
【符号の説明】
【００４２】
１，１０，２００，１００１，１２００，１３００，２００１  信号分離装置
【発明を実施するための最良の形態】
【００４３】
以下、本発明の実施の形態を図面を参照して説明する。
〔原理〕
まず、本発明の原理について説明する。
図１は、本発明の原理を備えた信号分離装置１の機能構成を例示したブロック図である。なお、後述のように、信号分離装置１は、例えば、公知のノイマン型のコンピュータに所定のプログラムを実行させることにより構成される。
信号分離装置１は、複数の信号源から発せられた源信号の混合からなる混合信号を当該源信号に分離する装置である。図１に例示するように、信号分離装置１は、周波数領域変換部２と、複素ベクトル生成部３と、正規化部４と、クラスタリング部５とを有している。ここで、正規化部４は、複素ベクトルが具備する特定の１つの要素を基準として当該複素ベクトルの各要素の偏角を正規化する第１正規化部４ａと、第１正規化部４ａで正規化された各要素の偏角を周波数に比例した値で除算する第２正規化部４ｂと、第２正規化部４ｂで正規化された各要素からなるベクトルのノルムを所定の値に正規化する第３正規化部４ｃとを有する。なお、第１正規化部４ａ及び第２正規化部４ｂにより、複素ベクトルの周波数依存性が排除される（周波数正規化）。【００４４】
信号分離装置１によって信号分離処理を行う場合、まず、複数のセンサで観測された混合信号（時間領域の信号）が周波数領域変換部２に入力される。周波数領域変換部２は、短時間離散フーリエ変換等によって、複数のセンサで観測された混合信号（時間領域の信号）を、周波数領域の混合信号に変換する。次に、複素ベクトル生成部３が、当該周波数領域の混合信号を用い、複素数の要素からなる複素ベクトルを生成する。次に、正規化部４が、当該複素ベクトルを正規化し、当該複素ベクトルの周波数依存性を排除した正規化ベクトルを生成する。【００４５】
図１の例の正規化では、まず、第１正規化部４ａが、時間周波数毎に、当該複素ベクトルが具備する特定の１つの要素を基準として当該複素ベクトルの各要素の偏角を正規化する。これにより、複素ベクトルの各要素の偏角は、源信号の位相や振幅に依存せず、各センサに対する信号源の相対位置や周波数のみに依存することになる（詳細は後述）。次に、第２正規化部４ｂが、第１正規化部４ａで正規化された各要素の偏角を周波数に比例した値で除算する。これにより、各複素ベクトルの各要素の周波数依存性が排除され、複素ベクトルは、各センサに対する各信号源の相対位置のみに依存したものに正規化される。さらに、第３正規化部４ｃが、第２正規化部４ｂで正規化された各要素からなるベクトルのノルムを所定の値に正規化する。【００４６】
次に、クラスタリング部５が、このような正規化が行われた正規化ベクトルのクラスタリングを行い、クラスタを生成する。これらのクラスタは、各センサに対する各信号源の相対位置のみに依存したものになる。分離信号生成部６は、これらのクラスタを利用して各種の信号分離処理を行い、周波数領域の分離信号を生成する。最後に、時間領域変換部が、周波数領域の分離信号を時間領域の分離信号に変換する。
上述のように、これらのクラスタを生成するために、各センサの厳密な位置情報を事前に取得しておく必要はない。また、これらのクラスタの生成には、全てのセンサでの観測信号の情報が用いられている。すなわち、本発明では、複数のセンサによる観測信号から得られた情報を容易かつ効率的に利用して信号分離処理を行うことができる。【００４７】
なお、ノルムを正規化しなくても、クラスタリングの方法を工夫することにより、各センサに対する各信号源の相対位置のみに依存したクラスタを生成することは可能である。しかし、クラスタリングの処理を簡略化するため、第３正規化部４ｃによるノルムの正規化を行うことが望ましい。
以下、本発明の各実施の形態を詳細に説明する。
〔第１の実施の形態（「第１の本発明」の例）〕
まず、本発明における第１の実施の形態について説明する。【００４８】
本形態は、上述の原理により、厳密なセンサ配置に係る情報の事前取得や煩雑な操作を必要とせずに正確にパーミュテーション問題を解決する形態である。なお、本形態では、後述する「基底ベクトル」が上述の「複素ベクトル」に相当する。
＜ハードウェア構成＞
図２は、第１の実施の形態における信号分離装置１０のハードウェア構成を例示したブロック図である。
図２に例示するように、この例の信号分離装置１０は、ＣＰＵ（Central Processing Unit）１０ａ、入力部１０ｂ、出力部１０ｃ、補助記憶装置１０ｆ、ＲＡＭ（Random Access Memory）１０ｄ、ＲＯＭ（Read Only Memory）１０ｅ及びバス１０ｇを有している。【００４９】
この例のＣＰＵ１０ａは、制御部１０ａａ、演算部１０ａｂ及びレジスタ１０ａｃ有し、レジスタ１０ａｃに読み込まれた各種プログラムに従って様々な演算処理を実行する。また、この例の入力部１０ｂは、データが入力される入力ポート、キーボード、マウス等であり、出力部１０ｃは、データを出力する出力ポート、ディスプレイ等である。補助記憶装置１０ｆは、例えば、ハードディスク、ＭＯ（Magneto-Optical disc）、半導体メモリ等であり、本形態の信号分離処理を実行するための信号分離プログラムを格納した信号分離プログラム領域１０ｆａ及びセンサで観測された時間領域の混合信号等の各種データが格納されるデータ領域１０ｆｂを有している。また、ＲＡＭ１０ｄは、例えば、ＳＲＡＭ （Static Random Access Memory）、ＤＲＡＭ （Dynamic Random Access Memory）等であり、信号分離プログラムが書き込まれる信号分離プログラム領域１０ｄａ及び各種データが書き込まれるデータ領域１０ｄｂを有している。また、この例のバス１０ｇは、ＣＰＵ１０ａ、入力部１０ｂ、出力部１０ｃ、補助記憶装置１０ｆ、ＲＡＭ１０ｄ及びＲＯＭ１０ｅを通信可能に接続している。【００５０】
＜ハードウェアとソフトウェアとの協働＞
この例のＣＰＵ１０ａは、読み込まれたＯＳ（Operating System）プログラムに従い、補助記憶装置１０ｆの信号分離プログラム領域１０ｆａに格納されている信号分離プログラムを、ＲＡＭ１０ｄの信号分離プログラム領域１０ｄａに書き込む。同様にＣＰＵ１０ａは、補助記憶装置１０ｆのデータ領域１０ｆｂに格納されている時間領域の混合信号等の各種データをＲＡＭ１０ｄのデータ領域１０ｄｂに書き込む。さらに、ＣＰＵ１０ａは、この信号分離プログラムや各種データが書き込まれたＲＡＭ１０ｄ上のアドレスをレジスタ１０ａｃに格納する。そして、ＣＰＵ１０ａの制御部１０ａａは、レジスタ１０ａｃに格納されたこれらのアドレスを順次読み出し、読み出したアドレスが示すＲＡＭ１０ｄ上の領域からプログラムやデータを読み出し、そのプログラムが示す演算を演算部１０ａｂに順次実行させ、その演算結果をレジスタ１０ａｃに格納していく。【００５１】
図３は、このようにＣＰＵ１０ａに信号分離プログラムが読み込まれることにより構成される信号分離装置１０のブロック図の例示である。また、図４Ａは、図３におけるパーミュテーション問題解決部１４０の詳細を例示したブロック図であり、図４Ｂは、図４Ａの基底ベクトル正規化部１４２の詳細を例示したブロック図である。
図３に例示するように、信号分離装置１０は、メモリ１００、周波数領域変換部１２０、分離行列算出部１３０、パーミュテーション問題解決部１４０、分離信号生成部１５０、時間領域変換部１６０及び制御部１７０を有している。また、この例のパーミュテーション問題解決部１４０は、逆行列算出部１４１（「複素ベクトル生成部」に対応）、基底ベクトル正規化部１４２（「正規化部」に対応）、クラスタリング部１４３、順列算出部１４４及び並び替え部１４５を有しており、基底ベクトル正規化部１４２は、周波数正規化部１４２ａ及びノルム正規化部１４２ｂを有している。また、周波数正規化部１４２ａは、第１正規化部１４２ａａと第２正規化部１４２ａｂとを有している。さらに、制御部１７０は一時メモリ１７１を有している。【００５２】
ここでメモリ１００及び一時メモリ１７１は、レジスタ１０ａｃ、補助記憶装置１０ｆのデータ領域１０ｆｂ或いはＲＡＭ１０ｄのデータ領域１０ｄｂ等に相当する。また、周波数領域変換部１２０、分離行列算出部１３０、パーミュテーション問題解決部１４０、分離信号生成部１５０、時間領域変換部１６０及び制御部１７０は、ＣＰＵ１０ａにＯＳプログラムや信号分離プログラムが読み込まれることにより構成されるものである。
なお、図３及び図４における破線の矢印は理論上の情報の流れを示し、実線の矢印は実際のデータの流れを示す。また、これらの図において制御部１７０に出入りするデータの流れに対応する矢印は省略してあり、図４における実際のデータの流れに対応する矢印も省略してある。【００５３】
＜処理＞
次に、本形態の信号分離装置１０の処理について説明する。なお、以下では、Ｎ個の源信号が混合され、Ｍ個のセンサで観測された状況を取り扱う。また、前処理において、各センサで観測された時間領域の混合信号xq（t）（q=1,...,M）がメモリ１００の記憶領域１０１に格納され、信号の伝達速度ｃ、Ｍ以下の自然数から選択された基準値Ｑ（Ｍ個のセンサから選択された１つの基準センサの添字番号）及び実数ｄの各パラメータが記憶領域１０７に格納されているものとする。【００５４】
図５は、第１の実施の形態における信号分離装置１０の処理の全体を説明するためのフローチャートである。以下、この図に沿って、本形態における信号分離装置１０の処理を説明していく。
［周波数領域変換部１２０の処理］
まず、周波数領域変換部１２０において、メモリ１００の記憶領域１０１から時間領域の混合信号xq（t）を読み出し、これらを短時間離散フーリエ変換等によって周波数毎の時系列の信号（「周波数領域の混合信号」と呼ぶ）Xq（f,τ）（q=1,...,M）に変換し、メモリ１００の記憶領域１０２に格納する（ステップＳ１）。【００５５】
［分離行列算出部１３０の処理］
次に、分離行列算出部１３０において、メモリ１００の記憶領域１０２から周波数領域の混合信号Xq（f,τ）を読み出す。周波数領域の混合信号Xq（f,τ）を読み出した分離行列算出部１３０は、これらからなる混合信号ベクトルX（f,τ）=[X1（f,τ）,...,XM（f,τ）]Tを用い、独立成分分析（ＩＣＡ）によって、周波数毎に第１の分離行列W（f）と分離信号ベクトルY（f,τ）=[Y1（f,τ）,...,YN（f,τ）]Tとを算出する。そして、算出された第１の分離行列W（f）はメモリ１００の記憶領域１０３に格納される（ステップＳ２）。【００５６】
ここで、分離行列算出部１３０において算出された第１の分離行列W（f）には、順序の任意性が含まれている。よって、以下のようにパーミュテーション問題解決部１４０において、第１の分離行列W（f）を順序の任意性を解決した第２の分離行列W’（f）に変更する。
［パーミュテーション問題解決部１４０の処理］
まず、逆行列算出部１４１において、メモリ１００の記憶領域１０３から第１の分離行列W（f）を読み出し、そのムーア・ペンローズ（Moore-Penrose）型一般化逆行列W+（f）=[A1（f）,...,AN（f）]（M=Nの場合は逆行列W-1（f）に一致）を計算し、これを構成する基底ベクトルAp（f）=[A1p（f）,...,AMp（f）]T （p=1,...,N）をメモリ１００の記憶領域１０４に格納する（ステップＳ３）。
【００５７】
次に、基底ベクトル正規化部１４２において、メモリ１００の記憶領域１０４から基底ベクトルAp（f）（p=1,...,N、f=0,fs/L,…,fs（L-1）/L）を読み出し、これらを正規化した正規化基底ベクトルAp''（f）を算出し、メモリ１００の記憶領域１０６に格納する（ステップＳ４）。なお、基底ベクトル正規化部１４２は、全ての基底ベクトルAp（f）（p=1,...,N、f=0,fs/L,…,fs（L-1）/L）を正規化し、それらを周波数に依存せず、信号源の位置のみに依存する正規化基底ベクトルAp''（f）に正規化する。その結果、それらをクラスタリングした際に、各クラスタが各信号源に対応するようになる。この正規化を適切に行わないとクラスタが形成されない。また、本形態における正規化は周波数正規化とノルム正規化との二段階からなる。周波数正規化は、周波数正規化部１４２ａ（図４Ｂ）において、基底ベクトルを周波数に依存しない周波数正規化ベクトルに正規化するものである。ノルム正規化は、ノルム正規化部１４２ｂにおいて、周波数正規化ベクトルをノルムが規定値（この例では１）をとる正規化基底ベクトルに正規化するものである。これらの正規化処理の詳細については後述する。【００５８】
次に、クラスタリング部１４３において、メモリ１００の記憶領域１０６から正規化基底ベクトルAp''（f）を読み出し、これらの正規化基底ベクトルAp''（f）をクラスタリングしてＮ個のクラスタＣｋ（k=1,...,N）を生成し、各クラスタＣｋとそれらのセントロイド（中心ベクトル）ηｋを特定する情報をメモリ１００の記憶領域１０８，１０９に格納する（ステップＳ５）。このクラスタリングは、例えば、各クラスタＣｋの要素（正規化基底ベクトルAv''（f））と各クラスタＣｋのセントロイドηｋとの間の二乗和Ｕｋの総和Ｕ【００５９】
【数１３】
0004406428.tif 000014
を最小値化することを基準に行われる。この最小値化は、例えば、非特許文献６などで解説されているk-meansクラスタリングを用いることによって効果的に行うことができる。なお、クラスタＣｋのセントロイドηｋは、【００６０】
【数１４】
0004406428.tif 000015
によって計算される。ここで｜Ck｜はクラスタCkにおける要素（正規化基底ベクトルAv'' （f））の数である。また、ここでは、距離としてユークリッド距離の自乗を用いているが、これを一般化したミンコフスキー距離などを用いてもよい。なお、正規化基底ベクトルAp'' （f）がクラスタを形成する理由については後述する。
次に、順列算出部１４４において、メモリ１００の記憶領域１０６から正規化基底ベクトルAp''（f）を、記憶領域１０９から各クラスタCｋのセントロイドηｋをそれぞれ読み出す。そして、順列算出部１４４は、これらを用い、第１の分離行列W（f）の要素を並び替えるための順列Πｆ（{1,2,…,N}から{1,2,…,N}への全単射な関数）を周波数ｆ毎に算出し、メモリ１００の記憶領域１１０に格納する（ステップＳ６）。なお、この順列Πｆは、【００６１】
【数１５】
0004406428.tif 000016
によって決定される。ここで式（１３）におけるargminΠ・は、・を最小値化するΠを意味する。なお、式（１３）に従って順列Πｆを決定する手順としては、例えば、採り得るすべての順列Π（Ｎ！通り）に対して
【００６２】
【数１６】
0004406428.tif 000017
を算出し、その最小値に対応するΠを順列Πｆとして決定する手順を例示できる。以下にこの具体例を示す。
［順列Πｆ決定の具体例１］
信号源の数Ｎが３であり、ある周波数ｆにおける正規化基底ベクトルA1''（f），A2''（f），A3''（f）と各セントロイドη１，η２，η３との距離の自乗が、以下の表に示されるものであったとする。【００６３】
【表１】
0004406428.tif 000018
この場合、式（１３）により得られる順列は、
Πｆ：［1,2,3］→[2,3,1]
となる。なぜなら、
‖η1- Aπ（1）''（f）‖2=‖η1- A2''（f）‖2=0.1
‖η2- Aπ（2）''（f）‖2=‖η2- A3''（f）‖2=0.2
‖η3- Aπ（3）''（f）‖2=‖η3- A1''（f）‖2=0.15
という組合せが、【００６４】
【数１７】
0004406428.tif 000019
を最小とするからである（順列Πｆ決定の具体例１の説明終わり）。
しかし、この手順はＮが大きくなると現実的ではない。そのため、その近似法として、順番に‖ηk-AΠ（k）''（f）‖2を最小値化するAΠ（k）''（f）を重複がないように選択していき、この選択されたAΠ（k）''（f）を正規化基底ベクトルAk''（f）に移す順列を順列Πｆとする手順等を用いてもよい。以下、上述の［順列Πｆ決定の具体例１］と同じ条件に対し、この近似法を適用して順列Πｆを決定する手順を示す。【００６５】
［順列Πｆ決定の具体例２］
まず、上記の表１の場合、距離の自乗の最小値は０．１（正規化基底ベクトルA2''（f）と各セントロイドη１との距離の自乗）であるため、Π（１）＝２を決定する。そして、正規化基底ベクトルA2''（f）と各セントロイドη１とに関連する行と列とを消すと以下のようになる。【００６６】
【表２】
0004406428.tif 000020
この表２の場合、距離の自乗の最小値は０．１５（正規化基底ベクトルA1''（f）と各セントロイドη3との距離の自乗）であるため、Π（３）＝１を決定する。そして、最後にΠ（２）に残りの３を割り当てる（［順列Πｆ決定の具体例２］の説明終わり）。
次に、並び替え部１４５において、メモリ１００の記憶領域１０３から第１の分離行列W（f）を、記憶領域１１０から順列Πｆを読み出す。そして、並び替え部１４５は、この第１の分離行列W（f）の行を順列Πｆに従って並び替えた第２の分離行列W’（f）を生成してメモリ１００の記憶領域１１１に格納する（ステップＳ７）。なお、第１の分離行列W（f）の行を順列Πｆに従って並び替えるとは、上述のムーア・ペンローズ型一般化逆行列W+（f）における要素AΠ（k）''（f）からAk''（f）への並び替えに対応する並び替えを、第１の分離行列W（f）に対して行うことを意味する。すなわち、第１の分離行列W（f）の順列Πｆ（ｋ）行目が第２の分離行列W’（f）のｋ行目になるように並び替える。上述の［順列Πｆ決定の具体例１，２］の場合、第１の分離行列W（f）の２，３，１行目が、それぞれ第２の分離行列W’（f）の１，２，３行目になる。【００６７】
［分離信号生成部１５０の処理］
その後、分離信号生成部１５０において、メモリ１００の記憶領域１０２から周波数領域の混合信号Xq（f,τ）を、記憶領域１１１から第２の分離行列W’（f）を読み出す。そして、分離信号生成部１５０は、周波数領域の混合信号Xq（f,τ）からなる混合信号ベクトルX（f,τ）=[X1（f,τ）,...,XM（f,τ）]Tと第２の分離行列W’（f）とを用い、分離信号ベクトル
Y（f,τ）=W’（f）・X（f,τ）
を算出し、これらの各要素である周波数領域の信号（「周波数領域の混合信号」と呼ぶ）Yp（f,τ）をメモリ１００の記憶領域１１２に格納する（ステップＳ８）。【００６８】
［時間領域変換部１６０の処理］
最後に時間領域変換部１６０において、メモリ１００の記憶領域１１２から周波数領域の分離信号Yp（f,τ）を読み出し、これを添え字ｐ毎に（Yp（f,τ）毎に）短時間逆フーリエ変換等によって時間領域の分離信号yp（t）に変換し、これらの時間領域の分離信号yp（t）をメモリ１００の記憶領域１１３に格納する（ステップＳ９）。
［正規化処理の詳細（ステップＳ４の詳細）］
次に、前述した基底ベクトル正規化部１４２における正規化処理の詳細（ステップＳ４の詳細）について説明する。【００６９】
図６は、この正規化処理の詳細を説明するためのフローチャートである。
まず、制御部１７０（図３）においてパラメータｐに１を代入し、これを一時メモリ１７１に格納する（ステップＳ１１）。また、制御部１７０においてパラメータｑに１を代入し、これを一時メモリ１７１に格納する（ステップＳ１２）。次に、周波数正規化部１４２ａ（図４）において、メモリ１００の記憶領域１０７から前述のパラメータｄ，ｃ，Ｑを読み出し、一時メモリ１７１からパラメータｐ，ｑを読み出し、基底ベクトルAp（f）の要素Aqp（f）に対し、【００７０】
【数１８】
0004406428.tif 000021
の演算を行い、当該演算結果Aqp’（f）を周波数正規化ベクトルAp’（f）の各要素Aqp’（f）としてメモリ１００の記憶領域１０５に格納する（ステップＳ１３）。なお、ａｒｇ［・］は偏角を意味し、ｊは虚数単位を意味する。
より詳細には、まず、周波数正規化部１４２ａの第１正規化部１４２ａａが、基底ベクトルAp（f）の特定の１つの要素AQp（f）を基準として当該基底ベクトルAp（f）の各要素Aqp（f）の偏角を、以下の演算によって正規化する。【００７１】
【数１９】
0004406428.tif 000022
次に、周波数正規化部１４２ａの第２正規化部１４２ａｂが、第１正規化部１４２ａａで正規化された各要素Aqp'''（f）の偏角を、以下のように周波数ｆに比例した値4fc-1dで除算する。【００７２】
【数２０】
0004406428.tif 000023
次に、制御部１７０において、一時メモリ１７１に格納されたパラメータｑがｑ＝Ｍを満たすか否かを判断する（ステップＳ１４）。ここでｑ＝Ｍでなければ、制御部１７０は、ｑ＋１の演算結果を新たなパラメータｑの値とし、これを一時メモリ１７１に格納し（ステップＳ１５）、処理をステップＳ１３へ戻す。一方、ｑ＝Ｍであれば、制御部１７０は、さらにｐ＝Ｎを満たすか否かを判断する（ステップＳ１６）。【００７３】
ここでｐ＝Ｎでなければ、制御部１７０において、ｐ＋１の演算結果を新たなパラメータｐの値とし、これを一時メモリ１７１に格納し（ステップＳ１７）、処理をステップＳ１２へ戻す。一方、ｐ＝Ｎであれば、制御部１７０においてパラメータｐに１を代入し、これを一時メモリ１７１に格納し（ステップＳ１８）、ノルム正規化部１４２ｂにおける処理を開始する。まず、ノルム正規化部１４２ｂにおいて、メモリ１００の記憶領域１０５から周波数正規化ベクトルAp’（f）の各要素Aqp’（f）を読み出し、【００７４】
【数２１】
0004406428.tif 000024
の演算を行って、周波数正規化ベクトルAp’（f）のノルム‖Ap’（f）‖を求め、周波数正規化ベクトルAp’（f）及びそのノルム‖Ap’（f）‖を一時メモリ１７１に格納する（ステップＳ１９）。
次に、ノルム正規化部１４２ｂにおいて、一時メモリ１７１から周波数正規化ベクトルAp’（f）及びそのノルム‖Ap’（f）‖を読み出し、
Ap'' （f）=Ap’（f）／‖Ap’（f）‖    …（18）
の演算を行って正規化基底ベクトルAp'' （f）を求め、これをメモリ１００の記憶領域１０６に格納する（ステップＳ２０）。【００７５】
その後、制御部１７０において、一時メモリ１７１に格納されたパラメータｐがｐ＝Ｎを満たすか否かを判断する（ステップＳ２１）。ここでｐ＝Ｎでなければ、制御部１７０は、ｐ＋１の演算結果を新たなパラメータｐの値とし、これを一時メモリ１７１に格納し（ステップＳ２２）、ステップＳ１９の処理に戻す。一方、ｐ＝Ｎであれば、制御部１７０は、ステップＳ４の処理を終了させる。
このように生成された正規化基底ベクトルAp''（f）は、周波数に依存せず、信号源の位置のみに依存するベクトルとなる。その結果、この正規化基底ベクトルAp''（f）はクラスタを形成することになる。以下にこの理由を説明する。【００７６】
［正規化基底ベクトルAp''（f）がクラスタを形成する理由］
基底ベクトルAp（f）の各要素Aqp（f）は、源信号ｐに対応する信号源ｋからセンサｑヘの周波数応答Ｈｑｋに比例した（ある複素数スカラが掛かった）ものになっている。これらの複素数スカラは離散時間に応じて変化するが、源信号ｐとセンサｑとに対応する複素数スカラと、源信号ｐとセンサＱとに対応する複素数スカラとの相対値は、離散時間が変化しても一定である（周波数ｆが同じであれば）。すなわち、周波数ｆが同じであれば、源信号ｐとセンサｑとに対応する複素数スカラの偏角と、源信号ｐとセンサＱとに対応する複素数スカラの偏角との相対値は一定である。
【００７７】
前述のように、周波数正規化部１４２ａの第１正規化部１４２ａａは、基底ベクトルAp（f）の特定の１つの要素AQp（f）を基準として当該基底ベクトルAp（f）の各要素Aqp（f）の偏角を正規化する。これにより、上述の複素数スカラの位相に起因する不確定性を取り除き、源信号ｐとセンサｑとに対応する基底ベクトルAp（f）の要素Aqp（f）の偏角を、源信号ｐと基準センサＱ（基準値Ｑに対応）とに対応する基底ベクトルAp（f）の要素AQp（f）の偏角に対する相対値として表現する。なお、この場合、要素AQp（f）の偏角に対応する当該相対値は０と表現される。信号源ｋからセンサqヘの周波数応答を、反射や残響の無い直接波モデルで近似して考える。すると、上記の第１正規化部１４２ａａにより正規化された偏角は、信号源ｋからセンサへの波の到達時間差と周波数ｆの双方に比例したものになる。ここでの到達時間差とは、信号源ｋからの波がセンサｑに到達するまでの時間と、当該波が基準センサＱに到達するまでの時間との時間差である。【００７８】
また、前述のように、第２正規化部１４２ａｂは、第１正規化部１４２ａａで正規化された各要素Aqp''' （f）の偏角を周波数ｆに比例した値で除算する。これにより、各要素Aqp''' （f）を、それらの偏角の周波数依存性を排除した各要素Aqp'（f）に正規化する。これにより、正規化された各要素Aqp'（f）は、直接波モデルに従えば、信号源ｋからセンサへの波の到達時間差のみに依存するものとなる。ここで、信号源ｋからセンサへの波の到達時間差は、信号源ｋ，センサｑ，基準センサＱの相対位置にのみ依存する。そのため、信号源ｋ，センサｑ，基準センサＱが同じであれば、周波数ｆが異なっても各要素Aqp'（f）の偏角は同一となる。従って、周波数正規化ベクトルAp' （f）は、周波数ｆには依存せず、信号源ｋの位置のみに依存する。【００７９】
そのため、周波数正規化ベクトルAp' （f）のノルムを正規化した正規化基底ベクトルAp''（f）のクラスタリングによって、同じ信号源毎に対応するクラスタが形成される。なお、実際の環境では、反射や残響などの影響により、直接波モデルは厳密には満たされないが、後述の実験結果に示すとおり十分に良い近似となっている。
次に、正規化基底ベクトルAp''（f）がクラスタを形成する理由を、モデルを用いて説明する。前述した式（１）に示されるインパルス応答hqk（r）を直接波（近距離場）混合モデルを用いて近似し、周波数領域で表すと【００８０】
【数２２】
0004406428.tif 000025
となる。ここで、ｄｑｋは信号源ｋとセンサｑとの間の距離であり、また、減衰1/dqｋは、距離ｄｑｋによって決まり、遅延（dqk-dQk）/cは、基準センサＱの位置で正規化された距離によって決まる。
また、独立成分分析（ICA）における順序の任意性とスケーリングの任意性とを考慮すると、基底ベクトルAp（f）と式（４）における信号源kから各センサへの周波数応答をまとめたベクトルHk（f）との間には、以下の関係が成り立つ。【００８１】
Ap（f）=εp・Hk（f），  Aqp（f）=εp・Hqk（f） …（20）
ここで、εpはスケーリングの任意性を表現する複素スカラ値であり、添え字ｋとｐが異なる可能性が順序の任意性を表現している。ここで、式（１６）（１８）（１９）（２０）より、
【００８２】
【数２３】
0004406428.tif 000026
となる。この式から分かるように、正規化基底ベクトルAp''（f）の各要素Aqp''（f）は、周波数ｆからは独立であり、信号源ｋとセンサｑの位置のみに依存する。従って、正規化基底ベクトルAp''（f）をクラスタリングすると、同じ信号源ごとに対応するクラスタが形成される。
同様なことは、信号の減衰を考慮していない近距離場混合モデルでモデル化した場合にも言える。すなわち、前述した式（１）に示される畳み込み混合モデルを、減衰を無視した近距離場混合モデルで近似し、周波数領域で表すと【００８３】
【数２４】
0004406428.tif 000027
となる。この場合、式（１６）（１８）（２０）（２２）より、【００８４】
【数２５】
0004406428.tif 000028
となる。この場合も、正規化基底ベクトルAp''（f）の各要素Aqp''（f）は、周波数ｆからは独立であり、信号源ｋとセンサｑの位置のみに依存している。
さらに、近距離場混合モデルだけではなく、遠距離場混合モデルでも同様なことがいえる。すなわち、前述した式（１）に示される畳み込み混合モデルを遠距離場混合モデルで近似し、周波数領域で表すと【００８５】
【数２６】
0004406428.tif 000029
となる。なお、SEq及びSEQはセンサｑ，Ｑの位置を示すベクトルである。またθkqQは、センサｑ，Ｑを結ぶ直線と、センサｑ，Ｑの中心点と信号源ｋとを結ぶ直線と、で成す角度を示す。この場合、式（１６）（１８）（２０）（２４）より、【００８６】
【数２７】
0004406428.tif 000030
となる。この場合も、正規化基底ベクトルAp''（f）の各要素Aqp''（f）は、周波数ｆからは独立であり、信号源ｋとセンサｑの位置のみに依存している。
また、式（２１）からパラメータｄの値は、ｄ＞ｄmax／２であることが望ましく（ｄmaxは要素AQp（f）に対応する基準センサＱと他のセンサとの最大距離を意味する。）、より好ましくはｄ≧ｄmaxであることが望ましく、さらにより好ましくはｄ＝ｄmaxであることが望ましいことが分かる。以下、この理由について説明する。
【００８７】
図７及び図８は、パラメータｄ毎の正規化基底ベクトルの要素Aqp''（f）とその偏角arg[Aqp''（f）]との関係を説明するための複素平面図である。なお、これらにおける横軸は実軸を縦軸は虚軸を示している。
図７Ａはｄmax／２≧dの場合における複素平面図である。ここで上述のdmaxの定義より、任意のｑ及びｋに対しdqk-dQkの絶対値はdmax以下となる。よってｄmax／２≧dの場合、（π/2）・（dqk-dQk）/d≦-π，（π/2）・（dqk-dQk）/d≧πと成り得る。その結果、式（２１）で表されるAqp''（f）の偏角arg[Aqp''（f）]は、2πを超えるα1≦arg[Aqp''（f）]≦α2（α1≦-π、α2≧π）の範囲に分布する可能性がある。そのため、異なる正規化基底ベクトルの要素Aqp''（f）の偏角が一致する可能性があり、前述したクラスタリングにおいて、異なる正規化基底ベクトルAp''（f）を同じクラスタにクラスタリングしてしまう可能性がある。よって、ｄ＞ｄmax／２であることが望ましい。しかし、この偏角の重複範囲に対応する正規化基底ベクトルAp''（f）のサンプルが存在しないならばｄmax／２≧ｄとしても問題はない。【００８８】
図７Ｂはｄmax／２＜ｄ＜ｄmaxの場合における複素平面図である。この場合、-π＜（π/2）・（dqk-dQk）/d＜-π/2，π/2＜（π/2）・（dqk-dQk）/d＜πと成り得る。その結果、式（２１）で表されるAqp''（f）の偏角arg[Aqp''（f）]は、β1≦arg[Aqp''（f）]≦β2（-π＜β1＜-π/2、π/2＜β2＜π）の範囲に分布する可能性がある。そのため、-π＜arg[Aqp''（f）]＜-π/2及びπ/2＜arg[Aqp''（f）]＜πの範囲において、異なる周波数正規化ベクトル要素間の偏角の差の増加に伴い、これらの周波数正規化ベクトル要素間の距離が単調増加しないこともありうる。これは、前述したクラスタリングの精度を低下させる可能性がある。よってｄ≧ｄmaxであることがより望ましい。【００８９】
図８Ａはｄ＝ｄmaxの場合における複素平面図であり、図８Ｂはｄ＞ｄmaxの場合における複素平面図である。ここでｄ＞ｄmaxの場合、-π/2＜（π/2）・（dqk-dQk）/d＜0，0＜（π/2）・（dqk-dQk）/d＜π/2と成り得る。その結果、式（２１）で表されるAqp''（f）の偏角arg[Aqp''（f）]は、図８Ｂに示すように、γ1≦arg[Aqp''（f）]≦γ2（-π/2＜γ1＜0、0＜γ2＜π/2）の範囲に分布する。そして、ｄが大きくなればなるほどその分布範囲は狭くなっていき、狭い範囲にクラスタが密集していく。これは、前述したクラスタリングの精度を低下させる。【００９０】
これに対し、ｄ＝ｄmaxである場合、-π/2≦（π/2）・（dqk-dQk）/d＜0，0＜（π/2）・（dqk-dQk）/d≦π/2と成り得る。その結果、式（２１）で表されるAqp''（f）の偏角arg[Aqp''（f）]は図８Ａに示すように-π/2≦arg[Aqp''（f）]≦π/2の範囲に分布する。この場合、周波数正規化ベクトルの要素間における偏角の差の増加に対して、それらの距離も単調増加するという関係を維持しつつ、できるだけ広い範囲にクラスタを分散させることができる。その結果、一般的にクラスタリングの精度を向上させることができる。
〔第２の実施の形態（「第１の本発明」の例）〕
次に本発明における第２の実施の形態について説明する。【００９１】
第１の実施の形態では、基底ベクトルから得られる情報によりパーミュテーション問題を解決したが、本形態では、これに特開２００４−１４５１７２号公報や「H. Sawada, R. Mukai, S. Araki, S. Makino, "A Robust and Precise Method for Solving the Permutation Problem of Frequency-Domain Blind Source Separation," IEEE Trans. Speech and Audio processing, vol. 12, no. 5, pp. 530-538, Sep. 2004.（以下「参考文献」と呼ぶ）」で示されているような分離信号のエンベロープの情報を統合することで、より精度よくパーミュテーション問題を解決する。なお、これらの文献では、基底ベクトルの代わりに信号源の方向に関する情報が利用されている。【００９２】
以下では、第１の実施の形態との相違点を中心に説明を行い、これと共通する事項については説明を省略する。
＜構成＞
図９は、第２の実施の形態における信号分離装置２００のブロック図の例示である。なお、本形態の信号分離装置２００もＣＰＵ１０ａ（図２）に信号分離プログラムが読み込まれることにより構成されるものである。また、図１０Ａは、図９におけるパーミュテーション問題解決部２４０の詳細を例示したブロック図であり、図１０Ｂは、図１０Ａの順列修正部２４７の詳細を例示したブロック図である。なお、図９及び図１０において第１の実施の形態と共通する部分については第１の実施の形態と同じ符号を付した。また、図９及び図１０における破線の矢印は理論上の情報の流れを示し、実線の矢印は実際のデータの流れを示す。また、これらの図において制御部１７０に出入りするデータの流れに対応する矢印は省略してあり、図１０における実際のデータの流れに対応する矢印も省略してある。【００９３】
本形態の第１の実施の形態のとの相違点は、主にパーミュテーション問題解決部２４０の構成である。すなわち、本形態のパーミュテーション問題解決部２４０は、第１の実施の形態のパーミュテーション問題解決部１４０に、順列評価部２４６と順列修正部２４７とを付加したものとなっている（図９，図１０Ａ）。なお、順列評価部２４６は、周波数毎に順列の信頼度を評価するものであり、順列修正部２４７は、順列の信頼度が低いと評価された周波数に対し、分離信号のエンベロープを用いて新たな順列を算出するものである。また、順列修正部２４７は、判定部２４７ａと分離信号生成部２４７ｂとエンベロープ算出部２４７ｃと順列再算出部２４７ｄと再並び替え部２４７ｅを有している（図１０Ｂ）。また、本形態では、順列算出部１４４と順列修正部２４７とによって請求項４の「順列算出部」を構成している。【００９４】
＜処理＞
図１１は、第２の実施の形態における信号分離装置２００の処理の全体を説明するためのフローチャートである。以下、この図に沿って、本形態における信号分離装置２００の処理を説明していく。
ステップＳ５１からＳ５７の処理は、第１の実施の形態におけるステップＳ１からＳ７までと同じであるため説明を省略する。本形態では、このステップＳ５７の処理の後に、順列評価部２４６において、周波数毎に順列Πfの信頼度を評価し、順列Πfの信頼度が低いと評価された周波数に対し、分離信号のエンベロープを用いて順列Πf’を算出し、この周波数のみについて第２の分離行列W’（f）の行を順列Πf’に従って並び替えて第３の分離行列W''（f）を生成し、メモリ１００の記憶領域１１０に格納する（ステップＳ５８）。なお、この処理の詳細については後述する。【００９５】
その後、分離信号生成部１５０において、メモリ１００の記憶領域１０２から周波数領域の混合信号Xq（f,τ）を、記憶領域１１１から第３の分離行列W''（f）を読み出す。そして、分離信号生成部１５０は、周波数領域の混合信号Xq（f,τ）からなる混合信号ベクトルX（f,τ）=[X1（f,τ）,...,XM（f,τ）]Tと第３の分離行列W''（f）とを用い、分離信号ベクトル
Y（f,τ）=W''（f）・X（f,τ）
を算出し、周波数領域の分離信号Yp（f,τ）をメモリ１００の記憶領域１１２に格納する（ステップＳ５９）。【００９６】
最後に時間領域変換部１６０において、メモリ１００の記憶領域１１２から周波数領域の分離信号Yp（f,τ）を読み出し、これを添え字ｐ毎に時間領域の分離信号yp（t）に変換し、これらの時間領域の分離信号yp（t）をメモリ１００の記憶領域１１３に格納する（ステップＳ６０）。
［ステップＳ５８の処理の詳細」
図１２及び図１３は、図１１におけるステップＳ５８の例を説明するためのフローチャートである。以下、このフローチャートに沿ってステップＳ５８の詳細を説明していく。【００９７】
まず、制御部１７０においてパラメータｆに０を代入し、集合Ｆを空集合とし、これらを示す情報を一時メモリ１７１に格納する（ステップＳ７１）。次に、順列評価部２４６において、メモリ１００の記憶領域１１０に格納された順列Πｆの信頼度を周波数毎に評価し、その評価結果trust（f）を一時メモリ１７１に格納する（ステップＳ７２）。ここで、順列Πｆの信頼性が高いとは、例えば、正規化基底ベクトルAp''（f）が、それぞれに対応するセントロイドηkに十分近いということである。また、正規化基底ベクトルAp''（f）が、それぞれに対応するセントロイドηkに十分近いかどうかは、例えば、正規化基底ベクトルAp''（f）とセントロイドηkとの距離がクラスタＣｋの分散Ｕｋ／｜Ｃｋ｜より小さいかどうか、すなわち、
Uk/|Ck|>‖ηk-AΠ（k）''（f）‖2    …（26）
により判定できる。よって、このステップでは、例えばまず、順列評価部２４６において、メモリ１００の記憶領域１０５から正規化基底ベクトルAp''（f）を、記憶領域１０９からセントロイドηｋを、記憶領域１１０から順列Πｆをそれぞれ読み出す。そして、順列評価部２４６は、周波数ｆ毎に、式（２６）を満たすか否かを判断し、満たす場合にはtrust（f）=1を、満たさない場合にはtrust（f）=0を出力し、一時メモリ１７１に格納する。【００９８】
次に、順列修正部２４７の判定部２４７ａが、一時メモリ１７１から周波数ｆ毎に評価結果trust（f）を読み出し、trust（f）=1であるか否かを判断する（ステップＳ７３）。ここで、trust（f）=0である場合、そのままステップＳ７６の処理に進む。一方、trust（f）=1である場合、制御部１７０において、集合Ｆと｛ｆ｝の和集合を新たな集合Ｆとして一時メモリ１７１に格納し（ステップＳ７４）、再並び替え部２４７ｅにおいて、この周波数ｆにおける第２の分離行列W'（f）を第３の分離行列W''（f）としてメモリ１００の記憶領域１１１に格納し（ステップＳ７５）、ステップＳ７６に進む。【００９９】
ステップＳ７６では、制御部１７０において、一時メモリ１７１に格納されたパラメータｆの値がf=（L-1）fs/Lを満たすか否かを判断し（ステップＳ７６）、これを満たさなければf+fs/Lの演算結果を新たなパラメータｆの値として（ステップＳ７７）一時メモリ１７１に格納し、ステップＳ７２に戻る。
一方、f=（L-1）fs/Lを満たす場合、分離信号生成部２４７ｂにおいて、集合Ｆに属しない周波数ｆを１つ選択する。そして、分離信号生成部２４７ｂは、この周波数ｆとその近傍であって集合Ｆに属する全ての周波数ｇ（g∈F，|g-f|≦δ，δは定数）とに対し、メモリ１００の記憶領域１０２から周波数領域の混合信号X（f,τ）=[X1（f,τ），…，XM（f,τ）]T, X（g,τ）=[X1（g,τ），…，XM（g,τ）]Tを読み出し、記憶領域１１１から第２の分離行列W’（f）,W’（g）を読み出し、
Y（f,τ）=W’（f）・X（f,τ）
Y（g,τ）=W’（g）・X（g,τ）
によって分離信号Y（f,τ）=[Y1（f,τ），…，YN（f,τ）]T, Y（g,τ）=[Y1（g,τ），…，YN（g,τ）]Tを算出して一時メモリ１７１に格納する（ステップＳ７８）。【０１００】
次に、エンベロープ算出部２４７ｃにおいて、一時メモリ１７１から全ての周波数領域の分離信号Yp（f,τ）, Yp（g,τ）を読み出し、これらのエンベロープ
vpf（τ）=|Yp（f,τ）|
vpg（τ）=|Yp（g,τ）|
を算出して一時メモリ１７１に格納する（ステップＳ７９）。
そして、順列再算出部２４７ｄにおいて、周波数の差δ以下の近傍における相関corの和の最大値【０１０１】
【数２８】
0004406428.tif 000031
を算出して一時メモリに格納する（ステップＳ８０）。ここで、Π’は周波数ｇに対し、既に決定された順列である。なお、この式における相関cor（Φ，Ψ）は、２つの信号Φ，Ψの相関を示し、以下のように定義される。
cor（Φ，Ψ）=（<Φ，Ψ>-<Φ>・<Ψ>）/（σΦ・σΨ）
ただし、<ζ>はζの時間平均であり、σΦはΦの標準偏差である。
【０１０２】
また、順列再算出部２４７ｄにおいて、この相関corの和を最大値化する順列【０１０３】
【数２９】
0004406428.tif 000032
を算出してメモリ１００の記憶領域１１０に格納する（ステップＳ８１）。ここで、Π’は周波数ｇに対し、既に決定された順列である。なお、argmaxΠνは、νを最大値化するΠを意味する。
次に、制御部１７０において集合Ｆと｛ζ｝（ζ=argmaxfRf）との和集合を新たな集合Ｆとして一時メモリ１７１に格納する（ステップＳ８２）。そして、再並び替え部２４７ｅにおいて、ｆ＝ζとし、第２の分離行列W’（f）の行を順列Πf’に従って並び替えた第３の分離行列W''（f）を生成し、これをメモリ１００の記憶領域１１１に格納する（ステップＳ８３）。【０１０４】
次に、制御部１７０において、一時メモリ１７１に格納された集合Ｆが全ての離散化された周波数f=0,fs/L,…,fs（L-1）/Lの要素を有するか否かを判断する（ステップＳ８４）。ここで、集合Ｆが全ての離散化された周波数f=0,fs/L,…,fs（L-1）/Lの要素を有しないのであれば、制御部１７０は、処理をステップＳ７８に戻す。一方、集合Ｆが全ての離散化された周波数f=0,fs/L,…,fs（L-1）/Lの要素を有するのであれば、制御部１７０は、ステップＳ５８の処理を終了させる。
なお、以上の方法を用いず、前述した特開２００４−１４５１７２号公報や「参考文献」に記載されたその他の方法を用い、ステップＳ５８の処理を実行してもよい。【０１０５】
＜実験結果＞
次に、上述した第１の実施の形態及び第２の実施の形態による音源分離実験結果を示す。
［実験結果１］
１つ目の実験は不規則なセンサ配置によるものである。実験条件は図１４Ａに示す通りである。ここでは４つの無指向性マイクを不規則に配置したものを用いたが、これらの配置情報としては、マイク間隔の上限が４ｃｍであるということだけを信号分離装置に与えた。音源の数は３つであり、３秒間の英語の音声をスピーカから流した。結果を図１４ＢにSIR（signal-to-interference ratio）によって示す。大きいSIRほど分離性能が良いことを示す。ここではパーミュテーション問題を解決する４種類の手法に関して結果を比較した。Envは分離信号のエンベロープ|Yp（f,τ）|の情報のみを用いたもの、Basisは正規化基底ベクトルAp''（f）のクラスタリングによるもの（第１の実施の形態の手法）、Basis+Envはそれら２種類の情報を統合してより精度良く問題を解決したもの（第２の実施の形態の手法）、Optimalは源信号skとインパルス応答hqk（r）を知って得た最適な順列によるものである。【０１０６】
結果を考察すると、Envだけでは性能にばらつきがあるが、第１の実施の形態のBasisはそれだけで十分な分離性能を得ている。また、双方を統合した第２の実施の形態のBasis+Envの結果はほぼOptimalに近い。このように不規則なセンサ配置を用いた場合でも、本発明により、周波数領域でのブラインド信号分離が高性能で達成された。
［実験結果２］
２つ目の実験は規則的なセンサ配置によるものである。実験条件は図１５Ａに示す通りである。３つの無指向性マイクを直線状に４ｃｍ間隔で配置した。音源の数は同様に３つであり、３秒問の英語の音声をスピーカから流した。結果を図１５Ｂに示す。この実験では、従来技術として示した信号源の位置に関する推定値をクラスタリングする方法も比較対象とし、６種類の手法に関して結果を比較した。DOAは信号の到来方向（DOA:direction-of-arrival）の推定値のみでパーミュテーション問題を解決したもの、DOA+Envはさらに分離信号のエンベロープ情報を統合したものある。【０１０７】
従来技術であるDOA及びDOA＋Envと、本発明であるBasis及びBasis+Envの結果を比較すると、従来手法が適用可能である規則的なセンサ配置に対しても、総じて本発明により結果が改善されていることがわかる。なお、演算量に関しては、従来技術とほぼ同等であった。
＜第１，２の実施の形態の特徴＞
以上より、上述した第１，２の実施の形態の特徴をまとめると以下のようになる。
（1）センサの配置情報を厳密に知る必要がなく、ある基準センサから他のセンサヘの距離の上限のみを知ればよいため、不規則なセンサ配置を採用でき、さらに位置をキャリブレーションする必要が無い。（2）基底ベクトルから得られる情報をすべて用いてクラスタリングを行うためパーミュテーション問題をより正確に解決することができ、信号分離の性能が向上する。【０１０８】
なお、本発明は上述の実施の形態に限定されるものではない。例えば、上述した各実施の形態では、一般化行列としてムーア・ペンローズ型一般化逆行列を用いたが他の一般化行列を利用する形態でもよい。
また、第１の実施の形態では、周波数正規化部１４２ａの第１正規化部１４２ａａが、基底ベクトルAp（f）の特定の１つの要素AQp（f）を基準とし、式（１５）によって当該基底ベクトルAp（f）の各要素Aqp（f）の偏角を正規化することとした。しかし、例えば、以下のような式に従って、第１正規化部１４２ａａが、基底ベクトルAp（f）の特定の１つの要素AQp（f）を基準として、当該基底ベクトルAp（f）の各要素Aqp（f）の偏角を正規化する構成であってもよい。【０１０９】
【数３０】
0004406428.tif 000033
ただし、・＊は、・の複素共役である。また、ψ｛・｝は関数であり、クラスタリング精度の観点から好ましくは単調増加関数であることが望ましい。
また、周波数正規化部１４２ａが式（１４）の代わりに、【０１１０】
【数３１】
0004406428.tif 000034
等の演算によって周波数正規化を行うこととしてもよい。ただし、ρは定数である（例えば、ρ＝１）。
また、上述の実施の形態では、ノルム正規化部１４２ｂにおいてノルムが１になるような正規化を行ったが、ノルムが１以外の規定値となるように正規化を行ってもよい。さらに、ノルム正規化部１４２ｂを設けず、ノルム正規化を行わない構成であってもよい。この場合、クラスタリング部１４３は、周波数正規化ベクトルAp’（f）のクラスタリングを行うことになる。しかし、周波数正規化ベクトルAp’（f）は、ノルムが統一されていない。そこで、この場合のクラスタリング基準は、ベクトルがノルムを含めて似ているかどうかではなく、ベクトルの方向のみが似ているかどうかになる。これは、類似度を用いた評価になる。類似度の１つとしてコサイン距離
cosθ=｜Ap’H（f）・ηk｜/（‖Ap’（f）‖・‖ηk‖）
を例示できる。ここでθは、周波数正規化ベクトルAp’（f）と、セントロイドηkのベクトルとがなす角度である。コサイン距離を用いる場合、クラスタリング部１４３は、コサイン距離の総和【０１１１】
【数３２】
0004406428.tif 000035
を最小値化するクラスタを生成する。なお、セントロイドηkは、各クラスタのメンバの平均として算出する。
さらに、上記の第２の実施の形態では、周波数毎に順列の信頼度を評価し、順列の信頼度が低いと評価された周波数に対し、分離信号のエンベロープを用いて新たな順列を算出することとした。しかし、分離信号のエンベロープとクラスタの中心ベクトルと正規化基底ベクトルとを一度に統合して順列を算出することとしてもよい。【０１１２】
また、まず分離信号のエンベロープを用いて順列を算出し、その順列の信頼度を周波数毎に評価し、信頼度が低いと評価された周波数に対して、第１の実施の形態の手法を適用して新たな順列を算出することとしてもよい。
さらに、上記の第２の実施の形態では、第２の分離行列W’（f）を用いて分離信号のエンベロープを算出したが、第１の分離行列W（f）或いはその行を並び替えた分離行列を用いてエンベロープを算出する構成としてもよい。
また、パラメータｄは、全てのセンサｑに対して同一であってもよいし、各センサｑに対応して複数の値を設定しておいてもい。例えば、基準センサとセンサｑとの距離を、それぞれのセンサｑに対応するパラメータｄの値としてもよい。【０１１３】
〔第３の実施の形態（「第２の本発明」の例）〕
次に、本発明における第３の実施の形態について説明する。
本形態は、前述の原理を用い、目的信号の方向に関する情報を知ることなく、複数の信号源から発せられた信号が混合された混合信号から目的信号を抽出する形態である。
＜構成＞
本形態の信号分離装置も、第１の実施の形態と同様、公知のノイマン型のコンピュータに信号分離プログラムが読み込まれることにより構成される。図１６は、第３の実施の形態における信号分離装置１００１の構成を例示したブロック図である。【０１１４】
この図に例示するように、信号分離装置１００１は、記憶領域１１０１〜１１１４を有するメモリ１１００、周波数領域変換部１１２０、信号分離部１１３０、目的信号選択部１１４０、時間周波数マスキング部１１５０（「分離信号生成部」に相当）、時間領域変換部１１６０、制御部１１７０及び一時メモリ１１８０を有する。ここでメモリ１１００及び一時メモリ１１８０は、例えば、レジスタ１０ａｃ、補助記憶装置１０ｆのデータ領域１０ｆｂ或いはＲＡＭ１０ｄのデータ領域１０ｄｂ（図２）の少なくとも一部に相当するが、特にこれに限定されない。また、周波数領域変換部１１２０、信号分離部１１３０、目的信号選択部１１４０、時間周波数マスキング部１１５０、時間領域変換部１１６０、制御部１１７０及び一時メモリ１１８０は、例えば、ＣＰＵ１０ａ（図２）にＯＳプログラムや信号分離プログラムが読み込まれることにより構成されるものである。【０１１５】
図１７Ａは、図１６における目的信号選択部１１４０の詳細構成を例示したブロック図であり、図１７Ｂは、図１７Ａにおける基底ベクトルクラスタリング部１１４２の詳細構成を例示したブロック図である。
これらの図に例示するように、目的信号選択部１１４０は、逆行列算出部１１４１（「複素ベクトル生成部」に対応）、基底ベクトルクラスタリング部１１４２及び選択部１１４３を有している。また、基底ベクトルクラスタリング部１１４２は、周波数正規化部１１４２ａ（「正規化部」を構成）、ノルム正規化部１１４２ｂ（「正規化部」を構成）、クラスタリング部１１４２ｃ及び分散判定部１１４２ｄを有している。また、周波数正規化部１１４２ａは、第１正規化部１１４２ａａ及び第２正規化部１１４２ａｂを有している。【０１１６】
図１８Ａは、図１６における時間周波数マスキング部１１５０の詳細構成を例示したブロック図であり、図１８Ｂは、図１８Ａにおけるマスク生成部１１５１の詳細を例示したブロック図である。
これらの図に例示するように、時間周波数マスキング部１１５０は、マスク生成部１１５１及びマスキング部１１５２を有している。また、マスク生成部１１５１は、白色化行列生成部１１５１ａ、白色化部１１５１ｂ、角度計算部１１５１ｃ及び関数演算部１１５１ｄを有している。【０１１７】
なお、図１６〜図１８における実線の矢印は実際のデータの流れを示し、破線の矢印は理論上の情報の流れを示すが、制御部１１７０及び一時メモリ１１８０に出入りするデータの流れは省略してある。また、信号分離装置１００１は制御部１１７０の制御のもと各処理を実行する。また、特に示さない限り、制御部１１７０は、必要なデータを一時メモリ１１８０に読み書きしながら各処理を実行していく。
＜処理＞
次に、本形態における信号分離装置１００１の処理について説明する。【０１１８】
空間内にＮ個の信号源ｋ（k∈{1,2,...,N}）が存在し、それらの信号sk（t）（ｔはサンプリング時刻）が混合され、Ｍ個のセンサｑ（q∈{1,2,...,M}）で混合信号xqとして観測されるとする。本形態では、混合信号x1（t）,...,xM（t）のみから、何れかの信号源から発せられた目的信号を抽出してその他の妨害信号を抑圧した信号y（t）を算出する。なお、信号源の数Ｎはセンサの数Ｍよりも多くても少なくてもよく、また同数であってもよい。さらに、信号源の数Ｎを事前に知る必要もなく、また信号源の数が数えられないような状況であっても構わない。【０１１９】
［処理の全体］
図１９は、第３の実施の形態における信号分離処理の全体を説明するためのフローチャートである。以下、この図を用いて本形態における信号分離処理の全体を説明する。
まず、前処理において、Ｍ個のセンサで観測された時間領域の混合信号xq（t）（q∈{1,...,M}）をメモリ１１００の記憶領域１１０１に格納しておく。信号分離処理が開始されると、まず、周波数領域変換部１１２０において、メモリ１１００の記憶領域１１０１から時間領域の混合信号xq（t）を読み出す。そして、周波数領域変換部１１２０は、これらを短時間フーリエ変換などにより周波数領域の混合信号Xq（f,τ）に変換し、メモリ１１００の記憶領域１１０２に格納する（ステップＳ１０１）。【０１２０】
次に、信号分離部１１３０において、メモリ１１００の記憶領域１１０２から周波数領域の混合信号Xq（f,τ）を読み出す。この例の信号分離部１１３０は、読み出した混合信号Xq（f,τ）を要素とする混合信号ベクトルX（f,τ）=[X1（f,τ）,…,XM（f,τ）]Tに、独立成分分析（ICA）を適用し、周波数ｆ毎にＭ行Ｍ列の分離行列W（f）=[W1（f）,...,WM（f）]H（＊Hは行列＊の複素共役転置行列）と分離信号ベクトル
Y（f,τ）=W（f）・X（f,τ）  …（30）
とを算出する（ステップＳ１０２）。なお、算出された分離行列W（f）はメモリ１１００の記憶領域１１０３に格納される。また、分離信号ベクトルY（f,τ）=[Y1（f,τ）,…,YM（f,τ）]Tを構成する分離信号Yp（f,τ）（p∈{1,..,M}）は、記憶領域１１０７に格納される。ステップＳ１０２の詳細については後述する。【０１２１】
次に、目的信号選択部１１４０において、メモリ１１００の記憶領域１１０３から分離行列W（f）を読み出し、その一般化逆行列の列である基底ベクトルを正規化し、その正規化された基底ベクトルをクラスタリングする。そして、目的信号選択部１１４０は、そのクラスタの分散を指標として、メモリ１１００の記憶領域１１０７の分離信号から目的信号を含む選択信号YI（f）（f,τ）とそれに対応する基底ベクトルAI（f）（f）とを周波数ｆ毎に選択し、メモリ１１００の記憶領域１１１１に格納する（ステップＳ１０３）。なお、本形態では、センサに近く、そのためセンサで観測されるパワーがその他の源信号と比べて支配的であり、さらに情報として有用である信号を目的信号として抽出する。このステップＳ１０３の詳細については後述する。【０１２２】
次に、時間周波数マスキング部１１５０において、メモリ１１００の記憶領域１１０２から周波数領域の混合信号Xq（f,τ）を読み出し、記憶領域１１０４から選択信号YI（f）（f,τ）に対応する基底ベクトルAI（f）（f）を読み出し、これらを用いて時間周波数マスクM（f,τ）を生成して記憶領域１１１２に格納する（ステップＳ１０４）。なお、ステップＳ１０４の処理（時間周波数マスキング部１１５０の処理）の詳細については後述する。
次に、時間周波数マスキング部１１５０において、メモリ１１００の記憶領域１１０７から目的信号選択部１１４０で選択された選択信号YI（f）（f,τ）を読み出し、記憶領域１１１２から時間周波数マスクM（f,τ）を読み出す。そして、時間周波数マスキング部１１５０は、この時間周波数マスクM（f,τ）を選択信号YI（f）（f,τ）に適用し、選択信号YI（f）（f,τ）に残存する妨害信号成分をより一層抑圧したマスク処理選択信号YI（f）'（f,τ）を生成し、メモリ１１００の記憶領域１１１３に格納する（ステップＳ１０５）。なお、このステップＳ１０５の処理の詳細（時間周波数マスキング部１１５０の処理）については後述する。【０１２３】
最後に、時間領域変換部１１６０において、メモリ１１００の記憶領域１１１３から、周波数領域における選択された分離信号YI（f）'（f,τ）を読み出し、それに短時間逆フーリエ変換などを施して時間領域の分離信号y（t）を生成し、メモリ１１００の記憶領域１１１４に格納する（ステップＳ１０６）。
［ステップＳ１０２の処理（信号分離部１１３０の処理）の詳細］
前述のように、この例の信号分離部１１３０は、独立成分分析（ICA）によって、混合信号ベクトルX（f,τ）=[X1（f,τ）,…,XM（f,τ）]TからＭ行Ｍ列の分離行列W（f）=[W1（f）,...,WM（f）]Hと分離信号ベクトルY（f,τ）=[Y1（f,τ）,…,YM（f,τ）]Tとを算出する（ステップＳ１０２）。独立成分分析（ICA）は、分離信号ベクトルY（f,τ）=[Y1（f,τ）,…,YM（f,τ）]Tの各要素が互いに独立になるように分離行列W（f）を算出する手法であり、そのアルゴリズムとしては、非特許文献４に記載されているものなど様々なものが紹介されている。なお、独立成分分析（ICA）は、パワーが小さく、ガウス性が高い妨害信号よりも、パワーが大きく、非ガウス性が高い本形態の目的信号をより有利に分離抽出できる手法である。【０１２４】
［ステップＳ１０３の処理（目的信号選択部１１４０の処理）の詳細］
独立成分分析（ICA）では、信号の独立性に着目してそれらの分離を行うため、得られる分離信号Yp（f,τ）には順序の任意性がある。その順序が入れ替わっても独立性が保たれるからである。従って、目的信号に対応する分離信号を周波数毎に選択しなければならない。目的信号選択部１１４０では、以下に示す処理により、この選択を行う。
図２０は、第３の実施の形態における目的信号選択部１１４０の処理の詳細を説明するためのフローチャートである。以下では、この図を用い、目的信号選択部１１４０の処理の詳細を説明する。【０１２５】
まず、逆行列算出部１１４１において、周波数毎にメモリ１１００の記憶領域１１０３からＭ行Ｍ列の分離行列W（f）を読み出し、それらの逆行列
W（f）-1=[A1（f）,...,AM（f）]（ただし、その列Ap（f）=[A1p（f）,...,AMp（f）]T）…（31）
を算出する。
ここで、式（３０）の両辺に式（３１）を掛け合わせることにより、周波数領域の混合信号X（f,τ）の分解【０１２６】
【数３３】
0004406428.tif 000036
を得ることができる。このAp（f）は基底ベクトルであり、それぞれが各周波数の分離信号Yp（f,τ）毎に対応する。式（３１）のように算出された基底ベクトルAp（f）はメモリ１１００の記憶領域１１０４に格納される（ステップＳ１１１）。
次に、基底ベクトルクラスタリング部１１４２において、全ての基底ベクトルAp（f）（p=1,...,M、f=0,fs/L,…,fs（L-1）/L）を正規化する。この正規化は、複数の信号源から発せられた信号の畳み込み混合を所定のモデル（例えば、近距離場モデル）に近似した場合に、正規化された基底ベクトルAp（f）が信号源の位置のみに依存したクラスタを形成するように行われる。この例では、第１の実施の形態と同様な周波数正規化とノルム正規化とを行う。【０１２７】
周波数正規化は、基底ベクトルクラスタリング部１１４２の周波数正規化部１１４２ａ（図１７Ｂ）において行われる。すなわち、周波数正規化部１１４２ａにおいて、メモリ１１００の記憶領域１１０４の基底ベクトルAp（f）（p=1,...,M、f=0,fs/L,…,fs（L-1）/L）を読み出し、これらを周波数に依存しない周波数正規化ベクトルAp'（f）に正規化してメモリ１１００の記憶領域１１０５に格納する（ステップＳ１１２）。なお、この正規化は基底ベクトルAp（f）の要素Aqp（f）毎に行われる（詳細は後述）。ノルム正規化は、基底ベクトルクラスタリング部１１４２のノルム正規化部１１４２ｂ（図１７Ｂ）において行われる。すなわち、ノルム正規化部１１４２ｂにおいて、メモリ１１００の記憶領域１１０５から周波数正規化ベクトルAp'（f）を読み出し、これらをノルムが規定値（この例では１）をとる正規化基底ベクトルAp''（f）に正規化してメモリ１１００の記憶領域１１０６に格納する（ステップＳ１１３）。なお、この正規化は周波数正規化ベクトルAp'（f）毎に行われる（詳細は後述）。【０１２８】
基底ベクトルの正規化が終了すると、次にクラスタリング部１１４２ｃ（図１７Ｂ）において、正規化基底ベクトルAp''（f）によって形成されるＭ個のクラスタCi（i∈{1,...,M}）を同定する。この例の場合、クラスタリング部１１４２ｃにおいて、メモリ１１００の記憶領域１１０６から正規化基底ベクトルAp''（f）を読み出し、これらをクラスタリングしてＭ個のクラスタＣi（i=1,...,M）を生成し、各クラスタＣiを特定する情報（例えば、各クラスタに属する正規化基底ベクトルAp''（f）を示す情報）と各クラスタＣiのセントロイド（中心ベクトル）ηiとをそれぞれメモリ１１００の記憶領域１１０９，１１１０に格納する（ステップＳ１１４）。なお、このクラスタリングは、例えば、各クラスタＣiの要素（正規化基底ベクトルAv''（f））と各クラスタＣiのセントロイドηiとの間の二乗和Ｕiの総和Ｕ【０１２９】
【数３４】
0004406428.tif 000037
を最小値化することを基準に行われる。この最小値化は、例えば、非特許文献６などで解説されているk-meansクラスタリングを用いることによって効果的に行うことができる。なお、クラスタＣiのセントロイドηiは、【０１３０】
【数３５】
0004406428.tif 000038
によって計算される。ここで｜Ci｜はクラスタCiにおける要素（正規化基底ベクトルAv'' （f））の数であり、‖＊‖はベクトル＊のノルムである。また、ここでは、距離としてユークリッド距離の自乗を用いているが、これを一般化したミンコフスキー距離などを用いてもよい。
Ｍ個のクラスタＣiが得られると、次に分散判定部１１４２ｄ（図１７Ｂ）において、目的信号に対応するクラスタを選択し、選択したクラスタを示す選択情報Ｉ（ｆ）をメモリ１１００の記憶領域１１１１に格納する（ステップＳ１１５）。なお本形態では、クラスタの分散Ui/|Ci|を指標に、目的信号を含む分離信号が選択される。すなわち、本形態における基底ベクトルの正規化は、複数の信号源から発せられた信号の畳み込み混合を所定のモデルに近似した場合に、正規化された基底ベクトルが信号源の位置のみに依存したクラスタを形成するように行われる。しかし、実環境では、このようなモデルに反映されない様々な因子が存在し、このような実環境とモデルとの乖離は、信号源からセンサまでの距離が大きいほど大きくなる。例えば、近距離場モデルでは、信号が壁などに反射することによる反射波成分が考慮されておらず、直接波成分に対する反射波成分の割合は、信号源がセンサから離れるに従って大きくなる。そのため、信号がセンサから離れるに従ってモデルが成り立ちにくくなる。その結果、センサに近い信号ほど実環境に近い条件で正規化ができ、実環境とモデルとの乖離に起因するクラスタの分散を小さくできる。本形態ではセンサに近い信号を目的信号とするため、分散の小さいクラスタを目的信号に対応するクラスタとして選択すればよいことになる。なお、この選択手順（ステップＳ１１５）の詳細については後述する。【０１３１】
各周波数ｆに対する選択情報I（f）が算出されると、次に、各周波数ｆにおける選択信号YI（f）（f,τ）とそれに対応する基底ベクトルAI（f）（f）とを選択する。すなわち、まず、選択部１１４３において、メモリ１１００の記憶領域１１１１から選択情報I（f）を読み出す。そして、選択部１１４３は、この選択情報I（f）に対応する分離信号を選択信号YI（f）（f,τ）として記憶領域１１０７から読み出し、それに対応する基底ベクトルAI（f）（f）を記憶領域１１０４から読み出し、これらを記憶領域１１１１に格納する（ステップＳ１１６）。【０１３２】
次に、ステップＳ１１２，１１３（図２０）の正規化の詳細を示す。
［ステップＳ１１２（周波数正規化）の詳細］
図２１Ａは、ステップＳ１１２の周波数正規化の詳細を説明するためのフローチャートである。
まず、制御部１１７０（図１６）においてパラメータｐに１を代入し、これを一時メモリ１１８０に格納する（ステップＳ１２１）。また、制御部１１７０においてパラメータｑに１を代入し、これを一時メモリ１１８０に格納する（ステップＳ１２２）。次に、周波数正規化部１１４２ａ（図１７Ｂ）において、メモリ１１００の記憶領域１１０８から前述のパラメータｄ，ｃ，Ｑを読み出し、記憶領域１１０４から基底ベクトルAp（f）の要素Aqp（f）を読み出し、一時メモリ１１８０からパラメータｐ，ｑを読み出す。そして、この周波数正規化部１１４２ａは、基底ベクトルAp（f）の要素Aqp（f）に対し、【０１３３】
【数３６】
0004406428.tif 000039
の演算を行い、当該演算結果Aqp’（f）を周波数正規化ベクトルAp’（f）の各要素Aqp’（f）としてメモリ１１００の記憶領域１１０５に格納する（ステップＳ１２３）。なお、ａｒｇ［・］は偏角を意味し、ｅｘｐはネピア数を意味し、ｊは虚数単位を意味する。より詳細には、前述の式（１５）（１６）の演算によって、この正規化を行う。
次に、制御部１１７０において、一時メモリ１１８０に格納されたパラメータｑがｑ＝Ｍを満たすか否かを判断する（ステップＳ１２４）。ここでｑ＝Ｍでなければ、制御部１１７０は、ｑ＋１の演算結果を新たなパラメータｑの値とし、これを一時メモリ１１８０に格納し（ステップＳ１２５）、処理をステップＳ１２３へ戻す。一方、ｑ＝Ｍであれば、制御部１１７０は、さらにｐ＝Ｍを満たすか否かを判断する（ステップＳ１２６）。【０１３４】
ここでｐ＝Ｍでなければ、制御部１１７０は、ｐ＋１の演算結果を新たなパラメータｐの値とし、これを一時メモリ１１８０に格納し（ステップＳ１２７）、処理をステップＳ１２２へ戻す。一方、ｐ＝Ｍであれば、制御部１１７０は、ステップＳ１１２の処理を終了させる（ステップＳ１１２（周波数正規化）の詳細の説明終わり）。
［ステップＳ１１３（ノルム正規化）の詳細］
図２１Ｂは、ステップＳ１１３のノルム正規化の詳細を説明するためのフローチャートである。【０１３５】
まず、制御部１１７０において、パラメータｐに１を代入し、これを一時メモリ１１８０に格納する（ステップＳ１３１）。次に、ノルム正規化部１１４２ｂにおいて、メモリ１１００の記憶領域１１０５から周波数正規化ベクトルAp’（f）の各要素Aqp’（f）を読み出し、【０１３６】
【数３７】
0004406428.tif 000040
の演算を行って、周波数正規化ベクトルAp’（f）のノルム‖Ap’（f）‖を求め、周波数正規化ベクトルAp’（f）及びそのノルム‖Ap’（f）‖を一時メモリ１１８０に格納する（ステップＳ１３２）。
次に、ノルム正規化部１１４２ｂにおいて、一時メモリ１１８０から周波数正規化ベクトルAp’（f）及びそのノルム‖Ap’（f）‖を読み出し、
Ap'' （f）=Ap’（f）／‖Ap’（f）‖    …（39）
の演算を行って正規化基底ベクトルAp'' （f）を求め、これをメモリ１１００の記憶領域１１０６に格納する（ステップＳ１３３）。その後、制御部１１７０において、一時メモリ１１８０に格納されたパラメータｐがｐ＝Ｍを満たすか否かを判断する（ステップＳ１３４）。ここでｐ＝Ｍでなければ、制御部１１７０は、ｐ＋１の演算結果を新たなパラメータｐの値とし、これを一時メモリ１１８０に格納し（ステップＳ１３５）、ステップＳ１３２の処理に戻す。一方、ｐ＝Ｍであれば、制御部１１７０は、ステップＳ１１３の処理を終了させる。なお、正規化基底ベクトルAp'' （f）がクラスタを形成する理由は、第１の実施の形態で示した通りである（ステップＳ１１３（ノルム正規化）の詳細の説明終わり）。【０１３７】
このように生成された正規化基底ベクトルAp''（f）は、第１の実施の形態で示した通り、周波数に依存せず、信号源の位置のみに依存するベクトルとなる。
［選択信号の選択手順（ステップＳ１１５）の詳細］
次に、前述した選択信号の選択手順（ステップＳ１１５）の詳細を例示する。
クラスタの選択手順１：
第１の例は、最も分散が小さいクラスタを目的信号に対応するクラスタとして選択するものである。図２２は、この第１の例を説明するためのフローチャートである。【０１３８】
この場合、まず分散判定部１１４２ｄ（図１７Ｂ）おいて、メモリ１１００の記憶領域１１０９からクラスタＣi（i∈{1,...,M}）を特定する情報を読み出し、さらに記憶領域１１０６，１１１０から、それぞれ正規化基底ベクトルAp''（f）∈Ｃi及びセントロイドηiを読み出す。そして、分散判定部１１４２ｄは、各ｉについて、式（３３）のようにＵｉを求め、Ｃｉに属する要素（正規化基底ベクトルAv'' （f））の数をカウントすることにより｜Ci｜を求め、クラスタＣiの分散Ui/|Ci|を算出して一時メモリ１１８０に格納する。そして、分散判定部１１４２ｄは、この一時メモリ１１８０に格納された分散Ui/|Ci|から最小なものを選択し、それに対応するクラスタを示す情報をクラスタ選択情報
ι=argmini Ui/|Ci| …（40）
として一時メモリ１１８０に格納する（ステップＳ１４１）。なお、式（４０）におけるargmini＊は、＊を最小にするｉを意味する。【０１３９】
次に、制御部１１７０（図１６）において、パラメータｆに０を代入し、これを一時メモリ１１８０に格納する（ステップＳ１４２）。
次に、分散判定部１１４２ｄは、一時メモリ１１８０からクラスタ選択情報ιを読み出し、メモリ１１００の記憶領域１１１０からクラスタ選択情報ιに対応するセントロイドηιを読み出す。また、分散判定部１１４２ｄは、メモリ１１００の記憶領域１１０６から正規化基底ベクトルAp'' （f）（p∈{1,...,M}）を読み出す。そして、分散判定部１１４２ｄは、周波数f毎に選択情報
I（f）=argminp‖Ap''（f）−ηι‖2  …（41）
を算出し、記憶領域１１１１に格納する（ステップＳ１４３）。【０１４０】
次に、制御部１１７０において、一時メモリ１１８０からパラメータｆを読み出し、f=（L-1）・fs/Lであるか否かを判断する（ステップＳ１４４）。ここでf=（L-1）・fs/Lでなければ、制御部１１７０は、このパラメータfにfs/Lを加算した値を新たなfとして一時メモリ１１８０に格納し（ステップＳ１４５）、処理をステップＳ１４３に戻す。一方、f=（L-1）・fs/Lであれば、制御部１１７０は、ステップＳ１１５の処理を終了させる。
クラスタの選択手順２：
第２の例は、予め定められたしきい値よりも分散が小さいクラスタを目的信号に対応するクラスタとして選択するものである。なお、このしきい値は、例えば実験結果又は経験に基づき定められる値であり、事前にメモリ１１００に格納しておくものである。【０１４１】
この場合、ステップＳ１４１（図２２）の処理の代わりに、分散判定部１１４２ｄにおいて、まず各クラスタの分散Ui/|Ci|を公知の並び替えアルゴリズムによって昇順或いは降順に並び替える。そして、分散判定部１１４２ｄは、メモリ１１００に格納されたしきい値を読み取り、分散Ui/|Ci|がしきい値よりも小さくなるクラスタを選択し、選択したクラスタに対応する添字ｉの集合をクラスタ選択情報ιとして一時メモリ１１８０に格納する。それ以外は、クラスタの選択手順１と同様である。
クラスタの選択手順３：
第３の例は、最も分散の小さいクラスタだけではなく、最も分散が小さいクラスタから所定番目に分散が小さいクラスタまで（例えば、最も分散の小さいクラスタから３番目に分散が小さいクラスタまで）を目的信号に対応するクラスタとして選択するものである。【０１４２】
この場合、ステップＳ１４１（図２２）の処理の代わりに、分散判定部１１４２ｄにおいて、まず各クラスタの分散Ui/|Ci|を、公知の並び替えアルゴリズムによって昇順或いは降順に並び替える。そして、分散判定部１１４２ｄは、最も分散Ui/|Ci|が小さいクラスタから所定番目に分散Ui/|Ci|が小さいクラスタまでを選択する。この選択後、分散判定部１１４２ｄは、選択したクラスタに対応する添字ｉの集合をクラスタ選択情報ιとして一時メモリ１１８０に格納する。それ以外は、クラスタの選択手順１と同様である。
なお、クラスタの選択手順１の代わりに分散の小ささが２番目以降のいずれかのクラスタを選択することとしてもよく、上述したクラスタの選択手順の一部を組み合わせた手順を適用してもよい（ステップＳ１１５の説明終わり／ステップＳ１０３の処理（目的信号選択部１１４０の処理）の詳細の説明終わり）。【０１４３】
［時間周波数マスキング部１１５０の処理（ステップＳ１０４，Ｓ１０５）の詳細］
次に、時間周波数マスキング部１１５０の処理について説明する。前述のように、時間周波数マスキング部１１５０は、目的信号選択部１１４０で選択された選択信号YI（f）（f,τ）に残存する妨害信号成分をさらに抑圧するものである。以下ではまず、なぜ選択信号YI（f）（f,τ）に妨害信号が残存するのかについて説明する。
前述した式（３０）を選択信号のみに関して書くと、
YI（f）（f,τ）=WI（f）H（f）・X（f,τ）  …（42）
となる。また、これに式（４）を代入して周波数ｆを省略して書くと、【０１４４】
【数３８】
0004406428.tif 000041
と表記できる。もしＮ≦Ｍであれば、独立成分分析（ICA）により、WIH・Hk=0,∀k∈{1,...,I-1,I+1,...,N}を満たすようなWIを設定することが可能である。そのため、式（４３）の第二項は０になる。しかし、より一般的に信号源数Ｎがセンサ数Ｍよりも多い場合には、WIH・Hk≠0,∀k∈κとなるκ⊆{1,...,I-1,I+1,...,N }が存在してしまう。その場合、選択信号YI（f）は不必要な残留成分（妨害信号の残留成分）【０１４５】
【数３９】
0004406428.tif 000042
を含んでしまう（以後ｆを省略せず）。
時間周波数マスキング部１１５０の目的は、この選択信号YI（f,τ）に含まれる不必要な残留成分を抑圧し、選択信号YI（f,τ）よりも妨害信号の残留成分が少ないマスク処理選択信号YI'（f,τ）を生成することである。そのために、時間周波数マスキング部１１５０は、まずマスク生成部１１５１（図１８）において、妨害信号の残留成分を有する時間周波数スロットに対して小さな値を採り、妨害信号の残留成分が少ない時間周波数スロットに対して大きな値を採る時間周波数マスク０≦M（f,τ）≦１を生成する。そして、マスキング部１１５２において、
YI（f）'（f,τ）=M（f,τ）・YI（f）（f,τ）  …（44）
に従ってマスキング処理を行い、マスク処理選択信号YI（f）'（f,τ）を出力する。以下、この詳細について説明する。【０１４６】
［ステップＳ１０４（マスク生成部１１５１）の処理の詳細］
図２３は、図１９におけるステップＳ１０４の詳細を説明するためのフローチャートである。以下、この図を用い、ステップＳ１０４（マスク生成部１１５１）の処理の詳細を説明する。
この例のマスク生成部１１５１は、周波数領域の混合信号ベクトルX（f,τ）を白色化する空間（白色化空間）において、混合信号ベクトルX（f,τ）と選択信号に対応する基底ベクトルAI（f）（f）とのなす角度θI（f）（f,τ）を求め、その角度θI（f）（f,τ）に基づき、時間周波数マスクを生成する。なお、白色化とは、混合信号ベクトルX（f,τ）を線形に変形し、その共分散行列が単位行列と等しくなるようにすることである。【０１４７】
そのためにまず、白色化行列生成部１１５１ａにおいて、周波数領域の混合信号Xq（f,τ）を用い、混合信号ベクトルX（f,τ）を白色化空間に移す白色化行列V（f）を生成する（ステップＳ１５１）。この例の場合、白色化行列生成部１１５１ａは、メモリ１１００の記憶領域１１０２から混合信号Xq（f,τ）を読み出し、R（f）=<X（f,τ）・X（f,τ）H>τとした場合におけるV（f）=R（f）−1/2を白色化行列V（f）として算出し、これを記憶領域１１１２に格納する。なお、<＊>τはベクトル＊を時間平均したベクトルを意味し、＊Hは＊の複素共役転置を意味し、R−1/2は、R-1/2・R・（R-1/2）H=I（Iは単位行列）を満たす行列を意味する。また、白色化行列V（f）を算出する代表的な方法として、Ｒ（f）を固有値分解R（f）=E（f）・D（f）・E（f）H（E（f）はユニタリ行列、D（f）は対角行列）し、V（f）=D（f）-1/2・E（f）Hの演算を行う方法を例示できる。なお、D（f）-1/2は、対角行列D（f）の各要素の（-1/2）乗した対角行列と等価となるため、対角行列D（f）の各要素の（-1/2）乗することによって算出できる。【０１４８】
次に、白色化部１１５１ｂにおいて、この白色化行列V（f）を用い、混合信号ベクトルX（f,τ）を白色化空間に写像した白色化混合信号ベクトルZ（f,τ）と、基底ベクトルAI（f）（f）を白色化空間に写像した白色化基底ベクトルBI（f）（f）とを算出する（ステップＳ１５２）。この例の場合、まず白色化部１１５１ｂにおいて、メモリ１１００の記憶領域１１０２から混合信号Xq（f,τ）を、記憶領域１１１１から選択信号YI（f）（f,τ）に対応する基底ベクトルAI（f）（f）を、記憶領域１１１２から白色化行列V（f）をそれぞれ読み出す。そして、白色化部１１５１ｂは、Z（f,τ）=V（f）・X（f,τ）の演算により白色化混合信号ベクトルZ（f,τ）を算出し、BI（f）（f）=V（f）・AI（f）（f）の演算により白色化基底ベクトルBI（f）（f）を算出し、これらをメモリ１１００の記憶領域１１１２に格納する。【０１４９】
次に、角度計算部１１５１ｃにおいて、白色化混合信号ベクトルZ（f,τ）と白色化基底ベクトルBI（f）（f）とがなす角度θI（f）（f,τ）を時間周波数毎に算出する（ステップＳ１５３）。この例の場合、まず角度計算部１１５１ｃにおいて、メモリ１１００の記憶領域１１１２から白色化混合信号ベクトルZ（f,τ）と白色化基底ベクトルBI（f）（f）とを読み出す。そして、角度計算部１１５１ｃは、
θI（f）（f,τ）=cos−1（|BI（f）H（f）・Z（f,τ）|／‖BI（f）（f）‖・‖Z（f,τ）‖）  …（45）
の演算により、各時間周波数スロットに対して角度θI（f）（f,τ）を算出し、記憶領域１１１２に格納する。なお、式（４５）における｜＊｜はベクトル＊の絶対値を意味し、‖＊‖はベクトル＊のノルムを意味する。【０１５０】
次に、関数演算部１１５１ｄにおいて、角度θI（f）（f,τ）を要素とする関数である時間周波数マスクM（f,τ）を生成する（ステップＳ１５４）。この例の場合、関数演算部１１５１ｄにおいて、まず、メモリ１１００の記憶領域１１０８から実数パラメータθＴ，ｇを、記憶領域１１１２から角度θI（f）（f,τ）を、それぞれ読み出す。そして、関数演算部１１５１ｄは、ロジスティック関数
M（θ（f,τ））=１／（1+eg・（θ（f,τ）−θT））   …（46）
を時間周波数マスクM（f,τ）として算出する（θＴはθＴを意味する）。なお、実数パラメータθＴ，ｇは、それぞれ時間周波数マスクM（f,τ）の転移点と傾き度を指定するパラメータであり、前処理によって記憶領域１１０８に格納しておいたものである。図２４Ａに、二種類の実数パラメータθＴ，ｇに対し、式（４６）に従って算出された時間周波数マスクM（f,τ）を例示する。この図に例示すように、実数パラメータθＴが小さいほど、時間周波数マスクM（f,τ）が大きな値（この例では１）をとる領域が狭くなる。これは実数パラメータθＴが小さいほど、マスク処理選択信号YI（f）'（f,τ）に残存する妨害信号成分が少なくなるが、その反面、マスク処理選択信号YI（f）'（f,τ）が不自然なものとなるといった傾向に現れる。例えば、音声信号を目的信号とする場合、実数パラメータθＴを小さくするほどミュージカルノイズ（musical noise）が大きくなる。また、実数パラメータｇが大きくなるほど時間周波数マスクM（f,τ）の波形（大きな値（この例では１）から小さな値（この例では０）への遷移）は急峻になる。ここで、できるだけマスク処理選択信号YI（f）'（f,τ）に残存する妨害信号成分を減らし、マスク処理選択信号YI（f）'（f,τ）を自然な信号とするためには実数パラメータｇを小さくし、時間周波数マスクM（f,τ）の波形を滑らかなものとすることが望ましい。
【０１５１】
なお、これらの実数パラメータθＴ，ｇを周波数毎に設定しておく構成としてもよい。また、さらに実数パラメータαを導入し、ロジスティック関数
M（θ（f,τ））=α／（1+eg・（θ（f,τ）−θT））   …（47）
を時間周波数マスクM（f,τ）としてもよい。さらに、角度θI（f）（f,τ）が０に近い領域において大きな値を採り、角度θI（f）（f,τ）が大きな領域において小さな値を採る０≦M（θ（f,τ））≦１の関数であれば、他の関数を時間周波数マスクM（f,τ）として用いてもよい（ステップＳ１０４（マスク生成部１１５１）の処理の詳細の説明終わり）。【０１５２】
［ステップＳ１０５（マスキング部１１５２）の処理の詳細］
マスキング部１１５２は、メモリ１１００の記憶領域１１１１から選択信号YI（f）（f,τ）を、記憶領域１１１２から時間周波数マスクM（f,τ）をそれぞれ読み出す。そして、マスキング部１１５２は、
YI（f）'（f,τ）=M（f,τ）・YI（f）（f,τ） …（48）
の演算によってマスク処理選択信号YI（f）'（f,τ）を求め、これをメモリ１１００の記憶領域１１１３に格納する（ステップＳ１０５（マスキング部１１５２）の処理の詳細の説明終わり）。【０１５３】
［時間周波数マスキングの効果］
次に、上述した時間周波数マスクM（f,τ）によるマスキングの効果について説明する。
信号源のスパース性が高く、信号源Sk（f,τ）がゼロに近くなる可能性が高いなら、式（４）は、【０１５４】
【数４０】
0004406428.tif 000043
として近似される。なお、式（４９）におけるｋは各信号源に対応する添字であり、時間周波数位置（f,τ）毎に決まる。よって、目的信号が唯一アクティブ或いはそれに近い時間周波数位置（f,τ）では、白色化混合信号ベクトルZ（f,τ）は、【０１５５】
【数４１】
0004406428.tif 000044
と近似できる。ここで、YI（f）（f,τ）はスカラ値である。また、前述のように白色化基底ベクトルBI（f）（f）は、
BI（f）（f）=V（f）・AI（f）（f） …（50）
である。以上より、目的信号が唯一アクティブ或いはそれに近い時間周波数位置（f,τ）では、白色化混合信号ベクトルZ（f,τ）と白色化基底ベクトルBI（f）（f）とがなす角度θI（f）（f,τ）は０に近くなることが分かる。そして、前述のように時間周波数マスクM（f,τ）は、角度θI（f）（f,τ）が０に近い領域において大きな値を採るのである。よって、時間周波数マスクM（f,τ）は、目的信号が唯一アクティブ或いはそれに近い時間周波数位置（f,τ）における選択信号YI（f）（f,τ）をマスク処理選択信号YI（f）'（f,τ）として抽出するものといえる（式（４８）参照）。【０１５６】
一方、仮にI（f）=1とおいた場合、目的信号が殆どアクティブでない時間周波数位置（f,τ）では、白色化混合信号ベクトルZ（f,τ）は、【０１５７】
【数４２】
0004406428.tif 000045
と近似できる。ここで、信号源のＮが、センサの数Ｍと同等かそれより少ないならば、白色化空間でのベクトルV（f）・H1（f）,...,V（f）・Hk（f）は互いに直交する。なお、式（５１）のSk（f,τ）はスカラ値である。これより、白色化混合信号ベクトルZ（f,τ）と白色化基底ベクトルBI（f）（f）とがなす角度θI（f）（f,τ）が大きくなることがわかる。またＮ＞Ｍでも、白色化基底ベクトルB1（f）（I（f）=1）は、目的信号以外のベクトルV（f）・H2（f）,...,V（f）・Hk（f）とは大きい角度を持つ傾向がある。以上より、目的信号が殆どアクティブでない時間周波数位置（f,τ）では、角度θI（f）（f,τ）が大きな値を採ることが分かる。そして、時間周波数マスクM（f,τ）は、角度θI（f）（f,τ）が０から離れた領域において小さな値を採るのであるから、目的信号が殆どアクティブでない時間周波数位置（f,τ）における選択信号YI（f）（f,τ）をマスク処理選択信号YI（f）'（f,τ）から排除するものといえる（式（２８）参照）。【０１５８】
以上より、時間周波数マスクM（f,τ）による時間周波数マスキング処理は、選択信号YI（f）（f,τ）に残存する妨害信号成分をより一層抑圧するものであることがわかる。
なお、この時間周波数マスキング処理は、音声や音楽などスパース性を持つ信号に対して特に有効である。スパース性が低い場合、目的信号がアクティブである時間周波数位置（f,τ）であっても他の妨害信号成分も多く含まれ、式（４９）の近似が成り立たなくなり、角度θI（f）（f,τ）が０から離れてしまうからである。すなわちスパースでない場合、例えば、図２４Ｂのように、ある時間周波数位置（f,τ）において目的信号（I（f）=1とする）に対応するベクトルV（f）・H1（f）の他に、妨害信号に対応するベクトルV（f）・H2（f），V（f）・H3（f）が共存することになる。この例の場合、白色化混合信号ベクトルZ（f,τ）は、【０１５９】
【数４３】
0004406428.tif 000046
となるため、白色化混合信号ベクトルZ（f,τ）と白色化基底ベクトルBI（f）（f）とがなす角度θI（f）（f,τ）も０から離れてしまう。これは目的信号がアクティブである時間周波数位置（f,τ）の信号をマスク処理選択信号YI（f）'（f,τ）から排除してしまう可能性があることを示している。
また、この時間周波数マスキング処理は、目的信号のパワーが妨害信号のパワーに比べて十分に大きい場合に特に有効である。すなわち、目的信号のパワーが妨害信号のパワーに比べ十分に大きい場合、たとえスパース性が低く、目的信号がアクティブである時間周波数位置（f,τ）に他の妨害信号成分も含まれる状況であっても、式（４９）の近似が比較的成り立ち、角度θI（f）（f,τ）が０に近づくからである。例えば、目的信号のパワーが妨害信号のパワーに比べ十分に大きい場合、式（５２）における妨害信号の寄与度は低くなり、白色化混合信号ベクトルZ（f,τ）と白色化基底ベクトルBI（f）（f）とがなす角度θI（f）（f,τ）も０に近づく。これは、目的信号がアクティブである時間周波数位置（f,τ）の信号が、マスク処理選択信号YI（f）'（f,τ）から排除されてしまう可能性を低くできることを示している。また、マスク処理選択信号YI（f）'（f,τ）に尚も残存する妨害信号成分を相対的に低くできることを意味している（ステップＳ１０５（マスキング部１１５２）の処理の詳細の説明終わり）。【０１６０】
〔第４の実施の形態（「第２の本発明」の例）〕
次に、本発明における第４の実施の形態について説明する。
本形態は、第３の実施の形態の変形例であり、時間周波数マスクによる処理を行わない点のみが第３の実施の形態と相違する。以下では、第３の実施の形態との相違点を中心に説明し、第３の実施の形態と共通する事項については説明を省略する。
＜構成＞
図２５は、第４の実施の形態における信号分離装置１２００のブロック図の例示である。【０１６１】
図２５に例示するように、本形態の信号分離装置１２００の信号分離装置１００１との相違点は、メモリ１１００に記憶領域１１１２，１１１３を有しない点、及び時間周波数マスキング部１１５０を有しない点である。
＜処理＞
次に、本形態における信号分離装置１２００の処理について説明する。
図２６は、第４の実施の形態における信号分離装置１２００の処理を説明するためのフローチャートである。以下では、第３の実施の形態との相違点を中心に説明する。【０１６２】
まず、第３の実施の形態と同様に、周波数領域変換部１１２０において、メモリ１１００の記憶領域１１０１から時間領域の混合信号xq（t）を読み出す。そして、周波数領域変換部１１２０は、これらを短時間フーリエ変換などにより周波数領域の混合信号Xq（f,τ）に変換し、メモリ１１００の記憶領域１１０２に格納する（ステップＳ１６１）。
次に、信号分離部１１３０において、メモリ１１００の記憶領域１１０２から周波数領域の混合信号Xq（f,τ）を読み出す。この例の信号分離部１１３０は、読み出した混合信号Xq（f,τ）を要素とする混合信号ベクトルX（f,τ）=[X1（f,τ）,…,XM（f,τ）]Tに、独立成分分析（ICA）を適用し、周波数ｆ毎にＭ行Ｍ列の分離行列W（f）=[W1（f）,...,WM（f）]H（＊Hは行列＊の複素共役転置行列）と分離信号ベクトルY（f,τ）=W（f）・X（f,τ）とを算出する（ステップＳ１６２）。なお、算出された分離行列W（f）はメモリ１１００の記憶領域１１０３に格納される。また、分離信号ベクトルY（f,τ）=[Y1（f,τ）,…,YM（f,τ）]Tを構成する分離信号Yp（f,τ）（p∈{1,..,M}）は記憶領域１１０７に格納される。【０１６３】
次に、目的信号選択部１１４０において、メモリ１１００の記憶領域１１０３から分離行列W（f）を読み出し、その一般化逆行列の列である基底ベクトルを正規化し、その正規化された基底ベクトルをクラスタリングする。そして、目的信号選択部１１４０は、そのクラスタの分散を指標として、メモリ１１００の記憶領域１１０７の分離信号から選択信号YI（f）（f,τ）を周波数ｆ毎に選択し、メモリ１１００の記憶領域１１１１に格納する（ステップＳ１６３）。
次に、時間領域変換部１１６０において、メモリ１１００の記憶領域１１１１から選択された分離信号YI（f）（f,τ）を読み出し、それに短時間逆フーリエ変換などを施して時間領域の分離信号y（t）を生成し、メモリ１１００の記憶領域１１１４に格納する（ステップＳ１６４）。【０１６４】
〔第５の実施の形態（「第２の本発明」の例）〕
次に、本発明における第５の実施の形態について説明する。
本形態は、第３の実施の形態の変形例であり、時間周波数マスクの生成方法のみが第３の実施の形態と相違する。以下では、第３の実施の形態との相違点を中心に説明し、第３の実施の形態と共通する事項については説明を省略する。
＜構成＞
図２７は、第５の実施の形態における信号分離装置１３００のブロック図の例示である。また、図２８Ａは、図２７における時間周波数マスキング部１３５０の詳細構成を示すブロック図であり、図２８Ｂは、図２８Ａのマスク生成部１３５１の詳細構成を示すブロック図である。なお、これらの図において第３の実施の形態と共通する部分については、第３の実施の形態の図と同じ符号を付した。【０１６５】
図２７に例示するように、本形態の信号分離装置１３００の信号分離装置１００１との相違点は、信号分離装置１３００が、時間周波数マスキング部１１５０の代わりに時間周波数マスキング部１３５０を有する点、メモリ１１００が記憶領域１１０８，１１１２の代わりに記憶領域１３０８，１３１２を有する点である。また、図２８Ａに示すように、時間周波数マスキング部１３５０は、マスク生成部１３５１及びマスキング部１１５２を有する。また、図２８Ｂに示すように、マスク生成部１３５１は、周波数正規化部１３５１ａ、ノルム正規化部１３５１ｂ、セントロイド抽出部１３５１ｃ、自乗距離算出部１３５１ｄ及び関数生成部１３５１ｅを有している。また、周波数正規化部１３５１ａは、第１正規化部１３５１ａａ及び第２正規化部１３５１ａｂを有しており、セントロイド抽出部１３５１ｃは、セントロイド選択部１３５１ｃａ及びノルム正規化部１３５１ｃｂを有している。【０１６６】
＜マスク生成処理＞
本形態と第３の実施の形態との相違点は、時間周波数マスクの生成処理（ステップＳ１０４）のみである。以下では、本形態の時間周波数マスクの生成処理のみを説明する。
図２９は、第５の実施の形態における時間周波数マスクの生成処理を説明するためのフローチャートである。また、図３０Ａは、図２９におけるステップＳ１７１の詳細を説明するためのフローチャートである。また図３０Ｂは、図２９におけるステップＳ１７２の詳細を説明するためのフローチャートである。以下、これらの図を用い、本形態の時間周波数マスクの生成処理を説明する。【０１６７】
まず、マスク生成部１３５１の周波数正規化部１３５１ａが、メモリ１１００の記憶領域１１０２に格納されている周波数領域の混合信号Xq（f,τ）から構成される混合信号ベクトルX（f，τ）を、周波数に依存しない周波数正規化ベクトルX’（f，τ）に正規化し（周波数正規化）、当該周波数正規化ベクトルX’（f，τ）の各要素Xq'（f,τ）をメモリ１１００の記憶領域１３１２に格納する（ステップＳ１７１）。
［周波数正規化処理（ステップＳ１７１）の詳細］
以下に、この周波数正規化処理（ステップＳ１７１）の詳細を説明する。【０１６８】
まず、制御部１１７０（図２７）においてパラメータｑに１を代入し、これを一時メモリ１１８０に格納する（ステップＳ１８１）。次に、周波数正規化部１３５１ａ（図２８Ｂ）において、メモリ１１００の記憶領域１３０８から前述のパラメータｄ，ｃ，Ｑを読み込み、記憶領域１１０２から各（f,τ）に対応する混合信号ベクトルX（f,τ）の要素Xq（f,τ）を読み込み、一時メモリ１１８０からパラメータｑを読み込む。そして、周波数正規化部１３５１ａは、【０１６９】
【数４４】
0004406428.tif 000047
の演算を行い、当該演算結果を周波数正規化ベクトルX'（f,τ）=[X1'（f,τ）,...,XM'（f,τ）]Tの各要素として、メモリ１１００の記憶領域１３１２に格納する（ステップＳ１８２）。なお、ａｒｇ［・］は偏角を意味し、ｊは虚数単位を意味する。
より詳細には、まず、周波数正規化部１３５１ａの第１正規化部１３５１ａａが、混合信号ベクトルX（f,τ）の特定の１つの要素Xq（f,τ）を基準として当該混合信号ベクトルX（f,τ）の各要素Xq（f,τ）の偏角を、以下の演算によって正規化する。【０１７０】
【数４５】
0004406428.tif 000048
次に、周波数正規化部１３５１ａの第２正規化部１３５１ａｂが、第１正規化部１３５１ａａで正規化された各要素Xq'''（f,τ）の偏角を、以下のように周波数ｆに比例した値4fc-1dで除算する。【０１７１】
【数４６】
0004406428.tif 000049
次に、制御部１１７０において、一時メモリ１１８０に格納されたパラメータｑがｑ＝Ｍを満たすか否かを判断する（ステップＳ１８３）。ここでｑ＝Ｍでなければ、制御部１１７０は、ｑ＋１の演算結果を新たなパラメータｑの値とし、これを一時メモリ１１８０に格納し（ステップＳ１８４）、処理をステップＳ１８２へ戻す。一方、ｑ＝Ｍであれば、制御部１１７０はステップＳ１７１の処理を終了させ、以下のステップＳ１７２の処理を実行させる（周波数正規化処理（ステップＳ１７１）の詳細終わり）。【０１７２】
次に、マスク生成部１３５１のノルム正規化部１３５１ｂが、メモリ１１００の記憶領域１３１２に格納されている各要素Xq'（f,τ）からなる周波数正規化ベクトルX’（f，τ）を、ノルムが所定の値（この例では１）になるノルム正規化ベクトルX’’（f，τ）に正規化し（ノルム正規化）、それら各要素Xq’’（f，τ）を記憶領域１３１２に格納する（ステップＳ１７２）。
［ノルム正規化処理（ステップＳ１７２）の詳細］
以下に、このノルム正規化処理（ステップＳ１７２）の詳細を説明する。【０１７３】
まず、ノルム正規化部１３５１ｂ（図２８Ｂ）において、メモリ１１００の記憶領域１３１２から、各（f,τ）に対応する周波数正規化ベクトルX'（f,τ）=[X1'（f,τ）,...,XM'（f,τ）]Tを読み込む。そして、ノルム正規化部１３５１ｂは、【０１７４】
【数４７】
0004406428.tif 000050
の演算を行ってそれらのノルム‖X'（f,τ）‖を求め、周波数正規化ベクトルX'（f,τ）及びノルム‖X'（f,τ）‖を一時メモリ１１８０に格納する（ステップＳ１８５）。
次に、ノルム正規化部１３５１ｂは、一時メモリ１１８０から各（f,τ）に対応する周波数正規化ベクトルX'（f,τ）及びノルム‖X'（f,τ）‖を読み出し、
X''（f,τ）=X'（f,τ）／‖X'（f,τ）‖
の演算を行ってノルム正規化ベクトルX''（f,τ）を求める（ステップＳ１８６）。【０１７５】
求められたノルム正規化ベクトルX''（f,τ）は、メモリ１１００の記憶領域１３１２に格納される。これによりステップＳ１７２の処理が終了する（ノルム正規化処理（ステップＳ１７２）の詳細の説明終わり）。
次に、セントロイド抽出部１３５１ｃのセントロイド選択部１３５１ｃａが、一時メモリ１１８０からクラスタ選択情報ιを読み込み（ステップＳ１４１参照）、このクラスタ選択情報ιに対応するセントロイドηιをメモリ１１００の記憶領域１１１０から読み込む（ステップＳ１７３）。次に、ノルム正規化部１３５１ｃｂが、セントロイド選択部１３５１ｃａが読み込んだセントロイドηιのノルムを所定の値（ステップＳ１７２の所定の値〔この例では１〕）に正規化する。このノルム正規化後のセントロイドηιをノルム正規化セントロイドηι’と呼ぶ（ステップＳ１７４）。なお、ノルム正規化の手順は、ステップＳ１８５，Ｓ１８６と同様である。また、ノルム正規化セントロイドηι’は、メモリ１１００の記憶領域１３１２に格納される。【０１７６】
次に、自乗距離算出部１３５１ｄが、メモリ１１００の記憶領域１３１２から、ノルム正規化ベクトルX''（f,τ）とノルム正規化セントロイドηι’とを読み込み、それらの距離の自乗
DS（f,τ）=‖ηι’-X''（f,τ）‖2
を算出し（ステップＳ１７５）、その距離の自乗DS（f,τ）をメモリ１１００の記憶領域１３１２に格納する。
次に、関数生成部１３５１ｅが、メモリ１１００の記憶領域１３１２から距離の自乗DS（f,τ）を読み込み、これを変数とする関数を用いて時間周波数マスクM（f,τ）を生成し、メモリ１１００の記憶領域１３１２に格納する（ステップＳ１７６）。具体的には、関数生成部１３５１ｅは、例えば、メモリ１１００の記憶領域１３０８から実数パラメータｇ，ＤＴを読み込み、以下のロジスティック関数を時間周波数マスクM（DS（f,τ））として生成する。なお、ＤＴは事前に記憶領域１３０８に格納しておいたものである。また、ｅはネピア数である。【０１７７】
【数４８】
0004406428.tif 000051
そして、以上のように生成された時間周波数マスクM（DS（f,τ））は、第３の実施の形態と同様、マスキング部１１５２でのマスキング処理に用いられる。  
〔実験結果〕
第３，第４の実施の形態の効果を示すために、マイクロホンに近い主要音声を強調して抽出する実験を行った。この実験では、図３１Ａに示す実験条件でインパルス応答hqk（r）を測定した。ここでは、カクテルパーティの状況を再現するようにスピーカを配置した。また、全てのスピーカの音量はほぼ等しくなるようにし、特別なスピーカの音量のみが著しく大きくならないようにした。マイクロホンでの混合は、８ｋＨｚでサンプリングした６秒間の英語音声を、測定したインパルス応答と畳み込むことにより作成した。マイクロホン配置は、図３１Ａに示すように３次元的にしたが、信号分離装置を実装したシステム（装置）には、基準マイクロホン（Mic.2）とその他のマイクロホン間の最大距離（3.5cm）のみを情報として与え、詳細なマイクロホン配置は与えていない。各実験において、マイクロホンに近い４つのスピーカ位置（a120,b120,c120,c170）の中から１つを目的音源として選び、残りの３つのスピーカは無音状態にした。マイクロホンから離れた６個のスピーカは、各実験で常に妨害音として作動させた。抽出結果はInput SIR−Output SIRであるSIR（signal−to−interference ratio）の改善量によって評価した。この値が大きいほど、より良く目的音声を抽出し、他の妨害音をより抑圧できたことを意味する。これら２種類のSIRは、【０１７８】
【数４９】
0004406428.tif 000052
はsk（t）からy1（t）へのインパルス応答である。
実験は、各目的音源位置に対して、７つの音声（目的音声１個と妨害音声６個）からなる組合せを１６種類作成して行った。図３１Ｂは、ICAのみの場合（第４の実施の形態）と、ICAと時間周波数マスキングの両方を組み合わせた場合（第３の実施の形態）とのSIRの平均改善量を示す表である。SIRの改善量は、目的音源の位置によって多少左右されるが、おおむね良好な結果を得ている。位置a120とb120とでの結果が良いのは、妨害音が異なった位置から到来するからである。位置c120とc170は、２次元的に見ると、同様の方向から多くの妨害音が到来するため困難な位置に思えるが、実際に位置c170での結果は大変良かった。これは、cl70の高さが妨害音と異なり、3次元状に配列されたマイクロホンによって、本システムが高さの違いを自動的に利用して抽出処理を行っているからである。また、図３１Ｂの表によると、時間周波数マスキングにより性能が改善されることが分かる。時間周波数マスクを決定する式（４６）に対して図３１Ａに示す３種類のパラメータを用いた。より小さいθTを用いることで、より大きいSIR改善量が達成されている。しかしながら、小さいθTによるいくつかの音は、不自然な音（musical noise）を伴った。実験の結果、多くの場合においてパラメータ（θT,g）=（0.333π,20）が、十分に妨害音の抑圧を行いながら、かつ自然な音を出力することを確認した。【０１７９】
〔変形例等〕
なお、本発明は上述の第３から第５の実施の形態に限定されるものではない。例えば、上述した各実施の形態では、信号分離部１１３０においてＭ行Ｍ列の分離行列W（f）を算出することとしたが、Ｎ行Ｍ列など正方行列以外の分離行列W（f）を算出する構成としてもよい。この場合、基底ベクトルは、分離行列W（f）の一般化逆行列W+（f）（例えばムーア・ペンローズ（Moore-Penrose）型一般化逆行列）の各列となる。
また、第３の実施の形態では、時間周波数マスクを用いて選択信号YI（f）（f,τ）から妨害信号成分をさらに抑制したマスク処理選択信号YI（f）'（f,τ）を生成することとしたが、他の方法によって妨害信号成分を抑制したマスク処理選択信号YI（f）'（f,τ）を生成することとしてもよい。例えば、信号源が２個だけである場合、抽出された分離信号Y1（f,τ），Y2（f,τ）の大きさを比較し、｜Y1（f,τ）｜＞｜Y2（f,τ）｜ならばY1（f,τ）をマスク処理選択信号YI（f）'（f,τ）として抽出し、｜Y1（f,τ）｜＜｜Y2（f,τ）｜ならばY2（f,τ）をマスク処理選択信号YI（f）'（f,τ）として抽出する時間周波数マスクを生成し、分離信号Y1（f,τ），Y2（f,τ）を要素とするベクトルに掛け合わせることとしてもよい。【０１８０】
さらに、上述の第３の実施の形態では、信号分離部１１３０が独立成分分析（ICA）を用いて分離行列と分離信号とを算出することとしたが、時間周波数マスク（時間周波数毎のマスクを意味し、例えば、１か０の値をとるバイナリマスク）によって観測信号から分離信号を抽出し（例えば、「O. Yilmaz and S. Rickard, "Blind separation of speech mixtures via time-frequency masking," IEEE Trans. an SP, vol. 52, no. 7, pp. 1830-1847, 2004.」参照）、その算出結果から分離行列を生成することとしてもよい。
また、上述の第３の実施の形態では、周波数正規化部１１４２ａの第１正規化部１１４２ａａが、基底ベクトルAp（f）の特定の１つの要素AQp（f）を基準とし、式（３５）の一部である式（１５）によって当該基底ベクトルAp（f）の各要素Aqp（f）の偏角を正規化することとした。しかし、例えば、前述した式（27-1）〜（27-3）等に従って、第１正規化部１１４２ａａが、基底ベクトルAp（f）の特定の１つの要素AQp（f）を基準として、当該基底ベクトルAp（f）の各要素Aqp（f）の偏角を正規化する構成であってもよい。【０１８１】
また、周波数正規化部１１４２ａが、式（３５）の代わりに、前述した式（28-1）〜（28-4）等の演算によって周波数正規化を行うこととしてもよい。
また、上述の第３の実施の形態では、ノルム正規化部１１４２ｂにおいてノルムが１になるような正規化を行ったが、ノルムが１以外の規定値となるように正規化を行ってもよい。さらに、ノルム正規化部１１４２ｂを設けず、ノルム正規化を行わない構成であってもよい。この場合のクラスタリングは、前述した通り、ベクトルの方向の類似度を基準にして行う。【０１８２】
また、パラメータｄは、全てのセンサｑに対して同一であってもよいし、各センサｑに対応して複数の値を設定しておいてもよい。例えば、基準センサとセンサｑとの距離を、それぞれのセンサｑに対応するパラメータｄの値としてもよい。
〔第６の実施の形態（「第３の本発明」の例）〕
次に、本発明における第６の実施の形態について説明する。
本形態は、前述の原理を用い、センサの厳密な配置情報を必要とすることなく、全ての観測信号から得られた情報を容易かつ効率的に利用して信号分離を行う形態である。本形態では、後述する「混合信号ベクトル」が前述の「複素ベクトル」に相当する。【０１８３】
＜構成＞
本形態の信号分離装置２００１も、第１の実施の形態と同様、公知のノイマン型のコンピュータに信号分離プログラムが読み込まれることにより構成される。図３２は、第６の実施の形態における信号分離装置２００１の構成を例示したブロック図である。また、図３３は、図３２における信号分離部２１２０の詳細を例示したブロック図である。なお、これらの図における実線の矢印は実際のデータの流れを示し、破線の矢印は理論的な情報の流れを示す。また、これらの図において制御部２１４０に出入りするデータの流れに対応する矢印は省略してある。【０１８４】
図３２，図３３に例示するように、信号分離装置２００１は、メモリ２１００、周波数領域変換部２１１０（「複素ベクトル生成部」に相当する機能を含む）、信号分離部２１２０、時間領域変換部２１３０及び制御部２１４０を有している。また、この信号分離部２１２０は、周波数正規化部２１２１（「正規化部」を構成）、ノルム正規化部２１２２（「正規化部」を構成）、クラスタリング部２１２３及び分離信号生成部２１２４を有している。また、周波数正規化部２１２１は、第１正規化部２１２１ａ及び第２正規化部２１２１ｂを有している。また、制御部２１４０は一時メモリ２１４１を有している。【０１８５】
ここでメモリ２１００及び一時メモリ２１４１は、レジスタ１０ａｃ、補助記憶装置１０ｆ及びＲＡＭ１０ｄ等に相当する。また、周波数領域変換部２１１０、信号分離部２１２０、時間領域変換部２１３０及び制御部２１４０は、ＣＰＵ１０ａにＯＳプログラムや信号分離プログラムが読み込まれ、ＣＰＵ１０ａがこれらを実行することにより構成されるものである。
＜処理＞
次に、本形態の信号分離装置２００１の処理について説明する。なお、以下では、Ｎ個の源信号が混合され、Ｍ個のセンサで観測された状況を取り扱う。また、前処理において、各センサで観測された時間領域の混合信号xq（t）（q=1,...,M）がメモリ２１００の記憶領域２１０１に格納され、信号の伝達速度ｃ、Ｍ以下の自然数から選択された基準値Ｑ，Ｑ’（それぞれＭ個のセンサから選択された基準センサの添字番号）及び実数ｄの各パラメータが記憶領域２１０５に格納されているものとする。【０１８６】
図３４は、本形態における信号分離装置２００１の処理の全体を説明するためのフローチャートである。以下、この図に沿って、本形態における信号分離装置２００１の処理を説明していく。
［処理の全体］
まず、周波数領域変換部２１１０において、メモリ２１００の記憶領域２１０１から時間領域の混合信号xq（t）を読み出し、これらを短時間離散フーリエ変換等によって周波数毎の時系列信号（「周波数領域の混合信号」と呼ぶ）Xq（f,τ）（q=1,...,M、f=0,fs/L,…,fs（L-1）/L、ｆｓはサンプリング周波数）に変換し、メモリ２１００の記憶領域２１０２に格納する（ステップＳ２０１）。【０１８７】
次に、信号分離部２１２０の周波数正規化部２１２１において、メモリ２１００の記憶領域２１０２から周波数領域の混合信号Xq（f,τ）を読み出す。周波数領域の混合信号Xq（f,τ）を読み出した周波数正規化部２１２１は、これらからなる混合信号ベクトルX（f,τ）=[X1（f,τ）,...,XM（f,τ）]Tを、周波数ｆに依存しない周波数正規化ベクトルX’（f,τ）に正規化する（ステップＳ２０２）。生成された各周波数正規化ベクトルX’（f,τ）はメモリ２１００の記憶領域２１０３に格納される。なお、ステップＳ２０２の処理の詳細については後述する。【０１８８】
次に、信号分離部２１２０のノルム正規化部２１２２において、メモリ２１００の記憶領域２１０３から各周波数正規化ベクトルX’（f,τ）を読み込み、これらをノルムが所定の値（例えば１）になるノルム正規化ベクトルX''（f,τ）に正規化する。そして、ノルム正規化部２１２２は、生成した各ノルム正規化ベクトルX''（f,τ）をメモリ２１００の記憶領域２１０４に格納する（ステップＳ２０３）。なお、この詳細については後述する。
次に、信号分離部２１２０のクラスタリング部２１２３において、メモリ２１００の記憶領域２１０４から各ノルム正規化ベクトルX''（f,τ）を読み込み、これらをクラスタリングし、クラスタを生成する。そして、クラスタリング部２１２３は、各クラスタを特定するクラスタ情報Ｃｋ（ｋ（ｋ=1,...,N）番目のクラスタのメンバX''（f,τ）を特定する情報）をメモリ２１００の記憶領域２１０６に格納する（ステップＳ２０４）。なお、この詳細については後述する。【０１８９】
次に、信号分離部２１２０の分離信号生成部２１２４において、メモリ２１００の記憶領域２１０６，２１０５からクラスタ情報Ｃｋ及び基準値Ｑ’を読み込む。そして、分離信号生成部２１２４は、クラスタ情報Ｃｋ及び基準値Ｑ’を用い、ｋ番目のクラスタに属するノルム正規化ベクトルX''（f,τ）に対応する混合信号ベクトルX（f,τ）からＱ’番目の要素XＱ’（f,τ）を記憶領域２１０２から抽出し、これをｋ番目の要素Yk（f,τ）とした分離信号ベクトルY（f,τ）を生成する。そして、分離信号生成部２１２４は、生成した分離信号ベクトルY（f,τ）をメモリ２１００の記憶領域２１０７に格納する（ステップＳ２０５）。なお、この詳細については後述する。【０１９０】
最後に、時間領域変換部２１３０において、メモリ２１００の記憶領域２１０７から分離信号ベクトルY（f,τ）を読み込み、その分離信号成分Yk（f,τ）を添字ｋ毎に短時間逆フーリエ変換等により、時間領域の分離信号ｙｋ（ｔ）に変換する。そして、時間領域変換部２１３０は、変換した時間領域の分離信号ｙｋ（ｔ）をメモリ２１００の記憶領域２１０８に格納する（ステップＳ２０６）。
次に、各処理の詳細について説明する。
［周波数正規化部２１２１・ノルム正規化部２１２２の処理の詳細］
周波数正規化部２１２１及びノルム正規化部２１２２は、全ての混合信号ベクトルX（f,τ）=[X1（f,τ）,...,XM（f,τ）]T（f=0,fs/L,…,fs（L-1）/L）を正規化し、それらを周波数に依存せず、信号源の位置のみに依存するノルム正規化ベクトルX''（f,τ）に正規化する。これにより、それらをステップＳ２０４でクラスタリングした際に、各クラスタが各信号源のみに対応するようになる。この正規化を適切に行わないとクラスタが形成されない。前述した通り、本形態における正規化は、周波数正規化とノルム正規化とからなる。周波数正規化は、周波数正規化部２１２１において、混合信号ベクトルX（f,τ）を周波数に依存しない周波数正規化ベクトルX’（f,τ）に正規化するものである。ノルム正規化は、ノルム正規化部２１２２において、周波数正規化ベクトルX’（f,τ）をノルムが規定値（この例では１）をとるノルム正規化ベクトルX''（f,τ）に正規化するものである。以下にこれらの正規化の詳細を説明する。【０１９１】
［周波数正規化部２１２１における処理の詳細（ステップＳ２０２の処理の詳細）］
図３５Ａは、図３４に示したステップＳ２０２の処理の詳細を説明するためのフローチャートである。以下、この図に沿ってステップＳ２０２の処理の詳細を説明する。
まず、制御部２１４０（図３２）においてパラメータｑに１を代入し、これを一時メモリ２１４１に格納する（ステップＳ２１１）。次に、周波数正規化部２１２１（図３２，図３３）において、メモリ２１００の記憶領域２１０５から前述のパラメータｄ，ｃ，Ｑを読み込み、記憶領域２１０２から各（f,τ）に対応する混合信号ベクトルX（f,τ）の要素Xq（f,τ）を読み込み、一時メモリ２１４１からパラメータｑを読み込む。そして、周波数正規化部２１２１は、【０１９２】
【数５０】
0004406428.tif 000053
の演算を行い、当該演算結果を周波数正規化ベクトルX'（f,τ）=[X1'（f,τ）,...,XM'（f,τ）]Tの各要素として、メモリ２１００の記憶領域２１０３に格納する（ステップＳ２１２）。なお、ａｒｇ［・］は偏角を意味し、ｊは虚数単位を意味する。
より詳細には、まず、周波数正規化部２１２１の第１正規化部２１２１ａが、混合信号ベクトルX（f,τ）の特定の１つの要素Xq（f,τ）を基準として当該混合信号ベクトルX（f,τ）の各要素Xq（f,τ）の偏角を、以下の演算によって正規化する。【０１９３】
【数５１】
0004406428.tif 000054
次に、周波数正規化部２１２１の第２正規化部２１２１ｂが、第１正規化部２１２１ａで正規化された各要素Xq'''（f,τ）の偏角を、以下のように周波数ｆに比例した値4fc-1dで除算する。【０１９４】
【数５２】
0004406428.tif 000055
次に、制御部２１４０において、一時メモリ２１４１に格納されたパラメータｑがｑ＝Ｍを満たすか否かを判断する（ステップＳ２１３）。ここでｑ＝Ｍでなければ、制御部２１４０は、ｑ＋１の演算結果を新たなパラメータｑの値とし、これを一時メモリ２１４１に格納し（ステップＳ２１４）、処理をステップＳ２１２へ戻す。一方、ｑ＝Ｍであれば、制御部２１４０はステップＳ２０２の処理を終了させ、以下のステップＳ２０３の処理を実行させる。【０１９５】
［ノルム正規化部２１２２における処理の詳細（ステップＳ２０３の詳細）］
図３５Ｂは、図３４に示したステップＳ２０３の処理の詳細を説明するためのフローチャートである。以下、この図に沿ってステップＳ２０３の処理の詳細を説明する。
ノルム正規化部２１２２（図３２，図３３）において、メモリ２１００の記憶領域２１０３から、各（f,τ）に対応する周波数正規化ベクトルX'（f,τ）=[X1'（f,τ）,...,XM'（f,τ）]Tを読み込む。そして、ノルム正規化部２１２２は、
【０１９６】
【数５３】
0004406428.tif 000056
の演算を行ってそれらのノルム‖X'（f,τ）‖を求め、周波数正規化ベクトルX'（f,τ）及びノルム‖X'（f,τ）‖を一時メモリ２１４１に格納する（ステップＳ２２１）。
次に、ノルム正規化部２１２２は、一時メモリ２１４１から各（f,τ）に対応する周波数正規化ベクトルX'（f,τ）及びノルム‖X'（f,τ）‖を読み出し、
X''（f,τ）=X'（f,τ）／‖X'（f,τ）‖    …（63）
の演算を行ってノルム正規化ベクトルX''（f,τ）を求める（ステップＳ２２２）。求められた各ノルム正規化ベクトルX''（f,τ）は、メモリ２１００の記憶領域２１０４に格納される。これによりステップＳ２０３の処理が終了する。【０１９７】
このように生成されたノルム正規化ベクトルX''（f,τ）は、周波数に依存せず、信号源の位置のみに依存するベクトルとなる。その結果、このノルム正規化ベクトルX''（f,τ）はクラスタを形成することになる。以下にこの理由を説明する。
［ノルム正規化ベクトルX''（f,τ）がクラスタを形成する理由］
本形態では源信号のスパース性を仮定しているので、混合信号ベクトルX（f,τ）の各要素Xｑ（f,τ）は、源信号ｐに対応する信号源ｋからセンサｑヘの周波数応答Ｈｑｋに比例した（ある複素数スカラである源信号Sk（f,τ）が掛かった）ものになっている（Xq（f,τ）=Hqk（f,τ）・Sk（f,τ））。【０１９８】
これらの源信号Sk（f,τ）は離散時間に応じて（すなわち位相に応じて）変化するが、当然ながら、周波数ｆが同じであれば、センサｑで観測される源信号Sk（f,τ）の偏角と、基準センサＱで観測される源信号Sk（f,τ）の偏角との相対値は離散時間が変化しても一定である。
前述のように、周波数正規化部２１２１の第１正規化部２１２１ａは、混合信号ベクトルX（f,τ）の特定の１つの要素XQ（f,τ）を基準として当該混合信号ベクトルX（f,τ）の各Xq（f,τ）の偏角を正規化する。【０１９９】
これにより、上述の源信号Sk（f,τ）の位相に起因する不確定性を取り除き、源信号ｐとセンサｑとに対応する混合信号ベクトルX（f,τ）の各要素Xq（f,τ）の偏角を、源信号ｐと基準センサＱ（基準値Ｑに対応）とに対応する混合信号ベクトルX（f,τ）の要素XQ（f,τ）の偏角に対する相対値として表現する。なお、この場合、要素XQ（f,τ）の偏角に対応する相対値は０と表現される。
信号源ｋからセンサqヘの周波数応答を、反射や残響の無い直接波モデルで近似して考える。すると、上記の第１正規化部２１２１ａにより正規化された偏角は、信号源ｋからセンサへの波の到達時間差と周波数ｆの双方に比例したものになる。ここでの到達時間差とは、信号源ｋからの波がセンサｑに到達するまで時間と、当該波が基準センサＱに到達するまでの時間との時間差である。【０２００】
また、前述のように、第２正規化部２１２１ｂは、第１正規化部２１２１ａで正規化された各要素Xq'''（f,τ）の偏角を周波数ｆに比例した値で除算する。これにより、各要素Xq''' （f,τ）を、それらの偏角の周波数依存性を排除した各要素Xq'（f,τ）に正規化する。これにより、正規化された各要素Xq'（f,τ）は、直接波モデルに従えば、信号源ｋからセンサへの波の到達時間差のみに依存するものとなる。ここで、信号源ｋからセンサへの波の到達時間差は、信号源ｋ，センサｑ，基準センサＱの相対位置にのみ依存する。そのため、信号源ｋ，センサｑ，基準センサＱが同じであれば、周波数ｆが異なっても各要素Xq'（f,τ）の偏角は同一となる。従って、周波数正規化ベクトルX'（f,τ）は、周波数ｆには依存せず、信号源ｋの位置のみに依存する。そのため、周波数正規化ベクトルX'（f,τ）のノルムを正規化したノルム正規化ベクトルX''（f,τ）のクラスタリングによって、同じ信号源毎に対応するクラスタが形成される。なお、実際の環境では、反射や残響などの影響により、直接波モデルは厳密には満たされないが、後述の実験結果に示すとおり十分に良い近似となっている。【０２０１】
次に、ノルム正規化ベクトルX''（f,τ）がクラスタを形成する理由を、モデルを用いて説明する。
前述した式（１）に示されるインパルス応答hqk（r）を直接波（近距離場）混合モデルを用いて近似し、周波数領域で表すと、【０２０２】
【数５４】
0004406428.tif 000057
となる。ここで、ｄｑｋは信号源ｋとセンサｑとの間の距離であり、γ（f）は周波数に依存した定数である。また、減衰γ（f）/dqｋは、距離ｄｑｋと定数γ（f）によって決まり、遅延（dqk-dQk）/cは、基準センサＱの位置で正規化された距離によって決まる。
また、信号のスパース性を仮定すると、各時間周波数（f,τ）において以下の関係が成り立つ。【０２０３】
Xq（f,τ）=Hqk（f,τ）・Sk（f,τ） …（65）
ここで、式（６２）（６３）（６４）（６５）より、【０２０４】
【数５５】
0004406428.tif 000058
となる。
この式から分かるように、ノルム正規化ベクトルX''（f,τ）の各要素Xq'' （f,τ）は、周波数ｆからは独立であり、信号源ｋとセンサｑの位置のみに依存する。従って、ノルム正規化ベクトルをクラスタリングすると、同じ信号源ごとに対応するクラスタが形成される。【０２０５】
また、同様なことは、信号の減衰を考慮していない近距離場混合モデルや遠距離場混合モデルでモデル化した場合にも言える（１の実施の形態と同様）。
また、第１の実施の形態と同様、式（６６）から、パラメータｄの値は、ｄ＞ｄmax／２であることが望ましく（ｄmaxは要素XQ'' （f,τ）に対応する基準センサと他のセンサとの最大距離を意味する。）、より好ましくはｄ≧ｄmaxであることが望ましく、さらにより好ましくはｄ＝ｄmaxであることが望ましいことが分かる。
図３７及び図３８は、パラメータｄ毎のノルム正規化ベクトルX''（f,τ）の要素Xq'' （f,τ）と、その偏角arg[Xq'' （f,τ）]との関係を説明するための複素平面図である。なお、これらにおける横軸は実軸を縦軸は虚軸を示している。【０２０６】
図３７Ａはｄmax／２≧dの場合における複素平面図である。ここで上述のdmaxの定義より、任意のｑ及びｋに対しdqk-dQkの絶対値はdmax以下となる。よってｄmax／２≧dの場合、（π/2）・（dqk-dQk）/d≦-π，（π/2）・（dqk-dQk）/d≧πと成り得る。その結果、式（６６）で表されるXq'' （f,τ）の偏角arg[Xq'' （f,τ）]は、2πを超えるα1≦arg[Xq'' （f,τ）]≦α2（α1≦-π、α2≧π）の範囲に分布する可能性がある。そのため、異なるノルム正規化ベクトルX''（f,τ）の要素Xq'' （f,τ）の偏角が一致する可能性があり、前述したクラスタリングにおいて、異なるノルム正規化ベクトルX''（f,τ）を同じクラスタにクラスタリングしてしまう可能性がある。よって、ｄ＞ｄmax／２であることが望ましい。しかし、この偏角の重複範囲に対応するノルム正規化ベクトルX''（f,τ）のサンプルが存在しないならばｄmax／２≧ｄとしても問題はない。【０２０７】
図３７Ｂはｄmax／２＜ｄ＜ｄmaxの場合における複素平面図である。この場合、-π＜（π/2）・（dqk-dQk）/d＜-π/2，π/2＜（π/2）・（dqk-dQk）/d＜πと成り得る。その結果、式（６６）で表されるXq'' （f,τ）の偏角arg[Xq'' （f,τ）]は、β1≦arg[Xq'' （f,τ）]≦β2（-π＜β1＜-π/2、π/2＜β2＜π）の範囲に分布する可能性がある。そのため、-π＜arg[Xq'' （f,τ）]＜-π/2及びπ/2＜arg[Xq'' （f,τ）]＜πの範囲において、異なるノルム正規化ベクトルX''（f,τ）の要素間における偏角の差の増加に伴い、これらの要素間の距離が単調増加しないこともありうる。これは、前述したクラスタリングの精度を低下させる可能性がある。よってｄ≧ｄmaxであることがより望ましい。【０２０８】
図３８Ａはｄ＝ｄmaxの場合における複素平面図であり、図３８Ｂはｄ＞ｄmaxの場合における複素平面図である。ここでｄ＞ｄmaxの場合、-π/2＜（π/2）・（dqk-dQk）/d＜0，0＜（π/2）・（dqk-dQk）/d＜π/2と成り得る。その結果、式（６６）で表されるXq'' （f,τ）の偏角arg[Xq'' （f,τ）]は、図３８Ｂに示すように、γ1≦arg[Xq'' （f,τ）]≦γ2（-π/2＜γ1＜0、0＜γ2＜π/2）の範囲に分布する。そして、ｄが大きくなればなるほどその分布範囲は狭くなっていき、狭い範囲にクラスタが密集していく。これは、前述したクラスタリングの精度を低下させる。【０２０９】
これに対し、ｄ＝ｄmaxである場合、-π/2≦（π/2）・（dqk-dQk）/d＜0，0＜（π/2）・（dqk-dQk）/d≦π/2と成り得る。その結果、式（６６）で表されるXq'' （f,τ）の偏角arg[Xq'' （f,τ）]は図３８Ａに示すように-π/2≦arg[Xq'' （f,τ）]≦π/2の範囲に分布する。この場合、ノルム正規化ベクトルX''（f,τ）の要素間における偏角の差の増加に対して、それらの距離も単調増加するという関係を維持しつつ、できるだけ広い範囲にクラスタを分散させることができる。その結果、一般的にクラスタリングの精度を向上させることができる（［周波数正規化部２１２１・ノルム正規化部２１２２の処理の詳細］の説明終わり）。【０２１０】
［クラスタリング部２１２３の処理の詳細（ステップＳ２０４の詳細）］
前述のようにクラスタリング部２１２３は、メモリ２１００の記憶領域２１０４からノルム正規化ベクトルX''（f,τ）を読み込み、これらをクラスタリングしてＭ個のクラスタを生成する。このクラスタリングは、例えば、各クラスタのメンバ（X''（f,τ）∈Ck）と各クラスタのセントロイドηｋとの間の二乗和Ｕｋの総和Ｕ【０２１１】
【数５６】
0004406428.tif 000059
を最小化することを基準に行われる。この最小化は、例えば、非特許文献６などで解説されているk-meansクラスタリングを用いることによって効果的に行うことができる。なお、クラスタ情報Ｃｋが示すクラスタのセントロイド（中心ベクトル）ηｋは、【０２１２】
【数５７】
0004406428.tif 000060
によって計算される。ここで｜Ck｜はクラスタ情報Ckが示すクラスタのメンバ（ノルム正規化ベクトルX''（f,τ））の数である。また、ここでは、距離としてユークリッド距離の自乗を用いているが、これを一般化したミンコフスキー距離などを用いてもよい（［クラスタリング部２１２３の処理の詳細］の説明終わり）。
［分離信号生成部２１２４の処理の詳細（ステップＳ２０５の詳細）］
図３６は、図３４に示したステップＳ２０５の処理の詳細を説明するためのフローチャートである。以下、この図に沿ってステップＳ２０５の処理の詳細を説明する。【０２１３】
まず、制御部２１４０（図３２）において、全てのパラメータｋ（ｋ＝１，...，Ｎ）及び時間周波数（ｆ，τ）（定義された範囲における全てのｆとτ）に対するＹｋ（ｆ，τ）の値を０に初期化し、これらをメモリ２１００の記憶領域２１０７に格納する（ステップＳ２３０）。
次に、制御部２１４０においてパラメータｋに１を代入し、これを一時メモリ２１４１に格納する（ステップＳ２３１）。次に、分離信号生成部２１２４（図３２，図３３）において、メモリ２１００の記憶領域２１０６のクラスタ情報Ｃｋを読み込み、これが示すｋ番目のクラスタのメンバ（ノルム正規化ベクトルX''（f,τ））を抽出し、これらを一時メモリ２１４１に格納する（ステップＳ２３２）。次に、分離信号生成部２１２４は、ステップＳ２３２で一時メモリ２１４１に格納された各ノルム正規化ベクトルX''（f,τ）を参照し、これらに対応する時間周波数（f,τ）の各混合信号ベクトルX（f,τ）をメモリ２１００の記憶領域２１０２から読み込み、これらを一時メモリ２１４１に格納する（ステップＳ２３３）。次に、分離信号生成部２１２４は、メモリ２１００の記憶領域２１０５から基準値Ｑ’を読み込み、一時メモリ２１４１から、ステップＳ２３３で格納された混合信号ベクトルX（f,τ）のＱ’番目の要素XQ’（f,τ）を抽出する（各時間周波数（f,τ）について）。そして、分離信号生成部２１２４は、抽出した要素XQ’（f,τ）を分離信号ベクトルY（f,τ）のｋ番目の要素Yk（f,τ）としてメモリ２１００の記憶領域２１０７の値を更新する（ステップＳ２３４）。すなわち、この例の分離信号生成部２１２４は、【０２１４】
【数５８】
0004406428.tif 000061
として要素Yk（f,τ）を抽出する。
次に、制御部２１４０において、一時メモリ２１４１に格納されたパラメータｋがｋ＝Ｎを満たすか否かを判断する（ステップＳ２３５）。ここでｋ＝Ｎでなければ、制御部２１４０は、ｋ＋１の演算結果を新たなパラメータｋの値とし、これを一時メモリ２１４１に格納し（ステップＳ２３６）、処理をステップＳ２３２へ戻す。一方、ｋ＝Ｎであれば、制御部２１４０はステップＳ２０５の処理を終了させる（［分離信号生成部２１２４の処理の詳細］の説明終わり）。【０２１５】
＜実験結果＞
次に、本形態による音源分離実験結果を示す。ここでは、本形態の効果を示すために二種類の信号分離の実験を行った。
１つ目の実験は、2つのセンサによる分離実験である。その実験条件は図３９Ａに示す通りである。信号源の数は３つであり、６秒間の英語の音声をスピーカから流した。また、この実験結果を図３９Ｂの表に示す。この表では、SIR（signal-to-interference ratio）の向上量を示している。大きい数字ほど分離性能が良いことを示す。ここでは、分離前のセンサでの観測結果［InputSIR］と、DOAによるクラスタリングによる結果［DOA（Previous）］と、本形態（正規化によるクラスタリング）による結果［Normalized obser. vector（Proposed）］を示している。この結果より、２つのセンサの場合、本形態の方法によって、DOAのみによるクラスタリング結果と同等の性能で信号分離が可能であることが確認できた。【０２１６】
２つ目の実験は、不規則なセンサ配置によるものである。実験条件は図４０Ａに示す通りである。この実験では、４つの無指向性マイク（センサ）を直線的ではない配置とした。また、これらの配置情報としては、マイク間隔の上限が４ｃmであるということだけを分離システムに与えた。信号源の数は４つであり、6秒問の英語の音声をスピーカから流した。このようなセンサや信号源の配置においてDOAを使おうとすると、それぞれのセンサペア毎にDOAを推定→センサペア毎にクラスタリング→全てのセンサペアについでのクラスタリング結果を統合、という繁雑な処理を行う必要がある。しかし、本形態の方法ではそのような繁雑な統合処理を行うことなく、図４０Ｂの表に示すような高い分離性能を得ることができる。また、2つ目の実験を図４１Ａのような条件で行った場合も同様に図４１Ｂの表に示すような高い分離性能を得ることができた。【０２１７】
＜第６の実施の形態の特徴＞
以上より、上述した実施の形態の特徴をまとめると以下のようになる。
（1）混合信号ベクトルから得られる情報をすべて用いてクラスタリングを行うため、全てのセンサの情報を有効に活用でき、信号分離の性能が向上する。
（2）センサの配置情報を厳密に知る必要がないため、不規則なセンサ配置を採用でき、さらにセンサ位置をキャリブレーションする必要が無い。
＜変形例等＞
なお、本発明は、上述の第６の実施の形態に限定されるものではない。例えば、第６の実施の形態では、周波数正規化部２１２１の第１正規化部２１２１ａが、混合信号ベクトルX（f,τ）の特定の１つの要素Xq（f,τ）を基準として当該混合信号ベクトルX（f,τ）の各要素Xq（f,τ）の偏角を、式（６１）の演算によって正規化することとした。しかし、例えば、以下のような式に従って、周波数正規化部２１２１の第１正規化部２１２１ａが、混合信号ベクトルX（f,τ）の特定の１つの要素Xq（f,τ）を基準として当該混合信号ベクトルX（f,τ）の各要素Xq（f,τ）の偏角を正規化する構成であってもよい。【０２１８】
【数５９】
0004406428.tif 000062
ただし、・＊は、・の複素共役である。また、ψ｛・｝は関数であり、クラスタリング精度の観点から好ましくは単調増加関数であることが望ましい。
また、周波数正規化部２１２１が式（６０）の代わりに、【０２１９】
【数６０】
0004406428.tif 000063
等の演算によって周波数正規化を行うこととしてもよい。ただし、ρは定数である（例えば、ρ＝１）。
また、上述の第６の実施の形態では、ノルム正規化部２１２２においてノルムが１になるような正規化を行ったが、ノルムが１以外の規定値となるように正規化を行ってもよい。さらに、ノルム正規化部２１２２を設けず、ノルム正規化を行わない構成であってもよい。この場合、クラスタリング部２１２３は、周波数正規化ベクトルX'（f,τ）のクラスタリングを行うことになる。しかし、周波数正規化ベクトルX'（f,τ）は、ノルムが統一されていない。そこで、この場合のクラスタリング基準は、ベクトルがノルムを含めて似ているかどうかではなく、ベクトルの方向のみが似ているかどうかになる。これは、類似度を用いた評価になる。類似度の１つとしてコサイン距離
cosθ=｜X'H（f,τ）・ηk｜/（‖X'（f,τ）‖・‖ηk‖）
を例示できる。ここでθは、周波数正規化ベクトルX'（f,τ）と、セントロイドηkのベクトルとがなす角度である。コサイン距離を用いる場合、クラスタリング部２１２３は、コサイン距離の総和【０２２０】
【数６１】
0004406428.tif 000064
を最小値化するクラスタを生成する。なお、セントロイドηkは、各クラスタのメンバの平均として算出する。
また、上述した基準値Ｑ，Ｑ’は互いに等しくてもよく、異なっていてもよい。
また、パラメータｄは、全てのセンサｑに対して同一であってもよいし、各センサｑに対応して複数の値を設定しておいてもよい。例えば、基準センサとセンサｑとの距離を、それぞれのセンサｑに対応するパラメータｄの値としてもよい。【０２２１】
また、分離信号生成部２１２４において、【０２２２】
【数６２】
0004406428.tif 000065
とする代わりに、バイナリマスク【０２２３】
【数６３】
0004406428.tif 000066
を作成し、
Yk（f,τ）=Mk（f,τ）XQ'（f,τ）
として分離信号ベクトルY（f,τ）のｋ番目の要素Yk（f,τ）を得ることとしてもよい。
さらに、上述の各実施の形態では、フーリエ変換及び逆フーリエ変換によって周波数領域と時間領域との間の変換を行うこととしたが、wavelet変換、DFTフィルタバンク、ポリフェイズフィルタバンクなどを用い、この変換を行うこととしてもよい（例えば、「R. E. Crochiere, L. R. Rabiner, "Multirate Digital Signal Processing." Eaglewood Cliffs, NJ: Prentice-Hall,1983 （ISBN 0-13-605162-6） ）。また、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。
【０２２４】
また、上述の各実施の形態の構成をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。
これらの処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよいが、具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto-Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。【０２２５】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
また、上述した実施形態とは別の実行形態として、コンピュータが可搬型記録媒体から直接このプログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。【０２２６】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【産業上の利用可能性】
【０２２７】
本技術により、様々な妨害信号が発生する実環境において、目的の信号を精度良く取り出すことが可能となる。音信号に対する応用例としては、例えば、音声認識器のフロントエンドとして働く音声分離システムなどが挙げられる。話者とマイクが離れた位置にあり、マイクが話者の音声以外を集音してしまうような状況でも、そのようなシステムを使うことで、話者の音声のみを取り出して正しく音声を認識することができる。
【図１】
0004406428.tif 000070
【図２】
0004406428.tif 000071
【図３】
0004406428.tif 000072
【図４】
0004406428.tif 000073
【図５】
0004406428.tif 000074
【図６】
0004406428.tif 000075
【図７】
0004406428.tif 000076
【図８】
0004406428.tif 000077
【図９】
0004406428.tif 000078
【図１０】
0004406428.tif 000079
【図１１】
0004406428.tif 000080
【図１２】
0004406428.tif 000081
【図１３】
0004406428.tif 000082
【図１４】
0004406428.tif 000083
【図１５】
0004406428.tif 000084
【図１６】
0004406428.tif 000085
【図１７】
0004406428.tif 000086
【図１８】
0004406428.tif 000087
【図１９】
0004406428.tif 000088
【図２０】
0004406428.tif 000089
【図２１】
0004406428.tif 000090
【図２２】
0004406428.tif 000091
【図２３】
0004406428.tif 000092
【図２４】
0004406428.tif 000093
【図２５】
0004406428.tif 000094
【図２６】
0004406428.tif 000095
【図２７】
0004406428.tif 000096
【図２８】
0004406428.tif 000097
【図２９】
0004406428.tif 000098
【図３０】
0004406428.tif 000099
【図３１】
0004406428.tif 000100
【図３２】
0004406428.tif 000101
【図３３】
0004406428.tif 000102
【図３４】
0004406428.tif 000103
【図３５】
0004406428.tif 000104
【図３６】
0004406428.tif 000105
【図３７】
0004406428.tif 000106
【図３８】
0004406428.tif 000107
【図３９】
0004406428.tif 000108
【図４０】
0004406428.tif 000109
【図４１】
0004406428.tif 000110