(57)【特許請求の範囲】
【請求項１】
１個以上の音源信号がノイズとともに混在している観測信号から、それぞれの音源に対応する分離信号を推定する音源分離装置であって、
複数の収音手段で観測された前記観測信号を周波数領域の信号に変換する周波数領域変換手段と、
前記観測信号と、各音源が各時間周波数において観測信号にどの程度寄与するかを表す音源占有度と、を用いて、空間パラメタを推定する空間パラメタ推定手段と、
前記観測信号と前記音源占有度とを用いて、前記音源信号の分散と空間相関行列とを含むモデルパラメタを推定し、前記音源信号の分散と空間相関行列に基づくウィーナーフィルタを生成するモデルパラメタ推定手段と、
前記観測信号と、前記空間パラメタまたは前記モデルパラメタと、を用いて、前記音源占有度を推定する音源占有度推定手段と、
前記ウィーナーフィルタを用いて、前記観測信号からソースイメージを推定し分離信号を生成する音源分離手段と、
前記分離信号を時間領域の信号に変換する時間領域変換手段と、を備える、
音源分離装置。
【請求項２】
請求項１記載の音源分離装置であって、
前記音源占有度を用いて、有効音源を推定する有効音源推定手段と、をさらに備え、
前記音源分離手段は、前記有効音源に対応する分離信号を生成する、
音源分離装置。
【請求項３】
１個以上の音源信号がノイズとともに混在している観測信号から、それぞれの音源に対応する分離信号を推定する音源分離装置であって、
複数の収音手段で観測された前記観測信号を周波数領域の信号に変換する周波数領域変換手段と、
前記観測信号と、各音源が各時間周波数において観測信号にどの程度寄与するかを表す音源占有度を用いて、前記音源信号の分散と空間相関行列とを含むモデルパラメタを推定し、前記音源信号の分散と空間相関行列に基づくウィーナーフィルタを生成するモデルパラメタ推定手段と、
前記観測信号と前記モデルパラメタとを用いて、前記音源占有度を推定する音源占有度推定手段と、
前記音源占有度を用いて、有効音源を推定する有効音源推定手段と、
前記ウィーナーフィルタを用いて、前記観測信号からソースイメージを推定し前記有効音源に対応する分離信号を生成する音源分離手段と、
前記分離信号を時間領域の信号に変換する時間領域変換手段と、を備える、
音源分離装置。
【請求項４】
１個以上の音源信号がノイズとともに混在している観測信号から、それぞれの音源に対応する分離信号を推定する音源分離方法であって、
複数の収音ステップで観測された前記観測信号を周波数領域の信号に変換する周波数領域変換ステップと、
前記観測信号と、各音源が各時間周波数において観測信号にどの程度寄与するかを表す音源占有度と、を用いて、空間パラメタを推定する空間パラメタ推定ステップと、
前記観測信号と、前記空間パラメタと、を用いて、前記音源占有度を推定する第一音源占有度推定ステップと、
前記観測信号と前記音源占有度とを用いて、前記音源信号の分散と空間相関行列とを含むモデルパラメタを推定し、前記音源信号の分散と空間相関行列に基づくウィーナーフィルタを生成するモデルパラメタ推定ステップと、
前記観測信号と、前記モデルパラメタと、を用いて、前記音源占有度を推定する第二音源占有度推定ステップと、
前記ウィーナーフィルタを用いて、前記観測信号からソースイメージを推定し分離信号を生成する音源分離ステップと、
前記分離信号を時間領域の信号に変換する時間領域変換ステップと、を備える、
音源分離方法。
【請求項５】
請求項４記載の音源分離方法であって、
前記音源占有度を用いて、有効音源を推定する有効音源推定ステップと、をさらに備え、
前記音源分離ステップは、前記有効音源に対応する分離信号を生成する、
音源分離方法。
【請求項６】
１個以上の音源信号がノイズとともに混在している観測信号から、それぞれの音源に対応する分離信号を推定する音源分離方法であって、
複数の収音ステップで観測された前記観測信号を周波数領域の信号に変換する周波数領域変換ステップと、
前記観測信号と、各音源が各時間周波数において観測信号にどの程度寄与するかを表す音源占有度を用いて、前記音源信号の分散と空間相関行列とを含むモデルパラメタを推定し、前記音源信号の分散と空間相関行列に基づくウィーナーフィルタを生成するモデルパラメタ推定ステップと、
前記観測信号と前記モデルパラメタとを用いて、前記音源占有度を推定する音源占有度推定ステップと、
前記音源占有度を用いて、有効音源を推定する有効音源推定ステップと、
前記ウィーナーフィルタを用いて、前記観測信号からソースイメージを推定し前記有効音源に対応する分離信号を生成する音源分離ステップと、
前記分離信号を時間領域の信号に変換する時間領域変換ステップと、を備える、
音源分離方法。
【請求項７】
請求項１から請求項３の何れかに記載の音源分離装置として、コンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は信号処理の技術分野に関し、特に１個以上の音源信号がノイズとともに混在している観測信号から、それぞれの音源に対応する分離信号を推定する音源分離技術に関する。
【背景技術】
【０００２】
マルチチャネルウィーナーフィルタを用いた音源分離装置が従来技術として知られている（図１参照）。この従来技術の詳細は、例えば、非特許文献１に記載されている。以下、概略を説明する。Ｋ個の音源１０（ｋ）（但し、ｋは音源のインデックス番号であり、ｋ＝１，…，Ｋ）それぞれから同時に発せられる音源信号ｓ（ｋ）（ｔ）を２個の収音手段（第一収音手段２１及び第二収音手段２２）で、ある収録時間内（例えば、５秒間）に観測する。この観測状況を状況Ｘとする。この状況Ｘの下、第一収音手段２１で収音された観測信号をｘ１（ｔ）とし、第二収音手段２２で収音された観測信号をｘ２（ｔ）とし、観測された観測信号をＸ（ｔ）＝［ｘ１（ｔ），ｘ２（ｔ）］Ｔとする。収音手段とは例えばマイクロホンのことであり、「Ｔ」は行列の転置を表し、ｔを離散時刻とし、ｔ＝１，…，Ｔとする。周波数領域変換部１１は、観測信号Ｘ（ｔ）を周波数領域に変換することで、周波数毎の時系列信号である観測信号ベクトルＸ（ｎ，ｆ）＝［ｘ１（ｎ，ｆ），ｘ２（ｎ，ｆ）］Ｔに変換する。以降、断りのない場合、観測信号とは、周波数毎の時系列信号である観測信号ベクトルＸ（ｎ，ｆ）＝［ｘ１（ｎ，ｆ），ｘ２（ｎ，ｆ）］Ｔを意味し、時間領域の観測信号の場合、それを明記する。周波数領域への変換は、例えば短時間フーリエ変換を用いれば良い。なお、フレーム数をＮとし、ｎはフレームに対応する時刻を表し、ｎ＝１，２，…，Ｎとする。また、ｆはサンプリング周波数ｆｓをＦ等分した離散点（周波数ビン）を表す。例えば、ｆ∈｛０，（１／Ｆ）ｆｓ，…，（（Ｆ−１）／Ｆ）ｆｓ｝である。
【０００３】
ここで、観測信号ベクトルは、以下の式（１）で表されると仮定する。
【０００４】
【数１】
0005406866.tif 000002
【０００５】
但し、ｃ（ｋ）（ｎ，ｆ）＝［ｃ（１，ｋ）（ｎ，ｆ），ｃ（２，ｋ）（ｎ，ｆ）］Ｔはｋ番目の音源信号ｓ（ｋ）（ｎ，ｆ）のソースイメージを表し、各成分ｃ（１，ｋ）（ｎ，ｆ），ｃ（２，ｋ）（ｎ，ｆ）は、それぞれ音源信号ｓ（ｋ）（ｎ，ｆ）の第一収音手段２１、第二収音手段２２における観測値を表す。つまり、ソースイメージｃ（ｋ）（ｎ，ｆ）は、各収音手段における観測値ｘ１（ｎ，ｆ），ｘ２（ｎ，ｆ）中の音源信号ｓ（ｋ）（ｎ，ｆ）に基づく信号である。また、ｈ（ｋ）＝［ｈ（１，ｋ），ｈ（２，ｋ）］Ｔであり、各成分ｈ（１，ｋ），ｈ（２，ｋ）は、それぞれ音源１０（ｋ）から第一収音手段２１、第二収音手段２２までの周波数応答を示し、ｓ（ｋ）（ｎ，ｆ）は音源信号ｓ（ｋ）（ｔ）を周波数領域で表現した信号であり、以降、断りのない場合、音源信号とは、周波数毎の時系列信号であるｓ（ｋ）（ｎ，ｆ）を意味し、時間領域の音源信号の場合、それを明記する。
【０００６】
音源信号の観測時間内においては、音源１０（ｋ）、第一収音手段２１及び第二収音手段２２は固定されており、また、Ｋ個の音源１０（ｋ）は全て、異なる位置に配置されているとする。すなわち、周波数応答ｈ（ｋ）は時刻ｎに依存せずに、ｋの値によって異なる値をとるものと仮定する。
【０００７】
従来技術では、ソースイメージｃ（ｋ）（ｎ，ｆ）の相関行列
【０００８】
【数２】
0005406866.tif 000003
【０００９】
が音源信号ｓ（ｋ）（ｎ，ｆ）の時間周波数毎の分散ｖ（ｋ）（ｎ，ｆ）と、時刻ｎによらない周波数毎のＭ行Ｍ列の空間相関行列Ｂ（ｋ）（ｆ）を用いて、
Ｒ（ｋ）（ｎ，ｆ）＝ｖ（ｋ）（ｎ，ｆ）Ｂ（ｋ）（ｆ）        （２）
Ｂ（ｋ）（ｆ）＝ｈ（ｋ）（ｆ）（ｈ（ｋ）（ｆ））Ｈ
とモデル化できると仮定する。但し、「Ｈ」はエルミート転置を表す。そして、空間相関行列Ｂ（ｋ）（ｆ）のクラスタリングを行なうことで、分散ｖ（ｋ）（ｎ，ｆ）と空間相関行列Ｂ（ｋ）（ｆ）を推定する（詳細は非特許文献１参照）。さらに、推定した分散ｖ（ｋ）（ｎ，ｆ）と空間相関行列Ｂ（ｋ）（ｆ）に基づくマルチチャネルウィーナーフィルタＷ（ｋ）（ｎ，ｆ）を用いて、音源分離を達成する。すなわち、空間相関行列Ｂ（ｋ）が、音源１０（ｋ）の位置毎に固有の値を持つことを利用して音源分離を行なう。具体的には、図１において、まずモデルパラメタ初期化部１２１において、分散ｖ（ｋ）（ｎ，ｆ）と空間相関行列Ｂ（ｋ）（ｆ）に適当な初期値を与える。
【００１０】
次に、モデルパラメタ推定部１２２において、分散ｖ（ｋ）（ｎ，ｆ）と空間相関行列Ｂ（ｋ）（ｆ）を周波数毎に推定する（詳細は非特許文献１参照）。この処理は周波数毎に独立に行なわれるため、各パラメタのインデックス（ｋ）と、実際の音源１０（ｋ）とが一致しない場合がある。そこで、パーミュテーション解決部１２３において、各パラメタのインデックスｋと実際の音源１０（ｋ）との対応を全ての周波数において一致させる。
【００１１】
音源分離部１３では、得られたパラメタを用いて、ウィーナーフィルタ
【００１２】
【数３】
0005406866.tif 000004
【００１３】
を計算する。但し、
【００１４】
【数４】
0005406866.tif 000005
【００１５】
とし、Ｒ（ｋ）（ｎ，ｆ）は式（２）により与えられる。そして音源分離部１３はさらにソースイメージｃ（ｋ）（ｎ，ｆ）の推定値である分離信号
ｃ＾（ｋ）（ｎ，ｆ）＝Ｗ（ｋ）（ｎ，ｆ）Ｘ（ｎ，ｆ）
を計算し出力する。最後に時間領域変換部１４にて周波数領域の分離信号ｃ＾（ｋ）（ｎ，ｆ）を時間領域の分離信号ｃ＾（ｋ）（ｔ）に変換する。
【先行技術文献】
【非特許文献】
【００１６】
【非特許文献１】N. Q. K. Duong, E. Vincent, and R. Gribonval, “Under-determined reverberant audio source separation using a full-rank spatial covariance model”, IEEE Transactions on Audio, Speech and Language Processing, 2010, vol. 18, no. 7, pp. 1830-1840.
【発明の概要】
【発明が解決しようとする課題】
【００１７】
従来技術は、よく響く（＝残響の多い）部屋において観測された信号も高精度・高品質で分離できることが報告されているが、モデルパラメタ推定のための繰り返し計算の収束が遅い。また、従来法は、音源数Ｋが既知である必要があり、音源数保持部１２４にて音源数Ｋの情報を与える必要がある。
【００１８】
本発明は、収束性の改善や音源数Ｋが未知の場合などを考慮した、高精度・高品質で音源分離可能な音源分離技術を提供することを目的とする。
【課題を解決するための手段】
【００１９】
上記の課題を解決するために、本発明の第一の態様に係る音源分離技術によれば、１個以上の音源信号がノイズとともに混在している観測信号から、それぞれの音源に対応する分離信号を推定する。複数の収音手段で観測された観測信号を周波数領域の信号に変換し、観測信号と、各音源が各時間周波数において観測信号にどの程度寄与するかを表す音源占有度と、を用いて、空間パラメタを推定し、観測信号と音源占有度とを用いて、音源信号の分散と空間相関行列とを含むモデルパラメタを推定し、音源信号の分散と空間相関行列に基づくウィーナーフィルタを生成し、観測信号と、空間パラメタまたはモデルパラメタと、を用いて、音源占有度を推定し、ウィーナーフィルタを用いて、観測信号からソースイメージを推定し分離信号を生成し、分離信号を時間領域の信号に変換する。
【００２０】
上記の課題を解決するために、本発明の第二の態様に係る音源分離技術によれば、１個以上の音源信号がノイズとともに混在している観測信号から、それぞれの音源に対応する分離信号を推定する。複数の収音手段で観測された観測信号を周波数領域の信号に変換し、観測信号と音源占有度とを用いて、音源信号の分散と空間相関行列とを含むモデルパラメタを推定し、音源信号の分散と空間相関行列に基づくウィーナーフィルタを生成し、観測信号と、モデルパラメタと、を用いて、音源占有度を推定し、音源占有度を用いて、有効音源を推定し、ウィーナーフィルタを用いて、観測信号からソースイメージを推定し有効音源に対応する分離信号を生成し、分離信号を時間領域の信号に変換する。
【発明の効果】
【００２１】
本発明は、音源信号ｓ（ｋ）（ｎ，ｆ）に関してスパース性の仮定を導入し、各音源が各時間周波数（ｎ，ｆ）において観測信号Ｘ（ｎ，ｆ）にどの程度寄与するかを表す音源占有度Ｍ（ｋ）（ｎ，ｆ）を導入する（音源占有度については参考文献１及び２参照）。音源占有度Ｍ（ｋ）（ｎ，ｆ）に基づく音源分離手法と従来のマルチチャネルウィーナーフィルタＷ（ｋ）（ｎ，ｆ）による手法とのハイブリッド構成とすることで、収束性の改善や音源数が未知の場合などを考慮した、高精度・高品質な音源分離を行なうことができるという効果を奏する。
【００２２】
［参考文献１］H. Sawada, S.  Araki, and S. Makino, “A two-stage frequency-domain blind source separation method for underdetermined convolutive mixtures”, in Proc. WASPAA2007, 2007, pp.139-142
［参考文献２］S. Araki, T. Nakatani, and H. Sawada, “Blind sparse source separation for unknown number of sources using Gaussian mixture model fitting with Dirichlet prior”, in Proc. ICASSP'09, 2009, pp.33-36
【図面の簡単な説明】
【００２３】
【図１】音源分離装置１の機能ブロック図。
【図２】音源分離装置２の構成例を示す機能ブロック図。
【図３】音源分離装置２の処理フローを示す図。
【図４】音源分離装置１及び２の収束性を示す図。
【図５】音源分離装置２の有効音源推定結果を示す図。
【図６】音源数が既知の音源分離装置１及び音源数が未知の音源分離装置２の音源分離性能を示す図。
【図７】音源分離装置３の構成例を示す機能ブロック図。
【図８】音源分離装置３の処理フローを示す図。
【図９】音源分離装置４の構成例を示す機能ブロック図。
【図１０】音源分離装置４の処理フローを示す図。
【発明を実施するための形態】
【００２４】
＜本実施形態のポイント＞
本実施形態では、空間相関行列Ｂ（ｋ）（ｆ）のクラスタリングを行なう際に、各音源信号ｓ（ｋ）（ｎ，ｆ）に関してスパース性の仮定を導入する。すなわち、音源信号ｓ（ｋ）（ｎ，ｆ）は稀にしか大きな値をとらず、ｓ（ｋ）（ｎ，ｆ）とｓ（ｋ’）（ｎ，ｆ）（但しｋ≠ｋ’とする）は各時間周波数（ｎ，ｆ）で互いに重ならない、ということを仮定する。これは、互いに異なる音声信号などで確認される性質である。これを仮定すると前記式（１）は、
X(n,f)=c(z(n,f))(n,f)+e(z(n,f))(n,f)=h(z(n,f))(f)s(z(n,f))(n,f)+e(z(n,f))(n,f)    （３）
となる。ここで、ｚ（ｎ，ｆ）は時間周波数（ｎ，ｆ）において占有的（支配的）な音源のインデックス番号を表し、ｅ（ｚ（ｎ，ｆ））（ｎ，ｆ）は、この占有的な音源以外の成分（占有的な音源のソースイメージｃ（ｚ（ｎ，ｆ））（ｎ，ｆ）から見ると雑音成分のソースイメージ）である。
【００２５】
さらに、観測信号Ｘ（ｎ，ｆ）が与えられた条件下で、ｚ（ｎ，ｆ）＝ｋ番目の音源信号が占有的であるという事後確率
【００２６】
【数５】
0005406866.tif 000006
【００２７】
を導入する。このＭ（ｋ）（ｎ，ｆ）を、音源占有度と呼ぶ。音源占有度の導入により、有効音源や有効音源数の推定が可能になったり、音源占有度Ｍ（ｋ）（ｎ，ｆ）の概念を利用する収束の早い音源分離手法と組み合わせることが可能となる。このため、音源数Ｋが未知であっても動作する、収束の速い音源分離手段を提供することが可能となる。
【００２８】
以下、本発明の実施形態について、説明する。
【００２９】
＜第一実施形態に係る音源分離装置２＞
図２及び図３を用いて第一実施形態に係る音源分離装置２を説明する。音源分離装置２は、周波数領域変換部２１と、空間パラメタ推定部２２１と、音源占有度推定部２５と、モデルパラメタ推定部２２２と、パーミュテーション解決部２２３と、音源分離部２３と、時間領域変換部２４と、有効音源推定部２６と、を備える。
【００３０】
Ｋ個の音源１０（ｋ）それぞれから同時に発せられる音源信号ｓ（ｋ）（ｔ）をＭ個の第ｍ収音手段２ｍで（但し、Ｋを１以上の整数とし、ｋ＝１，２，…，Ｋとし、Ｍを２以上の整数とし、ｍ＝１，２，…，Ｍとする）、ある収録時間内（例えば、５秒間）に観測する。この観測状況を状況Ｘとする。この状況Ｘの下、第ｍ収音手段２ｍで収音された観測信号をｘｍ（ｔ）とし、観測された観測信号をＸ（ｔ）＝［ｘ１（ｔ），ｘ２（ｔ），…，ｘＭ（ｔ）］Ｔとする。収音手段とは例えばマイクロホンのことである。音源分離装置２は、第ｍ収音手段２ｍを介して、Ｋ個の音源信号ｓ（ｋ）（ｔ）がノイズとともに混在している観測信号Ｘ（ｔ）を取得し、この観測信号から、それぞれの音源１０（ｋ）に対応する分離信号ｃ＾（ｋ）（ｔ）を推定する。なお、本実施形態では、説明を簡単にするために収音手段を２個（Ｍ＝２）とし、第一収音手段２１及び第二収音手段２２としているが、２個に限定されるものではない。
【００３１】
＜周波数領域変換部２１＞
周波数領域変換部２１は、Ｍ個の収音手段２ｍで観測されたＴ個の時間領域の観測信号Ｘ（ｔ）＝［ｘ１（ｔ），ｘ２（ｔ），…，ｘＭ（ｔ）］を受け取り（但し、ｔ＝１，２，…，Ｔ）、これを周波数毎の時系列信号である周波数領域の観測信号Ｘ（ｎ，ｆ）＝［ｘ１（ｎ，ｆ），ｘ２（ｎ，ｆ），…，ｘＭ（ｎ，ｆ）］Ｔに変換し（ｓ２１）、空間パラメタ推定部２２１、モデルパラメタ推定部２２２及び音源占有度推定部２５に出力する。以降、断りのない場合、観測信号とは、周波数毎の時系列信号である観測信号ベクトルＸ（ｎ，ｆ）を意味し、時間領域の観測信号の場合、それを明記する。周波数領域への変換は、例えば短時間フーリエ変換を用いれば良い。なお、本実施形態では、Ｎをフレームの個数とし、ｎをフレームに対応する時刻のインデックスとしてｎ＝１，２，…，Ｎとし、Ｆを周波数ビンの個数とし、ｆを周波数ビンに対応する周波数のインデックスとしてｆ＝１，２，…，Ｆとする。
【００３２】
＜空間パラメタ推定部２２１と音源占有度推定部２５＞
空間パラメタ推定部２２１は、観測信号Ｘ（ｎ，ｆ）を受け取り、観測信号Ｘ（ｎ，ｆ）と、各音源１０（ｋ）が各時間周波数（ｎ，ｋ）において観測信号Ｘ（ｎ，ｆ）にどの程度寄与するかを表す音源占有度Ｍ（ｋ）（ｎ，ｆ）と、を用いて、空間パラメタφを推定し（ｓ２２１）、音源占有度推定部２５に出力する。
【００３３】
音源占有度推定部２５は、観測信号Ｘ（ｎ，ｆ）と空間パラメタφとを受け取り、これらの値を用いて、音源占有度Ｍ（ｋ）（ｎ，ｆ）を推定し（ｓ２５１）、空間パラメタ推定部２２１へ出力する。例えば、以下のようにして空間パラメタφ、音源占有度Ｍ（ｋ）（ｎ，ｆ）を推定する。
【００３４】
観測信号ベクトルＸ（ｎ，ｆ）を空間パラメタφについてモデル化し、音源占有度Ｍ（ｋ）（ｎ，ｆ）を用いて観測信号ベクトルＸ（ｎ，ｆ）を分類する。なお、音源占有度Ｍ（ｋ）（ｎ，ｆ）を用いて、収束が早く、かつ、音源数Ｋが未知の場合にも動作する手法を用いることができることがポイントである。ここでの空間パラメタφとしては、例えば、ノルムで正規化された空間ベクトルｈ￣（ｋ）（ｆ）＝ｈ（ｋ）（ｆ）／‖ｈ（ｋ）（ｆ）‖（例えば、参考文献１参照）やマイク間位相差Ａ（ｎ，ｆ）＝ａｒｇ（ｘ２（ｎ，ｆ）／ｘ１（ｎ，ｆ））（例えば、参考文献２参照）を用いることができる。
【００３５】
本実施形態では、ノルムで正規化された観測信号Ｘ￣（ｎ，ｆ）＝Ｘ（ｎ，ｆ）／‖Ｘ（ｎ，ｆ）‖を、ノルムで正規化された空間ベクトルｈ￣（ｋ）（ｆ）＝ｈ（ｋ）（ｆ）／‖ｈ（ｋ）（ｆ）‖を用いて
【００３６】
【数６】
0005406866.tif 000007
【００３７】
とモデル化する。但し、σ（ｋ）（ｆ）はＸ￣（ｎ，ｆ）の分散を表し、本実施形態での空間パラメタは、φ＝｛ｈ￣（ｋ）（ｆ），σ（ｋ）（ｆ），α（ｋ）＝ｐ（ｚ（ｎ，ｆ）＝ｋ）である。なお、α（ｋ）は、ある音源のインデックスｋが占有的である音源のインデックスｚ（ｎ，ｆ）である確率を表す。また、式（４）において、音源数Ｋが既知の場合はＫ’＝Ｋとし、音源数Ｋが未知の場合は想定される音源数より十分大きな値をＫ’として用いる（例えばＫ’＝１０）。よって、以下において、ｋ＝１，２，…，Ｋ’とする。
【００３８】
空間パラメタφの推定は、対数尤度関数
【００３９】
【数７】
0005406866.tif 000008
【００４０】
を最大とする空間パラメタφを推定することで行なう。ここではＥＭアルゴリズムを用いて、空間パラメタφの推定を行なう。補助関数Ｑ（φ）は、
【００４１】
【数８】
0005406866.tif 000009
【００４２】
により与えられ、空間パラメタφは、以下のＥ−ｓｔｅｐとＭ−ｓｔｅｐを繰り返して計算することにより、高速に推定することができる。なお、φ’は現在の更新までで得られている空間パラメタである。
（空間パラメタφ及び音源占有度Ｍ（ｋ）（ｎ，ｆ）の初期化）
まず、空間パラメタφ及び音源占有度Ｍ（ｋ）（ｎ，ｆ）を初期化する（ｓ２２１１）。例えば、全てのｋ、ｎ、ｆに対して、以下のように初期値を与え、初期化する。
【００４３】
【数９】
0005406866.tif 000010
【００４４】
また、ｉ＝０とする。次に、更新回数ｉを１ずつ増やしながら（Ｓ２２１２）、収束するまで以下の計算を行なう。なお、二回目以降の処理においては、一回前の処理において生成された空間パラメタφと音源占有度推定部２５から受け取った最新の音源占有度Ｍ（ｋ）（ｎ，ｆ）を用いて以下の処理を行なう。
（Ｍ−ｓｔｅｐ（空間パラメタφの更新）（Ｓ２２１４））
空間パラメタ推定部２２１において、Ｍ−ｓｔｅｐ（空間パラメタφの更新）を行なう。以下の行列
【００４５】
【数１０】
0005406866.tif 000011
【００４６】
を求め、さらに、この行列Ｒの最大固有値に対応する固有ベクトルを求め、その固有ベクトルを新たなｈ￣（ｋ）（ｆ）とすることによって、ｈ￣（ｋ）（ｆ）を更新する。
【００４７】
また、Ｘ￣（ｎ，ｆ）の分散（σ（ｋ）（ｆ））２を
【００４８】
【数１１】
0005406866.tif 000012
【００４９】
として求め、σ（ｋ）（ｆ）を更新する。
【００５０】
また、ある音源のインデックスｋが占有的である音源のインデックスである確率α（ｋ）＝ｐ（ｚ（ｎ，ｆ）＝ｋ）を、ある音源の音源占有度Ｍ（ｋ）（ｎ，ｆ）の全ての時間周波数（ｎ，ｆ）に対する平均値として求める。つまり、α（ｋ）を
【００５１】
【数１２】
0005406866.tif 000013
【００５２】
として求め、α（ｋ）を更新する。
【００５３】
空間パラメタ推定部２２１は、推定し、更新した空間パラメタφを音源占有度推定部２５に出力する。
（Ｅ−ｓｔｅｐ（音源占有度Ｍ（ｋ）（ｎ，ｆ）の更新）（Ｓ２５１））
音源占有度推定部２５は、受け取った最新の空間パラメタφと、観測信号Ｘ（ｎ，ｆ）を用いて、
【００５４】
【数１３】
0005406866.tif 000014
【００５５】
を求め、音源占有度Ｍ（ｋ）（ｎ，ｆ）を更新する。
【００５６】
Ｍ−ｓｔｅｐとＥ−ｓｔｅｐを収束するまで（例えば、更新回数ｉが２０に達するまで、もしくは、Ｑ関数の値の変化量が十分小さくなるまで）繰り返す（ｓ２２１３）。
【００５７】
収束すると、空間パラメタ推定部２２１は、音源占有度推定部２５から受け取った最新の音源占有度Ｍ（ｋ）（ｎ，ｆ）をモデルパラメタ推定部２２２に出力する。
【００５８】
＜モデルパラメタ推定部２２２と音源占有度推定部２５＞
モデルパラメタ推定部２２２は、観測信号Ｘ（ｎ，ｆ）と音源占有度Ｍ（ｋ）（ｎ，ｆ）とを受け取り、これらの値を用いて、音源信号ｓ（ｋ）（ｎ，ｆ）の分散ｖ（ｋ）（ｎ，ｆ）と空間相関行列Ｂ（ｋ）（ｆ）とを含むモデルパラメタΘ＝｛θ（ｋ）｝ｋを推定し、分散ｖ（ｋ）（ｎ，ｆ）と空間相関行列Ｂ（ｋ）（ｆ）に基づくウィーナーフィルタＷ（ｋ）（ｎ，ｆ）を生成する（ｓ２２２）。なお、｛ａ（ｋ）｝ｋは全てのｋに関するａの集合を意味し、Θ＝｛θ（ｋ）｝ｋ＝｛θ（１），θ（２），…，θ（Ｋ’）｝である。
【００５９】
音源占有度推定部２５は、観測信号Ｘ（ｎ，ｆ）とモデルパラメタΘとを受け取り、これらの値を用いて、音源占有度を推定し（ｓ２５２）、モデルパラメタ推定部２２２に出力する。例えば、以下のようにしてモデルパラメタΘ、音源占有度Ｍ（ｋ）（ｎ，ｆ）を推定する。
【００６０】
まず、ソースイメージｃ（ｋ）（ｎ，ｆ）の相関行列Ｒ（ｋ）（ｎ，ｆ）＝ｃ（ｋ）（ｎ，ｆ）（ｃ（ｋ）（ｎ，ｆ））Ｈが音源信号ｓ（ｋ）（ｎ，ｆ）の分散ｖ（ｋ）（ｎ，ｆ）と、時刻ｎに依存しない空間相関行列Ｂ（ｋ）（ｆ）を用いて、
Ｒ（ｋ）（ｎ，ｆ）＝ｖ（ｋ）（ｎ，ｆ）Ｂ（ｋ）（ｆ）        （１２）
とモデル化できると仮定する。そして、この空間相関行列Ｂ（ｋ）（ｆ）のクラスタリングを行なうことで、モデルパラメタΘを推定する。このために、本実施形態では、以下のようにソースイメージｃ（ｋ）（ｎ，ｆ）と観測信号Ｘ（ｎ，ｆ）をモデル化する。
【００６１】
ソースイメージｃ（ｋ）（ｎ，ｆ）を、平均０、分散ｖ（ｋ）（ｎ，ｆ）Ｂ（ｋ）（ｆ）の複素正規分布Ｎｃを用いて、
p(c(k)(n,f);B(k)(f),v(k)(n,f))=Nc(c(k)(n,f);0,v(k)(n,f)B(k)(f))    (13)
とモデル化する。
【００６２】
また、式（３）の観測信号Ｘ（ｎ，ｆ）が、ほぼｚ（ｎ，ｆ）番目のソースイメージｃ（ｚ（ｎ，ｆ）））（ｎ，ｆ）のみから成り、それ以外の音源はほぼゼロとして観測されることを表現するために、観測信号Ｘ（ｎ，ｆ）を
【００６３】
【数１４】
0005406866.tif 000015
【００６４】
とモデル化する。なお、ｋ’はｚ（ｎ，ｆ）を除く音源のインデックス番号とし、ｋ’＝１，２，…，（ｚ（ｎ，ｆ）−１），（ｚ（ｎ，ｆ）＋１），…，Ｋ’であり、δはディラックのデルタ関数を表し、Ｕ（ｚ（ｎ，ｆ））（ｎ，ｆ）は占有的な音源のソースイメージｃ（ｚ（ｎ，ｆ））（ｎ，ｆ）に対する雑音成分ｅ（ｚ（ｎ，ｆ））（ｎ，ｆ）の相関行列であり、以下の式により表される。
【００６５】
【数１５】
0005406866.tif 000016
【００６６】
モデルパラメタ推定部２２２にて推定すべきモデルパラメタは、Θ＝｛θ（ｋ）｝ｋ＝（｛ｖ（ｋ）（ｎ，ｆ）｝ｋ，｛Ｂ（ｋ）（ｆ）｝ｋ）である。モデルパラメタΘの推定は、対数尤度関数
【００６７】
【数１６】
0005406866.tif 000017
【００６８】
を最大とするモデルパラメタΘを推定することで行なう。なお、Ｄは、全ての時間周波数（ｎ，ｆ）及び全ての音源のインデックスｋに対する、観測信号Ｘ（ｎ，ｆ）の集合と、占有的な音源のインデックスｚ（ｎ，ｆ）の集合と、ソースイメージｃ（ｋ）（ｎ，ｆ）の集合と、からなる集合を表し、∫ｃｄｃは、隠れ変数ｃ（ｋ）（ｎ，ｆ）についての周辺化を意味する。なお、式（１６）のＬ（ｎ，ｆ）は、
【００６９】
【数１７】
0005406866.tif 000018
【００７０】
である。本実施形態ではＥＭアルゴリズムを用いてモデルパラメタΘの推定を行なう。補助関数は、
【００７１】
【数１８】
0005406866.tif 000019
【００７２】
により与えられる。なお、Θ’は現在の更新までに得られているモデルパラメタである。
（モデルパラメタの初期化（ｓ２２２１））
モデルパラメタ推定部２２２は、空間パラメタ推定部２２１から音源占有度Ｍ（ｋ）（ｎ，ｆ）を受け取ると、まず、分散ｖ（ｋ）（ｎ，ｆ）と空間相関行列Ｂ（ｋ）（ｆ）を初期化する（ｓ２２２１）。例えば、空間パラメタ推定部２２１で受け取った音源占有度Ｍ（ｋ）（ｎ，ｆ）と観測信号Ｘ（ｎ，ｆ）の要素ｘ１（ｎ，ｆ）（第一収音手段２１の観測値）を用いて、
【００７３】
【数１９】
0005406866.tif 000020
【００７４】
として初期化する。
【００７５】
また、更新回数ｉ＝０とする。以下、ｉを増やしながら（Ｓ２２２２）、Ｅ−ｓｔｅｐとＭ−ｓｔｅｐを収束するまで繰り返し、モデルパラメタΘは、この繰り返しにより更新されながら推定される。
（Ｍ−ｓｔｅｐ（モデルパラメタΘの更新）（ｓ２２２４））
モデルパラメタ推定部２２２において、
【００７６】
【数２０】
0005406866.tif 000021
【００７７】
を計算する。なお、Ｔｒ（Ａ）は、行列Ａの対角成分の和を返す処理（トレース）を意味する。Ｅ−ｓｔｅｐとＭ−ｓｔｅｐとの繰り返し処理において、初めて式（２３）を計算する場合は、式（２３）におけるＭ（ｋ）（ｎ，ｆ）は、空間パラメタ推定部２２１の出力値である音源占有度Ｍ（ｋ）（ｎ，ｆ）を用い、２回目以降は音源占有度推定部２５の最新の出力値である音源占有度Ｍ（ｋ）（ｎ，ｆ）を用いる。ここで、
【００７８】
【数２１】
0005406866.tif 000022
【００７９】
であり、式（２８）におけるＷ（ｋ）（ｎ，ｆ）が、音源分離のためのマルチチャネルウィーナーフィルタＷ（ｋ）（ｎ，ｆ）であり、ソースイメージｃ（ｋ）（ｎ，ｆ）の推定値である分離信号ｃ＾（ｋ）（ｎ，ｆ）は、これを用いて式（２６）で計算される。
【００８０】
さらに、モデルパラメタ推定部２２２は、
【００８１】
【数２２】
0005406866.tif 000023
【００８２】
を求める。この値は、後述する式（３２）において用いる。
（Ｅ−ｓｔｅｐ（音源占有度Ｍ（ｋ）（ｎ，ｆ）の推定（ｓ２５２））
式（１９）のＱ関数の中のｐ（｛ｃ（ｋ）（ｎ，ｆ）｝ｋ，ｚ（ｎ，ｆ）｜Ｘ（ｎ，ｆ））の項は、
【００８３】
【数２３】
0005406866.tif 000024
【００８４】
と表すことができ、ｐ（ｚ（ｎ，ｆ）＝ｋ|Ｘ（ｎ，ｆ））が音源占有度Ｍ（ｋ）（ｎ，ｆ）に対応する。
【００８５】
音源占有度推定部２５は、音源占有度Ｍ（ｋ）（ｎ，ｆ）を以下の式により推定する。
【００８６】
【数２４】
0005406866.tif 000025
【００８７】
なお、ｃ＾（ｋ）（ｎ，ｆ）とｒ（ｋ）（ｎ，ｆ）は式（２６）および式（２７）で与えられる。また||Ａ||Ｂ＝ＡＨＢ−１Ａとする。
【００８８】
以上を収束するまで（例えば、更新回数ｉが２０に達するまで、もしくは、Ｑ関数の値の変化量が十分小さくなるまで）繰り返す（ｓ２２２３）。
【００８９】
収束すると、モデルパラメタ推定部２２２は、音源占有度推定部２５から受け取った最新の音源占有度Ｍ（ｋ）（ｎ，ｆ）と式（２８）により生成した最新のウィーナーフィルタＷ（ｋ）（ｎ，ｆ）を周波数毎に紐付けてパーミュテーション解決部２２３に出力する。
【００９０】
＜パーミュテーション解決部２２３＞
パーミュテーション解決部２２３は、音源占有度Ｍ（ｋ）（ｎ，ｆ）とこれに紐付けられたウィーナーフィルタＷ（ｋ）（ｎ，ｆ）を受け取り、音源毎にウィーナーフィルタＷ（ｋ）（ｎ，ｆ）をまとめる（ｓ２２３）。
【００９１】
空間パラメタφ及びモデルパラメタΘの推定は周波数毎に行なわれるため、各パラメタのインデックス番号ｋと、そのクラスタに対応する実際の音源１０（ｋ）とが一致しない場合がある。例えば、ある周波数ｆではｋ＝１が音源１０（１）に、ｋ＝２が音源１０（２）に対応するが、別の周波数ｆ’ではｋ＝１が音源１０（２）に、ｋ＝２が音源１０（１）に対応する、というように、周波数毎に対応関係がばらばらになってしまうことが一般的である。これをパーミュテーションの問題という。そこで、パーミュテーション解決部２２３において、全ての周波数ｆ＝１，２，…Ｆで各パラメタのインデックス（ｋ）と実際の音源１０（ｋ）とが完全に一対一対応するように整える。これは、例えば次のように行なわれる。
【００９２】
まず、各周波数ｆおよび各インデックスｋにおいて得られた各音源占有度Ｍ（ｋ）ｎ，ｆを、
γ（ｋ）（ｆ）＝［Ｍ（ｋ）（１，ｆ），・・・，Ｍ（ｋ）（Ｎ，ｆ）］
というベクトルとする。同じ音源であれば、音源占有度Ｍ（ｋ）（ｎ，ｆ）は、全ての周波数で同期する性質があることを利用し、異なる周波数間でのベクトルγ（ｋ）（ｆ）とγ（ｋ）（ｆ’）の相関が全ての周波数で最大となるように、インデックスの番号を入れ替える。すなわち、ベクトルａとｂの相関係数をρ（ａ，ｂ）とした場合に、
【００９３】
【数２５】
0005406866.tif 000026
【００９４】
を最大とするｋの配列Π（ｋ）（ｆ）を求める。ここで配列Π（ｋ）（ｆ）は、１，・・・，Ｋ’の整数が適切な順序で並んだ物であり、γ￣（ｋ）は、全ての周波数におけるインデックスｋに対応するγ（ｋ）（ｆ）の平均値である。上記Ｊの最大化は、例えば以下の繰り返し演算により行なうことができる。
【００９５】
【数２６】
0005406866.tif 000027
【００９６】
これにより全ての周波数で各パラメタのインデックスｋと音源１０（ｋ）との対応関係を揃えることができる。
【００９７】
パーミュテーション解決部２２３は、音源毎にまとめた音源占有度Ｍ（ｋ）（ｎ，ｆ）を有効音源推定部２６に出力する。さらに、音源毎にまとめた音源占有度Ｍ（ｋ）（ｎ，ｆ）に紐付けられたウィーナーフィルタＷ（ｋ）（ｎ，ｆ）を音源分離部２３に出力する。
【００９８】
＜有効音源推定部２６＞
有効音源推定部２６は、音源毎にまとめられた音源占有度Ｍ（ｋ）（ｎ，ｆ）を受け取り、この音源占有度Ｍ（ｋ）（ｎ，ｆ）を用いて、有効音源を推定し（ｓ２６）、有効音源を音源分離部２３に出力し、有効音源数Ｋを音源分離装置２の出力値として出力する。音源のインデックスｋ＝１，２，…，Ｋ’のうち、有効音源に対応するインデックスの集合を｛ｋ｝＊とし、その有効音源数Ｋとする。例えば、以下の方法で有効音源を推定する。本実施形態で利用する音源占有度Ｍ（ｋ）（ｎ，ｆ）は、各時間周波数（ｎ，ｆ）におけるｋ番目の音源１０（ｋ）の占有度を表しているため、パーミュテーション問題を解決した後の音源占有度Ｍ（ｋ）（ｎ，ｆ）の平均値を、各インデックスｋについて求めれば、占有度の高い音源を求めることが可能となる。よって、
【００９９】
【数２７】
0005406866.tif 000028
【０１００】
を計算し、ｐ（ｚ（ｎ，ｆ）＝ｋ）の値が予め設定した閾値ｔｈより大きいとき、インデックスｋを占有度の高い音源１０（ｋ）のインデックスとして判定し、そのｋの集合｛ｋ｝＊を出力する。また有効音源数Ｋも出力する。
【０１０１】
＜音源分離部２３＞
音源分離部２３は、ウィーナーフィルタＷ（ｋ）（ｎ，ｆ）と有効音源の集合｛ｋ｝＊を受け取り、有効音源の集合｛ｋ｝＊に対応するマルチチャネルウィーナーフィルタＷ（ｋ）（ｎ，ｆ）を生成する。さらに、音源分離部２３は、観測信号Ｘ（ｎ，ｆ）を受け取り、マルチチャネルウィーナーフィルタＷ（ｋ）（ｎ，ｆ）を用いて、式（２６）により観測信号Ｘ（ｎ，ｆ）からソースイメージｃ（ｋ）（ｎ，ｆ）を推定した、有効音源に対応する分離信号ｃ＾（ｋ）（ｎ，ｆ）を生成し（ｓ２３）、時間領域変換部２４に出力する。
【０１０２】
＜時間領域変換部２４＞
時間領域変換部２４は、周波数領域の分離信号ｃ＾（ｎ，ｆ）を受け取り、これを時間領域の分離信号ｃ＾（ｔ）に変換し（ｓ２４）、この値を音源分離装置２の出力値として出力する。なお、時間領域への変換は、周波数領域変換部２１で用いた変換方法に対応するものであればよい。
【０１０３】
＜効果＞
このような構成とすることで、各パラメタの収束の早く、音源数が未知の場合にも動作する、高精度・高品質な音源分離を行なうことができる。
【０１０４】
＜シミュレーション結果＞
第一実施形態の効果を調べるため、従来技術（非特許文献１参照）及び第一実施形態の音源分離装置で音源分離を行なった。実験にて、マイクロホン数は２、音源数は２または３とした。サンプリング周波数は８ｋＨｚ、マイクロホンの間隔は４ｃｍである。
【０１０５】
図４は、音源数Ｋが既知の場合（２または３）に、４通りの音声の組合せについて、信号全体の歪みの尺度ＳＤＲ（Signal to distortion ratio）を評価し、その平均を求めたものである。図４において、ＨＢ１は第一実施形態においてモデルパラメタΘの更新を１回のみにした場合の性能を、ＨＢ５０は従来技術においてモデルパラメタの更新を５０回にした場合の性能を示す。なお、図４及び図６において、実験時の残響時間は２５０ｍｓまたは４００ｍｓとしている。第一実施形態は、学習回数１回のみにもかかわらず、従来技術で、学習を５０回行なった場合よりも高い性能を示すことが分かる。これより、第一実施形態は、少ないモデルパラメタΘの更新で高い性能を示すことから、その収束性の早さが示された。
【０１０６】
図５は、音源数未知の場合に、Ｋ’＝８として実施例を用いた場合の重み係数ｐ（ｚ（ｎ，ｆ）＝ｋ）（式（３７）参照）をプロットしたものである。この結果より、有効音源推定部２６において有効音源及びその数の推定が可能であることが分かる。
【０１０７】
図６は、音源数Ｋが未知の場合に、４通りの音声組合せについてＳＤＲを評価し、その平均を求めた物である。図６において、ＨＢ（Ｋ  ｇｉｖｅｎ），ＨＢ（Ｋ  ｕｎｋｎｏｗｎ）はそれぞれ、正しい音源数Ｋを従来技術のシステムに与えた場合と、音源数未知の条件下でＫ’＝８として第一実施形態を適用した場合の性能を示している。第一実施形態は、音源数未知の場合でも、音源数既知の従来技術と同程度の分離性能を示すことがわかる。
【０１０８】
＜他の変形例＞
本実施形態においては、ソースイメージ及び観測信号のモデルとしてそれぞれ（１３）、（１４）を用いたが、それぞれ他の適切なモデルを用いることも可能である。
【０１０９】
本実施形態においては、有効音源の推定をモデルパラメタΘの推定の後に行なったが、これを空間パラメタφの推定後に行ない、モデルパラメタΘの推定は、推定された音源数Ｋの音源に対してのみ行なってもよい。この場合、空間パラメタφの推定後に、音源占有度Ｍ（ｋ）（ｎ，ｆ）とそれに紐付けられる空間パラメタφをパーミュテーション解決部２２３の入力とし、パーミュテーション解決部２２３は、音源毎に音源占有度Ｍ（ｋ）（ｎ，ｆ）（と空間パラメタφ）をまとめ、有効音源推定部２６に出力する。有効音源推定部２６は、上述の方法により、有効音源を推定し、有効音源に対応するインデックスの集合を｛ｋ｝＊と有効音源数Ｋと、それに対応する音源占有度Ｍ（ｋ）（ｎ，ｆ）と空間パラメタφをモデルパラメタ推定部２２２に出力する。
【０１１０】
本実施形態においては、有効音源推定部２６において、有効音源数Ｋを推定しているが、予め有効音源数Ｋが利用者等により与えられている場合には、有効音源推定部２６を設けなくともよい。その場合には、各部は、図示しない音源数保持部から記憶されている有効音源数Ｋを取得する。このような構成の場合には、各パラメタの収束の早く、かつ、高精度・高品質な音源分離を行なうことができる。
【０１１１】
なお、空間パラメタφ、モデルパラメタΘ及び音源占有度Ｍ（ｋ）（ｎ，ｆ）の初期値は上述した値以外の値であってもよい。例えば、各パラメタが取りうる値をランダムに設定してもよい。
【０１１２】
＜第二実施形態＞
図７及び図８を用いて第二実施形態に係る音源分離装置３を説明する。第一実施形態と異なる部分についてのみ説明する。音源分離装置３は、周波数領域変換部２１と、音源占有度推定部３５と、モデルパラメタ推定部３２２と、パーミュテーション解決部２２３と、音源分離部２３と、時間領域変換部２４と、有効音源推定部２６と、を備える。空間パラメタ推定部２２１を備えていない点、及び、音源占有度推定部３５とモデルパラメタ推定部３２２における処理内容が第一実施形態とは異なる。
【０１１３】
＜モデルパラメタ推定部３２２と音源占有度推定部３５＞
モデルパラメタ推定部３２２は、観測信号Ｘ（ｎ，ｆ）を受け取り、観測信号Ｘ（ｋ）（ｎ，ｆ）と音源占有度Ｍ（ｋ）（ｎ，ｆ）を用いて、音源信号の分散ｖ（ｋ）（ｎ，ｆ）と空間相関行列Ｂ（ｋ）（ｆ）とを含むモデルパラメタΘを推定し、音源信号の分散ｖ（ｋ）（ｎ，ｆ）と空間相関行列Ｂ（ｋ）（ｆ）に基づくウィーナーフィルタＷ（ｋ）（ｎ，ｆ）を生成する（ｓ３２２）。
【０１１４】
音源占有度推定部３５は、観測信号Ｘ（ｎ，ｆ）とモデルパラメタΘとを受け取り、これらの値を用いて、音源占有度Ｍ（ｋ）（ｎ，ｆ）を推定し（ｓ３５２）、モデルパラメタ推定部２２２に出力する。
【０１１５】
第一実施形態とは異なり、空間パラメタ推定部がないため、音源占有度推定部３５は、観測信号Ｘ（ｎ，ｆ）と空間パラメタφとを用いて、音源占有度Ｍ（ｋ）（ｎ，ｆ）を推定する必要がなく、音源占有度Ｍ（ｋ）（ｎ，ｆ）の初期値をモデルパラメタ推定部３２２において与える（ｓ３２２１）。例えば、全てのｋ、ｎ、ｆに対して、Ｍ（ｋ）（ｎ，ｆ）＝１として初期値を与える。
【０１１６】
他の処理については第一実施形態と同様である。
【０１１７】
＜効果＞
第二実施形態の音源分離装置３は、収束の高速化の効果はなくなるが、音源数が未知の場合にも動作し、高精度・高品質な音源分離を行なうことができる。
【０１１８】
＜第三実施形態＞
図９及び図１０を用いて第三実施形態に係る音源分離装置４を説明する。音源分離装置２と異なる部分についてのみ説明する。音源分離装置４は、周波数領域変換部２１と、空間パラメタ推定部４２１と、音源占有度推定部２５と、モデルパラメタ推定部４２２と、パーミュテーション解決部２２３と、音源分離部２３と、時間領域変換部２４と、有効音源推定部２６と、を備える。
【０１１９】
第一実施形態では、空間パラメタφの推定を十分収束するまで行なってから、モデルパラメタΘの推定を行なったが、本実施形態では、空間パラメタφとモデルパラメタΘの推定を、それぞれ１回ずつ更新しながら、全体としての最適化を行なう点が異なる。
【０１２０】
＜空間パラメタ推定部４２１＞
空間パラメタ推定部４２１は、観測信号Ｘ（ｎ，ｆ）を受け取り、観測信号Ｘ（ｎ，ｆ）と、各音源信号ｓ（ｋ）（ｎ，ｆ）に対応するソースイメージｃ（ｋ）（ｎ，ｆ）が各時間周波数（ｎ，ｋ）において観測信号Ｘ（ｎ，ｆ）にどの程度寄与するかを表す音源占有度Ｍ（ｋ）（ｎ，ｆ）と、を用いて、空間パラメタφを推定し（ｓ４２１）、音源占有度推定部４５に出力する。
【０１２１】
第一実施形態のｓ２２１１において説明した空間パラメタφ及び音源占有度Ｍ（ｋ）（ｎ，ｆ）の初期化とｓ２２２１において説明したモデルパラメタΘの初期化を繰り返し処理に先立ち行なう。
【０１２２】
以下、ｉを増やしながら（Ｓ２２２２）、空間パラメタ推定部４２１と音源占有度推定部２５とモデルパラメタ推定部４２２における処理を繰り返し、空間パラメタφ及びモデルパラメタΘは、この繰り返しにより更新されながら推定される。
【０１２３】
空間パラメタ推定部４２１において、空間パラメタφの推定（ｓ２２１４）は第一実施形態と同様の処理により行なわれ、空間パラメタφは音源占有度推定部２５に出力される。
【０１２４】
音源占有度推定部２５は、第一実施形態と同様の処理により、音源占有度Ｍ（ｋ）（ｎ，ｆ）を推定し（ｓ２５１）、空間パラメタ推定部４２１へ出力する。空間パラメタ推定部４２１は、受け取った音源占有度Ｍ（ｋ）（ｎ，ｆ）をモデルパラメタ推定部４２２に出力する。
【０１２５】
＜モデルパラメタ推定部４２２＞
モデルパラメタ推定部４２２は、第一実施形態と同様の処理により、モデルパラメタΘを推定し、これに基づくウィーナーフィルタＷ（ｋ）（ｎ，ｆ）を生成する（ｓ４２２）。なお、本実施形態では、空間パラメタ推定部４２１と音源占有度推定部２５とモデルパラメタ推定部４２２における処理をひとまとめとして処理として繰り返すため、モデルパラメタ推定部４２２においてｉを増やす必要はない。また、モデルパラメタ推定部４２２は、モデルパラメタΘを受け取るたびに収束しているか否かを判定し（ｓ４２２３）、収束すると、モデルパラメタ推定部２２２は、音源占有度推定部２５から受け取った最新の音源占有度Ｍ（ｋ）（ｎ，ｆ）と式（２８）により生成した最新のウィーナーフィルタＷ（ｋ）（ｎ，ｆ）を周波数毎に紐付けてパーミュテーション解決部２２３に出力する。収束していない場合には、第一実施形態と同様の処理により、モデルパラメタΘを推定し、音源占有度推定部２５に出力する。
【０１２６】
音源占有度推定部２５は、第一実施形態と同様の処理により音源占有度Ｍ（ｋ）（ｎ，ｆ）を推定し（ｓ２５２）、モデルパラメタ推定部４２２に出力する。さらに、モデルパラメタ推定部４２２は受け取った音源占有度Ｍ（ｋ）（ｎ，ｆ）を空間パラメタ推定部４２１に出力する。
【０１２７】
＜効果＞
このような構成とすることで第一実施形態と同様の効果を得ることができる。
【０１２８】
＜プログラム及び記録媒体＞
上述した音源分離装置は、コンピュータにより機能させることもできる。この場合はコンピュータに、目的とする装置（各種実施形態で図に示した機能構成をもつ装置）として機能させるためのプログラム、またはその処理手順（各実施形態で示したもの）の各過程をコンピュータに実行させるためのプログラムを、ＣＤ−ＲＯＭ、磁気ディスク、半導体記憶装置などの記録媒体から、あるいは通信回線を介してそのコンピュータ内にダウンロードし、そのプログラムを実行させればよい。
【０１２９】
＜その他の変形例＞
本発明は上記の実施形態及び変形例に限定されるものではない。例えば、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。また、各パラメタは、図示しない記憶部等に格納され、各部はこの記憶部から各パラメタを取得する構成としてもよい。
【符号の説明】
【０１３０】
２，３，４  音源分離装置
２ｍ  収音手段
１０  音源
２１  周波数領域変換部
２３  音源分離部
２４  時間領域変換部
２５，３５，４５  音源占有度推定部
２６  有効音源推定部
２２１，４２１  空間パラメタ推定部
２２２，３２２，４２２  モデルパラメタ推定部
２２３  パーミュテーション解決部
【図１】
0005406866.tif 000029
【図２】
0005406866.tif 000030
【図３】
0005406866.tif 000031
【図４】
0005406866.tif 000032
【図５】
0005406866.tif 000033
【図６】
0005406866.tif 000034
【図７】
0005406866.tif 000035
【図８】
0005406866.tif 000036
【図９】
0005406866.tif 000037
【図１０】
0005406866.tif 000038