"1個以上の音源信号がノイズとともに混在している観測信号から、それぞれの音源に対応する分離信号を推定する音源分離装置であって、複数の収音手段で観測された前記観測信号を周波数領域の信号に変換する周波数領域変換手段と、前記観測信号と、各音源が各時間周波数において観測信号にどの程度寄与するかを表す音源占有度と、を用いて、空間パラメタを推定する空間パラメタ推定手段と、前記観測信号と前記音源占有度とを用いて、前記音源信号の分散と空間相関行列とを含むモデルパラメタを推定し、前記音源信号の分散と空間相関行列に基づくウィーナーフィルタを生成するモデルパラメタ推定手段と、前記観測信号と、前記空間パラメタまたは前記モデルパラメタと、を用いて、前記音源占有度を推定する音源占有度推定手段と、前記ウィーナーフィルタを用いて、前記観測信号からソースイメージを推定し分離信号を生成する音源分離手段と、前記分離信号を時間領域の信号に変換する時間領域変換手段と、を備える、音源分離装置。請求項1記載の音源分離装置であって、前記音源占有度を用いて、有効音源を推定する有効音源推定手段と、をさらに備え、前記音源分離手段は、前記有効音源に対応する分離信号を生成する、音源分離装置。1個以上の音源信号がノイズとともに混在している観測信号から、それぞれの音源に対応する分離信号を推定する音源分離装置であって、複数の収音手段で観測された前記観測信号を周波数領域の信号に変換する周波数領域変換手段と、前記観測信号と、各音源が各時間周波数において観測信号にどの程度寄与するかを表す音源占有度を用いて、前記音源信号の分散と空間相関行列とを含むモデルパラメタを推定し、前記音源信号の分散と空間相関行列に基づくウィーナーフィルタを生成するモデルパラメタ推定手段と、前記観測信号と前記モデルパラメタとを用いて、前記音源占有度を推定する音源占有度推定手段と、前記音源占有度を用いて、有効音源を推定する有効音源推定手段と、前記ウィーナーフィルタを用いて、前記観測信号からソースイメージを推定し前記有効音源に対応する分離信号を生成する音源分離手段と、前記分離信号を時間領域の信号に変換する時間領域変換手段と、を備える、音源分離装置。1個以上の音源信号がノイズとともに混在している観測信号から、それぞれの音源に対応する分離信号を推定する音源分離方法であって、複数の収音ステップで観測された前記観測信号を周波数領域の信号に変換する周波数領域変換ステップと、前記観測信号と、各音源が各時間周波数において観測信号にどの程度寄与するかを表す音源占有度と、を用いて、空間パラメタを推定する空間パラメタ推定ステップと、前記観測信号と、前記空間パラメタと、を用いて、前記音源占有度を推定する第一音源占有度推定ステップと、前記観測信号と前記音源占有度とを用いて、前記音源信号の分散と空間相関行列とを含むモデルパラメタを推定し、前記音源信号の分散と空間相関行列に基づくウィーナーフィルタを生成するモデルパラメタ推定ステップと、前記観測信号と、前記モデルパラメタと、を用いて、前記音源占有度を推定する第二音源占有度推定ステップと、前記ウィーナーフィルタを用いて、前記観測信号からソースイメージを推定し分離信号を生成する音源分離ステップと、前記分離信号を時間領域の信号に変換する時間領域変換ステップと、を備える、音源分離方法。請求項4記載の音源分離方法であって、前記音源占有度を用いて、有効音源を推定する有効音源推定ステップと、をさらに備え、前記音源分離ステップは、前記有効音源に対応する分離信号を生成する、音源分離方法。1個以上の音源信号がノイズとともに混在している観測信号から、それぞれの音源に対応する分離信号を推定する音源分離方法であって、複数の収音ステップで観測された前記観測信号を周波数領域の信号に変換する周波数領域変換ステップと、前記観測信号と、各音源が各時間周波数において観測信号にどの程度寄与するかを表す音源占有度を用いて、前記音源信号の分散と空間相関行列とを含むモデルパラメタを推定し、前記音源信号の分散と空間相関行列に基づくウィーナーフィルタを生成するモデルパラメタ推定ステップと、前記観測信号と前記モデルパラメタとを用いて、前記音源占有度を推定する音源占有度推定ステップと、前記音源占有度を用いて、有効音源を推定する有効音源推定ステップと、前記ウィーナーフィルタを用いて、前記観測信号からソースイメージを推定し前記有効音源に対応する分離信号を生成する音源分離ステップと、前記分離信号を時間領域の信号に変換する時間領域変換ステップと、を備える、音源分離方法。請求項1から請求項3の何れかに記載の音源分離装置として、コンピュータを機能させるためのプログラム。従来技術は、よく響く=残響の多い部屋において観測された信号も高精度高品質で分離できることが報告されているが、モデルパラメタ推定のための繰り返し計算の収束が遅い。また、従来法は、音源数Kが既知である必要があり、音源数保持部124にて音源数Kの情報を与える必要がある。本発明は、収束性の改善や音源数Kが未知の場合などを考慮した、高精度高品質で音源分離可能な音源分離技術を提供することを目的とする。上記の課題を解決するために、本発明の第一の態様に係る音源分離技術によれば、1個以上の音源信号がノイズとともに混在している観測信号から、それぞれの音源に対応する分離信号を推定する。複数の収音手段で観測された観測信号を周波数領域の信号に変換し、観測信号と、各音源が各時間周波数において観測信号にどの程度寄与するかを表す音源占有度と、を用いて、空間パラメタを推定し、観測信号と音源占有度とを用いて、音源信号の分散と空間相関行列とを含むモデルパラメタを推定し、音源信号の分散と空間相関行列に基づくウィーナーフィルタを生成し、観測信号と、空間パラメタまたはモデルパラメタと、を用いて、音源占有度を推定し、ウィーナーフィルタを用いて、観測信号からソースイメージを推定し分離信号を生成し、分離信号を時間領域の信号に変換する。上記の課題を解決するために、本発明の第二の態様に係る音源分離技術によれば、1個以上の音源信号がノイズとともに混在している観測信号から、それぞれの音源に対応する分離信号を推定する。複数の収音手段で観測された観測信号を周波数領域の信号に変換し、観測信号と音源占有度とを用いて、音源信号の分散と空間相関行列とを含むモデルパラメタを推定し、音源信号の分散と空間相関行列に基づくウィーナーフィルタを生成し、観測信号と、モデルパラメタと、を用いて、音源占有度を推定し、音源占有度を用いて、有効音源を推定し、ウィーナーフィルタを用いて、観測信号からソースイメージを推定し有効音源に対応する分離信号を生成し、分離信号を時間領域の信号に変換する。本発明は、音源信号skn,fに関してスパース性の仮定を導入し、各音源が各時間周波数n,fにおいて観測信号Xn,fにどの程度寄与するかを表す音源占有度Mkn,fを導入する音源占有度については参考文献1及び2参照。音源占有度Mkn,fに基づく音源分離手法と従来のマルチチャネルウィーナーフィルタWkn,fによる手法とのハイブリッド構成とすることで、収束性の改善や音源数が未知の場合などを考慮した、高精度高品質な音源分離を行なうことができるという効果を奏する。[参考文献1]HSawadaSArakiandSMakinoAtwostagefrequencydomainblindsourceseparationmethodforunderdeterminedconvolutivemixturesinProcWASPAA20072007pp139142[参考文献2]SArakiTNakataniandHSawadaBlindsparsesourceseparationforunknownnumberofsourcesusingGaussianmixturemodelfittingwithDirichletpriorinProcICASSP092009pp3336<本実施形態のポイント>本実施形態では、空間相関行列Bkfのクラスタリングを行なう際に、各音源信号skn,fに関してスパース性の仮定を導入する。すなわち、音源信号skn,fは稀にしか大きな値をとらず、skn,fとsk’n,f但しk≠k’とするは各時間周波数n,fで互いに重ならない、ということを仮定する。これは、互いに異なる音声信号などで確認される性質である。これを仮定すると前記式1は、Xcehse3となる。ここで、zn,fは時間周波数n,fにおいて占有的支配的な音源のインデックス番号を表し、ezn,fn,fは、この占有的な音源以外の成分占有的な音源のソースイメージczn,fn,fから見ると雑音成分のソースイメージである。さらに、観測信号Xn,fが与えられた条件下で、zn,f=k番目の音源信号が占有的であるという事後確率を導入する。このMkn,fを、音源占有度と呼ぶ。音源占有度の導入により、有効音源や有効音源数の推定が可能になったり、音源占有度Mkn,fの概念を利用する収束の早い音源分離手法と組み合わせることが可能となる。このため、音源数Kが未知であっても動作する、収束の速い音源分離手段を提供することが可能となる。以下、本発明の実施形態について、説明する。<第一実施形態に係る音源分離装置2>図2及び図3を用いて第一実施形態に係る音源分離装置2を説明する。音源分離装置2は、周波数領域変換部21と、空間パラメタ推定部221と、音源占有度推定部25と、モデルパラメタ推定部222と、パーミュテーション解決部223と、音源分離部23と、時間領域変換部24と、有効音源推定部26と、を備える。K個の音源10kそれぞれから同時に発せられる音源信号sktをM個の第m収音手段2mで但し、Kを1以上の整数とし、k=1,2,...,Kとし、Mを2以上の整数とし、m=1,2,...,Mとする、ある収録時間内例えば、5秒間に観測する。この観測状況を状況Xとする。この状況Xの下、第m収音手段2mで収音された観測信号をxmtとし、観測された観測信号をXt=[x1t,x2t,...,xMt]Tとする。収音手段とは例えばマイクロホンのことである。音源分離装置2は、第m収音手段2mを介して、K個の音源信号sktがノイズとともに混在している観測信号Xtを取得し、この観測信号から、それぞれの音源10kに対応する分離信号c^ktを推定する。なお、本実施形態では、説明を簡単にするために収音手段を2個M=2とし、第一収音手段21及び第二収音手段22としているが、2個に限定されるものではない。<周波数領域変換部21>周波数領域変換部21は、M個の収音手段2mで観測されたT個の時間領域の観測信号Xt=[x1t,x2t,...,xMt]を受け取り但し、t=1,2,...,T、これを周波数毎の時系列信号である周波数領域の観測信号Xn,f=[x1n,f,x2n,f,...,xMn,f]Tに変換しs21、空間パラメタ推定部221、モデルパラメタ推定部222及び音源占有度推定部25に出力する。以降、断りのない場合、観測信号とは、周波数毎の時系列信号である観測信号ベクトルXn,fを意味し、時間領域の観測信号の場合、それを明記する。周波数領域への変換は、例えば短時間フーリエ変換を用いれば良い。なお、本実施形態では、Nをフレームの個数とし、nをフレームに対応する時刻のインデックスとしてn=1,2,...,Nとし、Fを周波数ビンの個数とし、fを周波数ビンに対応する周波数のインデックスとしてf=1,2,...,Fとする。<空間パラメタ推定部221と音源占有度推定部25>空間パラメタ推定部221は、観測信号Xn,fを受け取り、観測信号Xn,fと、各音源10kが各時間周波数n,kにおいて観測信号Xn,fにどの程度寄与するかを表す音源占有度Mkn,fと、を用いて、空間パラメタφを推定しs221、音源占有度推定部25に出力する。音源占有度推定部25は、観測信号Xn,fと空間パラメタφとを受け取り、これらの値を用いて、音源占有度Mkn,fを推定しs251、空間パラメタ推定部221へ出力する。例えば、以下のようにして空間パラメタφ、音源占有度Mkn,fを推定する。観測信号ベクトルXn,fを空間パラメタφについてモデル化し、音源占有度Mkn,fを用いて観測信号ベクトルXn,fを分類する。なお、音源占有度Mkn,fを用いて、収束が早く、かつ、音源数Kが未知の場合にも動作する手法を用いることができることがポイントである。ここでの空間パラメタφとしては、例えば、ノルムで正規化された空間ベクトルh ̄kf=hkf/∥hkf∥例えば、参考文献1参照やマイク間位相差An,f=argx2n,f/x1n,f例えば、参考文献2参照を用いることができる。本実施形態では、ノルムで正規化された観測信号X ̄n,f=Xn,f/∥Xn,f∥を、ノルムで正規化された空間ベクトルh ̄kf=hkf/∥hkf∥を用いてとモデル化する。但し、σkfはX ̄n,fの分散を表し、本実施形態での空間パラメタは、φ={h ̄kf,σkf,αk=pzn,f=kである。なお、αkは、ある音源のインデックスkが占有的である音源のインデックスzn,fである確率を表す。また、式4において、音源数Kが既知の場合はK’=Kとし、音源数Kが未知の場合は想定される音源数より十分大きな値をK’として用いる例えばK’=10。よって、以下において、k=1,2,...,K’とする。空間パラメタφの推定は、対数尤度関数を最大とする空間パラメタφを推定することで行なう。ここではEMアルゴリズムを用いて、空間パラメタφの推定を行なう。補助関数Qφは、により与えられ、空間パラメタφは、以下のE-stepとM-stepを繰り返して計算することにより、高速に推定することができる。なお、φ’は現在の更新までで得られている空間パラメタである。空間パラメタφ及び音源占有度Mkn,fの初期化まず、空間パラメタφ及び音源占有度Mkn,fを初期化するs2211。例えば、全てのk、n、fに対して、以下のように初期値を与え、初期化する。また、i=0とする。次に、更新回数iを1ずつ増やしながらS2212、収束するまで以下の計算を行なう。なお、二回目以降の処理においては、一回前の処理において生成された空間パラメタφと音源占有度推定部25から受け取った最新の音源占有度Mkn,fを用いて以下の処理を行なう。M-step空間パラメタφの更新S2214空間パラメタ推定部221において、M-step空間パラメタφの更新を行なう。以下の行列を求め、さらに、この行列Rの最大固有値に対応する固有ベクトルを求め、その固有ベクトルを新たなh ̄kfとすることによって、h ̄kfを更新する。また、X ̄n,fの分散σkf2をとして求め、σkfを更新する。また、ある音源のインデックスkが占有的である音源のインデックスである確率αk=pzn,f=kを、ある音源の音源占有度Mkn,fの全ての時間周波数n,fに対する平均値として求める。つまり、αkをとして求め、αkを更新する。空間パラメタ推定部221は、推定し、更新した空間パラメタφを音源占有度推定部25に出力する。E-step音源占有度Mkn,fの更新S251音源占有度推定部25は、受け取った最新の空間パラメタφと、観測信号Xn,fを用いて、を求め、音源占有度Mkn,fを更新する。M-stepとE-stepを収束するまで例えば、更新回数iが20に達するまで、もしくは、Q関数の値の変化量が十分小さくなるまで繰り返すs2213。収束すると、空間パラメタ推定部221は、音源占有度推定部25から受け取った最新の音源占有度Mkn,fをモデルパラメタ推定部222に出力する。<モデルパラメタ推定部222と音源占有度推定部25>モデルパラメタ推定部222は、観測信号Xn,fと音源占有度Mkn,fとを受け取り、これらの値を用いて、音源信号skn,fの分散vkn,fと空間相関行列Bkfとを含むモデルパラメタΘ={θk}kを推定し、分散vkn,fと空間相関行列Bkfに基づくウィーナーフィルタWkn,fを生成するs222。なお、{ak}kは全てのkに関するaの集合を意味し、Θ={θk}k={θ1,θ2,...,θK’}である。音源占有度推定部25は、観測信号Xn,fとモデルパラメタΘとを受け取り、これらの値を用いて、音源占有度を推定しs252、モデルパラメタ推定部222に出力する。例えば、以下のようにしてモデルパラメタΘ、音源占有度Mkn,fを推定する。まず、ソースイメージckn,fの相関行列Rkn,f=ckn,fckn,fHが音源信号skn,fの分散vkn,fと、時刻nに依存しない空間相関行列Bkfを用いて、Rkn,f=vkn,fBkf12とモデル化できると仮定する。そして、この空間相関行列Bkfのクラスタリングを行なうことで、モデルパラメタΘを推定する。このために、本実施形態では、以下のようにソースイメージckn,fと観測信号Xn,fをモデル化する。ソースイメージckn,fを、平均0、分散vkn,fBkfの複素正規分布Ncを用いて、pBvNc0vBとモデル化する。また、式3の観測信号Xn,fが、ほぼzn,f番目のソースイメージczn,fn,fのみから成り、それ以外の音源はほぼゼロとして観測されることを表現するために、観測信号Xn,fをとモデル化する。なお、k’はzn,fを除く音源のインデックス番号とし、k’=1,2,...,zn,f-1,zn,f1,...,K’であり、δはディラックのデルタ関数を表し、Uzn,fn,fは占有的な音源のソースイメージczn,fn,fに対する雑音成分ezn,fn,fの相関行列であり、以下の式により表される。モデルパラメタ推定部222にて推定すべきモデルパラメタは、Θ={θk}k={vkn,f}k,{Bkf}kである。モデルパラメタΘの推定は、対数尤度関数を最大とするモデルパラメタΘを推定することで行なう。なお、Dは、全ての時間周波数n,f及び全ての音源のインデックスkに対する、観測信号Xn,fの集合と、占有的な音源のインデックスzn,fの集合と、ソースイメージckn,fの集合と、からなる集合を表し、∫cdcは、隠れ変数ckn,fについての周辺化を意味する。なお、式16のLn,fは、である。本実施形態ではEMアルゴリズムを用いてモデルパラメタΘの推定を行なう。補助関数は、により与えられる。なお、Θ’は現在の更新までに得られているモデルパラメタである。モデルパラメタの初期化s2221モデルパラメタ推定部222は、空間パラメタ推定部221から音源占有度Mkn,fを受け取ると、まず、分散vkn,fと空間相関行列Bkfを初期化するs2221。例えば、空間パラメタ推定部221で受け取った音源占有度Mkn,fと観測信号Xn,fの要素x1n,f第一収音手段21の観測値を用いて、として初期化する。また、更新回数i=0とする。以下、iを増やしながらS2222、E-stepとM-stepを収束するまで繰り返し、モデルパラメタΘは、この繰り返しにより更新されながら推定される。M-stepモデルパラメタΘの更新s2224モデルパラメタ推定部222において、を計算する。なお、TrAは、行列Aの対角成分の和を返す処理トレースを意味する。E-stepとM-stepとの繰り返し処理において、初めて式23を計算する場合は、式23におけるMkn,fは、空間パラメタ推定部221の出力値である音源占有度Mkn,fを用い、2回目以降は音源占有度推定部25の最新の出力値である音源占有度Mkn,fを用いる。ここで、であり、式28におけるWkn,fが、音源分離のためのマルチチャネルウィーナーフィルタWkn,fであり、ソースイメージckn,fの推定値である分離信号c^kn,fは、これを用いて式26で計算される。さらに、モデルパラメタ推定部222は、を求める。この値は、後述する式32において用いる。E-step音源占有度Mkn,fの推定s252式19のQ関数の中のp{ckn,f}k,zn,f|Xn,fの項は、と表すことができ、pzn,f=kXn,fが音源占有度Mkn,fに対応する。音源占有度推定部25は、音源占有度Mkn,fを以下の式により推定する。なお、c^kn,fとrkn,fは式26および式27で与えられる。またAB=AHB-1Aとする。以上を収束するまで例えば、更新回数iが20に達するまで、もしくは、Q関数の値の変化量が十分小さくなるまで繰り返すs2223。収束すると、モデルパラメタ推定部222は、音源占有度推定部25から受け取った最新の音源占有度Mkn,fと式28により生成した最新のウィーナーフィルタWkn,fを周波数毎に紐付けてパーミュテーション解決部223に出力する。<パーミュテーション解決部223>パーミュテーション解決部223は、音源占有度Mkn,fとこれに紐付けられたウィーナーフィルタWkn,fを受け取り、音源毎にウィーナーフィルタWkn,fをまとめるs223。空間パラメタφ及びモデルパラメタΘの推定は周波数毎に行なわれるため、各パラメタのインデックス番号kと、そのクラスタに対応する実際の音源10kとが一致しない場合がある。例えば、ある周波数fではk=1が音源101に、k=2が音源102に対応するが、別の周波数f’ではk=1が音源102に、k=2が音源101に対応する、というように、周波数毎に対応関係がばらばらになってしまうことが一般的である。これをパーミュテーションの問題という。そこで、パーミュテーション解決部223において、全ての周波数f=1,2,...Fで各パラメタのインデックスkと実際の音源10kとが完全に一対一対応するように整える。これは、例えば次のように行なわれる。まず、各周波数fおよび各インデックスkにおいて得られた各音源占有度Mkn,fを、γkf=[Mk1,f,,MkN,f]というベクトルとする。同じ音源であれば、音源占有度Mkn,fは、全ての周波数で同期する性質があることを利用し、異なる周波数間でのベクトルγkfとγkf’の相関が全ての周波数で最大となるように、インデックスの番号を入れ替える。すなわち、ベクトルaとbの相関係数をρa,bとした場合に、を最大とするkの配列Πkfを求める。ここで配列Πkfは、1,,K’の整数が適切な順序で並んだ物であり、γ ̄kは、全ての周波数におけるインデックスkに対応するγkfの平均値である。上記Jの最大化は、例えば以下の繰り返し演算により行なうことができる。これにより全ての周波数で各パラメタのインデックスkと音源10kとの対応関係を揃えることができる。パーミュテーション解決部223は、音源毎にまとめた音源占有度Mkn,fを有効音源推定部26に出力する。さらに、音源毎にまとめた音源占有度Mkn,fに紐付けられたウィーナーフィルタWkn,fを音源分離部23に出力する。<有効音源推定部26>有効音源推定部26は、音源毎にまとめられた音源占有度Mkn,fを受け取り、この音源占有度Mkn,fを用いて、有効音源を推定しs26、有効音源を音源分離部23に出力し、有効音源数Kを音源分離装置2の出力値として出力する。音源のインデックスk=1,2,...,K’のうち、有効音源に対応するインデックスの集合を{k}とし、その有効音源数Kとする。例えば、以下の方法で有効音源を推定する。本実施形態で利用する音源占有度Mkn,fは、各時間周波数n,fにおけるk番目の音源10kの占有度を表しているため、パーミュテーション問題を解決した後の音源占有度Mkn,fの平均値を、各インデックスkについて求めれば、占有度の高い音源を求めることが可能となる。よって、を計算し、pzn,f=kの値が予め設定した閾値thより大きいとき、インデックスkを占有度の高い音源10kのインデックスとして判定し、そのkの集合{k}を出力する。また有効音源数Kも出力する。<音源分離部23>音源分離部23は、ウィーナーフィルタWkn,fと有効音源の集合{k}を受け取り、有効音源の集合{k}に対応するマルチチャネルウィーナーフィルタWkn,fを生成する。さらに、音源分離部23は、観測信号Xn,fを受け取り、マルチチャネルウィーナーフィルタWkn,fを用いて、式26により観測信号Xn,fからソースイメージckn,fを推定した、有効音源に対応する分離信号c^kn,fを生成しs23、時間領域変換部24に出力する。<時間領域変換部24>時間領域変換部24は、周波数領域の分離信号c^n,fを受け取り、これを時間領域の分離信号c^tに変換しs24、この値を音源分離装置2の出力値として出力する。なお、時間領域への変換は、周波数領域変換部21で用いた変換方法に対応するものであればよい。<効果>このような構成とすることで、各パラメタの収束の早く、音源数が未知の場合にも動作する、高精度高品質な音源分離を行なうことができる。<シミュレーション結果>第一実施形態の効果を調べるため、従来技術非特許文献1参照及び第一実施形態の音源分離装置で音源分離を行なった。実験にて、マイクロホン数は2、音源数は2または3とした。サンプリング周波数は8kHz、マイクロホンの間隔は4cmである。図4は、音源数Kが既知の場合2または3に、4通りの音声の組合せについて、信号全体の歪みの尺度SDRSignaltodistortionratioを評価し、その平均を求めたものである。図4において、HB1は第一実施形態においてモデルパラメタΘの更新を1回のみにした場合の性能を、HB50は従来技術においてモデルパラメタの更新を50回にした場合の性能を示す。なお、図4及び図6において、実験時の残響時間は250msまたは400msとしている。第一実施形態は、学習回数1回のみにもかかわらず、従来技術で、学習を50回行なった場合よりも高い性能を示すことが分かる。これより、第一実施形態は、少ないモデルパラメタΘの更新で高い性能を示すことから、その収束性の早さが示された。図5は、音源数未知の場合に、K’=8として実施例を用いた場合の重み係数pzn,f=k式37参照をプロットしたものである。この結果より、有効音源推定部26において有効音源及びその数の推定が可能であることが分かる。図6は、音源数Kが未知の場合に、4通りの音声組合せについてSDRを評価し、その平均を求めた物である。図6において、HBKgiven,HBKunknownはそれぞれ、正しい音源数Kを従来技術のシステムに与えた場合と、音源数未知の条件下でK’=8として第一実施形態を適用した場合の性能を示している。第一実施形態は、音源数未知の場合でも、音源数既知の従来技術と同程度の分離性能を示すことがわかる。<他の変形例>本実施形態においては、ソースイメージ及び観測信号のモデルとしてそれぞれ13、14を用いたが、それぞれ他の適切なモデルを用いることも可能である。本実施形態においては、有効音源の推定をモデルパラメタΘの推定の後に行なったが、これを空間パラメタφの推定後に行ない、モデルパラメタΘの推定は、推定された音源数Kの音源に対してのみ行なってもよい。この場合、空間パラメタφの推定後に、音源占有度Mkn,fとそれに紐付けられる空間パラメタφをパーミュテーション解決部223の入力とし、パーミュテーション解決部223は、音源毎に音源占有度Mkn,fと空間パラメタφをまとめ、有効音源推定部26に出力する。有効音源推定部26は、上述の方法により、有効音源を推定し、有効音源に対応するインデックスの集合を{k}と有効音源数Kと、それに対応する音源占有度Mkn,fと空間パラメタφをモデルパラメタ推定部222に出力する。本実施形態においては、有効音源推定部26において、有効音源数Kを推定しているが、予め有効音源数Kが利用者等により与えられている場合には、有効音源推定部26を設けなくともよい。その場合には、各部は、図示しない音源数保持部から記憶されている有効音源数Kを取得する。このような構成の場合には、各パラメタの収束の早く、かつ、高精度高品質な音源分離を行なうことができる。なお、空間パラメタφ、モデルパラメタΘ及び音源占有度Mkn,fの初期値は上述した値以外の値であってもよい。例えば、各パラメタが取りうる値をランダムに設定してもよい。<第二実施形態>図7及び図8を用いて第二実施形態に係る音源分離装置3を説明する。第一実施形態と異なる部分についてのみ説明する。音源分離装置3は、周波数領域変換部21と、音源占有度推定部35と、モデルパラメタ推定部322と、パーミュテーション解決部223と、音源分離部23と、時間領域変換部24と、有効音源推定部26と、を備える。空間パラメタ推定部221を備えていない点、及び、音源占有度推定部35とモデルパラメタ推定部322における処理内容が第一実施形態とは異なる。<モデルパラメタ推定部322と音源占有度推定部35>モデルパラメタ推定部322は、観測信号Xn,fを受け取り、観測信号Xkn,fと音源占有度Mkn,fを用いて、音源信号の分散vkn,fと空間相関行列Bkfとを含むモデルパラメタΘを推定し、音源信号の分散vkn,fと空間相関行列Bkfに基づくウィーナーフィルタWkn,fを生成するs322。音源占有度推定部35は、観測信号Xn,fとモデルパラメタΘとを受け取り、これらの値を用いて、音源占有度Mkn,fを推定しs352、モデルパラメタ推定部222に出力する。第一実施形態とは異なり、空間パラメタ推定部がないため、音源占有度推定部35は、観測信号Xn,fと空間パラメタφとを用いて、音源占有度Mkn,fを推定する必要がなく、音源占有度Mkn,fの初期値をモデルパラメタ推定部322において与えるs3221。例えば、全てのk、n、fに対して、Mkn,f=1として初期値を与える。他の処理については第一実施形態と同様である。<効果>第二実施形態の音源分離装置3は、収束の高速化の効果はなくなるが、音源数が未知の場合にも動作し、高精度高品質な音源分離を行なうことができる。<第三実施形態>図9及び図10を用いて第三実施形態に係る音源分離装置4を説明する。音源分離装置2と異なる部分についてのみ説明する。音源分離装置4は、周波数領域変換部21と、空間パラメタ推定部421と、音源占有度推定部25と、モデルパラメタ推定部422と、パーミュテーション解決部223と、音源分離部23と、時間領域変換部24と、有効音源推定部26と、を備える。第一実施形態では、空間パラメタφの推定を十分収束するまで行なってから、モデルパラメタΘの推定を行なったが、本実施形態では、空間パラメタφとモデルパラメタΘの推定を、それぞれ1回ずつ更新しながら、全体としての最適化を行なう点が異なる。<空間パラメタ推定部421>空間パラメタ推定部421は、観測信号Xn,fを受け取り、観測信号Xn,fと、各音源信号skn,fに対応するソースイメージckn,fが各時間周波数n,kにおいて観測信号Xn,fにどの程度寄与するかを表す音源占有度Mkn,fと、を用いて、空間パラメタφを推定しs421、音源占有度推定部45に出力する。第一実施形態のs2211において説明した空間パラメタφ及び音源占有度Mkn,fの初期化とs2221において説明したモデルパラメタΘの初期化を繰り返し処理に先立ち行なう。以下、iを増やしながらS2222、空間パラメタ推定部421と音源占有度推定部25とモデルパラメタ推定部422における処理を繰り返し、空間パラメタφ及びモデルパラメタΘは、この繰り返しにより更新されながら推定される。空間パラメタ推定部421において、空間パラメタφの推定s2214は第一実施形態と同様の処理により行なわれ、空間パラメタφは音源占有度推定部25に出力される。音源占有度推定部25は、第一実施形態と同様の処理により、音源占有度Mkn,fを推定しs251、空間パラメタ推定部421へ出力する。空間パラメタ推定部421は、受け取った音源占有度Mkn,fをモデルパラメタ推定部422に出力する。<モデルパラメタ推定部422>モデルパラメタ推定部422は、第一実施形態と同様の処理により、モデルパラメタΘを推定し、これに基づくウィーナーフィルタWkn,fを生成するs422。なお、本実施形態では、空間パラメタ推定部421と音源占有度推定部25とモデルパラメタ推定部422における処理をひとまとめとして処理として繰り返すため、モデルパラメタ推定部422においてiを増やす必要はない。また、モデルパラメタ推定部422は、モデルパラメタΘを受け取るたびに収束しているか否かを判定しs4223、収束すると、モデルパラメタ推定部222は、音源占有度推定部25から受け取った最新の音源占有度Mkn,fと式28により生成した最新のウィーナーフィルタWkn,fを周波数毎に紐付けてパーミュテーション解決部223に出力する。収束していない場合には、第一実施形態と同様の処理により、モデルパラメタΘを推定し、音源占有度推定部25に出力する。音源占有度推定部25は、第一実施形態と同様の処理により音源占有度Mkn,fを推定しs252、モデルパラメタ推定部422に出力する。さらに、モデルパラメタ推定部422は受け取った音源占有度Mkn,fを空間パラメタ推定部421に出力する。<効果>このような構成とすることで第一実施形態と同様の効果を得ることができる。<プログラム及び記録媒体>上述した音源分離装置は、コンピュータにより機能させることもできる。この場合はコンピュータに、目的とする装置各種実施形態で図に示した機能構成をもつ装置として機能させるためのプログラム、またはその処理手順各実施形態で示したものの各過程をコンピュータに実行させるためのプログラムを、CD-ROM、磁気ディスク、半導体記憶装置などの記録媒体から、あるいは通信回線を介してそのコンピュータ内にダウンロードし、そのプログラムを実行させればよい。<その他の変形例>本発明は上記の実施形態及び変形例に限定されるものではない。例えば、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。また、各パラメタは、図示しない記憶部等に格納され、各部はこの記憶部から各パラメタを取得する構成としてもよい。"