JP2013186383A - 音源分離装置、音源分離方法、およびプログラム - Google Patents

Abstract

【課題】空間的エイリアジングの問題が生じる場合にも、高い性能で音源分離を行う。
【解決手段】２個のマイクロホンで観測された観測信号の複素スペクトルを得、観測信号の複素スペクトルのマイクロホン間位相差を各時間周波数について計算する。１個のマイクロホンで観測された観測信号から得られた観測信号の複素スペクトル、特徴量作成部で得られたマイクロホン間位相差、および音源信号の複素スペクトルの分布を表す事前情報を用い、各時間周波数での各音源に対応する観測信号の複素スペクトルおよびマイクロホン間位相差の分布を表す確率モデルのパラメタを推定する。推定されたパラメタから得られる各時間周波数での各音源の観測信号の複素スペクトルおよびマイクロホン間位相差に対する寄与率を用い、各時間周波数での当該寄与率および観測信号の複素スペクトルから、分離信号の複素スペクトルを抽出し、それを時間領域の分離信号に変換する。
【選択図】図４

Description

本発明は信号処理の技術分野に属する。特に複数信号が混在している音響データから、それぞれの音源信号を推定する、音源分離技術に属する。その中でも特に、音源信号やそれらがどう混ざったかの情報を用いずに複数信号が混在している音響データのみから、それぞれの音源信号を推定する、ブラインド音源分離技術に属する。

音源がＮ_ｓ個（Ｎ_ｓは２以上の整数）あり、それらから同時に発せられた音源信号が混合し、２個のマイクロホンで観測され、観測信号ｘ（ｔ）＝［ｘ_１（ｔ），ｘ_２（ｔ）］^Ｔが得られたとする。ｔは離散的な時刻を表し、「α^Ｔ」はαの転置を表す。このようなＮ_ｓ個の音源信号が混ざった観測信号ｘ（ｔ）＝［ｘ_１（ｔ），ｘ_２（ｔ）］^Ｔから、それぞれの音源信号を推定する。このような技術としては、たとえば特許文献１に開示されたものがある。以下では特許文献１に開示された技術をもとに従来技術を説明する。

従来技術の構成を図１に例示する。
まず、周波数領域変換部１０１において、たとえば短時間フーリエ変換等により、時間領域の観測信号ｘ_Ｌ（ｔ）を所定の時間区間である時間フレームごとに周波数領域表現し、観測信号の複素スペクトルＸ_{ｎ，ｆ，Ｌ}を生成する。複素スペクトルＸ_{ｎ，ｆ，Ｌ}が以下のように表されると仮定する。

ここで、ｎは各時間フレームをそれぞれ代表する離散的な時刻（たとえば、各時間フレームの中心時刻）に対応する整数のインデックス、ｆは離散的な各周波数（０，ｆ_ｓ・（１／Ｎ_ｆ），...，ｆ_ｓ・（Ｎ_ｆ−１）／Ｎ_ｆ）に対応する整数のインデックス、Ｎ_ｆは周波数ｆの個数（Ｎ_ｆは１以上の整数、例えばＮ_ｆ≧２）、ｆ_ｓはサンプリング周波数、Ｌは各マイクロホンにそれぞれ対応するインデックス（Ｌ＝１，２）、ｍは各音源にそれぞれ対応するインデックス（ｍ＝１，...，Ｎ_ｓ）を表す。Ｈ_{ｆ，Ｌ，ｍ}は音源ｍからマイクロホンＬまでの周波数ｆでの周波数応答を表す。Ｓ_{ｎ，ｆ，ｍ}は、音源ｍから発せられた音源信号の、時刻ｎおよび周波数ｆでの複素スペクトルを表す。「音源ｍ」はインデックスｍに対応する音源、「時刻ｎ」はインデックスｎに対応する時刻、「周波数ｆ」はインデックスｆに対応する周波数、「マイクロホンＬ」はインデックスＬに対応するマイクロホンを意味する。時刻ｎおよび周波数ｆからなる組（時間周波数スロット（ｎ，ｆ））を「時間周波数（ｎ，ｆ）」と表現する。「・」は乗算演算子を表す。

音源分離を行うために、音源信号がスパースである、すなわち、各音源信号の複素スペクトルＳ_{ｎ，ｆ，ｍ}は稀にしか大きな値をとらず、Ｓ_{ｎ，ｆ，ｍ１}とＳ_{ｎ，ｆ，ｍ２}（ｍ１≠ｍ２、ｍ１，ｍ２∈｛１，...，Ｎ_ｓ｝）が各時間周波数（ｎ，ｆ）で互いに重ならない、ということを仮定する。これは、互いに異なる音源から発せられた音声信号などで確認される性質である。これを仮定すると、式（１）は以下のように近似できる。

ここでＳ_{ｎ，ｆ，ｍ}は、時間周波数（ｎ，ｆ）で支配的な音源信号を表す。

また、式（１）におけるＨ_{ｆ，Ｌ，ｍ}が、音源信号がマイクロホンに到達する間に部屋の残響の影響をうけないことを仮定する、以下の無響モデルで書けると仮定する。

ただし、ｅは自然対数の底、ｊは虚数単位、μ_ｍは音源ｍから発せられた音源信号がマイクロホン１，２に到達するまでの時間差である。

次に、位相差計算部１０２において、マイクロホン１と２における、観測信号の複素スペクトルＸ_{ｎ，ｆ，１}，Ｘ_{ｎ，ｆ，２}の位相差Ａ_ｎ，ｆ（以下「マイクロホン間位相差」と呼ぶ）を以下のように計算する。

ここで、式（２）と（３）より、音源ｍが支配的な時間周波数（ｎ，ｆ）では以下のようになることに注意。

ただし、ａｒｇ［α］はαの偏角を表す。

従来法では、次に位相差分類部１０３において、マイクロホン間位相差Ａ_ｎ，ｆをクラスタリングする。スパース性の仮定（式（２））より、音源ｍ’＝ｍ１だけが支配的な時間周波数（ｎ，ｆ）での時間差μ_ｍ１に対応するマイクロホン間位相差Ａ_ｎ，ｆ、および音源ｍ’＝ｍ２だけが支配的な時間周波数（ｎ，ｆ）での時間差μ_ｍ２に対応するマイクロホン間位相差Ａ_ｎ，ｆが求まっているので、マイクロホンからみた音源ｍ１の方向と音源ｍ２の方向とが異なっていれば、全時間周波数のマイクロホン間位相差Ａ_ｎ，ｆを分類できる。分類の方法としては、Ａ_ｎ，ｆに対してｋ−ｍｅａｎｓ法などを用いてクラスタリングを行なったり、Ａ_ｎ，ｆに対して混合正規分布をあてはめたりすることが一般的である。前者ではそれぞれのクラスタ、後者ではそれぞれの正規分布により、各音源信号の複素スペクトルを推定することができる。

