(57)【特許請求の範囲】
【請求項１】
複数の音源信号それぞれの音源パワー時系列全体の状態を表す音源パワーパラメータの事前確率密度関数と、その音源パワーパラメータが与えられた場合の各音源信号の各時間周波数点における音源パワー特徴量のモデルと、を記憶した音源モデル記憶部と、
上記複数の音源信号を複数のマイクロホンで収音した時間領域信号を時間周波数領域信号に変換した観測信号を入力として、各時間周波数点における音源位置特徴量と音源パワー特徴量を抽出する特徴抽出部と、
上記音源パワー特徴量と、上記観測信号が得られた下での占有的な音源の事後確率密度関数である音源占有度と、上記音源パワーパラメータの事前確率密度関数と上記各音源信号の音源パワー特徴量のモデルと、上記音源パワー特徴量のモデルと上記事前確率密度関数のそれぞれの挙動を制御するパラメータである音源モデルパラメータと、を入力として上記各音源の音源パワーパラメータを更新する音源パワーパラメータ更新部と、
上記音源位置特徴量と音源占有度を入力として、上記各音源の音源位置パラメータを更新する音源位置パラメータ更新部と、
上記音源パワー特徴量と上記音源パワーパラメータと上記音源占有度と、上記音源モデル記憶部に記憶された音源パワーパラメータの事前確率密度関数と音源パワー特徴量のモデルとを入力として上記音源モデルパラメータを更新する音源モデルパラメータ更新部と、
上記音源位置特徴量と音源パワー特徴量と各音源の更新された音源パワーパラメータと音源位置パラメータと、上記音源モデルパラメータと、上記音源モデル記憶部に記憶された音源パワーパラメータの事前確率密度関数と音源パワー特徴量のモデルとを入力として上記各音源の音源占有度を更新する音源占有度更新部と、
を具備する音源パラメータ推定装置。
【請求項２】
請求項１に記載した音源パラメータ推定装置において、
上記音源パワー特徴量の時系列は隠れマルコフモデルに従い、上記各音源信号のマルコフモデルの状態が既知の条件下における上記音源パワーパラメータの事後確率密度関数は共分散行列が対角行列で表されるガウス分布でモデル化されており、当該ガウス分布の平均と共分散行列の対角要素を上記音源モデルパラメータに含むことを特徴とする音源パラメータ推定装置。
【請求項３】
請求項１又は２に記載した音源パラメータ推定装置と、
その音源パラメータ推定装置が出力する音源パワー特徴量と、当該音源パラメータ推定装置が更新した音源占有度と音源パワー特徴量と音源パワーパラメータと音源モデルパラメータと、上記各音源信号の各時間周波数点における音源パワー特徴量のモデルと、を入力として複数の音源のそれぞれの音源分離信号を、最小自乗誤差推定により求める音源分離部と、
を具備する音源分離装置。
【請求項４】
複数の音源信号を複数のマイクロホンで収音した時間領域信号を時間周波数領域信号に変換した観測信号を入力として、各時間周波数点における音源位置特徴量と音源パワー特徴量を抽出する特徴抽出過程と、
上記音源パワー特徴量と、上記観測信号が得られた下での占有的な音源の事後確率密度関数である音源占有度と、音源モデル記憶部に記憶された音源パワーパラメータの事前確率密度関数と上記各音源信号の音源パワー特徴量のモデルと、上記音源パワー特徴量のモデルと上記事前確率密度関数のそれぞれの挙動を制御するパラメータである音源モデルパラメータと、を入力として上記各音源の音源パワーパラメータを更新する音源パワーパラメータ更新過程と、
上記音源位置特徴量と音源占有度を入力として、上記各音源の音源位置パラメータを更新する音源位置パラメータ更新過程と、
上記音源パワー特徴量と上記音源パワーパラメータと上記音源占有度と、上記音源モデル記憶部に記憶された音源パワーパラメータの事前確率密度関数と音源パワー特徴量のモデルとを入力として上記音源モデルパラメータを更新する音源モデルパラメータ更新過程と、
上記音源位置特徴量と音源パワー特徴量と各音源の更新された音源パワーパラメータと音源位置パラメータと、上記音源モデルパラメータと、上記音源モデル記憶部に記憶された音源パワーパラメータの事前確率密度関数と音源パワー特徴量のモデルとを入力として上記各音源の音源占有度を更新する音源占有度更新過程と、
を備える音源パラメータ推定方法。
【請求項５】
請求項４に記載した音源パラメータ推定方法において、
上記音源パワー特徴量の時系列は隠れマルコフモデルに従い、上記各音源信号のマルコフモデルの状態が既知の条件下における上記音源パワーパラメータの事後確率密度関数は共分散行列が対角行列で表されるガウス分布でモデル化されており、当該ガウス分布の平均と共分散行列の対角要素を上記音源モデルパラメータに含むことを特徴とする音源パラメータ推定方法。
【請求項６】
請求項４又は５に記載した音源パラメータ推定方法で抽出した音源パワー特徴量と、当該音源パラメータ推定方法で更新した音源占有度と音源パワー特徴量と音源パワーパラメータと音源モデルパラメータと、上記各音源信号の各時間周波数点における音源パワー特徴量のモデルと、を入力として複数の音源のそれぞれの音源分離信号を、最小自乗誤差推定により求める音源分離過程と、
を備える音源分離方法。
【請求項７】
請求項１乃至３の何れかに記載した音源パラメータ推定装置又は音源分離装置としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、複数の音源が同時に生成した音響信号が混ざって複数のマイクロホンで収音された観測信号から、各音源の音源パラメータを推定する音源パラメータ推定装置と、その音源パラメータに基づいて各音源を分離する音源分離装置と、それらの方法とプログラムに関する。
【背景技術】
【０００２】
図７に、非特許文献１に開示された従来の音源パラメータ推定装置９００の機能構成を示す。音源パラメータ推定装置９００は、音源モデル記憶部９２０と、特徴抽出部９１０と、音源パワーパラメータ更新部９３０と、音源位置パラメータ更新部９４０と、音源占有度更新部９５０と、を具備する。
【０００３】
音源モデル記憶部９２０は、複数の音源信号それぞれの音源パワー時系列全体の状態を表す音源パワーパラメータの事前確率密度関数p(q(l))（（l）はl番目の音源）と、その音源パワーパラメータq(l)が与えられた場合の各音源信号の各時間周波数点(n,k)における事後確率密度関数である音源パワー特徴量のモデルf=βq（l）,n,k(S)とを記憶する。特徴抽出部９１０は、複数の音源信号を複数（m個）のマイクロホンで収音した時間領域の信号を時間周波数領域信号に変換した観測信号x(m)n,kを入力として、各時間周波数点における音源位置特徴量An,kと音源パワー特徴量Xn,kを抽出する。
【０００４】
音源パワーパラメータ更新部９３０は、Ns個の音源ごとに音源パワー特徴量Xn,kと観測信号が得られた下での占有的な音源の事後確率密度関数である音源占有度^M(l)n,kと、音源パワー特徴量のモデルf=βq（l）,n,k(S)と音源パワーパラメータの事前確率密度関数p(q(l))とを入力として、音源パワーパラメータ^q(l)を更新する。音源位置パラメータ更新部９４０は、音源位置特徴量An,kと音源占有度^M(l)n,kを入力として、各音源の音源位置パラメータ^φ(l)を更新する。
【０００５】
音源占有度更新部９５０は、音源位置特徴量An,kと音源パワー特徴量Xn,kと各音源の更新された音源パワーパラメータ^q(l)と音源位置パラメータ^φ(l)と、音源パワー特徴量のモデルf=βq（l）,n,k(S)と音源パワーパラメータの事前確率密度関数p(q(l))とを基に各音源の音源占有度^M(l)n,kを更新する。
【０００６】
音源パラメータ推定装置９００の音源パラメータ推定技術を用いた音源分離装置は、音源パラメータ推定装置９００が出力する音源パワー特徴量と、更新した音源占有度と音源パワーパラメータと、各音源信号の各時間周波数点における音源パワー特徴量のモデルとを入力として複数の音源のそれぞれの音源分離信号を、最小自乗誤差推定により求める音源分離部を更に備える。
【先行技術文献】
【非特許文献】
【０００７】
【非特許文献１】Tomohiro Nakatani, Shoko Araki, Takuya Yoshioka, Masakiyo Fujimoto, “Multichannel source separation based on source location cue with log-spectral shaping by hidden Markov source model,” Proc. of Interspeech-2010,pp. 2766-2769, Sep., 2010.
【非特許文献２】中谷智広、荒木章子、吉岡卓也、藤本雅清、“DOAクラスタリングと音声の対数スペクトルHMMに基づく音源分離”日本音響学会2010年秋季研究発表会講演論文集、pp.577-580, 9月, 2010年.
【発明の概要】
【発明が解決しようとする課題】
【０００８】
従来の音源パラメータ推定装置は、音源モデル記憶部が記憶する音源パワーパラメータの事前確率密度関数と音源パワー特徴量のモデルのそれぞれの挙動を制御するパラメータである音源モデルパラメータが事前に与えられなければ音源パラメータの推定を行うことができなかった。
【０００９】
この発明はこの課題に鑑みてなされたものであり、音源モデルパラメータが事前に与えられていない場合でも、音源パラメータの一部としてその他の音源パラメータと一緒に音源モデルパラメータをも推定できる音源パラメータ推定装置と音源分離装置と、それらの方法とプログラムを提供することを目的とする。
【課題を解決するための手段】
【００１０】
この発明の音源パラメータ推定装置は、音源モデル記憶部と、特徴抽出部と、音源パワーパラメータ更新部と、音源位置パラメータ更新部と、音源モデルパラメータ更新部と、音源占有度更新部と、を具備する。音源モデル記憶部は、複数の音源信号それぞれの音源パワー時系列全体の状態を表す音源パワーパラメータの事前確率密度関数と、その音源パワーパラメータが与えられた場合の各音源信号の各時間周波数点における音源パワー特徴量のモデルとを記憶する。特徴抽出部は、複数の音源信号を複数のマイクロホンで収音した時間領域信号を時間周波数領域信号に変換した観測信号を入力として、各時間周波数点における音源位置特徴量と音源パワー特徴量を抽出する。音源パワーパラメータ更新部は、音源パワー特徴量と、観測信号が得られた下での占有的な音源の事後確率密度関数である音源占有度と、音源パワーパラメータの事前確率密度関数と各音源信号の音源パワー特徴量のモデルと、音源パワー特徴量のモデルと音源パワーパラメータの事前確率密度関数のそれぞれの挙動を制御するパラメータである音源モデルパラメータと、を入力として各音源の音源パワーパラメータを更新する。音源位置パラメータ更新部は、音源位置特徴量と音源占有度を入力として、各音源の音源位置パラメータを更新する。音源モデルパラメータ更新部は、音源パワー特徴量と音源パワーパラメータと音源占有度と、音源モデル記憶部に記憶された音源パワーパラメータの事前確率密度関数と音源パワー特徴量のモデルとを入力として音源モデルパラメータを更新する。音源占有度更新部は、音源位置特徴量と音源パワー特徴量と各音源の更新された音源パワーパラメータと音源位置パラメータと、音源モデルパラメータと、音源モデル記憶部に記憶された音源パワーパラメータの事前確率密度関数と音源パワー特徴量のモデルとを入力として各音源の音源占有度を更新する。
【発明の効果】
【００１１】
この発明の音源パラメータ推定装置によれば、音源モデルパラメータが事前に与えられていない場合でも、音源パラメータの一部としてその他の音源パラメータと一緒に音源モデルパラメータをも推定することができるので、事前に統計的性質を与えることができない多様な音源に対しても最適な音源パラメータを与えることができる。
【００１２】
その音源パラメータ推定装置を用いたこの発明の音源分離装置は、誤差の少ない音源分離信号を出力することが可能である。
【図面の簡単な説明】
【００１３】
【図１】この発明の音源パラメータ推定装置１００の機能構成例を示す図。
【図２】音源パラメータ推定装置１００の動作フローを示す図。
【図３】この発明の音源分離装置２００の機能構成例を示す図。
【図４】音源分離装置２００の動作フローを示す図。
【図５】実験結果を示す図であり、(ａ)は混合前の音声から推定した結果、（ｂ）は混合音から推定した結果を示す図である。
【図６】実験結果を示す図であり、分離後の信号のケプストラム歪みと混合音の数との関係を示す、（ａ）は混合条件１の混合音に対して処理した結果、（ｂ）と（ｃ）は残響のある別の環境で収録した混合音（混合条件２と３）に対して処理した結果を示す図である。
【図７】従来の音源パラメータ推定装置９００の機能構成を示す図。
【発明を実施するための形態】
【００１４】
以下、この発明の実施の形態を図面を参照して説明する。複数の図面中同一のものには同じ参照符号を付し、説明は繰り返さない。実施例の説明の前にこの発明の基本的な考えについて説明する。
【００１５】
〔この発明の基本的な考え〕
この発明では、予め音源モデルパラメータが与えられていなくとも音源パラメータの推定が行えるようにするために、ｌ番目の音源信号の音源パワー特徴量の時系列全体を{S(l)n,k}と表した時に、その同時確率密度関数を音源モデルパラメータψ(l)を含めてモデル化すると共に、複数の音が混在した信号の音源パワー特徴量Xn,kが与えられたときに、各音源の音源占有度、音源パワー特徴量のモデルと音源パワーパラメータの事前確率密度関数、および音源パワーパラメータに基づいて各音源の音源モデルパラメータを更新する音源モデルパラメータ更新部を有する点で新しい。
【００１６】
最初に説明に用いる記号について説明する。観測信号には、Ns個の音源信号が重畳しており、その音源信号をNm本のマイクロホンで収音する。m番目のマイクロホンから収音した収音信号を短時間フーリエ変換等を用いて周波数領域の信号に変換した観測信号をX(m)n,kと表記する。nはn番目の時間つまりフレーム番号、kはk番目の周波数つまりビン番号であり、n番目の時間及びk番目の周波数に対応する時間周波数点を参照する場合に、時間周波数点（n,k）と表記する。なお、記号＾の位置や添え字の表記とその位置は、式中の表記が正しい。
【００１７】
この発明では、ｌ番目の音源信号の音源パワー時系列全体{S(l)n,k}の同時確率密度関数を次式に示すようにモデル化する。
【００１８】
【数１】
0005438704.tif 000002
【００１９】
ここで、q(l)はｌ番目の音源の音源パワー時系列全体の状態を表す音源パワーパラメータを表す。以下では全ての音源のq(l)をまとめてq=[q(1),…，q(Ns)]とも表記する。ψ(l)はl番目の音源の音源モデルパラメータ全体を現す。全ての音源のψ(l)をまとめてψ=[ψ(1),…，ψ(Ns)]とも表記する。
【００２０】
また、βq(l),n,k,ψ(l)(S)は音源パワー特徴量のモデルであり、音源パワーパラメータq(l)と音源モデルパラメータψ(l)が与えられた下で各時間周波数点(n,k)の音源信号の音源パワーがS(l)n,kとなる確率密度関数である（式（３））。式（１）の総和演算は、q(l)が離散値ではなく連続値をとる場合にはq(l)に関する積分演算に置き換えて表現されるものとする。式（２）において、音源の状態が既知のもとでは、異なる時間周波数点における音源パワーS(l)n,kは相互に独立であるという仮定を導入している。
【００２１】
また、この発明では式（４）に示すように、各時間周波数点(n,k)において最も大きなエネルギーを持つ音源信号（以下、占有的な音源信号と称する）の音源パワーS(l)n,kは、観測信号の音源パワーと一致すると仮定する。
【００２２】
【数２】
0005438704.tif 000003
【００２３】
また、占有的ではない音源ｌに関しては、S(l)n,k≦Xn,kの関係を持つと仮定する。すると、各音源信号の状態が既知の条件の下で、観測信号の音源パワーXn,kの事後確率密度関数は次のように表現できることが知られている（参考文献：S. J. Rennie, J.R. Hershey, and P. A.01sen, “Hierarchical variational loopy belief propagation for multi-talker speech recognition,” Proc. ASRU-2009, pp. 176-181 2009.）。
【００２４】
【数３】
0005438704.tif 000004
【００２５】
この発明では、更に上式は次のように分解可能であると仮定している。
【００２６】
【数４】
0005438704.tif 000005
【００２７】
また、この発明では音源位置特徴量から音源位置パラメータφ＾(l)を推定するため、音源位置特徴量のモデルp(An,k;φ)を導入する。音源位置特徴量のモデルp(An,k;φ)は、各音源信号のエネルギーは異なる時間周波数点にわたり疎に分布していると仮定し、その時間周波数点において占有的な音源の音源位置のみに依存して決まると仮定する。
【００２８】
一般的に、全ての音源の音源位置パラメータφ(l)をまとめてφ=[ψ(1),…，ψ(Ns)]と表すと、音源位置特徴量のモデルp(An,k;φ)、つまり観測信号の音源位置特徴量の確率密度関数は、混合分布として式（８）に示すように展開することができる。
【００２９】
【数５】
0005438704.tif 000006
【００３０】
式（８）において、Zn,kは時間周波数点（n,k）において占有的な音源の番号を表す確率変数であり、Zn,k=ｌは、ｌ番目の音源が占有的な音源である場合を示す。また、p(Zn,k=ｌ)は、l番目の音源が時間周波数点（n,k）において占有的な音源になる事前確率密度関数を表している。更に、以降の説明では次の表記を用いることにする。
【００３１】
【数６】
0005438704.tif 000007
【００３２】
γφ(l),n,k(A)は、時間周波数点（n,k）において占有的な音源の番号がlの場合に、音源位置特徴量An,kが得られる確率密度関数を表す。これは、l番目の音源の音源位置パラメータφ(l)のみに依存するものとする。具体的なγφ(l),n,k(A)やφ(l)の定義については後述する。
【００３３】
式（８）のもと、γφ(l),n,k(A)が定義されている場合、音源位置パラメータφ(l)と占有的な音源の番号に関する事前確率密度関数p(Zn,k=l)が与えられれば、音源位置特徴量のモデルp(An,k;φ)は一意に定めることができる。逆に、音源位置特徴量An,kが観測された場合に、最尤推定などの方法に従い、音源位置パラメータと占有的な音源の番号に関する事前確率密度関数p(Zn,k=l)やその事後確率密度関数を推定することができる。
【００３４】
以上の定義に従うと、完全データの確率密度関数は式（１０）に示すように導出される。
【００３５】
【数７】
0005438704.tif 000008
【００３６】
式（１０）において、qが音源パワーパラメータ、ψが音源モデルパラメータ、φが音源位置パラメータであり、これらのパラメータがパラメータ推定の対象である。この発明では、次の対数尤度関数を最大化する値として、音源パワーパラメータqと音源モデルパラメータψと音源位置パラメータφを推定する。
【００３７】
【数８】
0005438704.tif 000009
【００３８】
式（１２）で、確率変数Zn,kは隠れ変数として扱われる。隠れ変数を含む対数尤度関数の最大化には、例えば、期待値最大化アルゴリズムなどを用いることができる。期待値最大化アルゴリズムでは、音源パワーパラメータの推定値^qと音源位置パラメータの推定値^φと音源モデルパラメータの推定値^ψに基づき、観測信号が得られた下での占有的な音源の番号の事後確率密度関数^M(l)n,k=p(Zn,k｜An,k,Xn,k,^q;^φ,^ψ)をも同時に推定する必要がある。この発明では、この関数の値を音源占有度と称し、この値も音源パラメータに含めて考える。
【００３９】
以上説明した考えで、音源パワー特徴量のモデルβq(l),n,k,ψ(l)(S)と、音源位置特徴量のモデルp(An,k;φ)の両者を考慮しながら最適な音源パラメータを推定することで音源パラメータの推定誤差を減らすことができる。また、音源位置特徴量のモデルp(An,k;φ)（式（８））と、音源パワー特徴量のモデル（式（７））に、占有的な音源の番号を表す変数Zn,kを共有化することで、２つの特徴量を考慮しながら音源パラメータ推定の計算を簡単にすることができる。
【００４０】
以上、述べたように、この発明の音源パラメータ推定方法によれば、ｌ番目の音源信号の音源パワー時系列全体{S(l)n,k}の同時確率密度関数を音源モデルパラメータψ(l)を含めてモデル化すると共に、複数の音が混在した信号の音源パワー特徴量Xn,kが与えられたときに、各音源の音源占有度、音源パワー特徴量のモデルと音源パワーパラメータの事前確率密度関数、および音源パワーパラメータに基づいて各音源の音源モデルパラメータを更新する音源モデルパラメータ更新部を有することで、予め音源モデルパラメータが与えられていなくとも音源パラメータの推定が行える。
【実施例１】
【００４１】
図１にこの発明の音源パラメータ推定装置１００の機能構成例を示す。その動作フローを図２に示す。音源パラメータ推定装置１００は、特徴抽出部９１０と、音源モデル記憶部９０と、音源の数に対応した数の音源パワーパラメータ更新部６０１〜６０Ｎｓと、音源パワーパラメータ更新部６０１〜６０Ｎｓと同じ数の音源位置パラメータ更新部９４０１〜９４０Ｎｓと、音源占有度更新部７０と、音源モデルパラメータ更新部８０と、を具備する。その各部の機能は、例えばＲＯＭ、ＲＡＭ、ＣＰＵ等で構成されるコンピュータに所定のプログラムが読み込まれて、ＣＰＵがそのプログラムを実行することで実現されるものである。
【００４２】
音源パラメータ推定装置１００は、従来技術で説明した音源パワーパラメータ推定装置９００に対して音源モデルパラメータ更新部８０を備える点で異なる。特徴抽出部９１０と音源位置パラメータ更新部９４０１〜９４０Ｎｓとは、参照符号から明らかなように音源パワーパラメータ推定装置９００と同じものである。
【００４３】
各機能部の動作を説明する。
【００４４】
〔特徴抽出部〕
特徴抽出部９１０は、複数の音源信号を複数（m本）のマイクロホンで収音した時間領域信号を時間周波数領域信号に変換した観測信号ｘ(m)n,kを入力として、各時間周波数点(n,k)における音源位置特徴量An,kと音源パワー特徴量Xn,kを抽出する。
【００４５】
音源パワー特徴量Xn,kは、例えば、１本目のマイクロホンが収音した信号の対数パワースペクトルを音源パワー特徴量として抽出する場合には式（１３）に示すように計算される。音源パワー特徴量Xn,kは、音源占有度更新部７０と音源パワーパラメータ更新部６０１〜６０Ｎｓに入力される。
【００４６】
【数９】
0005438704.tif 000010
【００４７】
音源位置特徴量An,kは、一般に各時間周波数点における異なるマイクロホン間での信号の位相差や強度比などに表れる。したがって、音源位置特徴量An,kは、信号の位相差や強度比を異なるマイクロホンペアごとにまとめて出来るベクトルであったり、そこから更に何らかの特徴抽出を行った結果の値として抽出される。例えば、２本のマイクロホンで収音した信号の位相差を音源位置特徴量An,kとして抽出する場合、式（１４）に示すように計算される。
【００４８】
【数１０】
0005438704.tif 000011
【００４９】
〔音源モデル記憶部〕
音源モデル記憶部９０は、複数の音源信号それぞれの音源パワー時系列全体の状態を表す音源パワーパラメータq(l)の事前確率密度関数p(q(l);ψ（ｌ）)と、その音源パワーパラメータq(l)が与えられた場合の各音源信号の各時間周波数点における音源パワー特徴量のモデルβq(l),n,k,ψ(l)(S)とを記憶する。(S)は音源パワー特徴量Xn,kを表す変数である。
【００５０】
〔音源占有度更新部〕
音源占有度更新部７０は、音源位置特徴量An,kと音源パワー特徴量Xn,kと各音源の更新された音源パワーパラメータ^q (l)と音源位置パラメータ^φ (l)と、音源モデルパラメータ^ψ (l)と、音源モデル記憶部９０に記憶された音源パワーパラメータの事前確率密度関数p(q(l);ψ（ｌ）)と音源パワー特徴量のモデルβq(l),n,k,ψ(l)(S)とを入力として上記各音源の音源占有度を更新する。
【００５１】
音源占有度更新部７０は、Σl^M(l)=1となるように、音源占有度^M (l)n,kを、例えば乱数で初期化する（ステップＳ７０）。若しくは、従来技術で説明した音源パラメータ推定装置９００と同じ方法を用いても良い。その後、音源パワーパラメータ更新部６０１〜６０Ｎｓと音源占有度更新部７０と音源位置パラメータ更新部９４０１〜９４０Ｎｓと音源モデルパラメータ更新部８０の各処理を収束するまで繰り返す。
【００５２】
〔音源パワーパラメータ更新部〕
音源パワーパラメータ更新部６０１〜６０Ｎｓは、音源パワー特徴量Xn,kと、各音源ｌごとに初期化された音源占有度^M (l)n,kと、音源モデル記憶部９０に記憶された音源パワーパラメータ^q (l)の事前確率密度関数p(q(l);ψ（ｌ）)と音源パワー特徴量のモデルβq(l),n,k,ψ(l)(S)と音源モデルパラメータ^ψ (l)とを入力として、音源パワーパラメータ^q (l)を式（１５）に示すように更新（M-step）する（ステップＳ６０）。
【００５３】
【数１１】
0005438704.tif 000012
【００５４】
〔音源位置パラメータ更新部〕
音源位置パラメータ更新部９４０１〜９４０Ｎｓは、各音源ｌごとに初期化された音源占有度^M (l)n,kと、音源位置特徴量An,kを入力として音源位置パラメータ^φ (l)を、式（１７）に示すように更新（M-step）する（ステップＳ９４０）。
【００５５】
【数１２】
0005438704.tif 000013
【００５６】
〔音源モデルパラメータ更新部〕
音源モデルパラメータ更新部８０は、音源パワー特徴量Xn,kと音源パワーパラメータ^q (l)と音源占有度^M (l)n,kと、音源モデル記憶部９０に記憶された音源パワーパラメータの事前確率密度関数p(q(l);ψ（ｌ）)と音源パワー特徴量のモデルβq(l),n,k,ψ(l)(S)とを入力として音源モデルパラメータ^ψ (l)を式（１５′）に示すように更新する（ステップＳ８１）。
【００５７】
【数１３】
0005438704.tif 000014
【００５８】
そして、音源占有度更新部７０は、各音源ｌごとに更新された音源パワーパラメータ^q(l)と音源モデルパラメータ^ψ(l)と音源位置特徴量An,kと音源パワー特徴量Xn,kと、音源モデル記憶部９０に記憶された音源パワーパラメータの事前確率密度関数p(q(l);ψ(l))と音源パワー特徴量のモデルβq(l),n,k,ψ(l)(S)と、を入力として、音源占有度^M (l)n,kを式（１８）に示すように更新（E-step）する（ステップＳ７１）。
【００５９】
【数１４】
0005438704.tif 000015
【００６０】
ステップＳ６０〜ステップＳ７１の処理は、収束が得られるまで繰り返される（ステップＳ７２のno）。より具体的な音源位置特徴量のモデル及び、音源パワー特徴量のモデルを用いた実施例２を次に説明する。
【実施例２】
【００６１】
先ず、特徴抽出部９１０は、式（１４）に基づきマイク間位相差を、音源位置特徴量An,kとして抽出する。また、各音源ｌに由来する観測信号のマイク間位相差は、各周波数ごとに異なる平均値μ(l)k、分散σ(l)kのガウス分布に従うと仮定する。すると式（９）は以下のように定義できる。
【００６２】
【数１５】
0005438704.tif 000016
【００６３】
但し、φ(l)k=[μ(l)k,σ(l)k]は、音源位置パラメータφ(l)のうち周波数kのみに関する部分を取り出したものであり、φ(l)は全ての周波数kについてφ(l)ｋを集めたφ(l)=[φ(l)1,…, φ(l)Nk]である。N（・）は、ガウス分布の確率密度関数を表す。
【００６４】
一方、特徴抽出部９１０は、式（１３）に基づき、どれか一つのマイクロホン信号の対数パワースペクトルを音源パワー特徴量Xn,kとして抽出するものとする。さらに、音源パワーパラメータq(l)は、q(l)=｛q(l)0,q(l)1,…}のように各時刻の状態を表す状態系列に分解され、一次のマルコフ過程に従い状態遷移が各時刻で起こると仮定する。
【００６５】
但し、音源パワーパラメータq(l)0は隠れマルコフモデルの初期状態を表す。更に、式（３）で定義される各時間周波数点（n,k）におけるS(l)n,kの事後確率密度関数は、その時刻の状態q(l)ｎのみに依存するガウス分布に従うと仮定する。これを数式で表すと次のようになる。
【００６６】
【数１６】
0005438704.tif 000017
【００６７】
ここで、π(l)i=p(q(l)0=i)は、隠れマルコフモデルの初期状態がｉである事前確率、α(l)i,j=p(q(l)n=j｜q(l)n-1=i)は、隠れマルコフモデルが状態ｉから状態ｊへ移る状態遷移確率、βi,n,k,ψ(l)(S)=p(S(l)n,k=S｜q(l)n=i；ψ(l))=N(S(l)n,k;μ(l)i,k,σ(l)i,k)は、隠れマルコフモデルの状態ｉにおける出力の確率密度関数であり、μ(l)i,k及びσ(l)i,kはその平均と分散である。
【００６８】
この定義に基づくと、音源モデルパラメータ^ψ(l)は、全てのi,j,k,lに対するπ(l)ｉ, α(l)i,j,μ(l)i,k,σ(l)i,kで構成される。この実施例では、全て若しくは一部の音源モデルパラメータが音源パラメータの一部として期待値最大化アルゴリズムにより推定される。
【００６９】
以上の仮定の下、図２で説明済みの期待値最大化アルゴリズムのM-step1は、各音源ｌごとに、音源パワーパラメータ更新部６０１〜６０Ｎｓが式（２２）を満たす状態時系列^q (l)=[^q (l)0,…,^q (l)Ns]を、Viterbiアルゴリズムを用いて更新する。
【００７０】
【数１７】
0005438704.tif 000018
【００７１】
また、M-step2は、各音源ｌごとに、音源位置パラメータ更新部９４０１〜９４０Ｎｓが、全ての周波数kで、φ(l)k=[μ(l)k,σ(l)k]を次のように更新する。
【００７２】
【数１８】
0005438704.tif 000019
【００７３】
また、M-step3は、音源ｌごとに、音源モデルパラメータ更新部８０が音源モデルパラメータ^ψ (l)を更新する。まず、π(l),α(l)を、i,jに関するπ(l)ｉ,α(l) ｉ,jの集合とすると、π(l),α(l)は以下のように更新される。
【００７４】
【数１９】
0005438704.tif 000020
【００７５】
上記の更新は、隠れマルコフモデルの学習のための既知の方法で容易に実現することが可能である。一方、μ(l),σ(l)を全てのi,kに関するμ(l)i,k,σ(l) i,kの集合とすると、μ(l),σ(l)の更新は、以下のように実現される。
【００７６】
【数２０】
0005438704.tif 000021
【００７７】
上記の更新を実現する一つの方法として、準ニュートン法、共役勾配法などに代表される逐次最大化アルゴリズムを上げることができる。これには、一般的に知られる多くのアルゴリズムを適用することができる。一方、上記の関数中に含まれるlog(βq(l)n,n,k,ψ(l)(Xn,k))は解析的な扱いが容易であるが、log(ρq(l)n,n,k,ψ(l)(Xn,k))は積分演算が含まれているため解析的な扱いが複雑になり、逐次最大化アルゴリズムに必要とされる計算を比較的複雑にしてしまう問題がある。これを回避する一つの方法は、log(ρq(l)n,n,k,ψ(l)(Xn,k))を解析的な扱いが比較的容易な関数で近似することである。以下では、これについて少し詳しく説明する。
【００７８】
まず、x=(Xn,k-μ(l)i,k)(σ(l) i,k)-1/2とおき、f(x)= log(ρq(l)n,n,k,ψ(l)(Xn,k))と等価な以下の正規表現に書き換えておく。
【００７９】
【数２１】
0005438704.tif 000022
【００８０】
すると、上記の関数f(x)の一次導関数df(x)/dxは、以下のどちらか一方の関数で近似できる。
【００８１】
【数２２】
0005438704.tif 000023
【００８２】
これらの関数は、どちらも解析的な扱いが容易であり、特にμ(l)i,k,σ(l) i,kの更新に用いることで、例えば、以下の２つの利点を得ることができる。
【００８３】
その１は、μ(l)i,kの更新に関してf(x)の一次導関数をg1(x)で近似することで、高速な収束が得られる繰り返し推定法の一つであるニュートン法に関しても大域的な収束が保証される。その２は、g2(x)は、log(ρq(l)n,n,k,ψ(l)(Xn,k))の一次導関数と同じ数学的な形式をしているため、μ(l)i,k,σ(l) i,kの更新においてf(x)の一次導関数を、g2(x)で近似することで、複雑な計算を行わずに逐次最大化アルゴリズムを適用することができる。
【００８４】
したがって、例えばg1(x)を用いると、ニュートン法によるμ(l)i,kの１回の更新は以下のように実現することができる。
【００８５】
【数２３】
0005438704.tif 000024
【００８６】
ここでν(l)iは、q(l)n=iを満たす時間nの集合を表す。ｇ1′(x)は、ｇ1(x)の一次導関数を表す。
【００８７】
一方、g2(x)を用いると、σ(l) i,kの更新は以下のように実現することができる。
【００８８】
【数２４】
0005438704.tif 000025
【００８９】
ただし、κ(l) i,kは以下の値をとるものとする。
【００９０】
【数２５】
0005438704.tif 000026
【００９１】
また、音源占有度更新部７０が行うE-stepは、音源占有度を式（３２）に示すように更新する。
【００９２】
【数２６】
0005438704.tif 000027
【００９３】
ここで、事前に学習した隠れマルコフモデルがない場合の音源モデルパラメータ^ψ (l)の初期化について他の方法を説明する。音源パワー特徴量の時系列を混合ガウスモデルでモデル化し、その結果得られた混合ガウス分布の分布パラメータψ′のうち混合比α′からπ(l),α(l)を、各ガウス分布の平均μ′i,kと分散σ′ i,kからμ(l)i,k，σ(l) i,kを定める。このとき、混合ガウスモデルは以下の形をとる。
【００９４】
【数２７】
0005438704.tif 000028
【００９５】
音源パワー特徴量Xn,kが与えられた条件下で、混合ガウス分布の分布パラメータを定めるためには、例えば、期待値最大化アルゴリズムなどのように、一般的に知られている方法を適用することができる。その結果得られた分布パラメータを元に、音源モデルパラメータ^ψ (l)を以下のように定めることができる。
【００９６】
【数２８】
0005438704.tif 000029
【００９７】
これにより、観測された音源パワー特徴量の分布を近似的に表現する音源モデルパラメータの初期化を行うことができる。
【００９８】
〔変形例１〕
次に、音源モデルパラメータ更新部８０が更新する音源モデルパラメータψ(l)の一つであるμ(l)i,kの更新方法の変形例を説明する。例えば、μ(l)i,kの事前確率密度関数p(μ(l)i,k)=N（μ(l)i,k;~μ (l)i,k，~σ (l) i,k）が与えられていると仮定する。この事前確率密度関数の分布パラメータとしては、上記して説明した例で述べた方法などに基づき初期化された音源モデルパラメータの値を用いて、~μ(l)i,k=^μ(l)i,k，~σ(l) i,k=^σ(l) i,kと定めることが効果的であることが実験により確認されている。
【００９９】
この事前確率密度関数を用いると^μ(l)i,kの更新式（２９）は、事後確率最大化基準に基づき以下のように修正される。
【０１００】
【数２９】
0005438704.tif 000030
【０１０１】
ここで、ρは事前確率密度関数の重みを調整するコントロールパラメータ（＞０）であり、事前確率密度関数を信頼する程度に基づいて自由に定めることができる。
【０１０２】
〔音源分離装置〕
図３に、この発明の音源分離装置２００の機能構成例を示す。その動作フローを図４に示す。音源分離装置２００は、上記した音源パラメータ推定装置１００と、音源分離部９５と、を具備する。
【０１０３】
音源分離部９５は、音源パラメータ推定装置１００が出力する音源パワー特徴量Xn,kと、更新した音源占有度^M (l)n,kと音源パワーパラメータ^q(l)nと音源モデルパラメータ^ψ(l)と、各音源信号の各時間周波数点における音源パワー特徴量のモデルβq(l),n,k,ψ(l)(S)と、を入力として複数の音源のそれぞれの音源分離信号^S(l)n,kを最小自乗誤差推定により求める。
【０１０４】
音源分離の方法は次式によって行う。
【０１０５】
【数３０】
0005438704.tif 000031
【０１０６】
〔確認実験〕
この発明の音源分離性能を評価する目的で確認実験を行った。実験条件を説明する。観測信号を３０組用意し、全ての観測信号において音源数はNs=2とした。各観測信号は、それぞれ２人の男性の発話、２人の女性の発話、若しくは１名の女性と１名の男性の発話の混合音で構成した。
【０１０７】
標本化周波数は16kHzとした。各観測信号に含まれる２つのマイクロホン信号は、各話者の発話に関するマイク間時間差がそれぞれ±1.5ミリ秒になるように、計算機上で信号を加算して合成した（混合条件１）。音源モデルパラメータの初期値は、上記した事前に学習した隠れマルコフモデルがない場合の初期化方法で初期化した値を用いた。そして、変形例１で説明した事前確率密度関数を利用した。各隠れマルコフモデルの状態数は４とした。
【０１０８】
実験結果を図５と図６に示す。図５の縦軸は、混合前の音声から推定した音源パワー特徴量に関する隠れマルコフモデルのそれぞれについての各状態の出力分布の平均のパワー(dB)、横軸は周波数(kHz)である。図５（ａ）は混合前の音声から推定した結果、図５（ｂ）は混合音から推定した結果を示す。
【０１０９】
混合音中から推定されたパラメータ（図５（ｂ））は、混合前の音声から推定されたパラメータと酷似しており、この発明のパラメータ推定精度の信頼性の高さを証明している。
【０１１０】
図６に、分離後の信号のケプストラム歪みと混合音の数との関係を示す。縦軸はケプストラム歪み(dB)、横軸は混合音の数（数が増えるほど観測信号が長くなる）である。太い実線（□）で示す特性は、この発明で音源モデルパラメータを混合前の各音声信号から学習した場合を示す。実線（○）で示す特性は、この発明で音源モデルパラメータを観測信号のみから推定した場合を示す。一点鎖線（△）で示す特性は、非特許文献１の方法を用いた場合を示す。
【０１１１】
図６（ａ）は、混合条件１の混合音に対して処理した結果、図６（ｂ）と（ｃ）は、残響のある別の環境で収録した混合音（混合条件２と３）に対して処理した結果を示している。全ての場合において、この発明の音源分離方法は、非特許文献１の方法よりもケプストラム歪みの小さな音源分離を実現している。また、音源モデルパラメータを事前学習していない場合でも、事前学習している場合に相当する性能が実現できていることが確認できる。
【０１１２】
上記装置における処理手段をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、各装置における処理手段がコンピュータ上で実現される。
【０１１３】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、DVD（Digital Versatile Disc）、DVD-RAM（Random Access Memory）、CD-ROM（Compact Disc Read Only Memory）、CD-R（Recordable）/RW（ReWritable）等を、光磁気記録媒体として、MO（Magneto Optical disc）等を、半導体メモリとしてEEP-ROM（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。
【０１１４】
また、このプログラムの流通は、例えば、そのプログラムを記録したDVD、CD-ROM等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記録装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【０１１５】
また、各手段は、コンピュータ上で所定のプログラムを実行させることにより構成することにしてもよいし、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【図１】
0005438704.tif 000032
【図２】
0005438704.tif 000033
【図３】
0005438704.tif 000034
【図４】
0005438704.tif 000035
【図５】
0005438704.tif 000036
【図６】
0005438704.tif 000037
【図７】
0005438704.tif 000038