(57)【特許請求の範囲】
【請求項１】
源信号の混合信号がＭ（Ｍ≧２）箇所のセンサでそれぞれ観測されて得られた観測信号xm(t)(m=1,...,M、tは時間）を、周波数領域の観測信号Xm(f,τ)（fは周波数、τは時間インデックス）に変換する周波数領域変換部と、
周波数領域の観測信号Xm(f,τ)を要素とする観測信号ベクトルX(f,τ)=[X1(f,τ),...,XM(f,τ)]Tを周波数f毎に独立にクラスタリングした場合に観測信号ベクトルX(f,τ)が属するクラスがCn(f) (n=1,...,N、N≧１)となる事象の事後確率P(Cn(f)|X(f,τ))を算出するクラス分類部と、
対応する周波数fが異なる上記事後確率P(Cn(f)|X(f,τ))間の類似度を指標として、上記事後確率P(Cn(f)|X(f,τ))と上記クラスCn(f)の番号nとの対応関係を並び替え、対応するクラスの番号が同一であって周波数が異なる事後確率間の類似度の総和が当該並び替え前よりも大きな事後確率P’(Ck(f)|X(f,τ)) (k=1,...,N)を生成するパーミュテーション問題解決部と、
上記パーミュテーション問題解決部で生成された上記事後確率P’(Ck(f)|X(f,τ))の大きさを指標とし、クラスCk(f)に属すると判定される上記周波数領域の観測信号Xm(f,τ)を周波数領域の分離信号Yn(f,τ)として抽出する分離部と、
を有することを特徴とする信号分離装置。
【請求項２】
請求項１に記載の信号分離装置であって、
上記クラス分類部は、
クラスCn(f)に属する観測信号ベクトルX(f,τ)の代表ベクトルをセントロイドan(f)とし、クラスCn(f)に属する観測信号ベクトルX(f,τ)の標準偏差をσn(f)とし、観測信号ベクトルX(f,τ)がクラスCn(f)に属する事象の確率密度関数をp(X(f,τ)|an(f),σn(f))とし、p(X(f,τ)|an(f),σn(f))の混合比をαn(f)とし、パラメータ集合θ(f)={a1(f),σ1(f),α1(f),...,aN(f),σN(f),αN(f)}とし、p(X(f,τ)|θ(f))=Σn=1Nαn(f)・p(X(f,τ)|an(f),σn(f))とした場合における、事後確率P(Cn(f)|X(f,τ),θ(f))=αn(f)・p(X(f,τ)|an(f),σn(f))／p(X(f,τ)|θ(f))を、上記パラメータ集合θ(f)を固定値として各観測信号ベクトルX(f,τ)について算出する事後確率計算部と、
各観測信号ベクトルX(f,τ)に対応するΣn=1N {P(Cn(f)|X(f,τ),θ(f))・logαn(f)・p(X(f,τ)|an(f),σn(f))}を周波数f毎に独立に加算したQ(f,θ(f))がそれぞれ最大となるパラメータ集合θ(f)を、事後確率P(Cn(f)|X(f,τ),θ(f))を固定値として算出するパラメータ推定部と、
所定の終了条件を満たすまで上記事後確率計算部の処理と上記パラメータ推定部の処理とを交互に実行させる第１演算制御部と、を有し、
上記事後確率P(Cn(f)|X(f,τ))は、
上記終了条件を満たした際に事後確率計算部で算出されていた最新の事後確率P(Cn(f)|X(f,τ),θ(f))である、
ことを特徴とする信号分離装置。
【請求項３】
請求項１又は２に記載の信号分離装置であって、
上記パーミュテーション問題解決部は、
対応する周波数fが異なる事後確率P(Cn(f)|X(f,τ))間の類似度を指標として、上記事後確率P(Cn(f)|X(f,τ))と上記クラスCn(f)の番号nとの対応関係を並び替える順列Πfを周波数f毎に生成するクラスタリング部と、
上記クラスタリング部で生成された順列Πfに従い、周波数f毎に上記事後確率P(Cn(f)|X(f,τ))と上記クラスCn(f)の番号nとの対応関係を並び替え、上記事後確率P’(Ck(f)|X(f,τ))を生成する並び替え部と、
を有することを特徴とする信号分離装置。
【請求項４】
請求項３に記載の信号分離装置であって、
上記クラスタリング部は、
取り扱う全周波数の集合をFとした場合における、すべての異なる周波数f∈Fの組合せに対応する事後確率P(Cn(f)|X(f,τ))間の類似度を指標とし、上記事後確率P(Cn(f)|X(f,τ))と上記クラスCn(f)の番号nとの対応関係を並び替える順列Πfを周波数f毎に生成する大域的最適化部を有する、
ことを特徴とする信号分離装置。
【請求項５】
請求項４に記載の信号分離装置であって、
上記大域的最適化部は、
上記順列Πfに従って周波数f毎に上記事後確率P(Cn(f)|X(f,τ))と上記クラスCn(f)の番号nとの対応関係を並び替えた事後確率P’(Ck(f)|X(f,τ))の周波数方向の平均値又は代表値をセントロイドck(τ)として算出するセントロイド算出部と、
上記事後確率P(Cn(f)|X(f,τ))とセントロイドcn(τ)との類似度を指標とし、上記順列Πfを周波数f毎に生成する順列生成部と、
所定の終了条件を満たすまでセントロイド算出部の処理と上記順列生成部の処理とを交互に実行させる第２演算制御部と、
を有することを特徴とする信号分離装置。
【請求項６】
請求項４又は５に記載の信号分離装置であって、
上記クラスタリング部は、
対応する周波数f∈Fが異なる事後確率P(Cn(f)|X(f,τ))間の類似度のうち、特定の周波数の組合せに対応する事後確率間の類似度のみを指標として用い、上記大域的最適化部で生成された順列Πfを更新し、新たな順列Πfを生成する局所的最適化部を更に有する、
ことを特徴とする信号分離装置。
【請求項７】
請求項６に記載の信号分離装置であって、
上記特定の周波数の組合せは、
周波数差が所定範囲内にある周波数の組合せ、及び／又は、倍音関係にある周波数の組合せである、
ことを特徴とする信号分離装置。
【請求項８】
周波数領域変換部が、源信号の混合信号がＭ（Ｍ≧２）箇所のセンサでそれぞれ観測されて得られた観測信号xm(t)(m=1,...,M、tは時間）を、周波数領域の観測信号Xm(f,τ)（fは周波数、τは時間インデックス）に変換する周波数領域変換過程と、
クラス分類部が、周波数領域の観測信号Xm(f,τ)を要素とする観測信号ベクトルX(f,τ)=[X1(f,τ),...,XM(f,τ)]Tを周波数f毎に独立にクラスタリングした場合に観測信号ベクトルX(f,τ)が属するクラスがCn(f) (n=1,...,N、Ｎ≧１)となる事象の事後確率P(Cn(f)|X(f,τ))を算出するクラス分類過程と、
パーミュテーション問題解決部が、対応する周波数fが異なる上記事後確率P(Cn(f)|X(f,τ))間の類似度を指標として、上記事後確率P(Cn(f)|X(f,τ))と上記クラスCn(f)の番号nとの対応関係を並び替え、対応するクラスの番号が同一であって周波数が異なる事後確率間の類似度の総和が当該並び替え前よりも大きな事後確率P’(Ck(f)|X(f,τ)) (k=1,...,N)を生成するパーミュテーション問題解決過程と、
分離部が、上記周波数領域の観測信号Xm(f,τ)と上記パーミュテーション問題解決部で生成された上記事後確率P’(Ck(f)|X(f,τ))とを用い、周波数領域の分離信号Yn(f,τ)を抽出する分離過程と、
を有することを特徴とする信号分離方法。
【請求項９】
請求項１から７の何れかに記載の信号分離装置としてコンピュータを機能させるためのプログラム。
【請求項１０】
請求項９に記載のプログラムを格納したコンピュータ読取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、信号処理の技術分野に属し、特に、複数の信号が空間内で混合されたものから、源信号をできるだけ正確に復元する信号分離技術に関する。
【背景技術】
【０００２】
［ブラインド信号分離］
まず、ブラインド信号分離の定式化を行う。扱う信号はあるサンプリング周波数fsでサンプリングされ、離散的に表現されるものとする。また、N個の信号が混合されてM個のセンサで観測されたとする。以下では、信号の発生源からセンサまでの距離により信号が減衰・遅延し、また壁や床などによる反射／残響が発生する状況を扱う。このような状況での混合は、源信号sn(t)(n=1,...,N)を発した信号源nからセンサm(m=1,...,M)へのインパルス応答hmn(r)による畳み込み混合
【０００３】
【数１】
0004769238.tif 000002
となる。ここでtはサンプリング時間を、rは掃引（時間シフトした信号のサンプル値それぞれに異なる係数を作用させる操作）のための変数を、それぞれ示している。一般的なインパルス応答hmn(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)を求めることにある。
【０００４】
［周波数領域における信号分離］
畳み込み混合の問題は扱いが繁雑である。よって、上述の式（１）に短時間離散フーリエ変換(DFT: Discrete Fourier Transform)を施して、信号を周波数領域に変換した上で分離の操作を行うことが有効である。上述の式（１）に短時間離散フーリエ変換を適用して周波数毎の時間系列を求めると以下のようになる。
【０００５】
【数２】
0004769238.tif 000003
ここでfは周波数であり、f=0, fs/L ,・・・, fs(L-1)/Lと離散化されている（fsはサンプリング周波数）。また、τは時間インデックスであり、jは虚数単位である。また、g(r)は窓関数である。ハニング窓g(r)=(1+cos(2π・r /L))/2などのg(0)にパワーの中心を持つ窓関数を用いることで、Xm(f,τ)は時間τを中心とする観測信号xm(t)の周波数特性を表現する。なお、Xm(f,τ)はＬサンプルにわたる情報を含んでいるため、すべての時間tを時間インデックスτとしてXm(f,τ)を求める必要はなく、適当な間隔で時間インデックスτを設定してXm(f,τ)を求める。
【０００６】
式(1)で示される時間領域での畳み込み混合を周波数領域での表現に変換すると、
【０００７】
【数３】
0004769238.tif 000004
と各周波数での単純混合に近似でき、分離の操作が単純になる。なお、Hmn(f)は信号源nからセンサmまでの周波数応答であり、Sn(f,τ)は式(2)と同様な式に従って源信号sn(t)に短時間離散フーリエ変換を施したものである。式(3)をベクトル表記すると、
【０００８】
【数４】
0004769238.tif 000005
となる。ここで、X(f,τ)=[X1(f,τ),...,XM(f,τ)]TはXm(f,τ)を要素とする観測信号ベクトルであり、Hn(f)=[H1n(f),...,HMn(f)]Tは信号源nからセンサmまでの周波数応答Hmn(f)を要素とするベクトルである。なお、[・]Tは[・]の転置を示す。
【０００９】
［スパース性に基づく信号分離］
ブラインド信号分離法の一つにスパース性に基づいて信号分離を行う方法がある（例えば、特許文献１等参照）。このスパース性に基づく信号分離の場合、信号源の数Nとセンサの数Mの関係にかかわらず（M≧2であればN>MでもN≦Mでも良い）、同一の仕組みにより分離の処理が可能である。これは独立成分分析（ICA: Independent Component Analysis）を用いたブラインド信号分離（例えば、特許文献2、非特許文献5等参照）と対比される。ICAを用いる場合には、信号源の数Nがセンサの数Mを超えない（N≦M）ことが強く望まれる。これに対して、スパース性に基づく信号分離では、センサの数に関する要求条件がより緩くなっており、より広い適用範囲が見込まれる。
【００１０】
ただし、スパース性による信号分離が有効に働くためには、対象となる源信号がスパース性を持つことが条件となる。スパース性とは、ほとんどの場合において信号の振幅が零に近く、大きな振幅となるのは稀であるという性質である。例えば、周波数領域での音声信号にはスパース性を十分に確認できる。スパース性を満たす源信号s1(t),...,sN(t)の場合、式(4)の混合過程は、さらに
X(f,τ)=Hp(f)・Sp(f,τ)     …(5)
と近似表現できる。ここで添字ｐは、時間周波数(f,τ)に依存したものとなる。ほとんどの場合において信号源の振幅が零に近いため、個々の時間周波数(f,τ)において最も振幅の大きい源信号Sp(f,τ)に関わる項だけで、式(4)が近似されている。
【００１１】
スパース性による信号分離では、各時間周波数(f,τ)において、どの源信号Sp(f,τ)の振幅が最も大きいかを推定する。言い替えると、観測信号ベクトルX(f,τ)をN個のクラスC1,...,CNに分類（クラスタリング）し、クラスCkには、源信号Sk(f,τ)が最も支配的な観測信号ベクトルX(f,τ)が属するようにする。ここで、サンプルXを観測した後における、サンプルXがクラスCkに属する事象の事後確率をP(Ck |X)で表記する。そのような事後確率をなんらかの方法で推定できれば、分離信号Yn(f,τ)は、例えば時間周波数マスキング
【００１２】
【数５】
0004769238.tif 000006
により構成できる。ここで、Jは分離信号を構成するために用いる基準センサの添字であり、1からMの範囲から選択される。
【００１３】
［クラス分類／事後確率の計算方法］
クラス分類或いは事後確率計算の方法として、信号源の方向や位置に相当する値を推定し、それに基づいてすべての時間周波数(f,τ)に関する観測信号ベクトルX(f,τ)を一気にクラス分類したり事後確率計算したりする方法が提案されている（例えば、特許文献１、特許文献3、非特許文献1、非特許文献2、非特許文献3等参照）。これらの方法では、1)信号源毎に推定された方向や位置に相当する値に従って、すべての時間周波数(f,τ)に関する観測信号ベクトルX(f,τ)をN個のクラスに分類したり事後確率を計算したりするプロセスと、2)分類された観測信号を元に、方向や位置に相当する値を信号源毎に再推定するプロセスとを行う。これらは、反射や残響の影響が比較的少ない場合には有効に働く。
【００１４】
また、観測信号ベクトルX(f,τ)をN個のクラスヘ分類する操作を、周波数毎に行う方法も提案されている。この場合は、ある周波数でのi番目のクラスと別の周波数でのi番目のクラスが、同じ信号源に対応するものかどうかが不明となる。従って、その後、同一信号源に対応するクラスを全周波数に渡って同定する必要がある。この問題は、ICAを用いたブラインド信号分離におけるパーミュテーション問題とほぼ同じである。これに対し、各周波数での分類結果から各信号源の方向や位置に相当する値を推定し、その推定結果に基づいてパーミュテーション問題を解決する方法（例えば、特許文献２、特許文献３、非特許文献４、非特許文献５等参照）や、周波数毎の分離信号エンベロープの相関係数の類似度に基づいてパーミュテーション問題を解決する方法 (例えば、非特許文献５、非特許文献６、非特許文献７、非特許文献８等）がこれまで用いられてきた。
【特許文献１】WO2005/024788
【特許文献２】WO2004/079388
【特許文献３】WO2006/085537
【非特許文献１】O. Yilmaz and S. Rickard, "Blind separation of speech mixtures via time-frequency masking," IEEE Trans. Signal Processing, vol. 52, no. 7, pp. 1830-1847, July 2004.
【非特許文献２】M. Mandel, D. Ellis, and T. Jebara, "An EM Algorithm for Localizing Multiple Sound Sources in Reverberant Environments," Advances in Neural Information Processing Systems, vol. 19, http://books.nips.cc/papers/files/nips19/NIPS2006_0202.pdf, 2006.
【非特許文献３】S. Araki, H. Sawada, R. Mukai and S. Makino, "Underdetermined Blind Sparse Source Separation for Arbitrarily Arranged Multiple Sensors," Signal Processing., vol. 87, no.8, 99. 1833-1847, 2007.
【非特許文献４】S. Winter, W. Kellermann, H. Sawada, and S. Makino, "MAP-Based Underdetermined Blind Source Separation of Convolutive Mixtures by Hierarchical Clustering and L1-Norm Minimization," EURASIP Journal on Advances in Signal Processing, 2007, Article ID 24717.
【非特許文献５】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.
【非特許文献６】R.K. Olsson and L.K. Hansen, "Blind Separation of More Sources than Sensors in Convolutive Mixtures," Proc. ICASSP 2006, May 2006, vol. V, pp. 657-660.
【非特許文献７】J. Anemuller, B. Kollmeier, "Amplitude Modulation Decorrelation for Convolutive Blind Source Separation," in Proc. ICA 2000, June 2000, pp. 215-220.
【非特許文献８】N. Murata, S. Ikeda, and A. Ziehe, "An Approach to Blind Source Separation Based on Temporal Structure of Speech Signals," Neurocomputing, vol. 41, pp. 1-24, Oct. 2001.
【発明の開示】
【発明が解決しようとする課題】
【００１５】
しかし、すべての時間周波数(f,τ)に関する観測信号ベクトルX(f,τ)を一気にクラス分類したり事後確率計算したりする方法の場合、反射や残響の影響が強い環境では信号源の方向や位置が正確に推定できず信号分離性能が劣化してしまう。
【００１６】
また、観測信号ベクトルX(f,τ)をN個のクラスヘ分類する操作を周波数毎に行い、各周波数での分類結果から各信号源の方向や位置に相当する値を推定し、その推定結果に基づいてパーミュテーション問題を解決する場合にも、反射や残響の影響が強い環境では信号源の方向や位置が正確に推定できず信号分離性能が劣化してしまう。
【００１７】
一方、観測信号ベクトルX(f,τ)をN個のクラスヘ分類する操作を周波数毎に行い、周波数毎の分離信号エンベロープの相関係数の類似度に基づいてパーミュテーション問題を解決する場合には、反射や残響の影響をそれほど受けることなく信号分離を行うことができる。しかし、この従来方法の場合、源信号が全周波数に渡って同じような振幅のエンベロープを持たない限り、周波数全体に渡って一貫性のあるパーミュテーション問題の解を得ることはできない。以下、このことを詳細に説明する。
【００１８】
この従来方法では、式(6)などに従って分離信号Yn(f,τ)を周波数f毎に計算した後、それらのエンベロープをvif(τ)=|Yn(f,τ)|として計算する。そして、それらの類似度を相関係数によって表現し、同一の源信号に対応するエンベロープ間の相関係数が最も大きくなると仮定してパーミュテーション問題を解決する。なお、系列長（時間インデックスτの数）がそれぞれTである２つの系列vif(τ)，vkg(τ)の相関係数は、
【００１９】
【数６】
0004769238.tif 000007
として計算される。ここで、
【００２０】
【数７】
0004769238.tif 000008
は、それぞれ、相関、平均、標準偏差である。また、相関係数は、−１から１までの値を取り、２つの系列が等しいときには１となる。また、
【００２１】
【数８】
0004769238.tif 000009
は、T個のτにそれぞれ対応するα(τ)の和を意味する。
【００２２】
図１８（ａ）に、信号源が３つ存在する場合の２つの周波数f=766Hz,g=906Hzにおける分離信号のエンベロープ系列v1f(τ),...,v3f(τ)，v1g(τ),...,v3g(τ)を例示する。なお、図１８（ａ）の横軸は時間（時間インデックスτ）を示し、縦軸はエンベロープを示す。また、エンベロープ系列v1f(τ),...,v3f(τ)，v1g(τ),...,v3g(τ)は、同じ信号源に対応するものが同じ添字となるようにパーミュテーションが揃えられている（例えば、v1f(τ)とv1g(τ)とは同じ信号源に対応する系列である）。そして、３つの信号源にそれぞれ対応する系列を濃い実線、薄い実線、破線で区別してある。
【００２３】
ここで、図１８（ａ）のエンベロープ系列間の相関係数を求めると以下のようになる。
【００２４】
【数９】
0004769238.tif 000010
【００２５】
一般に、周波数fとgが隣接や倍音の関係にあれば、同じ信号源に対応する分離信号のエンベロープ系列の相関係数の値は、異なる信号源に対応する分離信号のエンベロープ系列の相関係数の値よりも格段に大きくなる。しかし、式(8)の例では、周波数fとgが隣接や倍音の関係になっていないため、同じ信号源に対応する分離信号のエンベロープ系列の相関係数の値が、異なる信号源に対応する分離信号のエンベロープ系列の相関係数の値よりもさほど大きくなっていない。特に、相関係数ρ(v1f，v1g)とρ(v1f，v2g)との間では大小関係が逆転している。これでは周波数全体に渡って一貫性のあるパーミュテーション問題の解を得ることはできない。
【００２６】
本発明はこのような点に鑑みてなされたものであり、反響や残響の影響が強い環境であっても、また、信号源の全周波数に渡って同じような振幅のエンベロープを持たない場合であっても、高精度にブラインド信号分離を行うことができる技術を提供することを目的とする。
【課題を解決するための手段】
【００２７】
本発明では、まず、周波数領域変換部が、源信号の混合信号がＭ（Ｍ≧２）箇所のセンサでそれぞれ観測されて得られた観測信号xm(t)(m=1,...,M、tは時間）を、周波数領域の観測信号Xm(f,τ)（fは周波数、τは時間インデックス）に変換する。
【００２８】
次に、クラス分類部が、周波数領域の観測信号Xm(f,τ)を要素とする観測信号ベクトルX(f,τ)=[X1(f,τ),...,XM(f,τ)]Tを周波数f毎に独立にクラスタリングした場合に観測信号ベクトルX(f,τ)が属するクラスがCn(f)(n=1,...,N、N≧１)となる事象の事後確率P(Cn(f)|X(f,τ))を算出する。この処理は周波数f毎に独立に行われるため、反響や残響の影響が強い場合であっても処理の精度はさほど低下しない。また、当該クラス分類部で算出された事後確率P(Cn(f)|X(f,τ))に対応するクラスCn(f)の番号nは信号源に対応する。しかし、クラスCn(f)の番号nと信号源との対応関係は周波数f毎に相違する可能性が高い。
【００２９】
次に、パーミュテーション問題解決部が、対応する周波数fが異なる事後確率P(Cn(f)|X(f,τ))間の類似度を指標として、事後確率P(Cn(f)|X(f,τ))とクラスCn(f)の番号nとの対応関係を並び替え、対応するクラスの番号が同一であって周波数が異なる事後確率間の類似度の総和が当該並び替え前よりも大きな事後確率P’(Ck(f)|X(f,τ)) (k=1,...,N)を生成する。ここで、同じ信号源から発せられた信号が支配的な観測信号ベクトルX(f,τ)に対応する事後確率P(Cn(f)|X(f,τ))は、周波数fが相違する場合であっても類似度が大きい。この性質は全周波数中の多くの組合せについて妥当なものである。パーミュテーション問題解決部は、この性質を利用してパーミュテーション問題を解決する。
【００３０】
そして、分離部が、パーミュテーション問題解決部で生成された事後確率P’(Ck(f)|X(f,τ))の大きさを指標とし、クラスCk(f)に属すると判定される周波数領域の観測信号Xm(f,τ)を周波数領域の分離信号Yn(f,τ)として抽出する。
【発明の効果】
【００３１】
本発明では、反響や残響の影響が強い環境であっても、また、信号源の全周波数に渡って同じような振幅にエンベロープを持たない場合であっても、高精度にブラインド信号分離を行うことができる。
【発明を実施するための最良の形態】
【００３２】
以下、本発明を実施するための最良の形態を図面を参照して説明する。
【００３３】
〔原理〕
まず、本形態の信号分離の原理について説明する。
【００３４】
本形態では、周波数領域の観測信号Xm(f,τ)を要素とする観測信号ベクトルX(f,τ)=[X1(f,τ),...,XM(f,τ)]TをN個のクラスへ分類（クラスタリング）する操作を周波数f毎に行う。従ってパーミュテーション問題を解決する必要がある。本形態では、この分類操作において観測信号ベクトルX(f,τ)がクラスCn(f)に属する事後確率P(Cn(f)|X(f,τ))を明示的に計算しておき、この事後確率の系列を用いてパーミュテーション問題を解決する。このように周波数毎の分類操作において事後確率を明示的に計算し、その事後確率の系列を用いてパーミュテーション問題を高精度に解決する部分に本形態の特徴がある。
【００３５】
すなわち、本形態では、事後確率P(Cn(f)|X(f,τ))の系列（「アクティブ系列」と呼ぶ）
vnf(τ)= P(Cn(f)|X(f,τ))     …(9)
の類似度を求め、同一の源信号に対応するアクティブ系列間の類似度が最も大きくなると仮定してパーミュテーション問題を解決する。
【００３６】
図１８（ｂ）は、前述の図１８（ａ）と同じ観測信号に基づき生成されたアクティブ系列vnf(τ)= P(Cn(f)|X(f,τ))を示したグラフである。なお、図１８（ｂ）の横軸は時間（時間インデックスτ）を示し、縦軸は事後確率を示す。また、アクティブ系列v1f(τ),...,v3f(τ)，v1g(τ),...,v3g(τ)は、同じ信号源に対応するものが同じ添字となるようにパーミュテーションが揃えられている。そして、３つの信号源にそれぞれ対応する系列を濃い実線、薄い実線、破線で区別してある。
【００３７】
ここで、図１８（ｂ）と前述の図１８（ａ）と比較すれば分かるように、同じ信号源に対応するアクティブ系列の相関関係は、同じ信号源に対応するエンベロープ系列の相関関係よりも強いことが分かる。
【００３８】
例えば、図１８（ａ）に例示した周波数fのエンベロープ系列の時間３秒直後では、破線で示されたエンベロープ系列の振幅が格段に大きくなっており、濃い実線や薄い実線で示されたエンベロープ系列の振幅は零に近い。しかし、このようなエンベロープ系列の特徴は、図１８（ａ）の周波数ｇのエンベロープ系列の時間３秒直後には顕著に表れていない。これに対し、図１８（ｂ）に示す周波数fのアクティブ系列の時間３秒直後と、周波数ｇのアクティブ系列の時間３秒直後とは互いに類似した振幅を持つ。
【００３９】
また、図１８（ｂ）に示すアクティブ系列について式(7)に示した相関係数を求めると
【００４０】
【数１０】
0004769238.tif 000011
【００４１】
となる。このように、同じ信号源に対応するアクティブ系列の相関係数の値は、異なる信号源に対応するアクティブ系列の相関係数の値よりも格段に大きくなる。
【００４２】
なお、上記の例では二つの周波数f=766Hz及びg=906Hzの組み合わせを選択して、それらに関する相関係数を計算した結果を示した。しかし、同じ信号源に対応する系列間の相関係数の値が大きくなるという傾向は、多くの周波数の組み合わせに関し、エンベロープ系列よりもアクティブ系列のほうが顕著である。
【００４３】
以上より、本形態では、事後確率の系列であるアクティブ系列の類似度を用いることで、エンベロープ系列を用いていた従来技術よりも高精度にパーミュテーション問題を解決することができる。
【００４４】
〔第１実施形態〕
＜信号分離装置の構成＞
図１は、本形態の信号分離装置１０の機能構成の全体を例示したブロック図である。また、図２（ａ）は、図１に示したクラス分類部１２０の機能構成の詳細を例示したブロック図である。また、図２（ｂ）は、図１に示したパーミュテーション問題解決部１３０の機能構成の詳細を例示したブロック図である。また、図３は、図２（ｂ）のクラスタリング部１３２の機能構成の詳細を例示したブロック図である。また、図４は、本形態の信号分離装置１０を構成するハードウェアの構成を例示したブロック図である。なお、各図において、実線の矢印はデータの流れを示し、点線の矢印は論理的な情報の流れを示す。しかし、制御部１６０やメモリ１７０等、一部の構成に対するデータの流れの表記は省略する。
【００４５】
以下、これらの図を用い、本形態の信号分離装置の構成を説明する。
【００４６】
［ハードウェア構成］
図４に例示するように、この例の信号分離装置１０は、ＣＰＵ（Central Processing Unit）１０ａ、入力部１０ｂ、出力部１０ｃ、補助記憶装置１０ｆ、ＲＡＭ（Random Access Memory）１０ｄ、ＲＯＭ（Read Only Memory）１０ｅ及びバス１０ｇを有している。
【００４７】
この例のＣＰＵ１０ａは、制御部１０ａａ、演算部１０ａｂ及びレジスタ１０ａｃ有し、レジスタ１０ａｃに読み込まれた各種プログラムに従って様々な演算処理を実行する。また、補助記憶装置１０ｆは、本形態の信号分離処理を実行するための信号分離プログラムを格納した信号分離プログラム領域１０ｆａ及びセンサで観測された時間領域の混合信号等の各種データが格納されるデータ領域１０ｆｂを有している。また、ＲＡＭ１０ｄは、信号分離プログラムが書き込まれる信号分離プログラム領域１０ｄａ及び各種データが書き込まれるデータ領域１０ｄｂを有している。また、この例のバス１０ｇは、ＣＰＵ１０ａ、入力部１０ｂ、出力部１０ｃ、補助記憶装置１０ｆ、ＲＡＭ１０ｄ及びＲＯＭ１０ｅを通信可能に接続している。
【００４８】
［ハードウェアとソフトウェアとの協働］
この例のＣＰＵ１０ａは、読み込まれたＯＳ（Operating System）プログラムに従い、補助記憶装置１０ｆの信号分離プログラム領域１０ｆａに格納されている信号分離プログラムを、ＲＡＭ１０ｄの信号分離プログラム領域１０ｄａに書き込む。同様にＣＰＵ１０ａは、補助記憶装置１０ｆのデータ領域１０ｆｂに格納されている時間領域の混合信号等の各種データをＲＡＭ１０ｄのデータ領域１０ｄｂに書き込む。さらに、ＣＰＵ１０ａは、この信号分離プログラムや各種データが書き込まれたＲＡＭ１０ｄ上のアドレスをレジスタ１０ａｃに格納する。そして、ＣＰＵ１０ａの制御部１０ａａは、レジスタ１０ａｃに格納されたこれらのアドレスを順次読み出し、読み出したアドレスが示すＲＡＭ１０ｄ上の領域からプログラムやデータを読み出し、そのプログラムが示す演算を演算部１０ａｂに順次実行させ、その演算結果をレジスタ１０ａｃに格納していく。
【００４９】
このようにＣＰＵ１０ａに信号分離プログラムが読み込まれることにより図１から図３に例示する機能構成を具備する信号分離装置１０が構築される。
【００５０】
図１に例示するように、本形態の信号分離装置１０は、メモリ１００，１７０と、周波数領域変換部１１０と、クラス分類部１２０と、パーミュテーション問題解決部１３０と、分離部１４０と、時間領域変換部１５０と、信号分離装置１０全体を制御する制御部１６０とを有する。また、図２（ａ）に例示するように、この例のクラス分類部１２０は、ノルム正規化部１２１と、モデル化部１２２とを有し、モデル化部１２２は、初期パラメータ設定部１２２ａと、事後確率計算部１２２ｂと、パラメータ推定部１２２ｃと、演算制御部１２２ｄとを有する。また、図２（ｂ）に例示するように、パーミュテーション問題解決部１３０は、アクティブ系列生成部１３１と、クラスタリング部１３２と、並び替え部１３３とを有する。また、図３に例示するように、この例のクラスタリング部１３２は、大域的最適化部１３２ａと、局所的最適化部１３２ｂとを有する。また、大域的最適化部１３２ａは、初期パラメータ設定部１３２ａａと、順列生成部１３２ａｂと、セントロイド算出部１３２ａｃと、演算制御部１３２ａｄとを有し、局所的最適化部１３２ｂは、順列生成部１３２ｂａと、演算制御部１３２ｂｂとを有する。
【００５１】
ここでメモリ１００，１７０は、レジスタ１０ａｃ、補助記憶装置１０ｆのデータ領域１０ｆｂ或いはＲＡＭ１０ｄのデータ領域１０ｄｂ等に相当する。また、周波数領域変換部１１０、クラス分類部１２０、パーミュテーション問題解決部１３０、分離部１４０、時間領域変換部１５０及び制御部１６０は、ＣＰＵ１０ａにＯＳプログラムや信号分離プログラムが読み込まれることにより構成されるものである。
【００５２】
＜信号分離方法＞
図５は、本形態の信号分離方法の全体を説明するためのフローチャートである。また、図６は、図５のクラス分類過程（ステップＳ２）の詳細を説明するためのフローチャートである。また、図７（ａ）は、図５のパーミュテーション問題解決過程（ステップＳ３）の詳細を説明するためのフローチャートであり、図７（ｂ）は、図７（ａ）のクラスタリング過程（ステップＳ２２）の詳細を説明するためのフローチャートである。また、図８（ａ）は、図７（ｂ）の大域的最適化過程（ステップＳ３１）の詳細を説明するためのフローチャートであり、図８（ｂ）は、図７（ｂ）の局所的最適化過程（ステップＳ３２）の詳細を説明するためのフローチャートである。以下、これらの図を用い、本形態の信号分離方法を説明する。なお、各演算は、制御部１６０の制御のもと実行される。また、明示しないか限り、各演算過程で生成された演算結果は逐一メモリ１７０に格納され、必要に応じて読み出されて他の演算に用いられる。
【００５３】
［信号分離方法の全体］
まず、源信号の混合信号がＭ（Ｍ≧２）箇所のセンサでそれぞれ観測されて得られた観測信号xm(t)(m=1,...,M、tは時刻）が信号分離装置１０（図１）のメモリ１００の記憶領域１０１に格納される。なお、観測信号xm(t)は、サンプリング周波数fsでサンプリングされた離散値である。
【００５４】
次に、周波数領域変換部１１０が、記憶領域１０１から観測信号xm(t)を読み込み、それらを周波数領域の観測信号Xm(f,τ)（fは周波数、τは時間インデックス）に変換し、生成された周波数領域の観測信号Xm(f,τ)をｆ，τとの対応関係が特定可能な状態でメモリ１００の記憶領域１０２に格納する（周波数領域変換過程／ステップＳ１）。
【００５５】
次に、クラス分類部１２０が、記憶領域１０２から周波数領域の観測信号Xm(f,τ)を読み込み、これらを要素とする観測信号ベクトルX(f,τ)=[X1(f,τ),...,XM(f,τ)]Tをスパース性に基づき周波数f毎に独立にクラスタリングした場合に観測信号ベクトルX(f,τ)が属するクラスがCn(f) (n=1,...,N、N≧１)となる事象の事後確率P(Cn(f)|X(f,τ))を算出する。算出された事後確率P(Cn(f)|X(f,τ))は、ｎ，ｆ，τとの対応関係が特定可能な状態でメモリ１００の記憶領域１０３に格納される（クラス分類過程／ステップＳ２）。
【００５６】
次に、パーミュテーション問題解決部１３０が、記憶領域１０３から事後確率P(Cn(f)|X(f,τ))を読み込み、同じ信号源に対応するクラスがすべての周波数に渡って同じ添字（クラスCn(f)の番号n）を持つように、クラスCn(f)の番号nを付け替える。すなわち、パーミュテーション問題解決部１３０は、対応する周波数fが異なる事後確率P(Cn(f)|X(f,τ))間の類似度を指標として、事後確率P(Cn(f)|X(f,τ))とクラスCn(f)の番号nとの対応関係を並び替え、対応するクラスの番号が同一であって周波数が異なる事後確率間の類似度の総和が当該並び替え前よりも大きな事後確率P’(Ck(f)|X(f,τ)) (k=1,...,N)を生成する。生成された事後確率P’(Ck(f)|X(f,τ))は、ｋ，ｆ，τとの対応関係が特定可能な状態でメモリ１００の記憶領域１０４に格納される（パーミュテーション問題解決過程／ステップＳ３）。
【００５７】
次に、分離部１４０が、記憶領域１０２から周波数領域の観測信号Xm(f,τ)を読み込み、記憶領域１０４から事後確率P’(Ck(f)|X(f,τ))を読み込み、それを用い、周波数領域の分離信号Yn(f,τ)を抽出する。すなわち、分離部１４０が、事後確率P’(Ck(f)|X(f,τ))の大きさを指標とし、クラスCk(f)に属すると判定される周波数領域の観測信号Xm(f,τ)を周波数領域の分離信号Yk(f,τ)として抽出する。抽出された分離信号Yk(f,τ)は、ｆ，τとの対応関係が特定可能な状態でメモリ１００の記憶領域１０５に格納される（分離過程／ステップＳ４）。
【００５８】
最後に、時間領域変換部１５０が、記憶領域１０５から分離信号Yk(f,τ)を読み込み、分離信号Yk(f,τ)を時間領域の分離信号yk(t)に変換してメモリ１００の記憶領域１０６に格納する（時間領域変換過程／ステップＳ５）。
【００５９】
［周波数領域変換過程（ステップＳ１）の詳細］
周波数領域変換過程（ステップＳ１）は、例えば、式(2)に従った短時間フーリエ変換によって行う。
【００６０】
［クラス分類過程（ステップＳ２）の詳細］
本形態のクラス分類過程では、クラスCn(f)の代表ベクトルであるセントロイドan(f)と観測信号ベクトルX(f,τ)との距離に基づいて事後確率のモデルP(Cn(f)|X(f,τ), θ(f))（θ(f)はパラメータ集合）を生成し、事後確率P(Cn(f)|X(f,τ), θ(f))の計算とパラメータ集合θ(f)の推定とを、所定の終了条件を満たすまで交互に繰り返し、事後確率P(Cn(f)|X(f,τ))を求める。
【００６１】
以下では、混合ガウス分布で観測信号ベクトルX(f,τ)の集合をモデル化する方法を例示する。まず、ガウス分布の混合数を、仮定した源信号の数Nとし、観測信号ベクトルX(f,τ)がクラスCn(f)に属する事象の確率密度関数を、ガウス分布
【００６２】
【数１１】
0004769238.tif 000012
でモデル化する。なお、an(f)は、ｎ番目のクラスCn(f)に属する観測信号ベクトルX(f,τ)のセントロイドである。本形態では、an(f)のノルムが所定値（例えば１）に正規化されている。また、σn(f)は、ｎ番目のクラスCn(f)に属する観測信号ベクトルX(f,τ)の標準偏差であり、(σn(f))2は、ｎ番目のクラスCn(f)に属する観測信号ベクトルX(f,τ)の分散である。また、・Hは・の複素共役転置である。また、‖・‖は・のノルムを示す。
【００６３】
ここで、{ (an(f))H・X(f,τ)}・an(f)は、セントロイドan(f)が張る部分空間への観測信号ベクトルX(f,τ)の直交射影である。そのため、式(11)の‖X(f,τ)-{ (an(f))H・X(f,τ)}・an(f)‖は、観測信号ベクトルX(f,τ)とセントロイドan(f)が張る部分空間との最少距離を示し、これが小さいほど観測信号ベクトルX(f,τ)がクラスCn(f)に属する尤度が高くなる。
【００６４】
次に、各ガウス分布p(X(f,τ)|an(f),σn(f))の混合比をαn(f)（0<αn(f)<1とα1+…+αＮ =1とを満たす）とし、パラメータ集合をθ(f)={a1(f),σ1(f),α1(f),...,aN(f),σN(f),αN(f)}とすると、混合ガウス分布による密度関数は、
p(X(f,τ)|θ(f))=Σn=1Nαn(f)・p(X(f,τ)|an(f),σn(f))     …(12)
と表現される。
【００６５】
ここで、周波数f毎にＴ個の観測信号ベクトルX(f,τ)のサンプルが得られたとする。最尤推定の原理では、対数尤度の和
ΣτTlog p(X(f,τ)|θ(f))=ΣτTlogΣn=1Nαn(f)・p(X(f,τ)|an(f),σn(f))
を最大化するパラメータ集合θ(f)を求めるが、この形では、対数の中に確率密度関数p(X(f,τ)|an(f),σn(f))の和が含まれているため、計算が困難となる。そこで、本形態では、EMアルゴリズム（例えば、「汪金芳,田栗正章,手塚集,樺島祥介,上田修功,「計算統計I確率計算の新しい手法」,統計科学のフロンティア11,ISBN4-00-006851-2」等参照）を用いてパラメータ推定を行う。EMアルゴリズムでは、対数尤度の代わりにいわゆるQ関数
Q(f,θ(f))=ΣτTΣn=1N {P(Cn(f)|X(f,τ),θ(f))・logαn(f)・p(X(f,τ)|an(f),σn(f))}    …(13)
を最大化するパラメータ集合θ(f)を求める。ここで、P(Cn(f)|X(f,τ),θ(f))は、観測信号ベクトルX(f,τ)を周波数f毎に独立にクラスタリングした場合に観測信号ベクトルX(f,τ)が属するクラスがCn(f)となる事象の、観測信号ベクトルX(f,τ)を得た後における事後確率であり、ベイズの定理により、
P(Cn(f)|X(f,τ),θ(f))=αn(f)・p(X(f,τ)|an(f),σn(f))/p(X(f,τ)|θ(f))
…(14)
と書き下せる。
【００６６】
すなわち、本形態のクラス分類過程では、現在のパラメータ集合θ(f)を固定したまま、式 (11)(12)(14)を用いて、すべての観測信号ベクトルX(f,τ)とクラスCn(f)に対し、周波数f毎に事後確率P(Cn(f)|X(f,τ),θ(f))を計算する事後確率計算過程（E-step）と、事後確率P(Cn(f)|X(f,τ),θ(f))を固定したまま、式(13)のQ関数Q(f,θ(f))が最大となるパラメータ集合θ(f)を計算するパラメータ推定過程（M-step）とを、所定の終了条件を満たすまで繰り返し、終了条件を満たした時点の事後確率P(Cn(f)|X(f,τ),θ(f))を事後確率P(Cn(f)|X(f,τ))として出力する。
【００６７】
なお、この例のパラメータ推定過程でのセントロイドan(f)は、クラスCn(f)に関する相関行列
R=ΣτT P(Cn(f)|X(f,τ),θ(f))・X(f,τ)・X H(f,τ)
の最大固有値として算出される。また、分散(σn(f))2は、
【００６８】
【数１２】
0004769238.tif 000013
として計算される。また、混合比αn(f)は、
【００６９】
【数１３】
0004769238.tif 000014
として計算される。
【００７０】
また、事後確率計算過程とパラメータ推定過程とは、所定の終了条件を満たすまで繰り返されるが、その最終結果は初期値によって異なったものになることがあるため、初期値の設定は重要である。一般的には、狭い範囲に集中するセントロイドを初期値に設定したり、実際のサンプルからあまりにも乖離した初期値を設定したりすることは避けるべきである。適切にクラスタリングが行われない場合があるからである。好ましい初期値の設定方法には、特に制限はないが、例えば、以下のような方法を例示できる。
【００７１】
まず、セントロイドa1(f),...,aN(f)の初期値には、例えば、観測信号ベクトルX(f,τ)のＴ個のサンプルからランダムに選択したN個のサンプルを用いる。また、分散(σn(f))2の初期値には、例えば、観測信号ベクトルX(f,τ)毎にセントロイドan(f)との２乗距離の最小値（観測信号ベクトルX(f,τ)と何れかのセントロイドan(f)との２乗距離）を求め、それらを平均した値
【００７２】
【数１４】
0004769238.tif 000015
を用いる。また、混合比αn(f) の初期値は、例えば、αn(f)=1/Nと設定する。
【００７３】
また、上述のようにセントロイドan(f)と観測信号ベクトルX(f,τ)との距離に基づいて事後確率をモデル化する場合、その最尤推定に用いるサンプルである観測信号ベクトルX(f,τ)のノルム‖X(f,τ)‖は所定値（例えば１）に正規化されていることが望ましい。本来、セントロイドan(f)と観測信号ベクトルX(f,τ)との距離に基づいてクラス分類を行う場合、それらのノルムが各n、f、τにおいて一定値に正規化されていないと、セントロイドan(f)と観測信号ベクトルX(f,τ)との距離を、各n、f、τに渡って厳密に比較評価し、クラス分類を行うことはできない。すなわち、セントロイドan(f)と観測信号ベクトルX(f,τ)との距離に基づいて事後確率をモデル化する場合、サンプルである観測信号ベクトルX(f,τ)のノルムが大きいほどセントロイドan(f)との距離が大きくなり、適切な最尤推定を行うことができない場合がある。これは、事後確率P(Cn(f)|X(f,τ))の推定精度を低下させる。よって、本形態では、尤推定に用いるサンプルである観測信号ベクトルX(f,τ)のノルムを、例えば、
【００７４】
【数１５】
0004769238.tif 000016
により１に正規化する。なお、α←βは、βの値をαの値とすることを意味する。
【００７５】
以下、本形態のクラス分類過程（ステップＳ２）の詳細を図６のフローチャートに沿って説明する。なお、ここでは周波数fに関する処理のみを説明するが、クラス分類過程（ステップＳ２）の処理は周波数f毎に独立に行われる。
【００７６】
まず、クラス分類部１２０のノルム正規化部１２１（図２）にメモリ１００から読み込まれた観測信号ベクトルX(f,τ)が入力され、ノルム正規化部１２１が当該観測信号ベクトルX(f,τ)のノルムを所定値に正規化する。具体的には、ノルム正規化部１２１は、例えば、式(16)に従って観測信号ベクトルX(f,τ)のノルムを１に正規化する。また、正規化された測信号ベクトルX(f,τ)は、メモリ１７０に格納される（ノルム正規化過程／ステップＳ１１）。
【００７７】
次に初期パラメータ設定部１２２ａが、メモリ１７０からノルムが正規化された観測信号ベクトルX(f,τ)を読み込み、初期パラメータ設定部１２２ａは、これらを用いてパラメータ集合θ(f)={a1(f),σ1(f),α1(f),...,aN(f),σN(f),αN(f)}の初期値を設定する。この初期値の設定は、例えば、先に例示した初期値の設定方法に従って行う。また、設定されたパラメータ集合θ(f)の初期値は、メモリ１７０に格納される（初期パラメータ設定過程／ステップＳ１２）。なお、Nの値については、信号源数が既知なのであればその数をNとして用いてもよいし、厳密な信号源数が未知の場合には、経験則等から仮定できる信号源数をNとして用いてもよい。
【００７８】
次に、事後確率計算部１２２ｂは、メモリ１７０からノルムが正規化された観測信号ベクトルX(f,τ)と最新のパラメータ集合θ(f)とを読み込み、パラメータ集合θ(f)を固定値として、式 (11)(12)(14)を用いて、すべての観測信号ベクトルX(f,τ)とクラスCn(f)に対し、周波数f毎に事後確率P(Cn(f)|X(f,τ),θ(f))を算出する。算出された事後確率P(Cn(f)|X(f,τ),θ(f))は、メモリ１７０に格納される（事後確率計算過程／ステップＳ１３）。
【００７９】
次に、パラメータ推定部１２２ｃが、メモリ１７０から観測信号ベクトルX(f,τ)と最新の事後確率P(Cn(f)|X(f,τ),θ(f))を読み込み、読み込んだ事後確率P(Cn(f)|X(f,τ),θ(f))を固定値として、式(13)に示したQ関数Q(f,θ(f))〔「各観測信号ベクトルX(f,τ)に対応するΣn=1N {P(Cn(f)|X(f,τ),θ(f))・logαn(f)・p(X(f,τ)|an(f),σn(f))}を周波数f毎に独立に加算したQ(f,θ(f))」に相当〕がそれぞれ最大となるパラメータ集合θ(f)を算出する。この算出方法は、前述した通りである。算出されたパラメータ集合θ(f)は、メモリ１７０に格納される（パラメータ推定過程／ステップＳ１４）。
【００８０】
次に、演算制御部１２２ｄが、所定の終了条件を満たした否かを判定する（終了条件判定過程／ステップＳ１５）。なお「所定の終了条件」としては、例えば、以下を例示できる。
【００８１】
・ステップＳ１３で固定値として用いたパラメータ集合θ(f)からなるベクトルと、ステップＳ１４で新たに算出されたパラメータ集合θ(f)からなるベクトルとの距離が所定値以下（又は未満）であること。
【００８２】
・ステップＳ１３で算出された最新の事後確率P(Cn(f)|X(f,τ),θ(f))と、それよりも１つ前に算出された事後確率P(Cn(f)|X(f,τ),θ(f))（ステップＳ１２又は１つ前のループのステップＳ１３で作成された事後確率）との差の合計が所定値以下（又は未満）であること。
【００８３】
・ステップＳ１３とＳ１４の処理を所定回数繰り返したこと。
【００８４】
ここで、所定の終了条件を満たしていないと判定された場合、処理がステップＳ１３に戻される。一方、所定の終了条件を満たしたと判定された場合、すべてのCn(f)及びX(f,τ)にそれぞれ対応する最新の事後確率P(Cn(f)|X(f,τ),θ(f))がメモリ１７０からクラス分類部１２０に読み込まれ、これらが事後確率P(Cn(f)|X(f,τ))として出力される（事後確率出力過程／ステップＳ１６）。  
【００８５】
［パーミュテーション問題解決過程（ステップＳ３）の詳細］
次に、図７（ａ）を用い、パーミュテーション問題解決過程（ステップＳ３）の詳細を説明する。
【００８６】
まず、パーミュテーション問題解決部１３０（図２（ｂ））のアクティブ系列生成部１３１にメモリ１００から読み込まれた各事後確率P(Cn(f)|X(f,τ))が入力される。アクティブ系列生成部１３１は、式(9)に従い、各事後確率P(Cn(f)|X(f,τ))に対応するアクティブ系列vnf(τ)を生成して出力する（アクティブ系列生成過程／ステップＳ２１）。
【００８７】
各アクティブ系列vnf(τ)はクラスタリング部１３２に入力され、クラスタリング部１３２はそれらのクラスタリングを行う。ここでのクラスタリングは一般的なものとは少し異なり、クラスタリング結果は周波数毎の順列Πfとして表現される。より具体的には、クラスタリング部１３２は、対応する周波数fが異なる系列vnf(τ)間の類似度を指標として、事後確率P(Cn(f)|X(f,τ))とクラスCn(f)の番号nとの対応関係を並び替える順列Πfを周波数f毎に生成する。生成された順列Πfはメモリ１７０に格納される（クラスタリング過程／ステップＳ２２）。
【００８８】
次に、並び替え部１３３がメモリ１００から事後確率P(Cn(f)|X(f,τ))を読み込み、メモリ１７０から順列Πfを読み込む。そして、並び替え部１３３は、順列Πfに従い、周波数f毎に事後確率P(Cn(f)|X(f,τ))とクラスCn(f)の番号nとの対応関係を並び替え、事後確率P’(Ck(f)|X(f,τ))を生成する。生成された事後確率P’(Ck(f)|X(f,τ))はメモリ１７０に格納される（並び替え過程／ステップＳ２３）。
【００８９】
［クラスタリング過程（ステップＳ２２）の詳細］
次に、図７（ｂ）を用い、クラスタリング過程（ステップＳ２２）の詳細について説明する。
【００９０】
まず、クラスタリング部１３２の大域的最適化部１３２ａにアクティブ系列vnf(τ)が入力される。大域的最適化部１３２ａは、すべての異なる周波数f∈Fの組合せに対応するアクティブ系列vnf(τ)間の類似度を指標とし、事後確率P(Cn(f)|X(f,τ))とクラスCn(f)の番号nとの対応関係を並び替える順列Πfを周波数f毎に生成する。なお、Fは取り扱う全周波数ビンの集合を意味する。生成された順列Πfはメモリ１７０に格納される（大域的最適化過程／ステップＳ３１）。ここで、前述したように、アクティブ系列vnf(τ)間の類似度は、多くの周波数の組み合わせに関し、同じ信号源に対応する系列間の類似度が大きくなるという傾向をもつ。また、大域的最適化過程では、すべての異なる周波数f∈Fの組合せに対応する系列vnf(τ)間の類似度を指標として順列Πfを生成する。よって、この大域的最適化過程により、周波数全体に渡って一貫性のあるパーミュテーション問題の解となる順列Πfを求めることができる。
【００９１】
次に、局所的最適化部１３２ｂにアクティブ系列vnf(τ)と大域的最適化過程（ステップＳ３１）で生成された順列Πfとが入力される。局所的最適化部１３２ｂは、対応する周波数f∈Fが異なる系列vnf(τ)間の類似度のうち、特定の周波数の組合せに対応する系列間の類似度のみを指標として用い、大域的最適化部で生成された順列Πfを更新し、新たな順列Πfを生成する。生成された順列Πfはメモリ１７０に格納される（局所的最適化過程／ステップＳ３２）。ここで、「特定の周波数の組合せ」として、同じ信号源に対応するアクティブ系列vnf(τ)間の類似度が特に大きくなる周波数の組み合わせを選択することにより、大域的最適化過程で生成された順列Πfを、パーミュテーション問題をより高精度に解決できる順列Πfに補正することができる。なお、「特定の周波数の組合せ」としては、周波数差が所定範囲内にある周波数の組合せや、倍音関係にある周波数の組合せを例示できる。
【００９２】
このように、本形態のクラスタリング過程では、大域的最適化過程（ステップＳ３１）と局所的最適化過程（ステップＳ３２）を順に適用して、各周波数fでの順列Πfを算出することとしたため、周波数全体に渡って一貫性のある高精度な順列Πfを得ることができる。
【００９３】
［大域的最適化過程（ステップＳ３１）の詳細］
次に、大域的最適化過程（ステップＳ３１）の詳細を例示する。
【００９４】
この例では、異なる周波数f∈Fの組合せに対応するアクティブ系列vnf(τ)間の類似度を直接的に指標として用いるのではなく、各信号源に対応するセントロイドcn(τ)とアクティブ系列vnf(τ)との類似度を指標とすることで、異なる周波数f∈Fの組合せに対応するアクティブ系列vnf(τ)間の類似度を間接的に指標として用いる。これにより、大域的最適化過程の演算精度と演算効率が向上する。
【００９５】
具体的には、この例の大域的最適化過程では、信号源毎にセントロイドcn(τ)を推定し、目的関数
【００９６】
【数１６】
0004769238.tif 000017
を最大化する順列Πfを求める。この目的関数は、順列Πfによって並び替えられたアクティブ系列vkf(τ)｜k=Πf(n)とｎ番目の信号源に対応するセントロイドcn(τ)との類似度をすべての信号源及び全ての周波数で足し合わせたものである。
【００９７】
なお、ρ(vkf,cn)は、アクティブ系列vnf(τ)とセントロイドcn(τ)との類似度を示す関数値であり、例えば、アクティブ系列vnf(τ)とセントロイドcn(τ)との相関係数（式(7)参照）である。しかし、アクティブ系列vnf(τ)とセントロイドcn(τ)との類似度を示すのであれば、別の関数値をρ(vkf,cn)として用いてもよい。例えば、
【００９８】
【数１７】
0004769238.tif 000018
としてもよいし、その他のアクティブ系列vnf(τ)とセントロイドcn(τ)との距離D(vkf(τ),cn(τ))に対して単調減少の関係にある関数をρ(vkf,cn)としてもよい。また、系列長１のアクティブ系列vnf(τ)とセントロイドcn(τ)との類似度を示す関数値をρ(vkf,cn)として用いてもよい。例えば、
ρ(vkf,cn)= vkf(τ)・cn(τ)     …(19)
としてもよい。
【００９９】
式(17)の目的関数は、よく知られたk-means法（例えば、「R.O. Duda, P. E. hart, and D. G. Stork, Pattern Classification, Wiley Interscience, 2nd edition, 2000」等参照）と同じように、セントロイドcn(τ)と順列Πfとを交互に最適化することで最大化することができる。以下、図８（ａ）を用い、本形態の大域的最適化過程（ステップＳ３１）の詳細を例示する。
【０１００】
まず、クラスタリング部１３２（図３）の初期パラメータ設定部１３２ａａにメモリ１７０から読み込まれたアクティブ系列vnf(τ)が入力され、初期パラメータ設定部１３２ａａは、これらのアクティブ系列vnf(τ)を用いてセントロイドcn(τ)の初期値を設定する。初期パラメータ設定部１３２ａａは、例えば、読み込まれたアクティブ系列vnf(τ)から時間インデックスτ毎にN個のサンプルを選択し、それらをN個のセントロイドcn(τ)の初期値とする。生成されたセントロイドcn(τ)の初期値はメモリ１７０に格納される（初期パラメータ設定過程／ステップＳ４１）。
【０１０１】
次に、順列生成部１３２ａｂが、メモリ１７０からアクティブ系列vnf(τ)とセントロイドcn(τ)とを読み込み、すべての周波数f∈Fに対応するアクティブ系列vnf(τ)とセントロイドcn(τ)との類似度を指標とし、順列Πfを周波数f毎に生成する。本形態の順列生成部１３２ａｂは、周波数f毎に、アクティブ系列vnf(τ)とセントロイドcn(τ)との類似度を最大化させる順列Πfを以下ように決定する。なお、argmaxΠαは、αを最大にする順列Πを意味する。
【０１０２】
【数１８】
0004769238.tif 000019
そして、このように生成された順列Πfはメモリ１７０に格納される（順列生成過程／ステップＳ４２）。
【０１０３】
次に、セントロイド算出部１３２ａｃが、メモリ１７０からアクティブ系列vnf(τ)(n=1,...,N)と最新の順列Πfとを読み込み、順列Πfに従って周波数f毎にアクティブ系列vnf(τ)と番号nとの対応関係を並び替えたアクティブ系列vkf(τ)(k=1,...,N)の周波数方向の平均値又は代表値をセントロイドck(τ)として算出する。なお、並び替えたアクティブ系列vkf(τ)の周波数方向の平均値をセントロイドck(τ)とする場合、セントロイド算出部１３２ａｃは、例えば、
【０１０４】
【数１９】
0004769238.tif 000020
によってセントロイドck(τ)を算出する。ここで、|F|は取り扱う全周波数ビンの集合Fの要素数を意味する。また、並び替えたアクティブ系列vkf(τ)の周波数方向の代表値をセントロイドck(τ)とする場合、セントロイド算出部１３２ａｃは、例えば、kとτとが同じアクティブ系列vkf(τ)の集合毎にアクティブ系列vkf(τ)の振幅の頻度を求め、当該集合内で頻度が最大となる振幅を持つ何れかのサンプルを、k，τに対応するセントロイドck(τ)とする。
【０１０５】
そして、このように生成されたセントロイドck(τ)はメモリ１７０に格納される（セントロイド算出過程／ステップＳ４３）。
【０１０６】
次に、演算制御部１３２ａｄが、所定の終了条件を満たしたか否かを判定する（終了条件判定過程／ステップＳ４４）。なお「所定の終了条件」としては、例えば、以下を例示できる。
【０１０７】
・ステップＳ４３で生成されたセントロイドck(τ)と、前回のループで生成されたセントロイドck(τ)（前回のループで生成されたセントロイドck(τ)が存在しない場合はセントロイドck(τ)の初期値） との距離が所定値以下（又は未満）であること。
【０１０８】
・ステップＳ４２で生成された順列Πfと前回のループで生成された順列Πfとの相違箇所が所定個以下（又は未満）であること。
【０１０９】
・ステップＳ４２とＳ４３の処理を所定回数繰り返したこと。
【０１１０】
ここで、所定の終了条件を満たしていないと判定された場合、処理がステップＳ４２に戻される。一方、所定の終了条件を満たしていると判定された場合、大域的最適化部１３２ａは、メモリ１７０に格納されている最新の順列Πfを出力する（順列出力過程／ステップＳ４５）。
【０１１１】
［局所的最適化過程（ステップＳ３２）の詳細］
次に、図８（ｂ）を用い、局所的最適化過程（ステップＳ３２）の詳細を例示する。
【０１１２】
まず、局所的最適化部１３２ｂ（図３）の順列生成部１３２ｂaに、順列Πfとアクティブ系列vnf(τ)とが入力される。なお、１回目のループの際に入力される順列Πfは大域的最適化過程（ステップＳ３１）で出力された順列であり、２回目以降のループの際に入力される順列Πfは前回の順列生成過程（ステップＳ４６）で生成された順列である。
【０１１３】
順列生成部１３２ｂaは、入力された順列Πfの一部である順列Πg(g∈R(f))を固定し、順列Πgに従ってアクティブ系列vng(τ)と番号nとの対応関係を並び替えたアクティブ系列vk’g(τ)と、順列Πに従ってアクティブ 系列vnf(τ)と番号nとの対応関係を並び替えたアクティブ系列vkf(τ)との類似度の全信号源についての総和を最大にする、当該順列Πを、新たな順列Πfとして算出する（順列生成過程／ステップＳ４６）。この処理は各周波数fについて行われ、例えば、以下の式に従って行われる。 
【０１１４】
【数２０】
0004769238.tif 000021
なお、R(f)は、周波数fと特定の関係にある周波数gの集合であり、周波数fとgとの組み合わせが「特定の周波数の組合せ」に相当する。すなわち、R(f)は、アクティブ系列vkf(τ)とvk’g(τ)とが同じ信号源に相当するものであったときに、これらの類似度が特に大きくなるような周波数gを要素に持つ集合であることが望ましい。典型的には、周波数fの近傍の周波数の集合A(f)と、周波数fと倍音関係にある周波数の集合H(f)とに対し、R(f)= A(f)∪H(f)であることが望ましい。また、周波数fの近傍の周波数の集合A(f)は、例えば、
Α(f)={f-3Δf, f-2Δf, f-Δf, f+Δf, f+2Δf, f+3Δf}
として定義できる。ここで、Δf＝(1/L)fsは、隣り合う周波数ビン間の周波数の差である。  また、倍音関係にある周波数の集合Ηは、例えば、
Η(f)={round(f/2)-Δf, round(f/2), round(f/2)+Δf, 2f-Δf, 2f, 2f+Δf}
として定義できる。ここで、round(・)は、周波数の集合Fから・に最も近い周波数を意味する。その他、R(f)=A(f)又はR(f)=H(f)とする構成も可能である。
【０１１５】
次に、演算制御部１３２ｂｂが、所定の終了条件を満たしたか否かを判定する（終了条件判定過程／ステップＳ４７）。なお「所定の終了条件」としては、例えば、以下を例示できる。
【０１１６】
・ステップＳ４６で更新された順列Πfの更新箇所が所定数以下（又は未満）であること。
【０１１７】
・ステップＳ４６の処理を所定回数繰り返したこと。
【０１１８】
ここで、所定の終了条件を満たしていないと判定された場合、処理がステップＳ４６に戻される。一方、所定の終了条件を満たしていると判定された場合、局所的最適化部１３２ｂは最新の順列Πfを出力する（順列出力過程／ステップＳ４８）。
【０１１９】
［並び替え過程（ステップＳ２３）の詳細］
前述のように、並び替え部１３３は、順列Πfに従い、周波数f毎に事後確率P(Cn(f)|X(f,τ))とクラスCn(f)の番号nとの対応関係を並び替え、事後確率P’(Ck(f)|X(f,τ))を生成する。具体的には、例えば、
P’(Ck(f)|X(f,τ))←P(Cn(f)|X(f,τ))|n=Πf(k), ∀k,f,τ    …(24)
に従い、事後確率P’(Ck(f)|X(f,τ))を生成する（パーミュテーション問題解決過程（ステップS３）の説明終わり）。
【０１２０】
［分離過程（ステップＳ４）の詳細］
前述のように、分離過程では、分離部１４０（図１）が、事後確率P’(Ck(f)|X(f,τ))の大きさを指標とし、クラスCk(f)に属すると判定される周波数領域の観測信号Xm(f,τ)を周波数領域の分離信号Yk(f,τ)として抽出する。例えば、以下に例示する時間周波数マスキングによる方法が、比較的簡単な分離方法となる。なお、m’=1,...,Mであり、観測信号X m’(f,τ)はセンサm’に対応する観測信号である。
【０１２１】
【数２１】
0004769238.tif 000022
その他、L1ノルム最小化規範やその近似方法による分離方法を用いてもよい（例えば、特許文献１や非特許文献４等参照）。
【０１２２】
［時間領域変換過程（ステップＳ５）の詳細］
最後に、時間領域変換部１５０が、分離信号Yk(f,τ)を時間領域の分離信号yk(t)に変換する。この処理は、例えば、短時間逆フーリエ変換等によって行う。
【０１２３】
〔第２実施形態〕
次に、本発明の第２実施形態について説明する。
【０１２４】
本形態は第１実施形態のクラス分類部及びクラス分類過程の変形例である。以下では、クラス分類部及びクラス分類過程の相違点のみを説明する。
【０１２５】
＜クラス分類部及びクラス分類過程の変形例１＞
図９は、第１実施形態のクラス分類部１２０の変形例であるクラス分類部２２０の機能構成を示したブロック図である。
【０１２６】
クラス分類部１２０とクラス分類部２２０との構成上の相違点は、クラス分類部１２０がパラメータ集合θ(f)の初期値を設定する初期パラメータ設定部１２２ａを具備するモデル化部１２２を有していたのに対し、クラス分類部２２０が事後確率P(Cn(f)|X(f,τ))の初期値を設定する初期パラメータ設定部２２２ａを具備するモデル化部２２２を有する点である。
【０１２７】
また、クラス分類部１２０とクラス分類部２２０との処理上の相違点は、クラス分類部１２０では初期パラメータ設定部１２２ａがステップＳ１２（図６）でパラメータ集合θ(f)の初期値を生成していたのに対し、クラス分類部２２０では初期パラメータ設定部２２２ａが事後確率P(Cn(f)|X(f,τ))の初期値を生成する点と、ステップＳ１２とＳ１３との順序が逆になる点である。
【０１２８】
なお、事後確率P(Cn(f)|X(f,τ))の初期値の設定方法については特に制限はないが、広範囲に分布する観測信号ベクトルX(f,τ)が同一のクラスCn(f)に属するような初期値設定は好ましくない。適切なクラス分類がなされない場合があるからである。好ましい初期値の設定方法には、特に制限はないが、例えば、以下の(1)〜(4)のような手順を例示できる。
【０１２９】
(1)各観測信号ベクトルX(f,τ)のクラス分けパターンをランダムに複数パターン生成する。
【０１３０】
(2)各クラス分けパターンに対し、クラスCn(f)毎のセントロイドの組み合わせを算出する。
【０１３１】
(3)セントロイドの組み合わせに毎に、異なるクラスのセントロイド間の内積を求め、内積が最小となる組み合わせを選択する。
【０１３２】
(4)選択されたセントロイドの組み合わせに対応するクラス分けパターンに従い、事後確率P(Cn(f)|X(f,τ))の初期値を設定する。例えば、そのクラス分けパターンに従うと観測信号ベクトルX(f,τ)がクラスCn(f)に属することになる場合、P(Cn(f)|X(f,τ))=1とし、P(Cn'(f)|X(f,τ))=0  (n'≠n)として事後確率P(Cn(f)|X(f,τ))の初期値を設定する。
【０１３３】
＜クラス分類部及びクラス分類過程の変形例２＞
この変形例では、観測信号ベクトルX(f,τ)のノルムを正規化することなくモデル化部の処理を行う。第１実施形態で述べたように、セントロイドan(f)と観測信号ベクトルX(f,τ)との距離に基づいて事後確率をモデル化する場合、サンプルである観測信号ベクトルX(f,τ)のノルムが正規化されていないと事後確率P(Cn(f)|X(f,τ))の推定精度が低下してしまう。しかし、高い推定精度が要求されない用途に用いる場合や観測信号ベクトルX(f,τ)のノルムが安定している場合などには、クラス分類過程での観測信号ベクトルX(f,τ)のノルムの正規化を省略してもよい。
【０１３４】
また、セントロイドan(f)と観測信号ベクトルX(f,τ)の方向のみの類似度に基づいて事後確率をモデル化する場合には、観測信号ベクトルX(f,τ)のノルムを正規化する必要はない。例えば、セントロイドan(f)と観測信号ベクトルX(f,τ)とのコサイン距離
｜XH（f,τ）・an(f)｜/（‖X（f,τ）‖・‖an(f)‖）      …(26)
を用い、前述の式(11)の替わりに、
【０１３５】
【数２２】
0004769238.tif 000023
でモデル化してもよい。なお、式(26)の｜α｜はαの絶対値である。
【０１３６】
図１０（ａ）は、第１実施形態のクラス分類部１２０（図２（ａ））からノルム正規化部１２１を排除したクラス分類部３２０の機能構成を示したブロック図である。また、図１０（ｂ）は、図９に示したクラス分類部２２０からノルム正規化部１２１を排除したクラス分類部４２０の機能構成を示したブロック図である。このようにクラス分類部がノルム正規化部を具備しない構成であってもよい。
【０１３７】
＜クラス分類部及びクラス分類過程の変形例３＞
この変形例では、第１実施形態のモデル化部１２２のように最尤推定によって事後確率P(Cn(f)|X(f,τ))を推定するのではなく、よく知られたK-means法によって観測信号ベクトルX(f,τ)をクラスタリングして各観測信号ベクトルX(f,τ)が属するクラスCn(f)を１つずつ推定する。そして、観測信号ベクトルX(f,τ)がクラスCn(f)に属する場合の事後確率P(Cn(f)|X(f,τ))を１とし、観測信号ベクトルX(f,τ)がクラスCn(f)に属しない場合の事後確率P(Cn(f)|X(f,τ))を０とする。すなわち、この場合の事後確率P(Cn(f)|X(f,τ))は０と１のみをとる。
【０１３８】
図１１は、このようなクラス分類部５２０の変形例を示したブロック図であり、図１２は、このような場合のクラス分類過程を説明するためのフローチャートである。なお、これまで説明したのと同じ構成部分については同じ符号を付した。
【０１３９】
図１１に示すように、この例のクラス分類部５２０は、ノルム正規化部１２１と、初期パラメータ設定部５２２と、事後確率計算部５２３と、セントロイド算出部５２４と、演算制御部５２５とを有する。以下、この例のクラス分類部５２０によって行われるクラス分類過程を説明する。
【０１４０】
まず、ノルム正規化部１２１に観測信号ベクトルX(f,τ)が入力され、ノルム正規化部１２１は、各観測信号ベクトルX(f,τ)のノルムを所定値（例えば１）に正規化した各観測信号ベクトルX (f,τ)を生成してメモリ１７０（図１）に格納する（ノルム正規化過程／ステップＳ６１）。
【０１４１】
次に、初期パラメータ設定部６２２がメモリ１７０からノルムが正規化された観測信号ベクトルX(f,τ)を読み込み、これらを用い、周波数f毎に各クラスCi(f)(i=1,...,N)のセントロイドci(f)の初期値を設定してメモリ１７０に格納する（初期パラメータ設定過程／ステップＳ６２）。例えば、これらの観測信号ベクトルX(f,τ)からランダムにN個のサンプルを選択し、それらをセントロイドci(f)の初期値とする。
【０１４２】
次に、事後確率計算部５２３が、メモリ１７０からセントロイドci(f)とノルムが正規化された観測信号ベクトルX(f,τ)とを読み込み、各セントロイドci(f)を固定値として、観測信号ベクトルX(f,τ)毎に、 観測信号ベクトルX(f,τ)との距離が最も近いセントロイドcB(f)(B⊂{1,...,N})を選択し、 
P(Cn(f)|X(f,τ))=1  (if n=B)
P(Cn(f)|X(f,τ))=0  (if n≠B)
として、事後確率を算出してメモリ１７０に格納する（事後確率計算過程／ステップＳ６３）。
【０１４３】
次に、セントロイド算出部５２４が、メモリ１７０から事後確率P(Cn(f)|X(f,τ))と正規化された観測信号ベクトルX(f,τ)とを読み込み、読み込んだ事後確率P(Cn(f)|X(f,τ))を固定値として、各クラスCn(f)のセントロイドcn(f)を算出してメモリ１７０に格納する（セントロイド算出過程／ステップＳ６４）。例えば、クラスCn(f)に関する相関行列
R=ΣτT P(Cn(f)|X(f,τ))・X(f,τ)・X H(f,τ)
の最大固有値としてセントロイドcn(f)が算出される。
【０１４４】
次に、演算制御部５２５が、所定の終了条件を満たしたか否かを判定する（終了条件判定過程／ステップＳ６５）。なお「所定の終了条件」としては、例えば、以下を例示できる。
【０１４５】
・ステップＳ６４で生成されたセントロイドcn(f)と、前回のループで生成されたセントロイドcn(f)（前回のループで生成されたセントロイドcn(f)が存在しない場合はセントロイドcn(f)の初期値） との距離が所定値以下（又は未満）であること。
【０１４６】
・ステップＳ６３とＳ６４の処理を所定回数繰り返したこと。
【０１４７】
ここで、所定の終了条件を満たしていないと判定された場合、処理がステップＳ６３に戻される。一方、所定の終了条件を満たしたと判定された場合、すべてのCn(f)及びX(f,τ)にそれぞれ対応する最新の事後確率P(Cn(f)|X(f,τ))が出力される（事後確率出力過程／ステップＳ６６）。
【０１４８】
また、図１３（ａ）のクラス分類部６２０のように、初期パラメータ設定部６２２が、正規化された観測信号ベクトルX(f,τ)を用いて事後確率P(Cn(f)|X(f,τ))の初期値を生成し、セントロイド算出部５２４がステップＳ６４と同様に各クラスCn(f)のセントロイドcn(f)を算出し、事後確率計算部５２３がステップＳ６３と同様に事後確率P(Cn(f)|X(f,τ))を算出し、所定の終了条件を満たすまでセントロイド算出部５２４と事後確率計算部５２３との処理を繰り返すこととしてもよい。
【０１４９】
また、図１３（ｂ）のクラス分類部７２０のように、図１１のクラス分類部５２０からノルム正規化部１２１を排除し、観測信号ベクトルX(f,τ)のノルムを正規化せずに事後確率P(Cn(f)|X(f,τ))を算出してもよい。また、図１３（ｃ）のクラス分類部８２０のように、図１３（ａ）のクラス分類部６２０からノルム正規化部１２１を排除し、観測信号ベクトルX(f,τ)のノルムを正規化せずに事後確率P(Cn(f)|X(f,τ))を算出してもよい。
【０１５０】
この場合、事後確率計算部５２３は、観測信号ベクトルX(f,τ)との距離が最も近いセントロイドcB(f)(B⊂{1,...,N})を選択するのではなく、例えば、測信号ベクトルX(f,τ)毎に、観測信号ベクトルX(f,τ)とのコサイン距離
cosθ=｜XH（f,τ）・cB(f)｜/（‖X（f,τ）‖・‖cB(f)‖） …(28)
が最も近いセントロイドcB(f)を選択し、
P(Cn(f)|X(f,τ))=1  (if n=B)
P(Cn(f)|X(f,τ))=0  (if n≠B)
として、事後確率を算出してメモリ１７０に格納する。
【０１５１】
〔第３実施形態〕
次に、本発明の第３実施形態について説明する。
【０１５２】
本形態は、第１実施形態のパーミュテーション問題解決部が具備するクラスタリング部及びそのクラスタリング過程の変形例である。以下では、クラスタリング部及びクラスタリング過程の相違点のみを説明する。
【０１５３】
＜クラスタリング部及びクラスタリング過程の変形例１＞
図１４（ａ）は、第１実施形態のクラスタリング部１３２（図３）の変形例であるクラスタリング部２３２の機能構成を示したブロック図である。
【０１５４】
クラスタリング部１３２とクラスタリング部２３２との相違点は、クラスタリング部１３２が大域的最適化部１３２ａと局所的最適化部１３２ｂを有していたのに対し、クラスタリング部２３２が大域的最適化部１３２ａを有するが局所的最適化部１３２ｂを有しない点である。すなわち、クラスタリング部２３２は、大域的最適化部１３２ａによって算出された順列Πfをクラスタリング部２３２の出力とする。
【０１５５】
＜クラスタリング部及びクラスタリング過程の変形例１＞
図１４（ｂ）は、第１実施形態のクラスタリング部１３２（図３）の変形例であるクラスタリング部３３２の機能構成を示したブロック図である。また、図１５は、クラスタリング部３３２の局所的最適化部３３２ｂが行う局所的最適化過程を説明するためのフローチャートである。
【０１５６】
クラスタリング部１３２とクラスタリング部３３２との相違点は、クラスタリング部１３２が大域的最適化部１３２ａと局所的最適化部１３２ｂを有していたのに対し、クラスタリング部３３２が局所的最適化部１３２ｂを有するが大域的最適化部１３２ａを有しない点である。すなわち、クラスタリング部３３２は、局所的最適化部１３２ｂのみによって算出された順列Πfをクラスタリング部３３２の出力とする。その相違から、クラスタリング部３３２は、さらに初期順列設定部３３２ｂａを有する。
【０１５７】
以下、図１５を用い、クラスタリング部３３２の局所的最適化部１３２ｂが行う局所的最適化過程を説明する。
【０１５８】
まず、初期順列設定部３３２ｂａが初期順列Πfを生成してメモリ１７０（図１）に格納する（初期順列設定過程／ステップＳ７１）。次に、局所的最適化部３３２ｂ（図１４（ｂ））の順列生成部１３２ｂaに、順列Πfとアクティブ系列vnf(τ)とが入力される。なお、１回目のループの際に入力される順列Πfは初期順列設定過程（ステップＳ７１）で出力された初期順列Πfであり、２回目以降のループの際に入力される順列Πfは前回の順列生成過程（ステップＳ７２）で生成された順列である。
【０１５９】
順列生成部１３２ｂaは、第１実施形態の順列生成過程（ステップＳ４６）と同様に新たな順列Πfとして算出する（順列生成過程／ステップＳ７２）。
【０１６０】
次に、演算制御部１３２ｂｂが、第１実施形態の終了条件判定過程（ステップＳ７３）と同様に、所定の終了条件を満たしたか否かを判定する（終了条件判定過程／ステップＳ７３）。
【０１６１】
ここで、所定の終了条件を満たしていないと判定された場合、処理がステップＳ７２に戻される。一方、所定の終了条件を満たしていると判定された場合、局所的最適化部３３２ｂは最新の順列Πfを出力する（順列出力過程／ステップＳ７４）。
【０１６２】
＜クラスタリング部及びクラスタリング過程の変形例２＞
第１実施形態のクラスタリング部１３２は、各信号源に対応するセントロイドcn(τ)とアクティブ系列vnf(τ)との類似度を指標とすることで、異なる周波数f∈Fの組合せに対応するアクティブ系列vnf(τ)間の類似度を間接的に指標として用い、順列Πfを生成していた。
【０１６３】
しかし、異なる周波数f∈Fの組合せに対応するアクティブ系列vnf(τ)間の類似度を直接的に指標として用いて順列Πfを生成してもよい。この場合には、例えば、すべての周波数f∈Fに対応するアクティブ系列vnf(τ)間の類似度を指標として階層的クラスタリングを行い、 そのクラスタリング結果を用いて順列Πfを生成する。以下、このような手法を例示する。
【０１６４】
図１６（ａ）（ｂ）は、それぞれ、すべての周波数f∈Fに対応するアクティブ系列vnf(τ)間の類似度を指標として階層的クラスタリングを行い、そのクラスタリング結果を用いて順列Πfを生成する大域的最適化部４３２ａを具備するクラスタリング部４３２，５３２を示したブロック図である。
【０１６５】
図１６（ａ）（ｂ）に示すように、ラスタリング部４３２，５３２の相違点は、ラスタリング部４３２がさらに局所的最適化部１３２ｂを具備し、大域的最適化部４３２ａで生成された順列Πfを局所的最適化部１３２ｂで補正して出力するものであるのに対し、ラスタリング部５３２が大域的最適化部４３２ａで生成された順列Πfをそのまま出力するものである点である。また、大域的最適化部４３２ａは、階層的クラスタリング部４３２ａａと、演算制御部４３２ａｂと、クラス選択部４３２ａｃと、順列生成部４３２ａｄとを有する。なお、局所的最適化部１３２ｂについては第１実施形態で説明済であるため、以下では、局域的最適化部４３２ａの説明のみを行う。
【０１６６】
図１７は、大域的最適化部４３２ａが行う大域的最適化過程（ステップＳ３１）を説明するためのフローチャートである。以下、この図に従って大域的最適化部４３２ａが行う大域的最適化過程を説明する。
【０１６７】
まず、階層的クラスタリング部４３２ａａに各アクティブ系列vnf(τ)が入力され、階層的クラスタリング部４３２ａａは、すべての周波数f∈Fに対応するアクティブ系列vnf(τ)間の類似度を指標として凝集型の階層的クラスタリングを１階層分行い、各アクティブ系列vnf(τ)が属するクラスの情報CLq{vnf(τ)}を生成してメモリ１７０に格納する（階層的クラスタリング過程／ステップＳ８１） 。なお、CLq{vnf(τ)}は、アクティブ系列vnf(τ)がクラスCLqに属することを意味する。また、アクティブ系列vnf(τ)間の類似度としては、アクティブ系列vnf(τ)間の相関係数（式(7)参照）や、式(18)や式(19)のセントロイドcn(τ)をアクティブ系列に置換した関数等を用いることができる。
【０１６８】
次に、演算制御部４３２ａｂがメモリ１７０から最新のクラスの情報CLq{vnf(τ)}を読み込み、クラスの総数が仮定された信号源Ｎと定数βとの和以下であるか否かを判定する（終了条件判定過程／ステップＳ８２）。なお、（クラスの総数）≦Ｎを満たすか否かではなく、（クラスの総数）≦Ｎ＋βを満たすか否かを判定する理由は、生成されたクラスがすべて信号源に対応するとは限らないからである。すなわち、クラスタリング精度によっては、何れの信号源にも対応しないクラスが誤って生成される可能性がある。そのため、凝集型の階層的クラスタリングをクラスの総数がＮ以下になるまで実行すると、信号源に対応する正しいクラスが統合され、正しいクラスの数が信号源の数N未満となる可能性があるからである。よって、ここでは、何れの信号源にも対応しないクラスの想定数以上の定数βを設定しておき、信号源に対応する正しいクラスどうしが統合される事態を防止する。
【０１６９】
ここで、（クラスの総数）≦Ｎ＋βを満たさないと判定された場合には、処理がステップＳ７１に戻される。一方、（クラスの総数）≦Ｎ＋βを満たすと判定された場合には、クラス選択部４３２ａｃがメモリ１７０から最新のクラスの情報CLq{vnf(τ)}を読み込み、メンバーの多い方から順にN個のクラスを選択し、選択したクラスの情報CL’u{vnf(τ)}(u=1,...,N)をメモリ１７０に格納する（クラス選択過程／ステップＳ８３） 。この判定は、信号源に対応する正しいクラスのメンバー数は、何れの信号源にも対応しない誤ったクラスのメンバー数よりも大きい、という仮定に基づくものである。
【０１７０】
次に、順列生成部４３２ａｄがメモリ１７０から、クラス選択部４３２ａｃで選択されたクラスの情報CL’u{vnf(τ)}を読み込み、これらを用い、vnf(τ)をvuf(τ)の値として並び替える順列Πfを周波数f毎に生成してメモリに格納する（順列出力過程／ステップＳ８４） 。
【０１７１】
〔実験結果１〕
図１９及び図２０は、音声信号を対象として第１実施形態の信号分離を行った場合のスペクトログラムとクラス分類結果とを示した図である。なお、各図の横軸は時間であり、縦軸は周波数である。また、図１９(a)に示すN=3個の源信号Sn(f,τ)(n=1,2,3)が混ざり合い、その結果、図１９(b)に示すM=2個のセンサでの観測信号Xm(f,τ)(m=1,2)が得られた場合を例示する。
【０１７２】
まず、クラス分類部１２０によって周波数f毎にこのような観測信号ベクトルX(f,τ)をN=3個のクラスに分類することにより、図２０（ａ）に示す事後確率P(Cn(f)|X(f,τ))が得られる。なお、図２０（ａ）では、各(f,τ)に対応する事後確率P(Cn(f)|X(f,τ))の大きさ（0≦P(Cn(f)|X(f,τ))≦1）を色の濃度で表現している。色が薄いほど事後確率P(Cn(f)|X(f,τ))が小さく、色が濃いほど事後確率P(Cn(f)|X(f,τ))が大きい。
【０１７３】
次に、パーミュテーション解決部１３０により、事後確率とクラスとのの対応関係を並び替えると、図２０（ｂ）に示す事後確率P’(Ck(f)|X(f,τ)) (k=1,2,3)が得られる。図２０（ｂ）より、事後確率のパーミュテーション問題が解決されていることが分かる。
【０１７４】
その後、分離部１４０により、観測信号Xm(f,τ)と事後確率P’(Ck(f)|X(f,τ))とを用いて分離信号Yk(f,τ)を生成すると、図２０（ｃ）の結果が得られる。
【０１７５】
〔実験結果２〕
次に、第１実施形態の効果を示すために、図２１（ａ）に示す実験条件と図２１（ｂ）に示す３つのマイクロホンと４つのスピーカの配置を用いて実験を行った。４つの音を同時に鳴らした時の混合音を３つのマイクロホンで観測し、その観測信号のみからそれぞれの音に対応する分離信号を算出するという問題設定である。様々な音声信号の組合せで評価できるように、スピーカからマイクロホンまでのインパルス応答を測定し、音声信号をインパルス応答に畳み込んで混合することで観測信号を生成した。分離性能は、signal-to-interference ratio(SIR)の改善量で評価した。これは、各出力i毎に、出力SIRと入力SIRの差OutputSIRi-InputSIRiとして計算される。入力SIRと出力SIRは、それぞれ以下の式で計算される。
【０１７６】
【数２３】
0004769238.tif 000024
【０１７７】
ここで、J∈｛1,...,M｝はある選択された基準センサの番号を示す。また、源信号skのうち出力yi(t)に出てきた成分をyik(t)と表記する。この定義により、yi(t)=Σk=1Nyik(t)が満たされる。
【０１７８】
４つの音声の組合せを８通り用いて実験を行った。また、スピーカからマイクロホンまでの２種類の距離（60, 120cm）と６種類の残響時間（130, 200, 270, 320, 380, 450ms）を試した。
【０１７９】
図２１（ｃ）に、４種類のパーミュテーション解決法を用いた結果を、SIR改善量の全出力に関する平均値で示す。“TDOA”と“Envelope”は,それぞれ、従来の技術に相当するものであり、信号の方向や位置に相当する値（センサ間到達時間差、TDOA: Time Difference Of Arreival）を推定することに基づくものと、分離信号エンベロープの相関係数に基づくものである。“Posterior”は、第１実施形態の事後確率の系列を用いるものである。“Optima1”と書かれたものは、信号源に関する情報を用いて最適なパーミュテーションを算出したものである。現実的な状況では、そのような情報は得られないが、性能の上限を示す目的で掲載した。
【０１８０】
図２１（ｃ）の結果から以下のことが考察できる。“TDOA”は、スピーカからマイクロホンまでの距離が短い（60cm）場合や残響時間が短い場合（130ms）には適度に良い結果を出しているが、スピーカからマイクロホンまでの距離が長く(120cm)、残響の影響が大きい場合には性能が劣化している。“Envelope”は、多くの場合、それほど良い結果にはなっていない。“Posterior”は、“Optima1”以外の現実的な方法の中で最も良い性能を達成しており、第１実施形態の効果が確認できる。
【０１８１】
〔その他の変形例等〕
なお、本発明は上述の各実施形態に限定されるものではない。例えば、上述の各種初期パラメータの設定（ステップＳ４１等）には観測信号ベクトルを用いることとしたが、初期パラメータを固定値とする構成であってもよい。
【０１８２】
また、上述の実施形態では、得られた周波数領域の分離信号を時間領域に変換することとしたが、得られた周波数領域の分離信号をそのまま出力する構成であってもよい。
【０１８３】
さらに、上述の実施形態では、時間領域と周波数領域との変換に短時間フーリエ変換を用いることとしたが、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）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【０１８９】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【産業上の利用可能性】
【０１９０】
本技術により、様々な妨害信号が発生する実環境において、目的の信号を精度良く取り出すことが可能となる。音信号に対する応用例としては、音声認識器のフロントエンドとして働く音源分離システムなどが挙げられる。話者とマイクが離れた位置にあり、マイクが話者の音声以外を収音してしまうような状況でも、そのようなシステムを使うことで、話者の音声のみを取り出して正しく音声を認識することができる。
【図面の簡単な説明】
【０１９１】
【図１】図１は、第１実施形態の信号分離装置の機能構成の全体を例示したブロック図である。
【図２】図２（ａ）は、図１に示したクラス分類部の機能構成の詳細を例示したブロック図である。また、図２（ｂ）は、図１に示したパーミュテーション問題解決部の機能構成の詳細を例示したブロック図である。
【図３】図３は、図２（ｂ）のクラスタリング部の機能構成の詳細を例示したブロック図である。
【図４】図４は、第１実施形態の信号分離装置を構成するハードウェアの構成を例示したブロック図である。
【図５】図５は、第１実施形態の信号分離方法の全体を説明するためのフローチャートである。
【図６】図６は、図５のクラス分類過程（ステップＳ２）の詳細を説明するためのフローチャートである。
【図７】図７（ａ）は、図５のパーミュテーション問題解決過程（ステップＳ３）の詳細を説明するためのフローチャートである。また、図７（ｂ）は、図７（ａ）のクラスタリング過程（ステップＳ２２）の詳細を説明するためのフローチャートである。
【図８】図８（ａ）は、図７（ｂ）の大域的最適化過程（ステップＳ３１）の詳細を説明するためのフローチャートである。また、図８（ｂ）は、図７（ｂ）の局所的最適化過程（ステップＳ３２）の詳細を説明するためのフローチャートである。
【図９】図９は、第１実施形態のクラス分類部の変形例であるクラス分類部の機能構成を示したブロック図である。
【図１０】図１０（ａ）（ｂ）は、第１実施形態のクラス分類部（図２（ａ））からノルム正規化部を排除したクラス分類部の機能構成を示したブロック図である。
【図１１】図１１は、クラス分類部の変形例を示したブロック図である。
【図１２】図１２は、クラス分類過程の変形例を説明するためのフローチャートである。
【図１３】図１３（ａ）（ｂ）（ｃ）は、クラス分類部の変形例を示したブロック図である。
【図１４】図１４（ａ）（ｂ）は、第１実施形態のクラスタリング部（図３）の変形例を示したブロック図である。
【図１５】図１５は、局所的最適化過程の変形例を説明するためのフローチャートである。
【図１６】図１６（ａ）（ｂ）は、それぞれ、すべての周波数f∈Fに対応するアクティブ系列vnf(τ)間の類似度を指標として階層的クラスタリングを行い、そのクラスタリング結果を用いて順列Πfを生成する大域的最適化部を具備するクラスタリング部を示したブロック図である。
【図１７】図１７は、大域的最適化過程（ステップＳ３１）の変形例を説明するためのフローチャートである。
【図１８】図１８（ａ）は、信号源が３つ存在する場合の２つの周波数f=766Hz,g=906Hzにおける分離信号のエンベロープ系列v1f(τ),...,v3f(τ)，v1g(τ),...,v3g(τ)を例示したグラフである。図１８（ｂ）は、図１８（ａ）と同じ観測信号に基づき生成されたアクティブ系列vnf(τ)= P(Cn(f)|X(f,τ))を示したグラフである。
【図１９】図１９は、音声信号を対象として第１実施形態の信号分離を行った場合のスペクトログラムとクラス分類結果とを示した図である。
【図２０】図２０は、音声信号を対象として第１実施形態の信号分離を行った場合のスペクトログラムとクラス分類結果とを示した図である。
【図２１】図２１（ａ）は実験条件を示した表であり、図２１（ｂ）は、実験に用いた３つのマイクロホンと４つのスピーカの配置を示した図である。また、図２１（ｃ）は、実験結果を示したグラフである。
【符号の説明】
【０１９２】
１０  信号分離装置
【図１】
0004769238.tif 000025
【図２】
0004769238.tif 000026
【図３】
0004769238.tif 000027
【図４】
0004769238.tif 000028
【図５】
0004769238.tif 000029
【図６】
0004769238.tif 000030
【図７】
0004769238.tif 000031
【図８】
0004769238.tif 000032
【図９】
0004769238.tif 000033
【図１０】
0004769238.tif 000034
【図１１】
0004769238.tif 000035
【図１２】
0004769238.tif 000036
【図１３】
0004769238.tif 000037
【図１４】
0004769238.tif 000038
【図１５】
0004769238.tif 000039
【図１６】
0004769238.tif 000040
【図１７】
0004769238.tif 000041
【図１８】
0004769238.tif 000042
【図１９】
0004769238.tif 000043
【図２０】
0004769238.tif 000044
【図２１】
0004769238.tif 000045