よって次に音源分離部１０４において、位相差分類部１０３で得られた情報ξを用い、音源ｍに対応するクラスタを形成している時間周波数（ｎ，ｆ）では１に近い値を、それ以外の時間周波数では０に近い値をとるマスクＭ_{ｎ，ｆ，ｍ}を作る。Ｍ_{ｎ，ｆ，ｍ}としては、あるｍ番目のクラスタを形成している時間周波数（ｎ，ｆ）でのみ１をとり、他の時間周波数（ｎ，ｆ）では０をとるように設計することもできるし、上述の混合正規分布から計算される、各時間周波数（ｎ，ｆ）での各音源の寄与率を用いることもできる。このＭ_{ｎ，ｆ，ｍ}を全ての音源ｍについて作る。さらに、音源分離部１０４で、マスクＭ_{ｎ，ｆ，ｍ}を観測信号の１つ（ここでは複素スペクトルＸ_{ｎ，ｆ，１}）に乗算し、分離信号の複素スペクトルＹ_{ｎ，ｆ，ｍ}を得る。
Y_n,f,m=X_n,f,1・M_n,f,m (5)

最後に、時間領域変換部１０５において、得られた分離信号の複素スペクトルＹ_{ｎ，ｆ，ｍ}を時間領域信号に変換し、分離信号ｙ_ｍ（ｔ）を得る。

特開２０１０−１８７０６６号公報

上述したようなマイクロホン間位相差のクラスタリングに基づく音源分離では、式（４−２）における２π・ｆ・μ_ｍがとる値が、−πからπの範囲を越えることがある場合に、性能が不十分になる欠点がある。これは、式（４−２）の２π・ｆ・μ_ｍは−πからπの範囲を越える値を取ることができるにもかかわらず、式（４−１）のＡ_ｎ，ｆは

の関係を満たす

として、−πからπの範囲に丸められてしまうことに起因する。ここでｋは、μ_ｍの値が既知であれば一意に決まる整数である。これを「空間的エイリアジング」という。空間的エイリアジングは、２つのマイクロホン間隔が広いためにマイクロホン間時間差μ_ｍが大きい場合や周波数ｆが高い場合に、２π・ｆ・μ_ｍの絶対値が大きくなり生じる。

図２は、空間的エイリアジングの例を示す。この図の横軸は周波数ｆを表す、縦軸はマイクロホン間位相差Ａ_ｎ，ｆを表す。具体的には、３０度方向の音源１（ｍ＝１）および１２０度方向の音源２（ｍ＝２）から発せられた音源信号を４ｃｍ間隔の２個のマイクロホンで観測し、サンプリング周波数ｆ_ｓを１６ｋＨｚとした場合の２π・ｆ・μ_ｍ（ｍ＝１，２）の値（図の破線）と式（４−１）のＡ_ｎ，ｆの値（図の実線）をプロットしたものである。ただし、２個のマイクロホンを通る直線に直交するマイクロホン正面を９０度と定義している。

図２の実線に示すように、マイクロホン間位相差の値Ａ_ｎ，ｆがπから−πの範囲に丸められ、複数個の音源に対応するマイクロホン間位相差が同じ値をとる周波数が存在する（図２中の○で示した箇所）。このような周波数では、マイクロホン間位相差を手掛かりとした音源の分類ができないため、従来法では、このような周波数における音源分離が不可能となる。そのため、空間的エイリアジングが生じる場合の音源分離性能が不十分であった。

本発明はこのような点に鑑みてなされたものであり、空間的エイリアジングの問題が生じる場合にも、高い性能で音源分離を行うことが可能な技術を提供することを目的とする。

本発明では、２個のマイクロホンで観測された観測信号を、所定の時間区間ごとに周波数領域表現し、観測信号の複素スペクトルを得、観測信号の複素スペクトルのマイクロホン間位相差を各時間周波数について計算する。任意の１個のマイクロホンで観測された観測信号から得られた観測信号の複素スペクトル、特徴量作成部で得られたマイクロホン間位相差、および音源信号の複素スペクトルの分布を表す事前情報を用い、各時間周波数での各音源に対応する観測信号の複素スペクトルおよびマイクロホン間位相差の分布を表す確率モデルのパラメタを推定する。推定されたパラメタから得られる各時間周波数での各音源の観測信号の複素スペクトルおよびマイクロホン間位相差に対する寄与率を用い、各時間周波数での当該寄与率および観測信号の複素スペクトルから、分離信号の複素スペクトルを抽出し、それを時間領域の分離信号に変換する。

本発明では、マイクロホン間位相差だけでなく、音源信号の複素スペクトルの分布を表す事前情報も用いて音源の分類を行う。これにより、空間的エイリアジングの問題が生じる場合にマイクロホン間位相差が重なり分類ができなかった周波数についても、音源信号の複素スペクトルの分布を手掛りに音源を分類することができる。そのため、空間的エイリアジングの問題が生じる場合にも、高い性能で音源分離を行うことができる。

図１は従来の音源分離装置の構成を例示する図である。 図２は空間的エイリアジングの例を示す図である。 図３は実施形態の音源分離装置の全体構成を例示する図である。 図４はパラメタ推定部の詳細構成を例示するための図である。 図５は実施形態のパラメタ推定部の処理を説明するための図である。 図６Ａは従来の信号分離方法で得られた分離信号の複素スペクトルを例示した図である。図６Ｂは実施形態の信号分離方法で得られた分離信号の複素スペクトルを例示した図である。

以下、図面を参照して本発明の実施形態を説明する。なお、以下の説明において、テキスト中で使用する記号「^→」は、本来直前の文字の真上に記載されるべきものであるが、テキスト記法の制限により、当該文字の直後に記載する。式中においてはこれらの記号は本来の位置に記述している。

図３に例示するように、本形態の信号分離装置１０は、周波数領域変換部１１、特徴量作成部１２、パラメタ推定部１３、音源分離部１４、および時間領域変換部１５を有する。図４に例示するように、本形態のパラメタ推定部１３は、事後確率推定部１３１ａを含むＥステップ計算部１３１、平均パラメタ推定部１３２ａと分散パラメタ推定部１３２ｂと重みパラメタ推定部１３２ｃとスペクトルパラメタ推定部１３２ｄとを含むＭステップ計算部１３２、スペクトル事前分布情報保持部１３３、および判定部１３４を有する。

