"複数の原信号がノイズとともに混合され、二個のマイクで観測される状況で、観測信号からそれぞれの原信号を分離抽出する音源分離装置であって、前記原信号の音源のインデックスをkとし、周波数領域の前記観測信号xf,t=[xf,t,L,xf,t,R]Tから、前記二個のマイクへの前記原信号の到達時間差δf,kと雑音のパワースペクトルσf2と原信号のスペクトルsf,t,kと音源存在確率pk|θとを推定するパラメタ推定手段と、前記観測信号xf,tと推定されたパラメタθ={δf,k,σf2,sf,t,k,pk|θ}とから分離信号yf,t,kを生成する分離信号生成手段とを含み、前記パラメタ推定手段は、観測信号xf,tが音源kに帰属する期待値を示す事後確率をmf,t,kとし、前記二個のマイクの間隔をDとし、前記原信号の速度をcとし、φf=sinc2πfD/cとし、ξf,t,k=[xf,t,R-φfxt,t,L-sf,t,k]とし、前記スペクトルsf,t,k及び前記ξf,t,kの位相をそれぞれψsk及びψξkとし、前記到達時間差δf,kをとして推定する逆正接計算手段と、推定された前記到達時間差δf,kが内包する±πの不定性を補正する補正手段とを有する、音源分離装置。複数の原信号がノイズとともに混合され、二個のマイクで観測される状況で、観測信号からそれぞれの原信号を分離抽出する音源分離方法であって、前記原信号の音源のインデックスをkとし、周波数領域の前記観測信号xf,t=[xf,t,L,xf,t,R]Tから、前記二個のマイクへの前記原信号の到達時間差δf,kと雑音のパワースペクトルσf2と原信号のスペクトルsf,t,kと音源存在確率pk|θとを推定するパラメタ推定ステップと、前記観測信号xf,tと推定されたパラメタθ={δf,k,σf2,sf,t,k,pk|θ}とから分離信号yf,t,kを生成する分離信号生成ステップとを含み、前記パラメタ推定ステップは、観測信号xf,tが音源kに帰属する期待値を示す事後確率をmf,t,kとし、前記二個のマイクの間隔をDとし、前記原信号の速度をcとし、φf=sinc2πfD/cとし、ξf,t,k=[xf,t,R-φfxt,t,L-sf,t,k]とし、前記スペクトルsf,t,k及び前記ξf,t,kの位相をそれぞれψsk及びψξkとし、前記到達時間差δf,kをとして推定する逆正接計算ステップと、推定された前記到達時間差δf,kが内包する±πの不定性を補正する補正ステップとを有する、音源分離方法。コンピュータを請求項1記載の音源分離装置として機能させるためのプログラム。しかしながら、従来技術では、到達時間差δkの解析的な更新式は与えられていないため、多くの計算コストを要する全探索操作によって、到達時間差δkを推定する必要がある。本発明は、到達時間差δの性質に着目し、到達時間差δを効率的に推定する方法を与え、従来技術において必要であった全探索操作を不要とし、高速な音源分離技術を提供することを目的とする。上記の課題を解決するために、本発明の第一の態様によれば、複数の原信号がノイズとともに混合され、二個のマイクで観測される状況で、観測信号からそれぞれの原信号を分離抽出する。原信号の音源のインデックスをkとし、周波数領域の観測信号xf,t=[xf,t,L,xf,t,R]Tから、二個のマイクへの原信号の到達時間差δf,kと雑音のパワースペクトルσf2と原信号のスペクトルsf,t,kと音源存在確率pk|θとを推定する。観測信号xf,tと推定されたパラメタθ={δf,k,σf2,sf,t,k,pk|θ}とから分離信号yf,t,kを生成する。観測信号xf,tが音源kに帰属する期待値を示す事後確率をmf,t,kとし、二個のマイクの間隔をDとし、原信号の速度をcとし、φf=sinc2πfD/cとし、ξf,t,k=[xf,t,R-φfxt,t,L-sf,t,k]とし、スペクトルsf,t,k及びξf,t,kの位相をそれぞれψsk及びψξkとし、到達時間差δf,kをとして推定する。推定された到達時間差δf,kが内包する±πの不定性を補正する。本発明は、到達時間差δを効率的に推定する方法を与え、高速な音源分離ができるという効果を奏する。以下、本発明の実施形態について説明する。なお、以下の説明に用いる図面では、同じ機能を持つ構成部や同じ処理を行うステップには同一の符号を記し、重複説明を省略する。また、ベクトルや行列の各要素単位で行われる処理は、特に断りが無い限り、そのベクトルやその行列の全ての要素に対して適用されるものとする。<第一実施形態に係る音源分離装置1>図1は本実施形態に係る音源分離装置1の機能ブロック図を、図2はその処理フローを示す。音源分離装置1は周波数領域変換部11とパラメタ推定部12と分離信号生成部13と時間領域変換部14とを含み、さらにパラメタ推定部12はEステップ計算部121とMステップ計算部122とを含む。まず、観測信号について説明する。信号原以下音源とも言うがK個あり、kk=1,2,...,Kを音源のインデックスとし、音源kから発せられる信号以下原信号というをsktとする。複数の原信号s1t,,sKtがノイズとともに二個のマイクL,Rで観測される状況で、マイクLで観測される時間領域の観測信号をxLtとし、マイクRで観測される時間領域の観測信号をxRtとし、二個のマイクL,Rで観測される時間領域の観測信号をxt=[xLt,xRt]Tとする。Tは転置を表す。tはフレーム番号及びそのフレーム番号に対応する時刻を表す。Tを時間フレームの総数とすると、t=0,1,...,T-1である。ここで周波数領域の観測信号をxf,t=[xf,t,L,xf,t,R]Tと表記する。なお、fはサンプリング周波数fsをF等分した離散点であり、f∈{0,fs/F,...,F-1fs/F}である。以降、断りのない場合、観測信号とは、周波数領域の観測信号xf,t=[xf,t,L,xf,t,R]Tを指し、時間領域の観測信号の場合はそれを明記する。ここで、観測信号は、と表されると仮定する。ここで、sf,t,kは原信号sktのスペクトル以下音源スペクトルともいうを、nf,k=[nf,k,L,nf,k,R]はマイクL,Rにおける加算的雑音を表す。またbf,k=[bf,k,L,bf,k,R]Tは音源kに関するステアリングベクトル音源kの方向を特定するベクトルであり、方向ベクトルともいうであり、原信号sktのマイクL、Rへの到達時間差をδkとすると、である非特許文献1参照。本実施形態では、原信号の観測時間内においては、音源及びマイクは固定されており、またK個の音源は全て、異なる位置に配置されているとする。すなわち、ステアリングベクトルbf,kは時間tに依らず、kの値によって異なる値を取るものと仮定する。音源分離の目的は、観測信号xf,tのみを用いて、全ての音源スペクトルsf,t,kを推定することである。本実施形態に係る音源分離装置1は、時間領域の観測信号xt=[xLt,xRt]Tを入力とし、時間領域の分離信号推定された各原信号yktを出力する。以下、各部の処理内容を説明する。<周波数領域変換部11>まず、周波数領域変換部11は、マイクL、Rで収音した時間領域の観測信号xt=[xLt,xRt]Tを入力とし、これを短時間フーリエ変換等により周波数領域の観測信号xf,t=[xf,t,L,xf,t,R]Tに変換しs1、パラメタ推定部12及び分離信号生成部13に出力する。<パラメタ推定部12>次に、パラメタ推定部12は、観測信号xf,tから、音源分離のために必要なパラメタθを推定するs2。なおθ={δf,k,σf2,sf,t,k,pk|θ}であり、δf,kは原信号の、二個のマイクへの到達時間差を表す。σf2は雑音のパワースペクトルを、sf,t,kは音源スペクトルを、pk|θは音源存在確率混合信号中の音源kの寄与率を表す。本実施形態では、上記パラメタを推定するために、以下の2つの仮定を用いる。仮定1は、雑音nが、平均0、共分散行列σf2Vfの正規分布に従う定常雑音でモデル化できるという仮定である。ここでσf2は周波数fにおける雑音のパワーであり、Vfは例えば拡散性雑音の場合で与えられる非特許文献1参照。ここでcは原信号の速度音速等、Dは二個のマイクの間隔である。仮定2は、原信号がスパースな信号つまり、成分のうち0でないもの非零の成分がまばらである信号、言い換えると、ほとんどのまたは多くの成分が0である信号であるという仮定である。すなわち、ある時間周波数f,tにおいて、たかだか1つの原信号のみが支配的であると仮定する。仮定2によると、式2はと表記できる非特許文献1参照。上記仮定1、仮定2に基づけば、到達時間差δf,kとなる方向に存在する音源kから発せられる原信号が到来してxf,tが観測される尤度は、で与えられる非特許文献1参照。Hはエルミート転置を表す。ここでbf,kは式3で表される。これより、到達時間差δf,k、音源スペクトルsf,t,k及び雑音のパワースペクトルσf2は、対数尤度関数を最大化するパラメタとして最尤推定により求める非特許文献1参照。但し、上記において、pk|δf,k,σf2,sf,t,kは音源存在確率を表し、混合信号中の音源kの寄与率であり、である非特許文献1参照。具体的には、期待値最大化法以下EMアルゴリズムともいうを適用し、Q関数を最大とするパラメタθ={δf,k,σf2,sf,t,k,pk|θ}を、以下のEステップ及びMステップの繰り返しにより求める非特許文献1参照。ここでmf,t,kは後述する式11で、pxf,t|k,θ’は前述の式6で与えられ、θ’は、1回前の繰り返しで得られているパラメタを意味する。この繰り返し計算を、パラメタ推定部12にて行なう。以下、パラメタ推定部12の処理の詳細を説明する。図3はEステップ計算部121の機能ブロック図を、図4はMステップ計算部122の機能ブロック図を、図5はパラメタ推定部12の処理フローを示す。まず、パラメタ推定部12は、各パラメタを初期化するs20。パラメタのうちσf2、δf,k、及びpk|θを初期化する。例えば、σf2=|xf,t=0,L|2,δf,k=D/ccosαkcは原信号の速度音速等、Dはマイク間隔,αkは音源kの方向の初期値-π/2~π/2の間の適当な値,pk|θ=1/Kとする。さらに、初期化したパラメタδf,kとxf,t=0とを用いて、式3及び後述する式19に基づき、sf,t,kを初期化する。更新回数n=0とする。なお、最大更新回数N及び収束判定閾値Δは、当該装置の設計者や利用者等により予め設定されているものとする。Eステップ計算部121は事後確率推定部1211とQ関数計算部1212とを含む図3参照。Mステップ計算部122は時間差推定部1221と音源スペクトル推定部1222と雑音パワー推定部1223と音源存在確率推定部1224とを含み、時間差推定部1221はさらに逆正接計算部12211と時間補正部12212とを備える図4参照。事後確率推定部1211とQ関数計算部1212とにおける処理を併せてEステップと呼び、時間差推定部1221と音源スペクトル推定部1222と雑音パワー推定部1223と音源存在確率推定部1224とにおける処理を併せてMステップと呼ぶ。事後確率推定部1211Eステップ計算部121の事後確率推定部1211は、観測信号xf,tと、一回前の繰り返しで得られているパラメタθ’但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の事後確率推定においては、前述の初期化したパラメタとを入力とし、これらの値を用いて事後確率mf,t,k=pk|xf,t,θ’を以下の式11により求めs22、Q関数計算部1212とMステップ計算部122とに出力する。なお、pxf,t|k,θ’は式6により与えられる。なお、事後確率mf,t,kは観測信号xf,tが音源kに帰属する事後確率を表す。逆正接計算部12211Mステップ計算部122の時間差推定部1221の逆正接計算部12211は、観測信号xf,tと、事後確率mf,t,kと、一回前の繰り返しで得られているパラメタθ’より詳しく説明するとθ’のうちの音源スペクトルsf,t,kである。但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の逆正接計算においては、前述の初期化したパラメタとを入力とし、これらの値を用いて、到達時間差δf,kより詳しく言うと到達時間差δf,kに2πfを乗じた値2πfδf,kを以下の式13により推定しs25、時間補正部12212に出力する。ψsk但し添え字skはskを表す及びψξk但し添え字ξkはξkを表すは、それぞれsf,t,k及びξf,t,kの位相を表す。なお、式13の導出については後述する。なお、従来技術では式1に基づきδkを離散全探索によって推定するため多くの計算コストを要していたが、本実施形態では式13に基づきδf,kを算出するため全探索を要せず計算コストを小さくできる。また、式13に基づき周波数f毎に到達時間差δf,kを求める点が従来技術とは異なる。時間補正部12212Mステップ計算部122の時間差推定部1221の時間補正部12212は、観測信号xf,tと、事後確率mf,t,kと、一回前の繰り返しで得られているパラメタθ’但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の時間補正においては、前述の初期化したパラメタとを入力とし、これらの値を用いて、式13にて推定した到達時間差δf,kが内包する±πの不定性を補正するs26。この補正は、式13の左辺2πfδf,kが-πからπの値を取るのに対し、式13の右辺の逆正接が-π/2からπ/2の値しか返すことができないため、2πfδf,kが-π~-π/2及びπ/2~πの範囲を取る場合の値を正しく求めるために必要である。式13で得られた値を2πfδ’f,kと記載すると、補正は以下のように行なう。ここでξf,t,k、φfはそれぞれ式14、15で与えられ、ψsk但し添え字skはskを表す、ψξk但し添え字ξkはξkを表すはそれぞれsf,t,k及びξf,t,kの位相を表す。なお、式17、18の導出については後述する。さらに、時間補正部12212は、補正した値2πfδf,kを2πfで除算し、到達時間差δf,kを求め、音源スペクトル推定部1222と雑音パワー推定部1223とQ関数計算部1212とに出力する。音源スペクトル推定部1222Mステップ計算部122の音源スペクトル推定部1222は、到達時間差δf,kと観測信号xf,tを入力とし、これらの値を用いて、音源スペクトルsf,t,kを以下の式19により推定しs27、雑音パワー推定部1223とQ関数計算部1212とに出力する。ここで、bf,k、Vfはそれぞれ式3、4で表される。雑音パワー推定部1223Mステップ計算部122の雑音パワー推定部1223は、到達時間差δf,kと音源スペクトルsf,t,kと観測信号xf,tと事後確率mf,t,kとを入力とし、これらの値を用いて、雑音のパワースペクトルσf2を以下の式20により推定しs28、Q関数計算部1212に出力する。なお、Tは時間フレームの総数である。音源存在確率推定部1224Mステップ計算部122の音源存在確率推定部1224は、事後確率mf,t,kを入力とし、音源存在確率を以下の式21により推定しs29、Q関数計算部1212に出力する。Q関数計算部1212Eステップ計算部121のQ関数計算部1212は、観測信号xf,tと、パラメタθ={δf,k,σf2,sf,t,k,pk|θ}と、事後確率mf,t,kと、を入力とし、これらの値を用いてQ関数を上述の式10により求めるs30。s20~s30までの処理を終えると、パラメタ推定部12は、条件1Q関数の値の変化量|Qθ|θn-1-Qθ|θn|が所定の収束判定閾値Δより小さくなるか、または、条件2更新回数nが所定の最大更新回数N例えばN=20以上か否かを判定するs31。パラメタ推定部12は、条件1、条件2の何れかを満たしたときは、パラメタ推定部12はその時点で取得している最新の事後確率mf,t,kと到達時間差δf,kを分離信号生成部13に出力する。パラメタ推定部12は、条件1、条件2の何れも満たさないときは、EステップとMステップを繰り返すs31、s21。なお、図示しない記憶部にパラメタθとQ関数の値Qθ|θnとを記憶しておき、次の繰り返しの際に用いる。<分離信号生成部13>分離信号生成部13は、事後確率mf,t,kと到達時間差δf,kと観測信号xf,tとを入力とし、以下の式22により、分離信号yf,t,kを生成しs3、時間領域変換部14へ出力する。なお、音源スペクトルsf,t,kは、音源の発する原信号のスペクトルを推定したものであるが、この音源スペクトルを単純に時間領域の信号に変換した場合には、他の音源の発する原信号が残ることがある。上述の式22によって、音源kの発する原信号のみを抽出、分離することができる。<時間領域変換部14>時間領域変換部14は、分離信号yf,t,kを入力とし、周波数領域変換部11において行った周波数領域変換方法に対応する時間領域変換方法例えば短時間フーリエ逆変換で、分離信号yf,t,kを時間領域の分離信号yktに変換しs4、音源分離装置1の出力値として出力する。<本実施形態のポイント>以下、本実施形態のポイントを説明し、式13、17、18の導出方法を説明する。本実施形態では、性質1到達時間差δf,kがRチャネルとLチャネルの位相差に影響を与える値であること性質2位相差が周期的な値を取る量であることという2つの性質を利用して、到着時間差δf,kを推定する。ここで性質1は、式2にてノイズnf,tが十分に小さい場合を考えれば明らかである。性質2は、ある2つの位相を表わす量の差Θが、-π≦Θ<πの範囲の値だけではなく、Θ±2πMMは任意の整数という周期的な不定性を内包する値を取る性質を持つことを意味する。性質2について、さらに以下にて説明を行なう。式6の右辺における、expのカッコの中のベクトル及び行列x,b,Vを、それぞれの成分で表して整理すると、となるCはδf,kに依らない定数。式32は、位相や角度の分布のように周期的な値を取る変数に対する分布であるVonMises分布と同じ形をしていることが分かる参考文献1参照。[参考文献1]C.M.ビショップ著,元田ら訳,パターン認識と機械学習上,シュプリンガージャパン,2006.ここで-π<x≦π、μは分布の平均-π<μ≦π、κ>0は分布の集中度パラメタ正規分布での1/分散に相当、I0xは0次の第1種ベッセル関数である。すなわち、式33の変数xが式32のψξk-ψskに対応し、式33の平均μが式32の2πfδf,kに対応し、式33の集中度κが式32の|ξf,t,k||sf,t,k|/σf21-φf2に対応する。説明をより直感的にするため、雑音のパラメタがφf=0である場合を考えるこれは、雑音が式4に示す拡散性雑音ではなく、分散σf2のガウス雑音が観測xL,xRにそれぞれ乗ることを意味する。このとき、式14によりξf,t,k=xf,t,Rとなるため、式32において、ψξkはψxf,t,R但し、添え字xf,t,Rはxf,t,Rを表し、ψxf,t,RはマイクRの観測信号xf,t,Rの位相を表すとなる。また、式3、5により但し式5においてnf,t=0、sf,t,k=xf,t,Lとなるため、式32においてψskはψxf,t,L但し、添え字xf,t,Lはxf,t,Lを表し、ψxf,t,LはマイクLの観測信号xf,t,Lの位相を表すとなる。また、前述の通りξf,t,k=xf,t,Rであり、式3、5よりxf,t,R=ej2πfδk,fsf,t,k但し、eの添え字δk,fはδk,fを表すとなるため、式32において、|ξf,t,k|=|ej2πfδk,fsf,t,k|=|sf,t,k|となる。よって、式32はとなる。これとVonMises分布との解釈を合わせると、1RチャネルとLチャネルの位相差ψxf,t,R-ψxf,t,Lという周期的な値を取る変数の分布が、平均μ=2πfδf,kを取る、2VonMises分布の集中度κが、SN比|sf,k,t|2/σf2に対応するようになる。すなわち、SN比が低い条件が悪いと、RチャネルとLチャネルの位相差の値の集中度が下がる=分散が大きくなる。これは、SN比が低い条件においては位相差の測定値がばらつく現象と対応する、の2点が言える。本実施形態では、性質1、性質2を利用しており、実施の手続きとしては、RチャネルとLチャネルの位相差に関する量が、周期的な値を取る変数に対する分布であるVonMises分布で表現できることに着目し、VonMises分布のパラメタに対する最尤推定により、分布の平均値μ=2πfδf,kを推定することで、信号到達時間差パラメタδf,kを推定する。式31のpxf,t|k,θをQ関数の式10に代入し、を解くことで、式13が得られる。また時間補正部12212における関数Fの式17、18は、Q関数の2階微分である。時間補正部12212では、式13で得られた解δ’f,kがQ関数の極大値極小値のどちらを与えるかを、F2πfδ’f,kの値の正負にて調べ、δ’f,kが極小値を与える場合に、上述のπまたは-πの補正を行なう。<効果>従来法においては、時間差推定部において解析的な更新式が与えられていなかったため、多くの計算コストを要する全探索操作が必要であった。よって、時間差推定部の計算コストを削減し、高速な音源分離手段を提供することが課題である。本実施形態では、到着時間差δf,kが、マイクRとマイクLの位相差に影響を与える値であることと、位相差が周期的な値を取る性質を持つことに着目し、到着時間差δf,kを推定する。これにより、従来必要であった全探索操作が不要となるため、高速な音源分離手段を提供することが可能となる。<シミュレーション結果>発明の効果を示すため実験を行なった。図6に示す部屋において、音源数は2つ70度と150度又は3つ30,70,150度とした。音源は、英語話者音声を用い、音源数2及び3の場合それぞれにおいて10通りの音源組合せにて実験を行なった。雑音としては、平均0、共分散行列σf2VfVfは式4により与えられるのガウスノイズをSN比約25dBにて重畳した。部屋の残響時間は130ms、サンプリング周波数は8kHz,短時間フーリエ変換の窓長及びシフト長はそれぞれ64ms,16msとした。従来法では到達時間差δkを推定する際に、音源位置Θ図6参照を0度から180度まで1度きざみで変化させ、それに対応する181種類の到達時間差δkの値について全探索を行なった。図7に、SIR信号対雑音比、雑音には他話者音声含む、SDR信号対歪み比、及びパソコンIntelXeon登録商標X5650267GHz×2CPUにおける計算時間を示す。図7Aは音源が二つの場合70度と150度を、図7Bは音源が三つの場合30度,70度,150度を示す。SIRとSDRは大きな値であるほど良い性能であることを示す。数字は、10通りの音源組合せの平均値である。図7A、Bに示す通り、本実施形態は、1/10程度の計算時間で、従来法とほぼ同程度の分離性能を達成できることが見てとれる。<その他の変形例>本実施形態のポイントは上述の通り、到達時間差の推定方法である。従って、他の処理やパラメタの推定方法については、上記の実施形態に限定されるものではなく、他の従来技術を用いてもよい。本実施形態では、各部間で直接データを受け渡しているが、図示しない記憶部を介して、各データを読み書きしてもよい。また、本実施形態では、雑音パワー推定部1223は音源スペクトルsf,t,kを音源スペクトル推定部1222から取得しているが、到達時間差δf,kと観測信号xf,tとを用いて式19に基づき雑音パワー推定部1223で計算する構成としてもよい。本実施形態では、拡散性雑音の場合を想定しているが、他の特性を持つ雑音であってもよい。その場合、雑音の特性に応じて、式4や式15のsinc2πfD/cを適宜変更すればよい。分離信号生成部13は、事後確率mf,t,kと音源スペクトルsf,t,kとを入力とし、式22に代えて、以下の式により、分離信号yf,t,kを生成してもよい。この場合、パラメタ推定部12は、Q関数の値の変化量|Qθ|θn-1-Qθ|θn|が所定の収束判定閾値Δより小さくなったとき、または、更新回数nが所定の最大更新回数N以上になったときs31に取得している最新の事後確率mf,t,kと音源スペクトルsf,t,kを分離信号生成部13に出力する。このような構成により、式22と同様に分離信号yf,t,kを生成することができる式19参照。Mステップの計算順序は、本実施形態の計算順序に限らない。例えば、時間差推定部1221と音源スペクトル推定部1222の計算順序はどちらを先に行ってもよい。図8は音源スペクトル推定部1222の音源スペクトル推定処理s27を行った後に、時間差推定部1221の到達時間差推定処理s25、s26を行う場合の機能ブロック図を示す。以下、図8を用いて説明する。Mステップ計算部122の音源スペクトル推定部1222は、一回前の繰り返しで得られているパラメタθ’より詳しく説明するとθ’のうちの到達時間差δf,kである。但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の音源スペクトル計算においては、前述の初期化したパラメタと観測信号xf,tを入力とし、これらの値を用いて、音源スペクトルsf,t,kを式19により推定し、雑音パワー推定部1223とQ関数計算部1212と時間差推定部1221と時間補正部12212とに出力する。Mステップ計算部122の時間差推定部1221の逆正接計算部12211は、観測信号xf,tと、事後確率mf,t,kと、音源スペクトルsf,t,kとを入力とし、これらの値を用いて、到達時間差δf,kを式13により推定し、時間補正部12212に出力する。Mステップ計算部122の時間差推定部1221の時間補正部12212は、観測信号xf,tと、事後確率mf,t,kと、音源スペクトルsf,t,kと、一回前の繰り返しで得られているパラメタθ’より詳しく説明するとθ’のうちの雑音のパワースペクトルσf2である。但し、一回前の繰り返しで得られているパラメタθ’が存在しない場合、つまり、一回目の時間補正においては、前述の初期化したパラメタとを入力とし、これらの値を用いて、式13にて推定した到達時間差δf,kが内包する±πの不定性を式16に基づき補正する。さらに、時間補正部12212は、補正した値2πfδf,kを2πfで除算し、到達時間差δf,kを求め、雑音パワー推定部1223とQ関数計算部1212とに出力する。上述の構成とすることで、第一実施形態と同様の効果を得ることができる。なお、同様に、一回前の繰り返しで得られているパラメタθ’のうちの到達時間差δf,kに基づき雑音パワー推定部1223の雑音パワー推定処理s28を行った後に、時間差推定部1221の到達時間差推定処理s25、s26を行ってもよい。本発明は上記の実施形態及び変形例に限定されるものではない。例えば、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。<プログラム及び記録媒体>上述した音源分離装置は、コンピュータにより機能させることもできる。この場合はコンピュータに、目的とする装置各種実施例で図に示した機能構成をもつ装置として機能させるためのプログラム、またはその処理手順各実施例で示したものの各過程をコンピュータに実行させるためのプログラムを、CD-ROM、磁気ディスク、半導体記憶装置などの記録媒体から、あるいは通信回線を介してそのコンピュータ内にダウンロードし、そのプログラムを実行させればよい。"