音源パラメータ推定装置と音源分離装置とそれらの方法と、プログラムと記憶媒体
【課題】観測信号の音源位置特徴量と音源パワー特徴量の両者を考慮して音源パラメータを推定する音源パラメータ推定装置を提供する。
【解決手段】音源モデル記憶部は、複数の音源パワーパラメータの事前確率密度関数と、その音源パワーパラメータが与えられた場合の各音源信号の各時間周波数点における事後確率密度関数とを記憶する。音源パワーパラメータ更新部は、音源パワー特徴量と観測信号が得られた下での占有的な音源の事後確率密度関数である音源占有度を入力として、音源パワーパラメータの事前確率密度関数と各音源信号の事後確率密度関数に基づいて複数の音源の各音源の音源パワーパラメータを更新する。音源位置パラメータ更新部は、音源位置特徴量と音源占有度を入力として、各音源の音源位置パラメータを更新する。音源占有度更新部は、各音源の更新された音源パワーパラメータと音源位置パラメータを元に各音源の音源占有度を更新する。
【解決手段】音源モデル記憶部は、複数の音源パワーパラメータの事前確率密度関数と、その音源パワーパラメータが与えられた場合の各音源信号の各時間周波数点における事後確率密度関数とを記憶する。音源パワーパラメータ更新部は、音源パワー特徴量と観測信号が得られた下での占有的な音源の事後確率密度関数である音源占有度を入力として、音源パワーパラメータの事前確率密度関数と各音源信号の事後確率密度関数に基づいて複数の音源の各音源の音源パワーパラメータを更新する。音源位置パラメータ更新部は、音源位置特徴量と音源占有度を入力として、各音源の音源位置パラメータを更新する。音源占有度更新部は、各音源の更新された音源パワーパラメータと音源位置パラメータを元に各音源の音源占有度を更新する。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、複数の音源が同時に生成した音響信号が混ざって複数のマイクロホンで収音された観測信号から、各音源の音源パラメータを推定する音源パラメータ推定装置と、その音源パラメータに基づいて各音源を分離する音源分離装置とそれらの方法と、プログラムと記憶媒体に関する。
【背景技術】
【０００２】
従来の音源パラメータ推定装置９００の機能構成例を図６に示す（非特許文献１）。音源パラメータ推定装置９００は、特徴抽出部９０と、パラメータ推定部９１を備える。
【０００３】
特徴抽出部９０は、複数（Ns個）の音源が同時に生成した音響信号が混ざって複数（Nm本）のマイクロホンで収音された時間領域信号のそれぞれに対し、短時間フーリエ変換などを適用して得られる時間周波数領域信号x(m)n,k（nは時間番号、kは周波数番号、mはマイクロホン番号）を観測信号として受け取り、各時間周波数点(n,k)ごとに観測信号の音源位置に関する音源位置特徴量An,kを抽出する。
【０００４】
パラメータ推定部９１は、音源位置特徴量An,kを受け取り、各１番目の音源ごとに音源位置パラメータφ^(l)k、及び各時間周波数点において当該音源が占有的な音源である確率を表す音源占有度M~(l)n,kからなる二種類の音源パラメータを推定し、音源占有度を出力する。そのため、パラメータ推定部９１は、音源占有度更新部９１０と、音源の数に対応する数の音源位置パラメータ更新部９１１1〜９１１Nsと、パーミュテーション解決部９１２と、を備える。
【０００５】
先ず、音源占有度更新部９１０は、各時間周波数点における各音源の音源占有度を更新する処理を行う。そして、音源位置パラメータ更新部９１１1〜９１１Nsが各音源の音源位置パラメータを更新する。この更新する処理を交互に繰り返して各音源のパラメータの推定値を得る。そして、パーミュテーション解決部９１２が各音源パラメータの推定値が全周波数で同一の音源に関するものになるよう、音源間で音源位置パラメータ、及び音源占有度の入れ替えを行い音源占有度M~(l)n,kを得る。
【０００６】
また、従来の音源分離装置は、音源パラメータ推定装置９００に加えて図示しない音源分離部を備え、音源分離部が、各時間周波数点に対応する音源パワー特徴量Xn,kと各音源lの音源占有度M~ (l)n,kを入力として、各音源の分離信号S＾(l)n,kを推定して出力する。音源パワー特徴量Xn,kは特徴抽出部９０から入手する。
【先行技術文献】
【非特許文献】
【０００７】
【非特許文献１】H. Sawada, S. Araki, S. Makino, “A Two-Stage FrequencyDomain Blind Source Separation Method for Underdetermined Convolutive Mixtures,” IEEE Workshop on Applications of Signal Processing to Audio and Acoustics(WASPAA 2007), pp. 139-142, Oct. 2007.
【発明の概要】
【発明が解決しようとする課題】
【０００８】
従来の音源パラメータ推定装置は、パラメータ推定の計算が複雑になることを避ける目的で各時間周波数点における各音源の音源占有度を音源位置特徴量のみに基づいて求めていた。したがって、推定された音源占有度や音源位置パラメータは、観測信号の音源パワー特徴量に対して最適なものでは無かった。また、従来の音源分離装置は、その音源パワー特徴量に対して最適でない音源占有度と音源位置パラメータに基づいて分離信号を推定するので、分離信号も最適なものに成り得ない。この結果、従来の方法では、音源位置パラメータの推定値と分離信号の推定値の両方に、比較的大きな推定誤差が含まれる課題が残されていた。
【０００９】
この発明は、このような課題に鑑みてなされたものであり、観測信号の音源位置特徴量と音源パワー特徴量の両者を考慮しても比較的簡単な計算で最適な音源パラメータを推定することが出来る音源パラメータ推定装置と音源分離装置と、それらの方法とプログラムと記憶媒体を提供することを目的とする。
【課題を解決するための手段】
【００１０】
この発明の音源パラメータ推定装置は、音源モデル記憶部と、特徴抽出部と、音源パワーパラメータ更新部と、音源位置パラメータ更新部と、音源占有度更新部と、を具備する。音源モデル記憶部は、複数の音源信号それぞれの音源パワー時系列全体の状態を表す音源パワーパラメータの事前確率密度関数と、その音源パワーパラメータが与えられた場合の各音源信号の各時間周波数点における事後確率密度関数とを記憶する。特徴抽出部は、複数の音源信号を複数のマイクロホンで収音した時間領域信号を時間周波数領域信号に変換した観測信号を入力として、各時間周波数点における音源位置特徴量と音源パワー特徴量を抽出する。音源パワーパラメータ更新部は、音源パワー特徴量と観測信号が得られた下での占有的な音源の事後確率密度関数である音源占有度を入力として、音源パワーパラメータの事前確率密度関数と各音源信号の事後確率密度関数に基づいて各音源の音源パワーパラメータを更新する。音源位置パラメータ更新部は、音源位置特徴量と音源占有度を入力として、各音源の音源位置パラメータを更新する。音源占有度更新部は、音源位置特徴量と音源パワー特徴量と各音源の更新された音源パワーパラメータと音源位置パラメータを元に各音源の音源占有度を更新する。
【００１１】
また、この発明の音源分離装置は、この発明の音源パラメータ推定装置が出力する更新した音源占有度と、音源パワー特徴量と、更新した音源パワーパラメータと、上記各音源信号の各時間周波数点における事後確率密度関数と、を入力として複数の音源のそれぞれの音源分離信号を最小自乗誤差推定により求める音源分離部を更に備える。
【発明の効果】
【００１２】
この発明の音源パラメータ推定装置によれば、観測信号の音源位置特徴量と音源パワー特徴量の両者を考慮しながら最適な音源パラメータを推定することが出来る。その音源パラメータを用いるこの発明の音源分離装置は、誤差の少ない分離信号を出力することが可能になる。
【図面の簡単な説明】
【００１３】
【図１】この発明の音源パラメータ推定装置１００の機能構成例を示す図。
【図２】音源パラメータ推定装置１００の動作フローを示す図。
【図３】この発明の音源分離装置２００の機能構成例を示す図。
【図４】音源分離装置２００の動作フローを示す図。
【図５】確認実験の結果を示す図。
【図６】従来の音源パラメータ推定装置９００の機能構成を示す図。
【発明を実施するための形態】
【００１４】
以下、この発明の実施の形態を図面を参照して説明する。複数の図面中同一のものには
同じ参照符号を付し、説明は繰り返さない。実施例の説明の前にこの発明の基本的な考えについて説明する。
【００１５】
〔この発明の基本的な考え〕
この発明は、観測信号の音源位置特徴量と音源パワー特徴量の両者を考慮しながら最適な音源パラメータを推定する点で新しい。最初に説明に用いる記号について説明する。観測信号には、Ns個の音源信号が重畳しており、その音源信号をNm本のマイクロホンで収音する。m番目のマイクロホンからの収音した収音信号を短時間フーリエ変換等を用いて周波数領域の信号に変換した観測信号をx(m)n,kと表記する。nはn番目の時間つまりフレーム番号、kはk番目の周波数つまりビン番号であり、n番目の時間及びk番目の周波数に対応する時間周波数点を参照する場合に、時間周波数点(n,k)と表記する。なお、記号＾の位置や添え字の表記とその位置は、式中の表記が正しい。
【００１６】
この発明では、各音源信号に関する音源パワー特徴量のモデルを導入する。以下では、{xn,k}は、全時間周波数点(n,k)にわたる変数xn,kの集合を表すものとする。各音源信号の音源パワー時系列全体を{S(l)n,k}と表したときに、その同時確率密度関数は式（１）〜式（３）に示すようにモデル化されると仮定する。
【００１７】
【数１】
【００１８】
ここで、q(l)はl番目の音源の音源パワー時系列全体の状態を表す音源パワーパラメータを表す。なお、以下では、全ての音源のq(l)をまとめてq=[q(1),…,q(Ns)]とも表記することにする。
【００１９】
音源パワーパラメータq（l）が与えられた下で各時間周波数点(n,k)の音源信号の音源パワーがS(l)n,kとなる確率密度関数を音源パワー特徴量のモデルと称しβq(l),n,k(S)（式（３））と表記する。なお、式（１）の総和演算は、q(l)が離散値ではなく連続値をとる場合には、q(l)に関する積分演算に置換えて表現されるものとする。また、式（２）において、音源の状態が既知のもとでは、異なる時間周波数点における音源パワーS(l)n,kは相互に独立であるという仮定を導入している。
【００２０】
また、この発明では式（４）に示すように、各時間周波数点(n,k)において最も大きなエネルギーを持つ音源信号（以下、占有的な音源信号と呼ぶ）の音源パワーS(l)n,kは、観測信号の音源パワーと一致すると仮定する。
【００２１】
【数２】
【００２２】
また、占有的ではない音源lに関しては、S(l)n,k≦Xn,kの関係を持つと仮定する。すると、各音源信号の状態が既知の条件の下で、観測信号の音源パワーXn,kの事後確率密度関数は次のように表現できることが知られている（参考文献：S.J. Rennie, J.R. Hershey, and P.A. Olsen, “Hierarchical variational loopy belief propagation for multi-talker speech recognition,” Proc. ASRU-2009, pp. 176-181, 2009.）。
【００２３】
【数３】
【００２４】
この発明では、更に、上式は次のように分解可能であると仮定して用いる。
【００２５】
【数４】
【００２６】
また、この発明では音源位置特徴量から音源位置パラメータを推定するため、音源位置特徴量のモデルp(An,k;φ)を導入する。音源位置特徴量のモデルp(An,k;φ)は、各音源信号のエネルギーは異なる時間周波数点にわたり疎に分布していると仮定し、その時間周波数点において占有的な音源の音源位置のみに依存して決まると仮定する。
【００２７】
一般的に、全ての音源の音源位置パラメータφ(l)をまとめてφ=[φ(1)，…，φ(Ns)]と表すと、音源位置特徴量のモデルp(An,k;φ)、つまり観測信号の音源位置特徴量の確率密度関数は、混合分布として式（８）に示すように展開することが出来る。
【００２８】
【数５】
【００２９】
式（８）において、Zn,kは時間周波数点(n,k)において占有的な音源の番号を表す確率変数であり、Zn,k=lは、l番目の音源が占有的な音源である場合を示す。また、p(Zn,k=l)は、l番目の音源が時間周波数点(n,k)において占有的な音源になる事前確率密度関数を表している。更に、以降の説明では次の表記を用いることにする。
【００３０】
【数６】
【００３１】
γφ(l),n,k(A)は、時間周波数点(n,k)において占有的な音源の番号がlの場合に、音源位置特徴量Aが得られる確率密度関数を表す。これは、l番目の音源の音源位置パラメータφ(l)のみに依存するものとする。具体的なγφ(l),n,k(A)やφ(l)の定義については後述する。また、ψl,n,kは、占有的な音源の番号に関する事前確率密度関数を表す。
【００３２】
式（８）のもと、γφ(l),n,k(A)が定義されている場合、音源位置パラメータφ(l)と占有的な音源の番号に関する事前確率密度関数、ψl,n,kが与えられれば、音源位置特徴量のモデルp(An,k;φ)は一意に定めることが出来る。逆に、音源位置特徴量An,kが観測された場合に、最尤推定などの方法に従い、音源位置パラメータと占有的な音源の番号に関する事前確率密度関数ψl,n,kやその事後確率密度関数を推定することが出来る。
【００３３】
以上の定義に従うと、完全データの確率密度関数は式（１１）に示すように導出される。
【００３４】
【数７】
【００３５】
式（１１）において、ｑが音源パワーパラメータ、φが音源位置パラメータであり、このパラメータがパラメータ推定の対象である。この発明では、次の対数尤度関数を最大化する値として、音源パワーパラメータと音源位置パラメータを推定する。
【００３６】
【数８】
【００３７】
式（１３）で、Zn,kは隠れ変数として扱われる。隠れ変数を含む対数尤度関数の最大化には、例えば、期待値最大化アルゴリズムなどを用いることが出来る。期待値最大化アルゴリズムでは、音源パワーパラメータの推定値q＾と音源位置パラメータの推定値φ＾に基づき、観測信号が得られたもとでの占有的な音源の番号の事後確率密度関数M^(l)n,k=p(Zn,k|An,k,Xn,k,q^;φ^)をも同時に推定する必要がある。この発明では、この関数の値を音源占有度と称し、この値も音源パラメータに含めて考える。
【００３８】
以上述べた考えで、音源パワー特徴量のモデルβq(l),n,kと、音源位置特徴量のモデルp(An,k;φ)の両者を考慮しながら最適な音源パラメータを推定することで音源位置パラメータの推定値の誤差を減らすことが出来る。また、音源位置特徴量のモデルp(An,k;φ)（式（８））と、音源パワー特徴量のモデル（式（７））に、占有的な音源の番号を表す変数Zn,kを共有化することで、２つの特徴量を考慮しながら音源パラメータ推定の計算を簡単にすることが出来る。
【実施例１】
【００３９】
図１にこの発明の音源パラメータ推定装置１００の機能構成例を示す。その動作フローを図２に示す。音源パラメータ推定装置１００は、特徴抽出部１０と、音源モデル記憶部２０と、音源の数に対応した数の音源パワーパラメータ更新部３０1〜３０Nsと、音源パワーパラメータ更新部３０1〜３０Nsと同じ数の音源位置パラメータ更新部４０1〜４０Nsと、音源占有度更新部５０と、を具備する。その各部の機能は、例えばＲＯＭ、ＲＡＭ、ＣＰＵ等で構成されるコンピュータに所定のプログラムが読み込まれて、ＣＰＵがそのプログラムを実行することで実現されるものである。
【００４０】
特徴抽出部１０は、従来技術で説明した音源パラメータ推定装置９００の特徴抽出部９０と基本的に同じであり、音源位置特徴量An,kに加えて音源パワー特徴量Xn,kを出力する点で異なる。
【００４１】
音源パワー特徴量Xn,kは、例えば、１本目のマイクロホンが収音した信号の対数パワースペクトルを音源パワー特徴量として抽出する場合には式（１４）に示すように計算される。
【００４２】
【数９】
【００４３】
音源位置特徴量An,kは、一般に各時間周波数点における異なるマイクロホン間での信号の位相差や強度比などに表れる。したがって、音源位置特徴量An,kは、信号の位相差や強度比を異なるマイクロホンペアごとにまとめて出来るベクトルであったり、そこから更に何らかの特徴抽出を行った結果の値として抽出される。例えば、２本のマイクロホンで収音した信号の位相差を音源位置特徴量An,kとして抽出する場合、式（１５）に示すように計算される。
【００４４】
【数１０】
【００４５】
音源モデル記憶部２０は、各音源信号の状態を表す音源パワーパラメータq(l)の事前確率密度関数p(q(l))と、その音源パワーパラメータq(l)が与えられた場合の各音源信号の各時間周波数点における事後確率密度関数βq(l),n,k(S)を記憶する。q(l)は音源パワーパラメータを表し、Sは音源パワー特徴量Xn,kを表す。
【００４６】
音源パワー特徴量Xn,kは、音源占有度更新部５０と音源パラメータ更新部１〜Ns（以降３０１〜３０Ns）に入力される。音源位置特徴量An,kは、音源占有度更新部５０と音源位置パラメータ更新部１〜Ns（以降４０１〜４０Ns）に入力される。
【００４７】
音源占有度更新部５０は、ΣlM^(l)n,k=1となるよう、音源占有度M^(l)n,kを、例えば乱数で初期化する（ステップＳ５０）。その後、音源パラメータ更新部３０１〜３０Nsと音源占有度更新部５０と音源位置パラメータ更新部４０１〜４０Nsが、次の処理が収束するまで繰り返す。
【００４８】
音源パワーパラメータ更新部３０１〜３０Nsは、各音源lごとに初期化された音源占有度M^(l)n,kと、音源モデル記憶部２０に記憶された音源パワーパラメータの事前確率密度関数p(q(l))と、音源パワー特徴量のモデルβq(l),n,kを入力として音源パワーパラメータを、式（１６）に示すように更新(M-step)する（ステップＳ３０）。
【００４９】
【数１１】
【００５０】
音源位置パラメータ更新部４０１〜４０Nsは、各音源lごとに初期化された音源占有度M^(l)n,kと、音源位置特徴量An,kを入力として音源位置パラメータφ＾(l)を、式（１８）に示すように更新(M-step)する（ステップＳ４０）。
【００５１】
【数１２】
【００５２】
音源占有度更新部５０は、各音源lごとに更新された音源パワーパラメータと音源位置パラメータφ＾(l)と音源位置特徴量An,kと音源パワー特徴量Xn,kを入力として、音源占有度M^(l)n,kを式（１９）に示すように更新（E-step）する（ステップＳ５１）。
【００５３】
【数１３】
【００５４】
ステップＳ３０〜ステップＳ５１の処理は、収束が得られるまで繰り返される（ステップＳ５２のno）。より具体的な音源位置特徴量のモデル及び、音源パワー特徴量のモデルを用いた実施例２を次に説明する。
【実施例２】
【００５５】
先ず、特徴抽出部１０は、式（１５）に基づきマイク間位相差を、音源位置特徴量An,kとして抽出する。また、各音源lに由来する観測信号のマイク間位相差は、各周波数ごとに異なる平均値μ(l)k、分散σ(l)kのガウス分布に従うと仮定する。すると式（９）は以下のように定義できる。
【００５６】
【数１４】
【００５７】
但し、φ(l)k＝[μ(l)k, σ(l)k]は、音源位置パラメータφ(l)のうち周波数ｋのみに関する部分を取り出したものであり、φ(l)は全ての周波数kについてφ(l)kを集めたφ(l)＝[φ(ｌ)1，…，φ(l)Ｎｋ]である。Ｎ（・）は、ガウス分布の確率密度関数を表す。
【００５８】
一方、特徴抽出部１０は、式（１４）に基づき、どれか一つのマイクロホン信号の対数パワースペクトルを音源パワー特徴量として抽出するものとする。また、各音源の音源パワーの時系列は、隠れマルコフモデルに従うと仮定する。すなわち、音源パワーパラメータq(l)は、q(l)＝{q(l)0, q(l)１,…}のように各時刻の状態を表す状態系列に分解され、一次のマルコフ過程に従い状態遷移が各時刻で起こると過程する。
【００５９】
但し、q(l)0は隠れマルコフモデルの初期状態を表す。更に、式（３）で定義される各時間周波数点(n,k)におけるS(l)n,kの事後確率密度関数は、その時刻の状態q(l)nのみに依存するガウス分布に従うと仮定する。これを数式で表すと次のようになる。
【００６０】
【数１５】
【００６１】
ここで、πi=p(q(l)0=i)は、隠れマルコフモデルの初期状態がiである事前確率、αi,j=p(q(l)n=j|q(l)n-1=i)は、隠れマルコフモデルが状態iから状態jへ移る状態遷移確率、βi,n,k(S)=p(S(l)n,k＝Ｓ|q(l)n=i)=N(S(l)n,k;μi,k,σi,k)は、隠れマルコフモデルの状態iにおける出力の確率密度関数であり、μi.k及びσi,kはその平均と分散である。
【００６２】
この発明では、全ての状態i,j及び全ての周波数kに対する隠れマルコフモデルのパラメータπi,αi,j,μi,k,σi,kは、音源信号のデータベースなどを用いて事前に学習されているものとする。
【００６３】
以上の仮定の下、図２で説明済みの期待値最大化アルゴリズムのM-step1は、各音源lごとに、音源パワーパラメータ更新部３０1〜３０Nsが式（２３）を満たす状態時系列q^(l)=[q^(l)0,…，q^Ns]を、Viterbiアルゴリズムを用いて更新する。
【００６４】
【数１６】
【００６５】
また、M-step2は、各音源lごとに、音源位置パラメータ更新４０1〜４０Nsが、全ての周波数kで、φ^(l)k＝｛μ^(l)k,σ^(l)k}を次のように更新する。
【００６６】
【数１７】
【００６７】
また、音源占有度更新部５０が行うE-stepは、音源占有度を式（２６）に示すように更新する。
【００６８】
【数１８】
【実施例３】
【００６９】
音源パワー特徴量のモデルの一部である式（７）の定義を、少し変形して式（２７）に示すようにすることで音質が改善する。
【００７０】
【数１９】
【００７１】
式（２７）において、左辺はZn,k=lが既知のもとでのXn,kの事後確率密度関数となっているところが、式（７）との違いである。この定義に従い、期待値最大化アルゴリズムを用いて音源パラメータ推定の手続きを導出すると、音源占有度更新部５０′で行うE-stepは、音源占有度M^(l)n,kを次に示すように更新する。
【００７２】
【数２０】
【００７３】
その他の手続きは実施例２と同じである。
【実施例４】
【００７４】
音源占有度M^(l)n,kの初期化について他の方法を説明する。その一つの方法は、従来技術を用いてμ^(l)kとσ(l)kを求めた後に、これらの値を音源位置パラメータの初期値として、音源占有度の初期値を次式で求める。
【００７５】
【数２１】
【００７６】
より少ない計算量で初期値を求める方法を更に説明する。各時間ごとの音源位置特徴量を全周波数でまとめて出来るベクトルAn=[An,1,An,2,…,An,Nk]（音源位置特徴量ベクトルと称する）を用い、更に次の仮定を導入する。
【００７７】
仮定その１「混合音中でも、多くの時間において、何れか一つの音源だけがほとんどの周波数kにおいて占有的である。」仮定その２「その多くの時間において、音源位置特徴量ベクトルAnは、占有的な音源の音源位置特徴量とほぼ一致する。」
この仮定に従うと、音源位置パラメータの概算値はAnをクラスタリングすることにより得られることになる。従って、例えば次のようにすることで、音源占有度M^(l)n,kの初期化を行うことが出来る。
【００７８】
先ず、全てのnに関するAnに、例えばk-meansクラスタリングを適用して得られるNs個のセントロイドを、各音源lの音源位置パラメータμ(l)=[ μ(l)1,μ(l)2,…,μ(l)K]の最初の推定値とする。また、全てのl,kに対してσ(l)k=1とする。次に、音源占有度M^(l)n,kの初期値を式（３０）に基づいて求める。
【００７９】
〔音源分離装置〕
図３にこの発明の音源分離装置２００の機能構成例を示す。その動作フローを図４に示す。音源分離装置２００は、上記した音源パラメータ推定装置１００と、音源分離部６０と、を具備する。音源分離部６０は、音源パラメータ推定装置１００が出力する更新した音源占有度M^(l)n,kと、音源パワー特徴量Xn,kと、更新した音源パワーパラメータｑ^(l)nと、上記各音源信号の各時間周波数点における事後確率密度関数βq(l),n,k(S)と、を入力として複数の音源のそれぞれの音源分離信号S^(l)n,kを、最小自乗誤差推定により求める。
【００８０】
音源分離信号の推定は式（３１）で行う。
【００８１】
【数２２】
【００８２】
〔確認実験〕
この発明の音源分離性能を評価する目的で確認実験を行った。１０組の観測信号から成る３つのテストデータセットを用意した。観測信号は合計３０組であり、全ての観測信号において、音源数はNs=2とした。３つのデータセットは、それぞれ、２人の男性の発話、２人の女性の発話、１名の女性と１名の男性の発話の混合音で構成されている。
【００８３】
標本化周波数は１６kHzとした。各観測信号に含まれる２つのマイクロホン信号は、各話者の発話に関するマイク間時間差がそれぞれ±1.5ミリ秒になるように、計算機上で信号を加算して合成した。また、３つの話者非依存隠れマルコフモデルを事前に用意した。各隠れマルコフモデルの状態数rは、それぞれr=4,16,256とした。
【００８４】
実験結果を図６に示す。図６の横軸は状態数r、縦軸は分離前と分離後の信号のケプストラム歪みを(dB)で示す。比較のために、この発明の方法の他に、OracleCBとPDで音源分離を行ったときの結果を一緒に示している。
【００８５】
OracleCBでは、先ず各音源に関して音を混ぜる前のクリーンな音声から隠れマルコフモデルの状態系列を求め、次に式（２６）に基づき音源パワーパラメータのみを用いて、音源占有度M^(l)n,kを推定した後に、音源分離を式（３１）に基づいて実施した。PDは非特許文献１に記載された方法である。
【００８６】
何れの状態数でもこの発明のケプストラム歪みが最小であり、大幅にケプストラム歪みを減少させる。この結果は、この発明の方法が、入力信号に対してより高精度な音源分離を実現できることを示している。このように、この発明の音源パラメータ推定方法と音源分離方法の有効性を確認することが出来た。
なお、上記装置における処理手段をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、各装置における処理手段がコンピュータ上で実現される。
【００８７】
また、上記方法及び装置において説明した処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されるとしてもよい。
【００８８】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）/ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。
【００８９】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記録装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【００９０】
また、各手段は、コンピュータ上で所定のプログラムを実行させることにより構成することにしてもよいし、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【技術分野】
【０００１】
この発明は、複数の音源が同時に生成した音響信号が混ざって複数のマイクロホンで収音された観測信号から、各音源の音源パラメータを推定する音源パラメータ推定装置と、その音源パラメータに基づいて各音源を分離する音源分離装置とそれらの方法と、プログラムと記憶媒体に関する。
【背景技術】
【０００２】
従来の音源パラメータ推定装置９００の機能構成例を図６に示す（非特許文献１）。音源パラメータ推定装置９００は、特徴抽出部９０と、パラメータ推定部９１を備える。
【０００３】
特徴抽出部９０は、複数（Ns個）の音源が同時に生成した音響信号が混ざって複数（Nm本）のマイクロホンで収音された時間領域信号のそれぞれに対し、短時間フーリエ変換などを適用して得られる時間周波数領域信号x(m)n,k（nは時間番号、kは周波数番号、mはマイクロホン番号）を観測信号として受け取り、各時間周波数点(n,k)ごとに観測信号の音源位置に関する音源位置特徴量An,kを抽出する。
【０００４】
パラメータ推定部９１は、音源位置特徴量An,kを受け取り、各１番目の音源ごとに音源位置パラメータφ^(l)k、及び各時間周波数点において当該音源が占有的な音源である確率を表す音源占有度M~(l)n,kからなる二種類の音源パラメータを推定し、音源占有度を出力する。そのため、パラメータ推定部９１は、音源占有度更新部９１０と、音源の数に対応する数の音源位置パラメータ更新部９１１1〜９１１Nsと、パーミュテーション解決部９１２と、を備える。
【０００５】
先ず、音源占有度更新部９１０は、各時間周波数点における各音源の音源占有度を更新する処理を行う。そして、音源位置パラメータ更新部９１１1〜９１１Nsが各音源の音源位置パラメータを更新する。この更新する処理を交互に繰り返して各音源のパラメータの推定値を得る。そして、パーミュテーション解決部９１２が各音源パラメータの推定値が全周波数で同一の音源に関するものになるよう、音源間で音源位置パラメータ、及び音源占有度の入れ替えを行い音源占有度M~(l)n,kを得る。
【０００６】
また、従来の音源分離装置は、音源パラメータ推定装置９００に加えて図示しない音源分離部を備え、音源分離部が、各時間周波数点に対応する音源パワー特徴量Xn,kと各音源lの音源占有度M~ (l)n,kを入力として、各音源の分離信号S＾(l)n,kを推定して出力する。音源パワー特徴量Xn,kは特徴抽出部９０から入手する。
【先行技術文献】
【非特許文献】
【０００７】
【非特許文献１】H. Sawada, S. Araki, S. Makino, “A Two-Stage FrequencyDomain Blind Source Separation Method for Underdetermined Convolutive Mixtures,” IEEE Workshop on Applications of Signal Processing to Audio and Acoustics(WASPAA 2007), pp. 139-142, Oct. 2007.
【発明の概要】
【発明が解決しようとする課題】
【０００８】
従来の音源パラメータ推定装置は、パラメータ推定の計算が複雑になることを避ける目的で各時間周波数点における各音源の音源占有度を音源位置特徴量のみに基づいて求めていた。したがって、推定された音源占有度や音源位置パラメータは、観測信号の音源パワー特徴量に対して最適なものでは無かった。また、従来の音源分離装置は、その音源パワー特徴量に対して最適でない音源占有度と音源位置パラメータに基づいて分離信号を推定するので、分離信号も最適なものに成り得ない。この結果、従来の方法では、音源位置パラメータの推定値と分離信号の推定値の両方に、比較的大きな推定誤差が含まれる課題が残されていた。
【０００９】
この発明は、このような課題に鑑みてなされたものであり、観測信号の音源位置特徴量と音源パワー特徴量の両者を考慮しても比較的簡単な計算で最適な音源パラメータを推定することが出来る音源パラメータ推定装置と音源分離装置と、それらの方法とプログラムと記憶媒体を提供することを目的とする。
【課題を解決するための手段】
【００１０】
この発明の音源パラメータ推定装置は、音源モデル記憶部と、特徴抽出部と、音源パワーパラメータ更新部と、音源位置パラメータ更新部と、音源占有度更新部と、を具備する。音源モデル記憶部は、複数の音源信号それぞれの音源パワー時系列全体の状態を表す音源パワーパラメータの事前確率密度関数と、その音源パワーパラメータが与えられた場合の各音源信号の各時間周波数点における事後確率密度関数とを記憶する。特徴抽出部は、複数の音源信号を複数のマイクロホンで収音した時間領域信号を時間周波数領域信号に変換した観測信号を入力として、各時間周波数点における音源位置特徴量と音源パワー特徴量を抽出する。音源パワーパラメータ更新部は、音源パワー特徴量と観測信号が得られた下での占有的な音源の事後確率密度関数である音源占有度を入力として、音源パワーパラメータの事前確率密度関数と各音源信号の事後確率密度関数に基づいて各音源の音源パワーパラメータを更新する。音源位置パラメータ更新部は、音源位置特徴量と音源占有度を入力として、各音源の音源位置パラメータを更新する。音源占有度更新部は、音源位置特徴量と音源パワー特徴量と各音源の更新された音源パワーパラメータと音源位置パラメータを元に各音源の音源占有度を更新する。
【００１１】
また、この発明の音源分離装置は、この発明の音源パラメータ推定装置が出力する更新した音源占有度と、音源パワー特徴量と、更新した音源パワーパラメータと、上記各音源信号の各時間周波数点における事後確率密度関数と、を入力として複数の音源のそれぞれの音源分離信号を最小自乗誤差推定により求める音源分離部を更に備える。
【発明の効果】
【００１２】
この発明の音源パラメータ推定装置によれば、観測信号の音源位置特徴量と音源パワー特徴量の両者を考慮しながら最適な音源パラメータを推定することが出来る。その音源パラメータを用いるこの発明の音源分離装置は、誤差の少ない分離信号を出力することが可能になる。
【図面の簡単な説明】
【００１３】
【図１】この発明の音源パラメータ推定装置１００の機能構成例を示す図。
【図２】音源パラメータ推定装置１００の動作フローを示す図。
【図３】この発明の音源分離装置２００の機能構成例を示す図。
【図４】音源分離装置２００の動作フローを示す図。
【図５】確認実験の結果を示す図。
【図６】従来の音源パラメータ推定装置９００の機能構成を示す図。
【発明を実施するための形態】
【００１４】
以下、この発明の実施の形態を図面を参照して説明する。複数の図面中同一のものには
同じ参照符号を付し、説明は繰り返さない。実施例の説明の前にこの発明の基本的な考えについて説明する。
【００１５】
〔この発明の基本的な考え〕
この発明は、観測信号の音源位置特徴量と音源パワー特徴量の両者を考慮しながら最適な音源パラメータを推定する点で新しい。最初に説明に用いる記号について説明する。観測信号には、Ns個の音源信号が重畳しており、その音源信号をNm本のマイクロホンで収音する。m番目のマイクロホンからの収音した収音信号を短時間フーリエ変換等を用いて周波数領域の信号に変換した観測信号をx(m)n,kと表記する。nはn番目の時間つまりフレーム番号、kはk番目の周波数つまりビン番号であり、n番目の時間及びk番目の周波数に対応する時間周波数点を参照する場合に、時間周波数点(n,k)と表記する。なお、記号＾の位置や添え字の表記とその位置は、式中の表記が正しい。
【００１６】
この発明では、各音源信号に関する音源パワー特徴量のモデルを導入する。以下では、{xn,k}は、全時間周波数点(n,k)にわたる変数xn,kの集合を表すものとする。各音源信号の音源パワー時系列全体を{S(l)n,k}と表したときに、その同時確率密度関数は式（１）〜式（３）に示すようにモデル化されると仮定する。
【００１７】
【数１】
【００１８】
ここで、q(l)はl番目の音源の音源パワー時系列全体の状態を表す音源パワーパラメータを表す。なお、以下では、全ての音源のq(l)をまとめてq=[q(1),…,q(Ns)]とも表記することにする。
【００１９】
音源パワーパラメータq（l）が与えられた下で各時間周波数点(n,k)の音源信号の音源パワーがS(l)n,kとなる確率密度関数を音源パワー特徴量のモデルと称しβq(l),n,k(S)（式（３））と表記する。なお、式（１）の総和演算は、q(l)が離散値ではなく連続値をとる場合には、q(l)に関する積分演算に置換えて表現されるものとする。また、式（２）において、音源の状態が既知のもとでは、異なる時間周波数点における音源パワーS(l)n,kは相互に独立であるという仮定を導入している。
【００２０】
また、この発明では式（４）に示すように、各時間周波数点(n,k)において最も大きなエネルギーを持つ音源信号（以下、占有的な音源信号と呼ぶ）の音源パワーS(l)n,kは、観測信号の音源パワーと一致すると仮定する。
【００２１】
【数２】
【００２２】
また、占有的ではない音源lに関しては、S(l)n,k≦Xn,kの関係を持つと仮定する。すると、各音源信号の状態が既知の条件の下で、観測信号の音源パワーXn,kの事後確率密度関数は次のように表現できることが知られている（参考文献：S.J. Rennie, J.R. Hershey, and P.A. Olsen, “Hierarchical variational loopy belief propagation for multi-talker speech recognition,” Proc. ASRU-2009, pp. 176-181, 2009.）。
【００２３】
【数３】
【００２４】
この発明では、更に、上式は次のように分解可能であると仮定して用いる。
【００２５】
【数４】
【００２６】
また、この発明では音源位置特徴量から音源位置パラメータを推定するため、音源位置特徴量のモデルp(An,k;φ)を導入する。音源位置特徴量のモデルp(An,k;φ)は、各音源信号のエネルギーは異なる時間周波数点にわたり疎に分布していると仮定し、その時間周波数点において占有的な音源の音源位置のみに依存して決まると仮定する。
【００２７】
一般的に、全ての音源の音源位置パラメータφ(l)をまとめてφ=[φ(1)，…，φ(Ns)]と表すと、音源位置特徴量のモデルp(An,k;φ)、つまり観測信号の音源位置特徴量の確率密度関数は、混合分布として式（８）に示すように展開することが出来る。
【００２８】
【数５】
【００２９】
式（８）において、Zn,kは時間周波数点(n,k)において占有的な音源の番号を表す確率変数であり、Zn,k=lは、l番目の音源が占有的な音源である場合を示す。また、p(Zn,k=l)は、l番目の音源が時間周波数点(n,k)において占有的な音源になる事前確率密度関数を表している。更に、以降の説明では次の表記を用いることにする。
【００３０】
【数６】
【００３１】
γφ(l),n,k(A)は、時間周波数点(n,k)において占有的な音源の番号がlの場合に、音源位置特徴量Aが得られる確率密度関数を表す。これは、l番目の音源の音源位置パラメータφ(l)のみに依存するものとする。具体的なγφ(l),n,k(A)やφ(l)の定義については後述する。また、ψl,n,kは、占有的な音源の番号に関する事前確率密度関数を表す。
【００３２】
式（８）のもと、γφ(l),n,k(A)が定義されている場合、音源位置パラメータφ(l)と占有的な音源の番号に関する事前確率密度関数、ψl,n,kが与えられれば、音源位置特徴量のモデルp(An,k;φ)は一意に定めることが出来る。逆に、音源位置特徴量An,kが観測された場合に、最尤推定などの方法に従い、音源位置パラメータと占有的な音源の番号に関する事前確率密度関数ψl,n,kやその事後確率密度関数を推定することが出来る。
【００３３】
以上の定義に従うと、完全データの確率密度関数は式（１１）に示すように導出される。
【００３４】
【数７】
【００３５】
式（１１）において、ｑが音源パワーパラメータ、φが音源位置パラメータであり、このパラメータがパラメータ推定の対象である。この発明では、次の対数尤度関数を最大化する値として、音源パワーパラメータと音源位置パラメータを推定する。
【００３６】
【数８】
【００３７】
式（１３）で、Zn,kは隠れ変数として扱われる。隠れ変数を含む対数尤度関数の最大化には、例えば、期待値最大化アルゴリズムなどを用いることが出来る。期待値最大化アルゴリズムでは、音源パワーパラメータの推定値q＾と音源位置パラメータの推定値φ＾に基づき、観測信号が得られたもとでの占有的な音源の番号の事後確率密度関数M^(l)n,k=p(Zn,k|An,k,Xn,k,q^;φ^)をも同時に推定する必要がある。この発明では、この関数の値を音源占有度と称し、この値も音源パラメータに含めて考える。
【００３８】
以上述べた考えで、音源パワー特徴量のモデルβq(l),n,kと、音源位置特徴量のモデルp(An,k;φ)の両者を考慮しながら最適な音源パラメータを推定することで音源位置パラメータの推定値の誤差を減らすことが出来る。また、音源位置特徴量のモデルp(An,k;φ)（式（８））と、音源パワー特徴量のモデル（式（７））に、占有的な音源の番号を表す変数Zn,kを共有化することで、２つの特徴量を考慮しながら音源パラメータ推定の計算を簡単にすることが出来る。
【実施例１】
【００３９】
図１にこの発明の音源パラメータ推定装置１００の機能構成例を示す。その動作フローを図２に示す。音源パラメータ推定装置１００は、特徴抽出部１０と、音源モデル記憶部２０と、音源の数に対応した数の音源パワーパラメータ更新部３０1〜３０Nsと、音源パワーパラメータ更新部３０1〜３０Nsと同じ数の音源位置パラメータ更新部４０1〜４０Nsと、音源占有度更新部５０と、を具備する。その各部の機能は、例えばＲＯＭ、ＲＡＭ、ＣＰＵ等で構成されるコンピュータに所定のプログラムが読み込まれて、ＣＰＵがそのプログラムを実行することで実現されるものである。
【００４０】
特徴抽出部１０は、従来技術で説明した音源パラメータ推定装置９００の特徴抽出部９０と基本的に同じであり、音源位置特徴量An,kに加えて音源パワー特徴量Xn,kを出力する点で異なる。
【００４１】
音源パワー特徴量Xn,kは、例えば、１本目のマイクロホンが収音した信号の対数パワースペクトルを音源パワー特徴量として抽出する場合には式（１４）に示すように計算される。
【００４２】
【数９】
【００４３】
音源位置特徴量An,kは、一般に各時間周波数点における異なるマイクロホン間での信号の位相差や強度比などに表れる。したがって、音源位置特徴量An,kは、信号の位相差や強度比を異なるマイクロホンペアごとにまとめて出来るベクトルであったり、そこから更に何らかの特徴抽出を行った結果の値として抽出される。例えば、２本のマイクロホンで収音した信号の位相差を音源位置特徴量An,kとして抽出する場合、式（１５）に示すように計算される。
【００４４】
【数１０】
【００４５】
音源モデル記憶部２０は、各音源信号の状態を表す音源パワーパラメータq(l)の事前確率密度関数p(q(l))と、その音源パワーパラメータq(l)が与えられた場合の各音源信号の各時間周波数点における事後確率密度関数βq(l),n,k(S)を記憶する。q(l)は音源パワーパラメータを表し、Sは音源パワー特徴量Xn,kを表す。
【００４６】
音源パワー特徴量Xn,kは、音源占有度更新部５０と音源パラメータ更新部１〜Ns（以降３０１〜３０Ns）に入力される。音源位置特徴量An,kは、音源占有度更新部５０と音源位置パラメータ更新部１〜Ns（以降４０１〜４０Ns）に入力される。
【００４７】
音源占有度更新部５０は、ΣlM^(l)n,k=1となるよう、音源占有度M^(l)n,kを、例えば乱数で初期化する（ステップＳ５０）。その後、音源パラメータ更新部３０１〜３０Nsと音源占有度更新部５０と音源位置パラメータ更新部４０１〜４０Nsが、次の処理が収束するまで繰り返す。
【００４８】
音源パワーパラメータ更新部３０１〜３０Nsは、各音源lごとに初期化された音源占有度M^(l)n,kと、音源モデル記憶部２０に記憶された音源パワーパラメータの事前確率密度関数p(q(l))と、音源パワー特徴量のモデルβq(l),n,kを入力として音源パワーパラメータを、式（１６）に示すように更新(M-step)する（ステップＳ３０）。
【００４９】
【数１１】
【００５０】
音源位置パラメータ更新部４０１〜４０Nsは、各音源lごとに初期化された音源占有度M^(l)n,kと、音源位置特徴量An,kを入力として音源位置パラメータφ＾(l)を、式（１８）に示すように更新(M-step)する（ステップＳ４０）。
【００５１】
【数１２】
【００５２】
音源占有度更新部５０は、各音源lごとに更新された音源パワーパラメータと音源位置パラメータφ＾(l)と音源位置特徴量An,kと音源パワー特徴量Xn,kを入力として、音源占有度M^(l)n,kを式（１９）に示すように更新（E-step）する（ステップＳ５１）。
【００５３】
【数１３】
【００５４】
ステップＳ３０〜ステップＳ５１の処理は、収束が得られるまで繰り返される（ステップＳ５２のno）。より具体的な音源位置特徴量のモデル及び、音源パワー特徴量のモデルを用いた実施例２を次に説明する。
【実施例２】
【００５５】
先ず、特徴抽出部１０は、式（１５）に基づきマイク間位相差を、音源位置特徴量An,kとして抽出する。また、各音源lに由来する観測信号のマイク間位相差は、各周波数ごとに異なる平均値μ(l)k、分散σ(l)kのガウス分布に従うと仮定する。すると式（９）は以下のように定義できる。
【００５６】
【数１４】
【００５７】
但し、φ(l)k＝[μ(l)k, σ(l)k]は、音源位置パラメータφ(l)のうち周波数ｋのみに関する部分を取り出したものであり、φ(l)は全ての周波数kについてφ(l)kを集めたφ(l)＝[φ(ｌ)1，…，φ(l)Ｎｋ]である。Ｎ（・）は、ガウス分布の確率密度関数を表す。
【００５８】
一方、特徴抽出部１０は、式（１４）に基づき、どれか一つのマイクロホン信号の対数パワースペクトルを音源パワー特徴量として抽出するものとする。また、各音源の音源パワーの時系列は、隠れマルコフモデルに従うと仮定する。すなわち、音源パワーパラメータq(l)は、q(l)＝{q(l)0, q(l)１,…}のように各時刻の状態を表す状態系列に分解され、一次のマルコフ過程に従い状態遷移が各時刻で起こると過程する。
【００５９】
但し、q(l)0は隠れマルコフモデルの初期状態を表す。更に、式（３）で定義される各時間周波数点(n,k)におけるS(l)n,kの事後確率密度関数は、その時刻の状態q(l)nのみに依存するガウス分布に従うと仮定する。これを数式で表すと次のようになる。
【００６０】
【数１５】
【００６１】
ここで、πi=p(q(l)0=i)は、隠れマルコフモデルの初期状態がiである事前確率、αi,j=p(q(l)n=j|q(l)n-1=i)は、隠れマルコフモデルが状態iから状態jへ移る状態遷移確率、βi,n,k(S)=p(S(l)n,k＝Ｓ|q(l)n=i)=N(S(l)n,k;μi,k,σi,k)は、隠れマルコフモデルの状態iにおける出力の確率密度関数であり、μi.k及びσi,kはその平均と分散である。
【００６２】
この発明では、全ての状態i,j及び全ての周波数kに対する隠れマルコフモデルのパラメータπi,αi,j,μi,k,σi,kは、音源信号のデータベースなどを用いて事前に学習されているものとする。
【００６３】
以上の仮定の下、図２で説明済みの期待値最大化アルゴリズムのM-step1は、各音源lごとに、音源パワーパラメータ更新部３０1〜３０Nsが式（２３）を満たす状態時系列q^(l)=[q^(l)0,…，q^Ns]を、Viterbiアルゴリズムを用いて更新する。
【００６４】
【数１６】
【００６５】
また、M-step2は、各音源lごとに、音源位置パラメータ更新４０1〜４０Nsが、全ての周波数kで、φ^(l)k＝｛μ^(l)k,σ^(l)k}を次のように更新する。
【００６６】
【数１７】
【００６７】
また、音源占有度更新部５０が行うE-stepは、音源占有度を式（２６）に示すように更新する。
【００６８】
【数１８】
【実施例３】
【００６９】
音源パワー特徴量のモデルの一部である式（７）の定義を、少し変形して式（２７）に示すようにすることで音質が改善する。
【００７０】
【数１９】
【００７１】
式（２７）において、左辺はZn,k=lが既知のもとでのXn,kの事後確率密度関数となっているところが、式（７）との違いである。この定義に従い、期待値最大化アルゴリズムを用いて音源パラメータ推定の手続きを導出すると、音源占有度更新部５０′で行うE-stepは、音源占有度M^(l)n,kを次に示すように更新する。
【００７２】
【数２０】
【００７３】
その他の手続きは実施例２と同じである。
【実施例４】
【００７４】
音源占有度M^(l)n,kの初期化について他の方法を説明する。その一つの方法は、従来技術を用いてμ^(l)kとσ(l)kを求めた後に、これらの値を音源位置パラメータの初期値として、音源占有度の初期値を次式で求める。
【００７５】
【数２１】
【００７６】
より少ない計算量で初期値を求める方法を更に説明する。各時間ごとの音源位置特徴量を全周波数でまとめて出来るベクトルAn=[An,1,An,2,…,An,Nk]（音源位置特徴量ベクトルと称する）を用い、更に次の仮定を導入する。
【００７７】
仮定その１「混合音中でも、多くの時間において、何れか一つの音源だけがほとんどの周波数kにおいて占有的である。」仮定その２「その多くの時間において、音源位置特徴量ベクトルAnは、占有的な音源の音源位置特徴量とほぼ一致する。」
この仮定に従うと、音源位置パラメータの概算値はAnをクラスタリングすることにより得られることになる。従って、例えば次のようにすることで、音源占有度M^(l)n,kの初期化を行うことが出来る。
【００７８】
先ず、全てのnに関するAnに、例えばk-meansクラスタリングを適用して得られるNs個のセントロイドを、各音源lの音源位置パラメータμ(l)=[ μ(l)1,μ(l)2,…,μ(l)K]の最初の推定値とする。また、全てのl,kに対してσ(l)k=1とする。次に、音源占有度M^(l)n,kの初期値を式（３０）に基づいて求める。
【００７９】
〔音源分離装置〕
図３にこの発明の音源分離装置２００の機能構成例を示す。その動作フローを図４に示す。音源分離装置２００は、上記した音源パラメータ推定装置１００と、音源分離部６０と、を具備する。音源分離部６０は、音源パラメータ推定装置１００が出力する更新した音源占有度M^(l)n,kと、音源パワー特徴量Xn,kと、更新した音源パワーパラメータｑ^(l)nと、上記各音源信号の各時間周波数点における事後確率密度関数βq(l),n,k(S)と、を入力として複数の音源のそれぞれの音源分離信号S^(l)n,kを、最小自乗誤差推定により求める。
【００８０】
音源分離信号の推定は式（３１）で行う。
【００８１】
【数２２】
【００８２】
〔確認実験〕
この発明の音源分離性能を評価する目的で確認実験を行った。１０組の観測信号から成る３つのテストデータセットを用意した。観測信号は合計３０組であり、全ての観測信号において、音源数はNs=2とした。３つのデータセットは、それぞれ、２人の男性の発話、２人の女性の発話、１名の女性と１名の男性の発話の混合音で構成されている。
【００８３】
標本化周波数は１６kHzとした。各観測信号に含まれる２つのマイクロホン信号は、各話者の発話に関するマイク間時間差がそれぞれ±1.5ミリ秒になるように、計算機上で信号を加算して合成した。また、３つの話者非依存隠れマルコフモデルを事前に用意した。各隠れマルコフモデルの状態数rは、それぞれr=4,16,256とした。
【００８４】
実験結果を図６に示す。図６の横軸は状態数r、縦軸は分離前と分離後の信号のケプストラム歪みを(dB)で示す。比較のために、この発明の方法の他に、OracleCBとPDで音源分離を行ったときの結果を一緒に示している。
【００８５】
OracleCBでは、先ず各音源に関して音を混ぜる前のクリーンな音声から隠れマルコフモデルの状態系列を求め、次に式（２６）に基づき音源パワーパラメータのみを用いて、音源占有度M^(l)n,kを推定した後に、音源分離を式（３１）に基づいて実施した。PDは非特許文献１に記載された方法である。
【００８６】
何れの状態数でもこの発明のケプストラム歪みが最小であり、大幅にケプストラム歪みを減少させる。この結果は、この発明の方法が、入力信号に対してより高精度な音源分離を実現できることを示している。このように、この発明の音源パラメータ推定方法と音源分離方法の有効性を確認することが出来た。
なお、上記装置における処理手段をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、各装置における処理手段がコンピュータ上で実現される。
【００８７】
また、上記方法及び装置において説明した処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されるとしてもよい。
【００８８】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）/ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。
【００８９】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記録装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【００９０】
また、各手段は、コンピュータ上で所定のプログラムを実行させることにより構成することにしてもよいし、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【特許請求の範囲】
【請求項１】
複数の音源信号それぞれの音源パワー時系列全体の状態を表す音源パワーパラメータの事前確率密度関数と、その音源パワーパラメータが与えられた場合の各音源信号の各時間周波数点における事後確率密度関数とを記憶した音源モデル記憶部と、
上記複数の音源信号を複数のマイクロホンで収音した時間領域信号を時間周波数領域信号に変換した観測信号を入力として、各時間周波数点における音源位置特徴量と音源パワー特徴量を抽出する特徴抽出部と、
上記音源パワー特徴量と、上記観測信号が得られた下での占有的な音源の事後確率密度関数である音源占有度を入力として、上記音源パワーパラメータの事前確率密度関数と上記各音源信号の事後確率密度関数に基づいて上記各音源の音源パワーパラメータを更新する音源パワーパラメータ更新部と、
上記音源位置特徴量と音源占有度を入力として、上記各音源の音源位置パラメータを更新する音源位置パラメータ更新部と、
上記音源位置特徴量と音源パワー特徴量と各音源の更新された音源パワーパラメータと音源位置パラメータを元に上記各音源の音源占有度を更新する音源占有度更新部と、
を具備する音源パラメータ推定装置。
【請求項２】
請求項１に記載した音源パラメータ推定装置において、
上記音源位置特徴量はマイク間位相差であり、上記音源パワー特徴量は上記観測信号の対数パワースペクトルであり、その音源パワー特徴量の時系列は隠れマルコフモデルに従うこと、
を特徴とする音源パラメータ推定装置。
【請求項３】
請求項１又は２に記載した音源パラメータ推定装置において、
各時間ごとの音源位置特徴量を全周波数でまとめて出来る音源位置特徴量ベクトルにクラスタリングを適用して得られる中央値を、各音源の音源位置パラメータの最初の推定値とし、その結果得られる音源位置特徴量が得られた下での占有的な音源の事後確率密度関数を、上記音源占有度の初期値とすることを特徴とする音源パラメータ推定装置。
【請求項４】
請求項１乃至３の何れかに記載した音源パラメータ推定装置と、
上記音源パラメータ推定装置が出力する更新した音源占有度と、音源パワー特徴量と、更新した音源パワーパラメータと、上記各音源信号の各時間周波数点における事後確率密度関数と、を入力として複数の音源のそれぞれの音源分離信号を、最小自乗誤差推定により求める音源分離部と、
を具備する音源分離装置。
【請求項５】
複数の音源信号を複数のマイクロホンで収音した時間領域信号を時間周波数領域信号に変換した観測信号を入力として、各時間周波数点における音源位置特徴量と音源パワー特徴量を抽出する特徴抽出過程と、
上記音源パワー特徴量と上記観測信号が得られた下での占有的な音源の事後確率密度関数である音源占有度を入力として、音源モデル記憶部に記憶された複数の音源信号それぞれの音源パワー時系列全体の状態を表す音源パワーパラメータの事前確率密度関数とその音源パワーパラメータが与えられた場合の各音源信号の事後確率密度関数に基づいて上記各音源の音源パワーパラメータを更新する音源パワーパラメータ更新過程と、
上記音源位置特徴量と音源占有度を入力として、上記各音源の音源位置パラメータを更新する音源位置パラメータ更新過程と、
上記音源位置特徴量と音源パワー特徴量と各音源の更新された音源パワーパラメータと音源位置パラメータを元に上記各音源の音源占有度を更新する音源占有度更新過程と、
を含む音源パラメータ推定方法。
【請求項６】
請求項５に記載した音源パラメータ推定方法において、
上記音源位置特徴量はマイク間位相差であり、上記音源パワー特徴量は上記観測信号の対数パワースペクトルでありその音源パワー特徴量の時系列は隠れマルコフモデルに従うこと、
を特徴とする音源パラメータ推定方法。
【請求項７】
請求項５又は６に記載した音源パラメータ推定方法において、
各時間ごとの音源位置特徴量を全周波数でまとめて出来る音源位置特徴量ベクトルにクラスタリングを適用して得られる中央値を各音源の音源位置パラメータの最初の推定値とし、その結果得られる音源位置特徴量が得られた下での占有的な音源の事後確率密度関数を、上記音源占有度の初期値とすることを特徴とする音源パラメータ推定方法。
【請求項８】
請求項５乃至７の何れかに記載した音源パラメータ推定方法と、
上記音源パラメータ推定方法で更新した音源占有度と、音源パワー特徴量と、更新した音源パワーパラメータと、上記各音源信号の各時間周波数点における事後確率密度関数と、を入力として複数の音源のそれぞれの音源分離信号を、最小自乗誤差推定により求める音源分離過程と、
を含む音源分離方法。
【請求項９】
請求項１乃至５の何れかに記載した音源パラメータ推定装置又は音源分離装置としてコンピュータを機能させるためのプログラム。
【請求項１０】
請求項９に記載した何れかのプログラムを記憶したコンピュータで読み取り可能な記憶媒体。
【請求項１】
複数の音源信号それぞれの音源パワー時系列全体の状態を表す音源パワーパラメータの事前確率密度関数と、その音源パワーパラメータが与えられた場合の各音源信号の各時間周波数点における事後確率密度関数とを記憶した音源モデル記憶部と、
上記複数の音源信号を複数のマイクロホンで収音した時間領域信号を時間周波数領域信号に変換した観測信号を入力として、各時間周波数点における音源位置特徴量と音源パワー特徴量を抽出する特徴抽出部と、
上記音源パワー特徴量と、上記観測信号が得られた下での占有的な音源の事後確率密度関数である音源占有度を入力として、上記音源パワーパラメータの事前確率密度関数と上記各音源信号の事後確率密度関数に基づいて上記各音源の音源パワーパラメータを更新する音源パワーパラメータ更新部と、
上記音源位置特徴量と音源占有度を入力として、上記各音源の音源位置パラメータを更新する音源位置パラメータ更新部と、
上記音源位置特徴量と音源パワー特徴量と各音源の更新された音源パワーパラメータと音源位置パラメータを元に上記各音源の音源占有度を更新する音源占有度更新部と、
を具備する音源パラメータ推定装置。
【請求項２】
請求項１に記載した音源パラメータ推定装置において、
上記音源位置特徴量はマイク間位相差であり、上記音源パワー特徴量は上記観測信号の対数パワースペクトルであり、その音源パワー特徴量の時系列は隠れマルコフモデルに従うこと、
を特徴とする音源パラメータ推定装置。
【請求項３】
請求項１又は２に記載した音源パラメータ推定装置において、
各時間ごとの音源位置特徴量を全周波数でまとめて出来る音源位置特徴量ベクトルにクラスタリングを適用して得られる中央値を、各音源の音源位置パラメータの最初の推定値とし、その結果得られる音源位置特徴量が得られた下での占有的な音源の事後確率密度関数を、上記音源占有度の初期値とすることを特徴とする音源パラメータ推定装置。
【請求項４】
請求項１乃至３の何れかに記載した音源パラメータ推定装置と、
上記音源パラメータ推定装置が出力する更新した音源占有度と、音源パワー特徴量と、更新した音源パワーパラメータと、上記各音源信号の各時間周波数点における事後確率密度関数と、を入力として複数の音源のそれぞれの音源分離信号を、最小自乗誤差推定により求める音源分離部と、
を具備する音源分離装置。
【請求項５】
複数の音源信号を複数のマイクロホンで収音した時間領域信号を時間周波数領域信号に変換した観測信号を入力として、各時間周波数点における音源位置特徴量と音源パワー特徴量を抽出する特徴抽出過程と、
上記音源パワー特徴量と上記観測信号が得られた下での占有的な音源の事後確率密度関数である音源占有度を入力として、音源モデル記憶部に記憶された複数の音源信号それぞれの音源パワー時系列全体の状態を表す音源パワーパラメータの事前確率密度関数とその音源パワーパラメータが与えられた場合の各音源信号の事後確率密度関数に基づいて上記各音源の音源パワーパラメータを更新する音源パワーパラメータ更新過程と、
上記音源位置特徴量と音源占有度を入力として、上記各音源の音源位置パラメータを更新する音源位置パラメータ更新過程と、
上記音源位置特徴量と音源パワー特徴量と各音源の更新された音源パワーパラメータと音源位置パラメータを元に上記各音源の音源占有度を更新する音源占有度更新過程と、
を含む音源パラメータ推定方法。
【請求項６】
請求項５に記載した音源パラメータ推定方法において、
上記音源位置特徴量はマイク間位相差であり、上記音源パワー特徴量は上記観測信号の対数パワースペクトルでありその音源パワー特徴量の時系列は隠れマルコフモデルに従うこと、
を特徴とする音源パラメータ推定方法。
【請求項７】
請求項５又は６に記載した音源パラメータ推定方法において、
各時間ごとの音源位置特徴量を全周波数でまとめて出来る音源位置特徴量ベクトルにクラスタリングを適用して得られる中央値を各音源の音源位置パラメータの最初の推定値とし、その結果得られる音源位置特徴量が得られた下での占有的な音源の事後確率密度関数を、上記音源占有度の初期値とすることを特徴とする音源パラメータ推定方法。
【請求項８】
請求項５乃至７の何れかに記載した音源パラメータ推定方法と、
上記音源パラメータ推定方法で更新した音源占有度と、音源パワー特徴量と、更新した音源パワーパラメータと、上記各音源信号の各時間周波数点における事後確率密度関数と、を入力として複数の音源のそれぞれの音源分離信号を、最小自乗誤差推定により求める音源分離過程と、
を含む音源分離方法。
【請求項９】
請求項１乃至５の何れかに記載した音源パラメータ推定装置又は音源分離装置としてコンピュータを機能させるためのプログラム。
【請求項１０】
請求項９に記載した何れかのプログラムを記憶したコンピュータで読み取り可能な記憶媒体。
【図１】
【図２】
【図３】
【図４】
【図５】
【図６】
【図２】
【図３】
【図４】
【図５】
【図６】
【公開番号】特開２０１２−４２６６４（Ｐ２０１２−４２６６４Ａ）
【公開日】平成２４年３月１日（２０１２．３．１）
【国際特許分類】
【出願番号】特願２０１０−１８３１３０（Ｐ２０１０−１８３１３０）
【出願日】平成２２年８月１８日（２０１０．８．１８）
【出願人】（０００００４２２６）日本電信電話株式会社 (13,992)
【公開日】平成２４年３月１日（２０１２．３．１）
【国際特許分類】
【出願日】平成２２年８月１８日（２０１０．８．１８）
【出願人】（０００００４２２６）日本電信電話株式会社 (13,992)
[ Back to top ]