本形態の信号分離装置１０は、ＣＰＵ（central processing unit）やＲＡＭ（random-access memory）などを含む公知または専用のコンピュータに所定のプログラムが読み込まれることで構成される装置である。信号分離装置１０への入力情報や各部での計算結果は、図示していない一時メモリに記憶され、必要に応じて読み出されて他の処理に利用される。信号分離装置１０は、図示していない制御部を有し、信号分離装置１０の処理は制御部の制御のもとで実行される。

複数の音源から発せられた音源信号が空間で混合され、図示していない２個のマイクロホン１，２で観測される。２個のマイクロホンで観測された観測信号はデジタル化され、時間領域の観測信号ｘ（ｔ）＝［ｘ_１（ｔ），ｘ_２（ｔ）］^Ｔが得られる。ここで、ｘ_Ｌ（ｔ）はマイクロホンＬ（Ｌ＝１，２）で観測された信号に対応する時間領域の観測信号である。

＜周波数領域変換部１１＞
時間領域の観測信号ｘ_Ｌ（ｔ）（Ｌ＝１，２）は周波数領域変換部１１に入力される。周波数領域変換部１１は、観測信号ｘ_Ｌ（ｔ）（Ｌ＝１，２）を、時間フレーム（所定の時間区間）ごとに周波数領域表現し、観測信号の複素スペクトルＸ_{ｎ，ｆ，Ｌ}（Ｌ＝１，２）を生成して出力する。例えば周波数領域変換部１１は、観測信号ｘ_Ｌ（ｔ）（Ｌ＝１，２）を、短時間フーリエ変換により、周波数毎の時系列信号である複素スペクトルＸ_{ｎ，ｆ，Ｌ}（Ｌ＝１，２）に変換して出力する。

＜特徴量作成部１２＞
複素スペクトルＸ_{ｎ，ｆ，Ｌ}（Ｌ＝１，２）は特徴量作成部１２に入力される。特徴量作成部１２は、任意の１個のマイクロホン（本形態の例ではマイクロホン１）で観測された観測信号から得られた複素スペクトルＸ_{ｎ，ｆ，１}を複素スペクトルＸ_ｎ，ｆとして出力する。また特徴量作成部１２は、各時間周波数（ｎ，ｆ）について、複素スペクトルＸ_{ｎ，ｆ，Ｌ}（Ｌ＝１，２）のマイクロホン間位相差Ａ_ｎ，ｆを式（４−１）に従って計算して出力する。

＜パラメタ推定部１３＞
本形態では、上記の複素スペクトルＸ_ｎ，ｆとマイクロホン間位相差Ａ_ｎ，ｆを、（モデル）パラメタθ（後述）による確率モデルで表し、この確率モデルのパラメタθと各時間周波数（ｎ，ｆ）で各音源ｍから発せられる音源信号の存在確率ｐ（ｍ｜Ｘ_ｎ，ｆ，Ａ_ｎ，ｆ，θ）とを推定することで音源分離を行なう。そのため、次にパラメタ推定部１３において、パラメタθの推定を行なう。以降、パラメタ推定部１３で用いる確率モデルを説明し、その後、パラメタ推定部１３の各構成部の処理を説明する。

≪本形態における確率モデル≫
［観測信号の生成モデル］
確率モデルのパラメタをまとめてθと表現し（詳しくは後述）、全ての時間周波数（ｎ，ｆ）におけるＸ_ｎ，ｆおよびＡ_ｎ，ｆの集合をそれぞれ｛Ｘ_ｎ，ｆ｝および｛Ａ_ｎ，ｆ｝とすると、全ての時間周波数（ｎ，ｆ）での各音源ｍに対応する観測信号の複素スペクトルＸ_ｎ，ｆおよびマイクロホン間位相差Ａ_ｎ，ｆの分布を表す確率モデルは以下のようになる。

ここで、Ｍはモデル数であり、ここでは既知の音源数Ｎ_ｓをモデル数Ｍとする（Ｍ＝Ｎ_ｓ）。ｐ（θ）は、モデルパラメタの事前分布であり、ｐ（ｍ）は各音源ｍが観測信号中にそれぞれどんな割合で含まれているかを示す混合重みである（Σ_ｍ＝１ ^Ｍｐ（ｍ）＝１）。

ここで、複素スペクトルＸ_ｎ，ｆとマイクロホン間位相差Ａ_ｎ，ｆは互いに独立であると仮定し、以下のように上記の確率モデルを表現する。
p(X_n,f,A_n,f|m,θ_X,θ_A)=p(X_n,f|m,θ_X)・p(A_n,f|m,θ_A) (8)
ここでｐ（Ｘ_ｎ，ｆ｜ｍ，θ_Ｘ）とｐ（Ａ_ｎ，ｆ｜ｍ，θ_Ａ）はそれぞれ、観測信号の複素スペクトルの確率モデルとマイクロホン間位相差の確率モデルであり、それぞれ以下のように定義する。

[マイクロホン間位相差の確率モデルｐ（Ａ_ｎ，ｆ｜ｍ，θ_Ａ）]
支配的な音源ｍについて、式（６）における整数ｋが既知であると仮定すると、マイクロホン間位相差Ａ_ｎ，ｆは、平均がμ_ｍ、分散がσ_ｍ ^２の以下の正規分布Ｎでモデル化できる。

すなわちマイクロホン間位相差の確率モデルのパラメタは、全ての音源ｍにおけるμ_ｍ，σ_ｍ ^２の集合θ_Ａ＝｛μ_ｍ，σ_ｍ ^２｝である。

ただし実際は、整数ｋは未知であり、整数ｋの推定が必要であるため、マイクロホン間位相差の確率モデルｐ（Ａ_ｎ，ｆ｜ｍ，θ_Ａ）は、上記の式（９−１）を整数ｋで周辺化した以下の式で与えられる。

ここで、ｐ（ｋ）は各ｋに対応する定数である。ｋは０を含む所定範囲内（−Ｋ≦ｋ≦Ｋ）に含まれる各整数である。Ｋはｋの絶対値の最大値であり、周波数ｆによって異なる値をとりうる。Ｋの例は以下の通りである。

はαよりも大きなαに最も近い整数、Ｄはマイクロホン間隔、ｃは音速を表す。また全ての周波数ｆで同じＫを用いても良い。ただし、この場合のＫは十分大きな値（例えば、サンプリング周波数ｆ_ｓの１／２の値をｆとして式（９−３）から得られるＫなど）であることが望ましい。

［観測信号の複素スペクトルの確率モデルｐ（Ｘ_ｎ，ｆ｜ｍ，θ_Ｘ）］
本形態でも音源のスパース性（式（２））を仮定する。ただし本形態では、式（２）においてＨ_{ｆ，Ｌ，ｍ}＝１を仮定する。この仮定はスパース性の一般性を損なわない。この仮定の下では、スパース性とは、ある時間周波数（ｎ，ｆ）において、ｍ番目の音源の音源信号のみが観測され（当該時間周波数（ｎ，ｆ）での複素スペクトルＳ_{ｎ，ｆ，ｍ}をＸ_{ｎ，ｆ，ｍ}と表す）、その他の音源ｍ’の音源信号は０として観測される（すなわち、Ｓ_{ｎ，ｆ，ｍ’}＝０，ｍ’≠ｍ，ｍ’∈｛１，...，Ｍ｝）、と解釈できる。これより、観測信号の複素スペクトルの確率モデルｐ（Ｘ_ｎ，ｆ｜ｍ，θ_Ｘ）を、以下のようにモデル化する。

ここでｐ（Ｓ_{ｎ，ｆ，ｍ}）は、音源信号ｍの複素スペクトルＳ_{ｎ，ｆ，ｍ}の確率モデルである。本形態では、平均０、分散γ_{ｎ，ｆ，ｍ} ^２の以下の複素正規分布Ｎ_ｃでｐ（Ｓ_{ｎ，ｆ，ｍ}）をモデル化する。
p(s_n,f,m) = N_c(s_n,f,m; 0, γ_n,f,m ²) (11)
ここで、分散γ_{ｎ，ｆ，ｍ} ^２は、各音源信号の複素スペクトルＳ_{ｎ，ｆ，ｍ}のパワーの期待値Ｅ［｜ｓ_{ｎ，ｆ，ｍ}｜^２］に相当する。

本形態では、分散γ_{ｎ，ｆ，ｍ} ^２の対数値であるスペクトルパラメタρ_{ｎ，ｆ，ｍ}について、事前分布を与える。
ρ_n,f,m=log(γ_n,f,m ²) (12)
ただし、「ｌｏｇ」は自然対数を表す。

具体的には、まず各音源ｍについて、各時刻ｎでの全ての周波数ｆにおけるスペクトルパラメタρ_{ｎ，ｆ，ｍ}を要素とするスペクトルベクトルρ^→ _ｎ，ｍ＝［ρ_{ｎ，１，ｍ}，...，ρ_{ｎ，Ｎｆ，ｍ}］^Ｔを考える。ただしｆ∈｛１，...，Ｎ_ｆ｝であり、添え字の「Ｎｆ」はＮ_ｆを意味する。次に、スペクトルベクトルρ^→ _ｎ，ｍに対する事前分布ｐ（ρ^→ _ｎ，ｍ）を、以下の混合正規分布でモデル化する（混合正規分布モデル：Gaussian Mixture Model（GMM））。

ここで、Ｉは上記の混合正規分布に用いられる正規分布の個数であり、ｉ（ｍ）は音源ｍのスペクトルベクトルρ^→ _ｎ，ｍの混合正規分布に含まれる各正規分布を表すインデックスである。ν^→ _ｉ（ｍ）＝［ν_{１，ｉ（ｍ）}，...，ν_{Ｎｆ，ｉ（ｍ）}］^Ｔ、Ｖ_ｉ（ｍ）、およびｐ（ｉ（ｍ））は、それぞれ、当該ＧＭＭの平均、共分散行列（対角要素がυ_{ｆ，ｉ（ｍ）}の対角行列）、および混合重みである。これらのパラメタは、学習用データを用いた事前学習によって得られる（学習法は後述）。

この事前分布の情報を用いて音源分離を行うことにより、音源信号の分布に近い分布を持つ分離信号を得ることができる。すなわち、空間的エイリアジングのために特定の周波数で大きく歪んだ周波数特性を持つ信号が分離信号として得られることの無いよう、制御することが可能となる。また上記より、観測信号の複素スペクトルの確率モデルのパラメタは、全ての時間周波数（ｎ，ｆ）および音源ｍにおけるρ_{ｎ，ｆ，ｍ}の集合θ_Ｘ＝｛ρ_{ｎ，ｆ，ｍ}｝である。なお、分散γ_{ｎ，ｆ，ｍ} ^２の対数値に事前分布を与えてモデル化するのではなく、分散γ_{ｎ，ｆ，ｍ} ^２そのものに事前分布を与えてモデル化してもよい。すなわち、式（１２）に代えてρ_{ｎ，ｆ，ｍ}＝γ_{ｎ，ｆ，ｍ} ^２としてもよい。

≪パラメタ推定部１３の詳細≫
パラメタ推定部１３（図４）は、パラメタθ＝｛θ_Ａ，ｐ（ｍ），θ_Ｘ｝＝｛μ_ｍ，σ_ｍ ^２，ｐ（ｍ），ρ_{ｎ，ｆ，ｍ}｝を繰り返し演算によって求める。以下θ^ｔ＝｛μ_ｍ ^ｔ，（σ_ｍ ^２）^ｔ，ｐ^ｔ（ｍ），ρ_{ｎ，ｆ，ｍ} ^ｔ｝をｔ回目の反復更新で得られたパラメタとし、θ_Ａ ^ｔ＝｛μ_ｍ ^ｔ，（σ_ｍ ^２）^ｔ｝およびθ_Ｘ ^ｔ＝｛ρ_{ｎ，ｆ，ｍ} ^ｔ｝とする。

図５に例示するように、はじめにステップＳ０において、パラメタθの初期値θ^０＝｛μ_ｍ ^０，（σ_ｍ ^２）^０，ｐ^０（ｍ），ρ^０ _{ｎ，ｆ，ｍ}｝、パラメタ更新回数を表すインデックスｔの初期値ｔ＝０、モデル数Ｍ、ｋの範囲を決めるＫ、パラメタ更新回数の最大値Ｔ、収束判定の閾値Δを設定する。パラメタθ^０の初期値としては、例えば、マイクロホン間隔Ｄで決まるμ_ｍのとりうる範囲（−Ｄ／ｃ〜Ｄ／ｃ）をＭ等分した値をμ_ｍ ^０、（σ_ｍ ^２）^０＝π／３０（ｍ＝１，...，Ｍ）、ｐ^０（ｍ）＝１／Ｍ（ｍ＝１，...，Ｍ）、ρ^０ _{ｎ，ｆ，ｍ}＝ｌｏｇ（Ｘ_ｎ，ｆ ^２）（ｍ＝１，...，Ｍ）とする。また、モデル数Ｍは既知の音源数Ｎ_ｓであり、Ｋは例えば前述の式（９−３）から得られる値であり、最大値Ｔは例えばＴ＝２０であり、閾値Δは例えば十分小さな正の値である。

また、予め学習しておいたスペクトルベクトルρ^→ _ｎ，ｍに対する事前分布ｐ（ρ^→ _ｎ，ｍ）のパラメタ｛ν_{ｆ，ｉ（ｍ）}，υ_{ｆ，ｉ（ｍ）}，ｐ（ｉ（ｍ））｝を、スペクトル事前分布情報保持部１３３に保持しておく。パラメタ｛ν_{ｆ，ｉ（ｍ）}，υ_{ｆ，ｉ（ｍ）}，ｐ（ｉ（ｍ））｝は、音源信号の複素スペクトルＳ_{ｎ，ｆ，ｍ}の分布を表す事前情報である。

事前分布ｐ（ρ^→ _ｎ，ｍ）の学習方法としては、以下の方法を例示できる。音源ｍから予め得ていた学習用の複素スペクトルをＳ_{ｎ，ｆ，ｍ}’としてρ_{ｎ，ｆ，ｍ}＝ｌｏｇ（(Ｓ_{ｎ，ｆ，ｍ}’)^２）を計算し、それらを用いてスペクトルベクトルρ^→ _ｎ，ｍ＝［ρ_{ｎ，１，ｍ}，...，ρ_{ｎ，Ｎｆ，ｍ}］^Ｔを構成する。このスペクトルベクトルρ^→ _ｎ，ｍに（１３）の混合正規分布をあてはめて、パラメタ｛ν_{ｆ，ｉ（ｍ）}，υ_{ｆ，ｉ（ｍ）}，ｐ（ｉ（ｍ））｝を計算する。混合正規分布のあてはめには、例えばＥＭアルゴリズム（参考文献１「汪ほか、「計算統計Ｉ−確率計算の新しい手法」、岩波書店」）を用いることができる。この場合には、参考文献１のＰ１６１の式（７６）のｘをスペクトルベクトルρ^→ _ｎ，ｍに、μ_ｉをν^→ _ｉ（ｍ）に、Σ_ｉをＶ_ｉ（ｍ）に、α_ｉをｐ（ｉ（ｍ））にそれぞれ置き換え、参考文献１の式（７８）（７９）（８０）にて、それぞれのパラメタを計算し、それらを｛ν_{ｆ，ｉ（ｍ）}，υ_{ｆ，ｉ（ｍ）}，ｐ（ｉ（ｍ））｝とする。

次にステップＳ１において、制御部（図示せず）がｔを１だけインクリメントする（ｔ＋１を新たなｔとする）。

［Ｅステップ計算部１３１］
次にステップＳ２（Ｅステップ）において、Ｅステップ計算部１３１の事後確率推定部１３１ａが、入力された複素スペクトルＸ_ｎ，ｆおよびマイクロホン間位相差Ａ_ｎ，ｆ、スペクトル事前分布情報保持部１３３に保存されているパラメタ｛ν_{ｆ，ｉ（ｍ）}，υ_{ｆ，ｉ（ｍ）}，ｐ（ｉ（ｍ））｝、ならびに一時メモリ（図示せず）に保存されているパラメタθ^ｔ−１＝｛μ_ｍ ^ｔ−１，（σ_ｍ ^２）^ｔ−１，ｐ^ｔ−１（ｍ），ρ^ｔ−１ _{ｎ，ｆ，ｍ}｝を用い、ｐ（Ｘ_ｎ，ｆ｜ｍ，θ_Ｘ ^ｔ−１）（式（１０）参照）、ｐ（Ａ_ｎ，ｆ｜ｍ，ｋ，θ_Ａ ^ｔ−１）（式（９−１）参照）、ｐ（ρ^→ _ｎ，ｍ ^ｔ）およびｐ（ρ^→ _ｎ，ｍ ^ｔ｜ｉ（ｍ））・ｐ（ｉ（ｍ））（式（１３）（１４）参照）を計算した上で、以下の事後確率を計算する。
ここで、Ｇ_{ｎ，ｆ，ｍ}は時間周波数（ｎ，ｆ）での音源ｍの複素スペクトルＸ_ｎ，ｆに対する寄与率を表し、Ｇ_{ｎ，ｆ，ｍ} ^ｔはｔ回目の反復更新で得られる寄与率Ｇ_{ｎ，ｆ，ｍ}を表す。ｚ_{ｎ，ｉ（ｍ）}は時刻ｎでのスペクトルベクトルρ^→ _ｎ，ｍがｉ（ｍ）番目の混合分布をとる事後確率を表す。ρ^→ _ｎ，ｍ ^ｔ＝［ρ_{ｎ，１，ｍ} ^ｔ，...，ρ_{ｎ，Ｎｆ，ｍ} ^ｔ］^Ｔはｔ回目の反復更新で得られるスペクトルパラメタρ_{ｎ，ｆ，ｍ} ^ｔを要素とするスペクトルベクトルを表す。また（１６）以降、Σ_ｋα（ｋ）はｋ＝−ＫからＫまでについてのα（ｋ）の和Σ_ｋ＝−Ｋ ^Ｋα（ｋ）、Σ_ｍα（ｍ）はｍ＝１からＭまでについてのα（ｍ）の和Σ_ｍ＝１ ^Ｍα（ｍ）、Σ_ｎα（ｎ）はｎ＝１からＮ_ｎまでについてのα（ｎ）の和Σ_ｎ＝１ ^Ｎｎα（ｎ）、Σ_ｆα（ｆ）はｆ＝１からＮ_ｆまでについてのα（ｆ）の和Σ_ｆ＝１ ^Ｎｆα（ｆ）である。ここで、Ｎ_ｆは周波数の個数であり、Ｎ_ｎはデータ全体での時間フレームの個数である。

［Ｍステップ計算部１３２］
次にステップＳ３（Ｍステップ）において、Ｍステップ計算部１３２が、入力された複素スペクトルＸ_ｎ，ｆおよびマイクロホン間位相差Ａ_ｎ，ｆ、ステップＳ２で得られた事後確率ｇ_{ｎ，ｆ，ｍ，ｋ} ^ｔ，ｚ_{ｎ，ｉ（ｍ）} ^ｔ、ならびにスペクトル事前分布情報保持部１３３に保存されている事前分布ｐ（ρ^→ _ｎ，ｍ）のパラメタ｛ν_{ｆ，ｉ（ｍ）}，υ_{ｆ，ｉ（ｍ）}，ｐ（ｉ（ｍ））｝を用い、以下のようにパラメタθ^ｔ＝｛μ_ｍ ^ｔ，（σ_ｍ ^２）^ｔ，ｐ^ｔ（ｍ），ρ_{ｎ，ｆ，ｍ} ^ｔ｝を得て出力する（更新式の導出は後述）。

［平均パラメタ推定部１３２ａ］
平均パラメタ推定部１３２ａは、マイクロホン間位相差Ａ_ｎ，ｆおよび式（１６）で得られた事後確率ｇ_{ｎ，ｆ，ｍ，ｋ} ^ｔを入力とし、以下のように平均μ_ｍ ^ｔ（各音源ｍに対応するマイクロホン間位相差Ａ_ｎ，ｆの平均）を計算して出力する。

すなわち、平均μ_ｍ ^ｔはマイクロホン間位相差Ａ_ｎ，ｆに２πｋを加算して得られる値の集合と複素スペクトルＸ_ｎ，ｆとから得られる値である。

［分散パラメタ推定部１３２ｂ］
分散パラメタ推定部１３２ｂは、マイクロホン間位相差Ａ_ｎ，ｆ、式（１６）で得られた事後確率ｇ_{ｎ，ｆ，ｍ，ｋ} ^ｔ、および式（１８）で得られた平均μ_ｍ ^ｔを入力とし、以下のように分散（σ_ｍ ^２）^ｔ（各音源ｍに対応するマイクロホン間位相差Ａ_ｎ，ｆの分散）を計算して出力する。

すなわち、分散（σ_ｍ ^２）^ｔはマイクロホン間位相差Ａ_ｎ，ｆに２πｋを加算して得られる値の集合と複素スペクトルＸ_ｎ，ｆとから得られる値である。

［スペクトルパラメタ推定部１３２ｄ］
スペクトルパラメタ推定部１３２ｄは、複素スペクトルＸ_ｎ，ｆ、式（１５）で得られた寄与率Ｇ_{ｎ，ｆ，ｍ} ^ｔ、式（１７）で得られた事後確率ｚ_{ｎ，ｉ（ｍ）} ^ｔ、およびスペクトル事前分布情報保持部１３３に保存されているパラメタν_{ｆ，ｉ（ｍ）}，υ_{ｆ，ｉ（ｍ）}を用い、以下の式（２０）を満たすスペクトルパラメタρ_{ｎ，ｆ，ｍ} ^ｔ（各音源ｍに対応する複素スペクトルＸ_ｎ，ｆの分散を表すスペクトルパラメタ）を計算して出力する。

ただし、式（２０）はρ_{ｎ，ｆ，ｍ} ^ｔについて解析的に解けないため、ニュートン法などを用いて、この式を満たすρ_{ｎ，ｆ，ｍ} ^ｔを計算する（具体的な方法は後述）。

［重みパラメタ推定部１３２ｃ］
重みパラメタ推定部１３２ｃは、式（１５）で得られた寄与率Ｇ_{ｎ，ｆ，ｍ} ^ｔを用い、以下のように混合分布の重みｐ^ｔ（ｍ）を計算して出力する。

［判定部１３４］
判定部１３４は、ｔ＝Ｔとなるか、もしくは、以下の収束条件｜Ｑ（θ｜θ^ｔ）−Ｑ（θ｜θ^ｔ−１）｜＜Δを満たすかを判定する（ステップＳ４）。
ただし、Ｑ（θ｜θ^ｔ）は、ｔ回目の反復更新で得られた事後確率ｇ_{ｎ，ｆ，ｍ，ｋ} ^ｔ，ｚ_{ｎ，ｉ（ｍ）} ^ｔおよびパラメタθ^ｔ＝｛μ_ｍ ^ｔ，（σ_ｍ ^２）^ｔ，ｐ^ｔ（ｍ），ρ_{ｎ，ｆ，ｍ} ^ｔ｝を用いて計算した以下のＱ関数を表す。

Ｑ（θ｜θ^ｔ−１）は、ｔ−１回目の反復更新で得られた事後確率ｇ_{ｎ，ｆ，ｍ，ｋ} ^ｔ−１，ｚ_{ｎ，ｉ（ｍ）} ^ｔ−１およびパラメタθ^ｔ−１＝｛μ_ｍ ^ｔ−１，（σ_ｍ ^２）^ｔ−１，ｐ^ｔ−１（ｍ），ρ_{ｎ，ｆ，ｍ} ^ｔ−１｝を用いて計算したＱ関数である。すなわち、収束条件｜Ｑ（θ｜θ^ｔ）−Ｑ（θ｜θ^ｔ−１）｜＜Δは、ｔ回目の反復更新で得られた事後確率ｇ_{ｎ，ｆ，ｍ，ｋ} ^ｔ，ｚ_{ｎ，ｉ（ｍ）} ^ｔおよびパラメタθ^ｔ＝｛μ_ｍ ^ｔ，（σ_ｍ ^２）^ｔ，ｐ^ｔ（ｍ），ρ_{ｎ，ｆ，ｍ} ^ｔ｝から得られるＱ（θ｜θ^ｔ）と、ｔ−１回目の反復更新で得られた事後確率ｇ_{ｎ，ｆ，ｍ，ｋ} ^ｔ−１，ｚ_{ｎ，ｉ（ｍ）} ^ｔ−１およびパラメタθ^ｔ−１＝｛μ_ｍ ^ｔ−１，（σ_ｍ ^２）^ｔ−１，ｐ^ｔ−１（ｍ），ρ_{ｎ，ｆ，ｍ} ^ｔ−１｝から得られるＱ（θ｜θ^ｔ−１）との差がΔ未満であるかの条件である。Σ_{｛ｉ（ｍ）｝}α（ｉ（ｍ））は、Σ_{｛ｉ（1）｝}α（ｉ（１））・Σ_{｛ｉ（2）｝}α（ｉ（２））…・Σ_{｛ｉ（Ｍ）｝}α（ｉ（Ｍ））を表す。

ステップＳ４でｎｏであればステップＳ１に戻る。一方、ｙｅｓであればパラメタ推定部１３が、パラメタθ^ｔ＝｛μ_ｍ ^ｔ，（σ_ｍ ^２）^ｔ，ｐ^ｔ（ｍ），ρ_{ｎ，ｆ，ｍ} ^ｔ｝をθ＝｛μ_ｍ，（σ_ｍ ^２），ｐ（ｍ），ρ_{ｎ，ｆ，ｍ}｝とし、パラメタθ^ｔから得られる寄与率Ｇ_{ｎ，ｆ，ｍ} ^ｔを寄与率Ｇ_{ｎ，ｆ，ｍ}として得て出力する。

次に、音源分離部１４において音源分離を行なう。最終的に得られたパラメタθから得られる事後確率Ｇ_{ｎ，ｆ，ｍ}（式（１５））の値は、各時間周波数（ｎ，ｆ）での各音源の観測信号の複素スペクトルＸ_{ｎ，ｆ，Ｌ}（Ｌ＝１，２）に対する寄与率に相当する。そのため、本形態の音源分離部１４は、各時間周波数（ｎ，ｆ）での当該寄与率Ｇ_{ｎ，ｆ，ｍ}および観測信号の複素スペクトルＸ_{ｎ，ｆ，Ｌ}（Ｌ＝１，２）の１つ（ここではＸ_ｎ，ｆ＝Ｘ_{ｎ，ｆ，１}）を入力とし、これらから、各音源ｍに対応する分離信号の複素スペクトルＹ_{ｎ，ｆ，ｍ}を抽出する。例えば、音源分離部１４は、寄与率Ｇ_{ｎ，ｆ，ｍ}を観測信号の複素スペクトルＸ_{ｎ，ｆ，１}に以下のように乗算し、分離信号の複素スペクトルＹ_{ｎ，ｆ，ｍ}を得る。
Y_n,f,m=G_n,f,m・X_n,f,1(22)

最後に時間領域変換部１５が、分離信号の複素スペクトルＹ_{ｎ，ｆ，ｍ}を入力とし、例えば逆フーリエ変換等によって、複素スペクトルＹ_{ｎ，ｆ，ｍ}を時間領域の分離信号ｙ_ｍ（ｔ）に変換して出力する。

＜Ｍステップの各計算式の導出＞
ここでは、Ｍステップ計算部１３２（ステップＳ３）が計算する各パラメタの導出法について述べる。本形態では、観測信号の複素スペクトルＸ_ｎ，ｆとマイクロホン間位相差Ａ_ｎ，ｆについてのパラメタθ＝｛μ_ｍ，σ_ｍ ^２，ｐ（ｍ），ρ_{ｎ，ｆ，ｍ}｝、ＭＡＰ推定により推定する。｛α_ｎ，ｆ｝を全ての時間周波数（ｎ，ｆ）におけるα_ｎ，ｆの集合、｛α_ｍ｝を全ての音源ｍにおけるα_ｍの集合とすると、コスト関数は、次のように与えられる。

このコスト関数Ｌ（θ）を最大とするパラメタθをＥＭアルゴリズムにより求める。ＥＭアルゴリズムのためのＱ関数は以下のようになる。

ここでｐ（Ｘ_ｎ，ｆ，Ａ_ｎ，ｆ｜ｍ，ｋ，θ）は式（８）、ｇ_{ｎ，ｆ，ｍ，ｋ}とｚ_{ｎ，ｉ（ｍ）}はそれぞれ（１６）、（１７）で与えられ、Ｑ関数の最後の項のｐ（ρ^→ _ｎ，ｍ｜ｉ（ｍ））・ｐ（ｉ（ｍ））はｐ（ｉ（ｍ））・Π_ｆＮ（ρ_{ｎ，ｆ，ｍ}；ν_{ｆ，ｉ（ｍ）}，υ_{ｆ，ｉ（ｍ）}）である（式（１４））。

Ｑ関数の平均μ_ｍ、分散σ_ｍ ^２、スペクトルパラメタρ_{ｎ，ｆ，ｍ}での偏微分結果をそれぞれゼロと置いて解くことで、式（１８）、（１９）、（２０）が得られる。また、式（２１）の重みのパラメタは、

の拘束のもとで解くことで得られる（導出終）

＜スペクトルパラメタ推定部１３２ｄでの具体的な計算＞
スペクトルパラメタ推定部１３２ｄでは、式（２０）を満たすスペクトルパラメタρ_{ｎ，ｆ，ｍ}を計算するが、この式は、ρ_{ｎ，ｆ，ｍ}について解析的に解けないため、ここではニュートン法を用いてρ_{ｎ，ｆ，ｍ}を計算する方法を示す。

式（２０）を整理して以下のように表現する。

ここで、ａ，ｂ，ｃはそれぞれ以下である。

さらに、

となる。

この式にＮｅｗｔｏｎ法を適用する。まずＹを０と１の間の任意の値などで初期化した後、

の計算を数回（５回程度）反復して行なう。ただし、「α←β」は、βを新たなαとすることを意味し、Ｆ’（Ｙ）は関数Ｆ（Ｙ）をＹについて微分することを意味する。これにより得た解Ｙを用いて、
ρ_n,f,m=log(b/a)-Y (26)
を計算し、スペクトルパラメタρ_{ｎ，ｆ，ｍ}を得る。

［本形態の効果］
本形態の効果を説明する。実験条件としては、２個または３個の音源が同時に鳴っている環境において２つのマイクロホンで収録を行い、従来法および提案法で音源分離を行なった。サンプリング周波数は１６ｋＨｚ、マイク間隔は２０ｃｍであり、８５０Ｈｚ以上で空間的エイリアジングが起こる条件である。

ここでの従来法は、特許文献１に示されている方法である。また提案法におけるスペクトルベクトルの事前分布の学習には、音源分離対象と同じ話者が話した１４３文（分離対象の文章は含まない）を用いた。音源分離性能を、信号対妨害音比（Signal to interference ratio: SIR）および、信号対歪み比（Signal to distortion ratio: SDR）にて評価した。６通りの音声組合せにて評価し、その平均を求めた。

図６Ａおよび図６Ｂは、それぞれ、従来法および提案法によって分離された信号のスペクトル（音源数３）を例示している。図６Ａ（従来法）では、音源（２）と音源（３）のスペクトルにおいて、横縞様の成分が多く見られる。例えば１秒を過ぎた付近と、２．５秒付近は、音源（２）と音源（３）は無音であるにもかかわらず横縞様の成分が見える。これは、空間的エイリアジングのためマイク間位相差が重なった周波数において、音源分離が十分に行なわれず、音源（１）の信号成分が、音源（２）と音源（３）の分離信号に残ってしまったためである。一方、図６Ｂ（提案法）では、横縞様の成分は減少していることが分かる（四角で囲んだ部分など）。これより、提案法は、空間的エイリアジングのためマイクロホン間位相差が重なった周波数においても、高精度で音源分離を可能にすることが分かる。

表１は、ＳＩＲとＳＤＲの値を示している。これより、発明法は従来法より高い分離性能を得られることが分かる。

なお、本発明は上述の実施の形態に限定されるものではない。たとえば、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。

上述の構成をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。このプログラムをコンピュータで実行することにより、上記処理機能がコンピュータ上で実現される。この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体の例は、非一時的な（non-transitory）記録媒体である。このような記録媒体の例は、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等である。

このプログラムの流通は、たとえば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。

このようなプログラムを実行するコンピュータは、たとえば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。処理の実行時、このコンピュータは、自己の記録装置に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。

上記実施形態では、コンピュータ上で所定のプログラムを実行させて本装置の処理機能が実現されたが、これらの処理機能の少なくとも一部がハードウェアで実現されてもよい。

１０ 音源分離装置
１２ 特徴量作成部
１３ パラメタ推定部
１４ 音源分離部
１３１ Ｅステップ計算部
１３１ａ 事後確率推定部
１３２ａ 平均パラメタ推定部
１３２ｂ 分散パラメタ推定部
１３２ｃ 重みパラメタ推定部
１３２ｄ スペクトルパラメタ推定部
１３３ スペクトル事前分布情報保持部

Claims (7)

1. 複数の音源が発した音源信号が混合され、２個のマイクロホンで観測される状況で、上記音源信号それぞれの推定信号である分離信号を分離抽出する音源分離装置であって、
   ２個のマイクロホンで観測された観測信号を、所定の時間区間ごとに周波数領域表現し、観測信号の複素スペクトルを得る周波数領域変換部と、
   上記観測信号の複素スペクトルのマイクロホン間位相差を各時間周波数について計算する特徴量作成部と、
   任意の１個のマイクロホンで観測された観測信号から得られた上記観測信号の複素スペクトル、上記特徴量作成部で得られた上記マイクロホン間位相差、および上記音源信号の複素スペクトルの分布を表す事前情報を用い、上記各音源に対応する観測信号の複素スペクトルおよびマイクロホン間位相差の分布を表す確率モデルのパラメタを推定するパラメタ推定部と、
   上記パラメタから得られる各時間周波数での上記各音源の上記観測信号の複素スペクトルに対する寄与率を用い、各時間周波数での当該寄与率および上記観測信号の複素スペクトルから、分離信号の複素スペクトルを抽出する音源分離部と、
   上記分離信号の複素スペクトルを時間領域の分離信号に変換する時間領域変換部と、
   を有する音源分離装置。
2. 請求項１に記載の音源分離装置であって、
   上記パラメタ推定部が、
   上記各音源に対応する上記マイクロホン間位相差の平均を全ての音源について計算する平均パラメタ推定部と、
   上記各音源に対応する上記マイクロホン間位相差の分散を全ての音源について計算する分散パラメタ推定部と、
   上記各音源に対応する上記観測信号の複素スペクトルの分散を表すスペクトルパラメタを全ての音源について計算するスペクトルパラメタ推定部と、
   上記音源信号の複素スペクトルの分布を表す事前情報を保持しているスペクトル事前分布情報保持部と、
   混合分布の重みを計算する重みパラメタ推定部と、を有し、
   上記確率モデルは、混合分布モデルであり、
   上記パラメタは、上記マイクロホン間位相差の平均、上記マイクロホン間位相差の分散、上記スペクトルパラメタ、および上記混合分布の重みを含み、
   上記マイクロホン間位相差の平均、上記マイクロホン間位相差の分散、および上記混合分布の重みは、上記観測信号の複素スペクトルと上記マイクロホン間位相差とから得られる値であり、
   上記観測信号の複素スペクトルの分散は、上記観測信号の複素スペクトルと上記マイクロホン間位相差と上記音源信号の複素スペクトルの分布を表す事前情報とから得られる値である、
   ことを特徴とする音源分離装置。
3. 請求項２に記載の音源分離装置であって、
   上記マイクロホン間位相差の平均、および上記マイクロホン間位相差の分散は、上記マイクロホン間位相差に２πｋを加算して得られる値の集合と上記観測信号の複素スペクトルとから得られる値であり、ｋは０を含む所定範囲内に含まれる各整数である、
   ことを特徴とする音源分離装置。
4. 複数の音源が発した音源信号が混合され、２個のマイクロホンで観測される状況で、上記音源信号それぞれの推定信号である分離信号を分離抽出する音源分離方法であって、
   ２個のマイクロホンで観測された観測信号を、所定の時間区間ごとに周波数領域表現し、観測信号の複素スペクトルを得る周波数領域変換ステップと、
   上記観測信号の複素スペクトルのマイクロホン間位相差を各時間周波数について計算する特徴量作成ステップと、
   任意の１個のマイクロホンで観測された観測信号から得られた上記観測信号の複素スペクトル、上記特徴量作成部で得られた上記マイクロホン間位相差、および上記音源信号の複素スペクトルの分布を表す事前情報を用い、上記各音源に対応する観測信号の複素スペクトルおよびマイクロホン間位相差の分布を表す確率モデルのパラメタを推定するパラメタ推定ステップと、
   上記パラメタから得られる各時間周波数での上記各音源の上記観測信号の複素スペクトルに対する寄与率を用い、各時間周波数での当該寄与率および上記観測信号の複素スペクトルから、分離信号の複素スペクトルを抽出する音源分離ステップと、
   上記分離信号の複素スペクトルを時間領域の分離信号に変換する時間領域変換ステップと、
   を有する音源分離方法。
5. 請求項４に記載の音源分離方法であって、
   上記パラメタ推定ステップが、
   上記各音源に対応する上記マイクロホン間位相差の平均を全ての音源について計算する平均パラメタ推定ステップと、
   上記各音源に対応する上記マイクロホン間位相差の分散を全ての音源について計算する分散パラメタ推定ステップと、
   上記各音源に対応する上記観測信号の複素スペクトルの分散を全ての音源について表すスペクトルパラメタを計算するスペクトルパラメタ推定ステップと、
   上記音源信号の複素スペクトルの分布を表す事前情報を保持しているスペクトル事前分布情報保持ステップと、
   混合分布の重みを計算する重みパラメタ推定ステップと、を有し、
   上記確率モデルは、混合分布モデルであり、
   上記パラメタは、上記マイクロホン間位相差の平均、上記マイクロホン間位相差の分散、上記スペクトルパラメタ、および上記混合分布の重みを含み、
   上記マイクロホン間位相差の平均、上記マイクロホン間位相差の分散、および上記混合分布の重みは、上記観測信号の複素スペクトルと上記マイクロホン間位相差とから得られる値であり、
   上記観測信号の複素スペクトルの分散は、上記観測信号の複素スペクトルと上記マイクロホン間位相差と上記音源信号の複素スペクトルの分布を表す事前情報とから得られる値である、
   ことを特徴とする音源分離方法。
6. 請求項５に記載の音源分離方法であって、
   上記マイクロホン間位相差の平均、および上記マイクロホン間位相差の分散は、上記マイクロホン間位相差に２πｋを加算して得られる値の集合と上記観測信号の複素スペクトルとから得られる値であり、ｋは０を含む所定範囲内に含まれる各整数である、
   ことを特徴とする音源分離方法。
7. 請求項１から３の何れかの音源分離装置としてコンピュータを機能させるためのプログラム。