JP6059072B2 - モデル推定装置、音源分離装置、モデル推定方法、音源分離方法及びプログラム - Google Patents

Description

この発明は、複数の音が混合した混合信号をもとに、混合前の原音（音源信号）を復元する音源分離技術に関し、特に、音響系に関する事前情報を用いずに音源分離を行うブラインド音源分離技術のうちクラスタリングに基づく音源分離を行う技術に関する。

まず、ブラインド音源分離技術全体における、クラスタリングに基づくブラインド音源分離技術の位置づけを説明する。ブラインド音源分離技術には、単一のマイクロフォンを用いる方法と複数のマイクロフォンを備えるマイクロフォンアレイを用いる方法とがある。後者の方法は、音源分離の手がかりとして音源位置の違いを利用できるため、一般に前者の方法よりも高い音源分離性能を達成することができる。複数のマイクロフォンを用いた音源分離技術としては、独立成分分析による方法とクラスタリングに基づく方法とが一般的である。後者の方法は、音源の数がマイクロフォンの数よりも多い場合にも適用できるという利点がある。

クラスタリングに基づく音源分離技術は、混合信号の各時間周波数成分に寄与する音源信号は高々一つであるという仮定（W-Disjoint Orthogonality: WDO）に基づく（非特許文献１参照）。この仮定は、例えば音声の場合のように、各音源信号が、
・ごく少数の時間周波数成分のみが大きく、
・その他の時間周波数成分はすべて０に近い
というスパース性を満たす場合に、特に精度よく成立することが知られている。

ここで、記法を定義する。m番目のマイクロフォンで観測された混合信号の時間周波数変換をy_mτωで表す。ここで、τはフレーム番号を表し、ωは角周波数を表す。また、各マイクロフォンで観測された混合信号の時間周波数変換を並べたベクトルを式(1)で表す。

ここで、Mはマイクロフォンの個数であり、・^T（上付き文字のT）はベクトル・の転置を表す。また、仮にk番目の音源のみが存在するとしたときに、１番目のマイクロフォンにおいて観測される信号の時間周波数変換をs_kτωで表す。また、s_kτωからy_τωへの伝達関数をh_kωで表す。s_kτωを１番目のマイクロフォンにおいて定義したため、h_kωの第１成分は１に等しいことに注意する。

上述の記法を用いて、混合信号のモデルを記述する。まず、WDOを仮定せず、y_τωにすべての音源信号が寄与する一般の場合を考えると、混合信号のモデルは式(2)で与えられる。

ここで、Kは音源数（既知と仮定する）である。

一方、WDOを仮定する場合、y_τωに寄与する音源信号の番号をd(τ,ω)で表すと、

であるから、式(2)は式(4)のように単純になる。

WDOの仮定のもとでは、各音源s_kτωを推定する音源分離の問題は、d(τ,ω)を推定する問題と等価である。以下でその理由を説明する。まず、式(4)の第一要素を抜き出すと式(5)を得る。

式(5)と式(3)とにより、各音源信号は次のように表せる。

ただし、m_kτωは時間周波数マスクと呼ばれ、式(8)で定義される。

以上より、d(τ,ω)を推定できれば、これを用いて式(8)によりマスクm_kτωを計算でき、このマスクm_kτωを式(7)のように混合信号y_1τωに乗算することによりs_kτωを計算することができる。

d(τ,ω)の推定は、通常、各時間周波数成分に対して、音源の方向や位置などの特徴量を抽出して、この特徴量に基づいて時間周波数成分を各音源に対応するクラスにクラスタリングすることにより行われる。このクラスタリングの精度は、クラスタリングに基づく音源分離の性能に大きく影響する。

クラスタリングの従来技術として非特許文献２に記載されている方法が挙げられる。この方法では、混合信号の各時間周波数成分に対して、y_τωをノルムが１になるように正規化した、式(9)で定義される特徴量を計算し、この特徴量に基づいて周波数ビンごとに時間周波数成分をクラスタリングする。

ここで、||・||はベクトルのユークリッドノルムを表す。

この方法では、ある周波数ビンにおけるi番目のクラスと別の周波数ビンにおけるi番目のクラスとが、必ずしも同じ音源に対応するとは限らないというパーミュテーション問題が生じる。したがって、この方法では、周波数ビンごとのクラスタリングを行った後で、すべての周波数ビンにおけるi番目のクラスが同じ音源に対応するように各周波数ビンにおけるクラス番号を付け替えるパーミュテーション解決を行うというように、二段階の処理が必要であった。

A. Jourjine, S. Rickard, and O. Yilmaz, "Blind separation of disjoint orthogonal signals: demixing N sources from 2 mixtures", in Proceedings of ICASSP, vol. 5, pp. 2985-2988, Jun. 2000. H. Sawada, S. Araki, and S. Makino, "Underdetermined convolutive blind source separation via frequency bin-wise clustering and permutation alignment", IEEE Transactions ASLP, vol. 19, no. 3, pp. 516-527, Mar. 2011.

非特許文献２に記載されたようなクラスタリングに基づく音源分離の従来技術においては、周波数ビンごとのクラスタリングを行った後で、パーミュテーション解決を行うという、二段階の処理が必要であった。

この発明の目的は、パーミュテーション問題を生じず、二段階の処理を必要としない音源分離手法を提供することである。

上記の課題を解決するために、この発明の一態様によるモデル推定装置は、kを音源のインデックスとし、τを時間フレームのインデックスとし、ωを角周波数とし、K個の音源からの信号が混合された混合信号をM個のマイクロフォンで観測した観測信号から、各観測信号の時間周波数成分からなる観測信号ベクトルy_τωに対応する特徴ベクトルx_τωを抽出する特徴抽出部と、特徴ベクトルx_τωを所定の確率モデルにあてはめ、各音源の尤度の時系列が周波数ビン間で同期しているほど高い評価値を与える評価関数を用いて、確率モデルのモデルパラメータを推定するモデル推定部と、を含む。確率モデルは、各音源に関する特徴ベクトルx_τωの分布の重み付き和で表される混合モデルであり、確率モデルの混合重みは、時間フレームτに依存し、角周波数ωに依存しない重みであり、確率モデルのモデルパラメータは、混合重みと、各音源に関する特徴ベクトルx_τωの分布のパラメータである。

この発明の一態様による音源分離装置は、kを音源のインデックスとし、τを時間フレームのインデックスとし、ωを角周波数とし、K個の音源からの信号が混合された混合信号をM個のマイクロフォンで観測した観測信号から、各観測信号の時間周波数成分からなる観測信号ベクトルy_τωに対応する特徴ベクトルx_τωを抽出する特徴抽出部と、特徴ベクトルx_τωを所定の確率モデルにあてはめ、各音源の尤度の時系列が周波数ビン間で同期しているほど高い評価値を与える評価関数を用いて、確率モデルのモデルパラメータを推定するモデル推定部と、モデルパラメータを用いて音源kに対応するマスクm_kτωを求めるマスク作成部と、観測信号ベクトルy_τωからマスクm_kτωを用いて時間周波数領域の分離音を計算する分離音作成部と、を含む。確率モデルは、各音源に関する特徴ベクトルx_τωの分布の重み付き和で表される混合モデルであり、確率モデルの混合重みは、時間フレームτに依存し、角周波数ωに依存しない重みであり、確率モデルのモデルパラメータは、混合重みと、各音源に関する特徴ベクトルx_τωの分布のパラメータである。

この発明によれば、パーミュテーション問題を生じず、二段階の処理を必要としない音源分離技術を実現することができる。これにより、例えば、音源位置などが時間的に変化する時変の環境での音声強調のためのオンライン音源分離を容易に実現することが可能となる。

モデル推定装置の機能構成を例示する図。 モデル推定装置の処理フローを例示する図。 音源分離装置の機能構成を例示する図。 音源分離装置の処理フローを例示する図。 実験環境を説明するための図。 実験結果を示す図。 実験結果を示す図。 実験結果を示す図。 実験結果を示す図。

以下、この発明の実施の形態について詳細に説明する。なお、図面中において同じ機能を有する構成部には同じ番号を付し、重複説明を省略する。

［発明のポイント］
詳細は後述するが、この発明の音源分離技術のポイントは以下の３点である。
１．一般に音源信号は、「音源信号の時間周波数変換の振幅値の時系列{|s_kτω|}_τが、周波数ビン間で類似する」という共通振幅変調の性質をもつ（例えば、非特許文献「G. J. Brown, “Computational Auditory Scene Analysis: A Representational Approach”, Ph.D. thesis, University of Sheffield, 1992.」を参照）。この発明では、この共通振幅変調の性質を、パーミュテーション問題を回避するための手掛かりとして利用できることに着目した。各音源信号に対する、周波数ビン間での振幅変調の共通性を、時変・周波数非依存の混合重みとしてモデル化することにより、パーミュテーションを引き起こさずに、クラスタリングできる。この時変・周波数非依存の混合重みは、従来用いられてきた時不変の混合重みとは一線を画す画期的なアイディアである。
２．混合重みの事前分布としてディリクレ分布を導入する。このディリクレ分布のハイパーパラメータを調整することにより、時変・周波数非依存の混合重みが極端に変動しすぎないように制御し、音源分離性能の低下を防ぐことができる。
３．パラメータ更新の反復ごとに、目的関数が最大となるように平均方向と密度パラメータを音源間で置換する。これにより、混合重みを時変・周波数非依存にすることによって生じる目的関数の局所解への収束を回避し、適切にモデルパラメータを推定することができる。

［第一実施形態］
この発明の第一実施形態は、複数の音源からの信号を複数個のマイクロフォンで観測し、モデルパラメータを推定するモデル推定装置である。

図１を参照して、第一実施形態のモデル推定装置Ａの機能構成例を説明する。モデル推定装置Ａは、周波数領域変換部１、特徴抽出部２及びモデル推定部３を含む。モデル推定部３は、事後確率計算部３１、パラメータ更新部３２及びパラメータ保持部３３を含む。パラメータ更新部３２は、混合重み更新手段３２１、相関行列更新手段３２２、平均方向更新手段３２３、密度パラメータ更新手段３２４及びパーミュテーション解決手段３２５を含む。

図２を参照して、モデル推定装置Ａの動作例を手続きの順に従って説明する。

周波数領域変換部１へM個のマイクロフォンにより観測された時間領域の混合信号~y_tが入力される。時間領域の混合信号~y_tは式(10)で定義される。

ここで、tは時間インデックスを表し、・^T（上付き文字のT）はベクトル・の転置を表し、~y_mtはm（1≦m≦M）番目のマイクロフォンで観測された時間領域の混合信号を表す。

周波数領域変換部１は入力された時間領域の混合信号~y_tから短時間フーリエ変換などにより時間周波数領域の観測信号ベクトルy_τωを生成し出力する（ステップＳ１）。時間周波数領域の観測信号ベクトルy_τωは式(11)で定義される。

ここで、τは時間フレームのインデックスを表し、ωは角周波数を表し、y_mτωは混合信号~y_mtの時間周波数領域での表現である。

特徴抽出部２は、周波数領域変換部１の出力する時間周波数領域の観測信号ベクトルy_τωを入力とし、特徴ベクトルx_τωを計算し出力する（ステップＳ２）。特徴ベクトルx_τωの計算は、時間周波数領域の観測信号ベクトルy_τωを正規化することにより行ってもよいし、時間周波数領域の観測信号ベクトルy_τωを白色化した後に正規化することによって行ってもよいし、時間周波数領域の観測信号ベクトルy_τωを正規化した後に白色化して再度正規化することにより行ってもよい。例えば、時間周波数領域の観測信号ベクトルy_τωを正規化することにより特徴ベクトルx_τωを計算する場合は式(12)により計算すればよい。

また、例えば、時間周波数領域の観測信号ベクトルy_τωを白色化した後に正規化する場合には、以下のように特徴ベクトルx_τωを計算すればよい。まず、時間周波数領域の観測信号ベクトルy_τωを用いて、式(13)により時間周波数領域の観測信号ベクトルy_τωの標本相関行列R_ω ^ｙを計算する。

ここで、Tはフレームの個数であり、・^H（上付き文字のH）はエルミート（Hermite）転置である。

次に、標本相関行列R_ω ^ｙの固有値と固有ベクトルを計算する。計算した標本相関行列R_ω ^ｙの固有値を、大きい順に並べたものをσ_ω1,σ_ω2,…,σ_ωMと表す。したがって、式(14)の関係が成り立つ。

ここで、標本相関行列R_ω ^ｙはエルミート行列であるから、固有値σ_ω1,σ_ω2,…,σ_ωMはすべて実数であることに注意する。また、固有値σ_ω1,σ_ω2,…,σ_ωMに対応し、正規直交系をなす標本相関行列R_ω ^ｙの固有ベクトルをu_ω1,u_ω2,…,u_ωMで表す。ここで、標本相関行列R_ω ^ｙはエルミート行列であるから、このような固有ベクトルが存在することに注意する。

次に、行列Σ_ωを式(15)により求め、行列U_ωを式(16)により求める。

次に、行列U_ω及び行列Σ_ωを用いて、時間周波数領域の観測信号ベクトルy_τωを白色化したベクトルy’_τωを式(17)により計算する。

最後に、次式のようにベクトルy’_τωをそのノルムで正規化することにより、特徴ベクトルx_τωを計算する。

モデル推定部３は、特徴ベクトルx_τωを特徴ベクトルの分布を表す確率モデルに当てはめ、確率モデルを評価する所定の評価関数を用いて、信号抽出に適した確率モデルのモデルパラメータを計算する。

特徴ベクトルの分布を表す確率モデルは、例えば、以下のようにモデル化されたものである。音源の位置が固定の場合、各周波数ビンにおいて、特徴ベクトルx_τωは理想的にはそれぞれの音源kごとに固有の値をとる。ただし、実際には、雑音・残響の影響、モデル化誤差などに起因する変動が存在するため、特徴ベクトルx_τωは音源kごとにある値を中心に分布する。そのため、この発明では、音源kに関する特徴ベクトルx_τωの分布を、例えば、以下のようにワトソン（Watson）分布でモデル化する。

ここで、a_kωは音源kに関する特徴ベクトルの分布の中心を表し、平均方向（mean orientation）と呼ばれ、κ_kωは音源kに関する特徴ベクトルの分布の広がりの小ささを表し、密度パラメータ（concentration parameter）と呼ばれる。M(a,b,x)はクンマー（Kummer）関数である。クンマー関数についての詳細は「S. Sra and D. Karp, “The multivariate Watson distribution: maximum-likelihood estimation and other aspects”, arXiv: 1104.4422v2, 2012.（参考文献１）」を参照されたい。ここで、特徴ベクトルの分布が周波数ビンごとに定義されることに注意する。

一般に音源信号は、「音源信号の時間周波数変換の振幅値の時系列{|s_kτω|}_τが、周波数ビン間で類似する」という共通振幅変調の性質をもつ（例えば、「G. J. Brown, “Computational Auditory Scene Analysis: A Representational Approach”, Ph.D. thesis, University of Sheffield, 1992.」を参照）。この発明では、この共通振幅変調の性質を、パーミュテーション問題を回避するための手掛かりとして利用できることに着目した。上述のWDO性の仮定にもとづき、この共通振幅変調の性質を、クラスタリングの枠組みにおいて利用しやすい表現で言い換えると、「観測信号に寄与する音源インデックスの時系列{d(τ,ω)}_τは、周波数ビン間で類似する」と言える。この発明では、この周波数ビン間での{d(τ,ω)}_τの類似性を、「d(τ,ω)の事前分布P(d(τ,ω)=k)が、フレームτに依存（時変）し、周波数ビン（角周波数ω）にはよらない（周波数非依存）」とモデル化する。このような各音源信号に対する、周波数ビン間での振幅変調の共通性を利用することにより、パーミュテーションを引き起こさずに、クラスタリングできる。これが上述で示したこの発明の３つのポイントのうちの第１のポイントである。この事前確率をα_kτにより表す。なお、α_kτはΣ_k=1 ^Kα_kτ=1を満たす。

この事前確率は、１個の時間フレームごとに変化すると仮定してもよいし、複数の時間フレームからなるブロックごとに変化すると仮定してもよい。事前確率が１個の時間フレームごとに変化すると仮定する場合、任意の音源kと任意の時間フレームτに対して、α_kτは独立変数であり、推定すべきパラメータである。

一方、事前確率が数個の時間フレームからなるブロックごとに変化すると仮定する場合、Bをブロックの総数とし、ブロック番号をb=1,2,…,Bとし、Jを各ブロック内における時間フレームの総数とし、各ブロック内における時間フレームの番号をj=1,2,…,Jとすると、τ=(b−1)×J+jと表せ、α_{k,(b−1)×J+j}（j=1,2,…,J）は等しくなるから、推定すべきパラメータである混合重みは~α_kb=α_{k,(b−1)×J+1}により定義される~α_kbである。以下では、特に断りのない限り、事前確率が１個の時間フレームごとに変化すると仮定する場合について説明する。

以上より、特徴ベクトルx_τωの尤度関数は、式(20)で表す混合モデルで与えられる。

ここで、Θは、式(21)に示すパラメータ集合である。

ここで、{α_kτ}_kτは式(22)により定義される。

他の同様の記法もこれにならって定義される。以降では、α_kτを混合重みと呼ぶ。混合重みα_kτが急激に変動して音源分離性能が低下することを防ぐため、混合重みα_kτの事前分布として式(23)に示すディリクレ（Dirichlet）分布を用いる。これがこの発明の第２のポイントである。

ここで、Γはガンマ関数であり、φはハイパーパラメータと呼ばれる。φの値を十分大きく定めることにより、混合重みα_kτの変動を抑えることができる。φの値を微調整する必要はないが、例えば、φ=1,10,100,1000などの値を用いることができる。

混合重みα_kτ以外のパラメータについては一様な事前分布を仮定する。したがって、p(Θ)=Π_τp({α_kτ}_k)である。

モデル推定部３では、特徴ベクトルx_τωを以上のようにモデル化された確率モデルに当てはめ、確率モデルを評価する所定の評価関数を用いて、事後確率及び信号抽出に適したパラメータ集合Θを求める。

以下、モデル推定部３の各部の処理を詳細に説明する。モデル推定部３は、図１に示すとおり、事後確率計算部３１、パラメータ更新部３２及びパラメータ保持部３３を含む。モデル推定部３での処理に先立ち、パラメータ集合Θの初期値をパラメータ保持部３３に用意しておく（ステップＳ０）。この初期値は、例えば、α_kτ=1/K、κ_kω=20とし、a_kωは{x_τω}_τωから無作為に選ぶことにより設定することができる。

事後確率計算部３１は、パラメータ保持部３３に記憶されたパラメータ集合Θから事後確率γ_kτω、すなわち特徴ベクトルx_τωが与えられたもとでd(τ,ω)=kとなる条件付き確率を式(24)により計算する（ステップＳ３１）。

パラメータ更新部３２は、図１に示すとおり、混合重み更新手段３２１、相関行列更新手段３２２、平均方向更新手段３２３、密度パラメータ更新手段３２４及びパーミュテーション解決手段３２５を含み、現在のパラメータ集合Θを更新して新たなパラメータ集合Θ’を生成する（ステップＳ３２）。

混合重み更新手段３２１は、事後確率γ_kτωを用いて、式(25)を計算することにより、混合重みα_kτを新しい値α’_kτに更新する。

ここで、Fは周波数ビンの個数を表す。φ=1のとき、α’_kτは全周波数ビンにわたる事後確率γ_kτωの平均値となることがわかる。φの増加とともに、α’_kτは定数1/Kに近づく。

相関行列更新手段３２２は、特徴ベクトルx_τωと事後確率γ_kτωを用いて、式(26)を計算することにより、各音源kに対する相関行列R_kωを新しい値R’_kωに更新する。

平均方向更新手段３２３は、相関行列R_kωの正規化された主成分ベクトルとして、平均方向a_kωを新しい値a’_kωに更新する。

密度パラメータ更新手段３２４は、相関行列R_kωの最大固有値λ_kωを用いて、密度パラメータк_kωを式(27)により新しい値к’_kωに更新する。

パーミュテーション解決手段３２５は、式(28)〜(30)に示すように、各周波数ビンにおいて、平均方向a’_kωと密度パラメータк’_kωを、事後確率p(Θ’|{x_τω}_τω)が最大になるように音源間で置換する（ステップＳ３２５）。これがこの発明の第３のポイントである。

ここで、Π:{1,2,…,K}→{1,2,…,K}である。

なお、以上では、混合重みが１個の時間フレームごとに変化する場合の処理について説明したが、混合重みが複数の時間フレームからなるブロックごとに変化する場合は、混合重み更新手段３２１における混合重みα_kτの更新式(25)において、分子の事後確率γ_kτωの時間フレームτに含まれるF個の時間周波数成分についての和を事後確率γ_kτωの時間フレームτを含むブロックbに含まれるF×J個の時間周波数成分の和に置き換え、分母のFをF×Jで置き換えればよい。一方、相関行列更新手段３２２、平均方向更新手段３２３、密度パラメータ更新手段３２４及びパーミュテーション解決手段３２５においては、混合重みが１個の時間フレームごとに変化する場合の処理と同一の処理を行えばよい。

以下、パラメータ更新部３２における各更新式の導出根拠を説明する。パラメータ更新はEM（Expectation-Maximization）アルゴリズムを導入して、それに基づき行う。なお、{d(τ,ω)}_τωは、EMアルゴリズムにおける隠れ変数として扱う。

まず、MAP（Maximum a posteriori）推定のためのコスト関数L(Θ)は、式(31)〜(33)により与えられる。

ここで、{x_τω}_τωは互いに独立であると仮定し、Θに依存しない定数項を無視した。この目的関数を式(34)に示す制約条件のもとで最大化する。

目的関数L(Θ)は、パーミュテーション問題がない場合に大きい値を取るため、L(Θ)の最大化によりパーミュテーション問題が回避できる。実際、式(33)の第一項から分かるように、目的関数L(Θ)が大きくなるのは、混合重みα_kτが大きい値をとるk、τに対し、音源kに対する尤度（もっともらしさ）p(x_τω|d(τ,ω)=k,a_kω,κ_kω)が大きい場合である。したがって、L(Θ)の最大化により、音源kに対する尤度の時系列{p(x_τω|d(τ,ω)=k,a_kω,κ_kω)}_τが周波数ビン間で同期する。このことと、上述の「観測信号に寄与する音源インデックスの時系列{d(τ,ω)}_τは、周波数ビン間で類似する」という性質を考え合わせると、L(Θ)はパーミュテーション問題がない場合に大きい値を取ることがわかる。EMアルゴリズムで用いる評価関数（Q関数）は式(35)(36)により与えられる。

更新後のパラメータ集合Θ'は次式により定義され、

Q関数を式(34)の制約のもとで最大にするものとして導かれる。すなわち、混合重みα_kτの新たな値α'_kτを求める式(25)は、ラグランジュ（Lagrange）の未定乗数法によって、式(37)(38)により導出される。

ここで、μはラグランジュの未定乗数である。

平均方向の算出方法は、クーラン・フィッシャー（Courant-Fischer）の定理を式(39)に適用することにより導出される。

また、密度パラメータの更新式(27)については、まず∂Q/∂κ_kω=0より式(40)を得る。

ここで、

であり、λ_kωは相関行列R_kωの最大固有値である。上式は、近似的に次のように解くことができる（詳細は上述の参考文献１参照）。

パラメータ保持部３３は、パラメータ更新部３２での更新処理により得られたパラメータ集合Θ’を記憶する（ステップＳ３３）。また、事後確率計算部３１での次回の処理の際には、記憶したパラメータ集合Θ’をパラメータ集合Θとして提供する。

ステップＳ３１からステップＳ３３までの処理は、事前に設定した最大反復回数max_iterに達するまで、またはパラメータ更新部３２における各パラメータの更新による変動幅が収束判定の閾値Δよりも小さくなるまで、反復して行う（ステップＳ９１）。最大反復回数max_iter及び閾値Δの具体的な値は、例えば、max_iter=100、Δ=10^-10とすることができる。

ステップＳ９１において、モデル推定部３における処理が最大反復回数max_iterに達した場合、または各パラメータの更新による変動幅が閾値Δよりも小さくなった場合、モデル推定部３は反復終了後の事後確率γ^o _kτωを出力する。

［第二実施形態］
この発明の第二実施形態は、第一実施形態のモデル推定装置Ａを用いて音源分離装置として構成した実施形態である。

図３を参照して、第二実施形態の音源分離装置Ｂの機能構成例を説明する。音源分離装置Ｂは、第一実施形態のモデル推定装置Ａの各部に加えて、音源分離部４及び時間領域変換部５を含む。音源分離部４は、マスク作成部４１及び分離音作成部４２を含む。

図４を参照して、音源分離装置Ｂの動作例を手続きの順に従って説明する。ステップＳ０からステップＳ９１までの処理は第一実施形態のモデル推定装置Ａの動作例と同様であるので詳細な説明は省略する。

音源分離部４は、周波数領域変換部１の出力する混合音の時間周波数変換y_τωと、事後確率計算部３１の出力する反復終了後の事後確率γ^o _kτωとを用いて、分離音の時間周波数変換^s_kτωを推定する。

マスク作成部４１は、反復終了後の事後確率γ^o _kτωを用いて、混合音に含まれる各音源を抽出するマスクm_kτωを作成する（ステップＳ４１）。まず、マスク作成部４１は、反復終了後の事後確率γ^o _kτωを用いて、式(42)によりd(τ,ω)の推定値^d(τ,ω)を計算する。

次に、マスク作成部４１は、マスクm_kτωを式(43)により計算する。

なお、マスク作成部４１は、マスクm_kτωを式(44)により求めてもよい。

分離音作成部４２は、式(45)により、マスクm_kτωを混合音の時間周波数変換y_1τωに乗算し、分離音の時間周波数変換^s_kτωを計算する。

時間領域変換部５は、音源kごとに、時間周波数領域の分離信号^s_kτωを時間領域の分離信号~^s_ktに変換して出力する（ステップＳ５）。

［実験結果］
この発明の効果を確認するために二通りの実験を行った。
＜実験１＞
１番目の実験は、以下の二点を実証するための実験である。
・本発明の方法により、パーミュテーション解決の後処理を行わなくても、音源分離が可能であること。
・混合重みのディリクレ事前分布が、音源分離性能の低下を防ぐ効果を持つこと。

そのために、以下に示す３つの異なるクラスタリング方法に基づく音源分離手法により特徴ベクトルx_τωを音源分離し、その性能を比較した。
１．本発明の方法１（ディリクレ事前分布あり）：時変・周波数非依存の混合重みを持つ混合ワトソンモデルによるクラスタリング方法に基づく音源分離手法である。混合重みの事前分布はディリクレ分布であると仮定し、ディリクレ分布のハイパーパラメータをφ=10,10²,10³として構成した。
２．本発明の方法２（ディリクレ事前分布なし）：本発明の方法１において、ディリクレ分布のハイパーパラメータをφ=1として、ディリクレ事前分布の効果をなくすように構成した。
３．比較手法：時不変・周波数非依存の混合重み0.5を持つ混合ワトソンモデルによる、周波数ビンごとのクラスタリングの後、非特許文献２に記載されているパーミュテーション解決技術を適用するように構成した。混合重みの事前分布は一様であると仮定した。

特徴ベクトルx_τωとしては、観測信号ベクトルy_τωを正規化したのち、白色化し、再び正規化したものを用いた。

図５を参照して実験環境を説明する。直方体の空間の中央付近に正三角形を描くように３つのマイクロフォンが配置され、それらのマイクロフォンを取り囲む円を描くように４つの音源が配置される。４つの音源は図５の紙面下方向を0°として反時計回りに70°、150°、245°、315°の位置に配置されている。３つのマイクロフォンの描く正三角形の中心から４つの音源までの距離は1.2メートルとした。３つのマイクロフォンと４つの音源は水平に設置され、その床面からの高さは1.2メートルとした。実験では、図５中に太字で描いた紙面右側の２つの音源（70°及び150°の位置）と、紙面左側の２つのマイクロフォン（m=1,2）を用いた。

表１にその他の条件を含めた実験条件をまとめる。

音源分離性能の評価尺度としては、SDR（signal-to-distortion ratio）を、２つの音源と８通りの話者組み合わせとに関して平均したものを用いた。以下では、このSDRの平均値を単にSDRと呼ぶ。SDRが大きいほど音源分離性能が高いことを意味する。SDRについての詳細は、「E. Vincent, H. Sawada, P. Bofill, S. Makino, and J. Rosca, “First stereo audio source separation evaluation campaign: data, algorithms and results”, in Proceedings of ICA, pp. 552-559, 2007.（参考文献２）」を参照されたい。

表２に、異なる残響時間（RT60）それぞれに対する各手法のSDRを示す。

本発明の各方法及び比較手法の実験結果により同等のSDRが得られたことがわかる。特に、本発明の方法１においてφ=10³とすると、残響時間が130,200,250,300,370ミリ秒である場合には、比較手法と同じか0.1〜0.4dB大きいSDRが得られた。残響時間が440ミリ秒である場合には比較手法のSDRの方が高かったが、その差はわずか0.1dBであった。この結果から、本発明の方法により、パーミュテーション解決の後処理を行わなくても良好に音源分離ができていることがわかる。

また、本発明の方法１と本発明の方法２とを比較すると、前者の方が高いSDRを与えた。特に、本発明の方法１において、φ=10³とすると、本発明の方法２と比べて、0.8〜1.3dB大きいSDRが得られており、比較手法と同等以上のSDRが得られた。このことから、混合重みのディリクレ事前分布が、音源分離性能の低下を防ぐ効果を持つことがわかる。

＜実験２＞
２番目の実験は、パラメータ更新部におけるパーミュテーション解決手段の有効性を検証するための実験である。これを検証するために、次の二つの方法によるパラメータ推定を比較した。
１．本発明の方法１：実験１で用いた本発明の方法１において、ディリクレ分布のハイパーパラメータをφ=10³として構成した。
２．本発明の方法３：本発明の方法１において、パーミュテーション解決手段３２５における処理を行わない。

平均方向の初期値としては、1〜2kHzの範囲において、真の値を音源間で交換することで得られる、パーミュテーションが生じている初期値を用いた。ここで、真の値は、仮にそれぞれの音源を別個に観測できるとして、それを用いて計算した。他の実験条件は、実験１におけるものと同一である。

本発明の各方法と比較手法とによる平均方向の推定値を比較するために、次式に示す平均方向のマイクロフォン間位相差をプロットした。

ここで、[・]_mは、ベクトルの第m要素を表す。

実験２の実験結果を図６〜９に示す。実験結果は、位相差δ_kω（縦軸、単位は度（°））を周波数ω/(2π)（横軸、単位はキロヘルツ（kHz））の関数としてプロットした。各クラスタkに対するプロットを異なるマーカで示しており、クラスタk=1は「＋」でプロットし、クラスタk=2は「×」でプロットした。

図６は、本発明の方法３の反復終了後の値である。図７は、本発明の方法１の反復終了後の値である。図８は、初期値である。図９は、真の値である。

図６に示す本発明の方法３の反復終了後のマイクロフォン間位相差は、図８に示す初期値の場合と非常に近い値をとることがわかる。EMアルゴリズムは、局所解への収束性が保証されたアルゴリズムであることから、上述の結果は、パーミュテーションの生じている初期値が目的関数の局所解となっていることを示している。

一方、図７に示す本発明の方法１の反復終了後のマイクロフォン間位相差は、図９に示す真の値の場合と非常に近い値をとることがわかる。

この結果より、本発明のパーミュテーション解決手段はパーミュテーション問題を回避するのに有効であり、本発明の方法においてパーミュテーション解決手段を用いないとパーミュテーション問題の生じている局所解に陥ってしまうことがわかる。

以上に説明したとおり、実験１の結果から、この発明によりパーミュテーション解決の後処理を行わなくても音源分離が可能であること、及び、混合重みのディリクレ事前分布が音源分離性能の低下を防ぐ効果を持つことが確認された。また、実験２の結果から、パラメータ更新部におけるパーミュテーション解決手段の有効性が確認された。

［プログラム、記録媒体］
この発明は上述の実施形態に限定されるものではなく、この発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。上記実施例において説明した各種の処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。

また、上記実施形態で説明した各装置における各種の処理機能をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記各装置における各種の処理機能がコンピュータ上で実現される。

この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。

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

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

また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。

Ａ モデル推定装置
Ｂ 音源分離装置
１ 周波数領域変換部
２ 特徴抽出部
３ モデル推定部
３１ 事後確率計算部
３２ パラメータ更新部
３２１ 混合重み更新手段
３２２ 相関行列更新手段
３２３ 平均方向更新手段
３２４ 密度パラメータ更新手段
３２５ パーミュテーション解決手段
３３ パラメータ保持部
４ 音源分離部
４１ マスク作成部
４２ 分離音作成部
５ 時間領域変換部

Claims (11)

1. kを音源のインデックスとし、τを時間フレームのインデックスとし、ωを角周波数とし、
   K個の音源からの信号が混合された混合信号をM個のマイクロフォンで観測した観測信号から、各観測信号の時間周波数成分からなる観測信号ベクトルy_τωに対応する特徴ベクトルx_τωを抽出する特徴抽出部と、
   上記特徴ベクトルx_τωを所定の確率モデルにあてはめ、各音源の尤度の時系列が周波数ビン間で同期しているほど高い評価値を与える評価関数を用いて、上記確率モデルのモデルパラメータを推定するモデル推定部と、
   を含み、
   上記確率モデルは、各音源に関する特徴ベクトルx_τωの分布の重み付き和で表される混合モデルであり、
   上記確率モデルの混合重みは、時間フレームτに依存し、角周波数ωに依存しない重みであり、
   上記確率モデルのモデルパラメータは、上記混合重みと、各音源に関する上記特徴ベクトルx_τωの分布のパラメータである
   モデル推定装置。
2. 請求項１に記載のモデル推定装置であって、
   音源kに関する上記特徴ベクトルx_τωの分布は、平均方向をa_kωとし、密度パラメータをκ_kωとするワトソン分布であり、
   音源kに関する上記特徴ベクトルx_τωの分布のパラメータは、上記平均方向a_kωと上記密度パラメータκ_kωである
   モデル推定装置。
3. 請求項１に記載のモデル推定装置であって、
   上記混合重みの事前分布は、音源kに依存しないハイパーパラメータφを各混合重みの指数とする上記混合重みについてのディリクレ分布である
   モデル推定装置。
4. 請求項２に記載のモデル推定装置であって、
   上記混合重みの事前分布は、音源kに依存しないハイパーパラメータφを各混合重みの指数とする上記混合重みについてのディリクレ分布である
   モデル推定装置。
5. 請求項４に記載のモデル推定装置であって、
   上記モデル推定部は、音源kに関する上記特徴ベクトルx_τωの分布と音源kの時間フレームτにおける混合重みとの積に基づいて、特徴ベクトルx_τωが与えられたもとで、x_τωに対応する上記観測信号ベクトルy_τωが音源kに対応するクラスタに属する条件付き確率を計算する事後確率計算部と、
   上記条件付き確率と上記ハイパーパラメータφとに基づいて、上記混合重みを更新する混合重み更新手段と、
   上記条件付き確率と上記特徴ベクトルx_τωとに基づいて、音源kに対する相関行列R_kωを計算する相関行列更新手段と、
   上記相関行列R_kωの正規化された主成分ベクトルを新たな値として上記平均方向a_kωを更新する平均方向更新手段と、
   上記相関行列R_kωの最大固有値に基づいて、上記密度パラメータκ_kωを更新する密度パラメータ更新手段と、
   周波数ビンごとに、上記評価関数が最大になるように、上記平均方向a_kωと上記密度パラメータκ_kωとをそれぞれ音源間で並べ替えるパーミュテーション解決手段と、
   を含むモデル推定装置。
6. 請求項５に記載のモデル推定装置であって、
   γ_kτωを上記条件付き確率とし、α_kτを上記混合重みとし、d(τ,ω)を上記観測信号ベクトルy_τωに寄与する音源番号とし、Fを周波数ビンの数とし、・^Hを・のエルミート転置とし、λ_kωを上記相関行列R_kωの最大固有値とし、
   上記事後確率計算部は、次式により上記条件付き確率を計算し、
   

   

   
   上記混合重み更新手段は、次式により求めたα’_kτを新たな値として上記混合重みを更新し、
   

   

   
   上記相関行列更新手段は、次式により求めたR’_kωを新たな値として上記相関行列R_kωを更新し、
   

   

   
   上記密度パラメータ更新手段は、次式により求めたκ’_kωを新たな値として上記密度パラメータκ_kωを更新する
   

   

   
   モデル推定装置。
7. kを音源のインデックスとし、τを時間フレームのインデックスとし、ωを角周波数とし、
   K個の音源からの信号が混合された混合信号をM個のマイクロフォンで観測した観測信号から、各観測信号の時間周波数成分からなる観測信号ベクトルy_τωに対応する特徴ベクトルx_τωを抽出する特徴抽出部と、
   上記特徴ベクトルx_τωを所定の確率モデルにあてはめ、各音源の尤度の時系列が周波数ビン間で同期しているほど高い評価値を与える評価関数を用いて、上記確率モデルのモデルパラメータを推定するモデル推定部と、
   上記モデルパラメータを用いて音源kに対応するマスクm_kτωを求めるマスク作成部と、
   上記観測信号ベクトルy_τωから上記マスクm_kτωを用いて時間周波数領域の分離音を計算する分離音作成部と、
   を含み、
   上記確率モデルは、各音源に関する特徴ベクトルx_τωの分布の重み付き和で表される混合モデルであり、
   上記確率モデルの混合重みは、時間フレームτに依存し、角周波数ωに依存しない重みであり、
   上記確率モデルのモデルパラメータは、上記混合重みと、各音源に関する上記特徴ベクトルx_τωの分布のパラメータである
   音源分離装置。
8. kを音源のインデックスとし、τを時間フレームのインデックスとし、ωを角周波数とし、
   特徴抽出部が、K個の音源からの信号が混合された混合信号をM個のマイクロフォンで観測した観測信号から、各観測信号の時間周波数成分からなる観測信号ベクトルy_τωに対応する特徴ベクトルx_τωを抽出する特徴抽出ステップと、
   モデル推定部が、上記特徴ベクトルx_τωを所定の確率モデルにあてはめ、各音源の尤度の時系列が周波数ビン間で同期しているほど高い評価値を与える評価関数を用いて、上記確率モデルのモデルパラメータを推定するモデル推定ステップと、
   を含み、
   上記確率モデルは、各音源に関する特徴ベクトルx_τωの分布の重み付き和で表される混合モデルであり、
   上記確率モデルの混合重みは、時間フレームτに依存し、角周波数ωに依存しない重みであり、
   上記確率モデルのモデルパラメータは、上記混合重みと、各音源に関する上記特徴ベクトルx_τωの分布のパラメータである
   モデル推定方法。
9. kを音源のインデックスとし、τを時間フレームのインデックスとし、ωを角周波数とし、
   特徴抽出部が、K個の音源からの信号が混合された混合信号をM個のマイクロフォンで観測した観測信号から、各観測信号の時間周波数成分からなる観測信号ベクトルy_τωに対応する特徴ベクトルx_τωを抽出する特徴抽出ステップと、
   モデル推定部が、上記特徴ベクトルx_τωを所定の確率モデルにあてはめ、各音源の尤度の時系列が周波数ビン間で同期しているほど高い評価値を与える評価関数を用いて、上記確率モデルのモデルパラメータを推定するモデル推定ステップと、
   マスク作成部が、上記モデルパラメータを用いて音源kに対応するマスクm_kτωを求めるマスク作成ステップと、
   分離音作成部が、上記観測信号ベクトルy_τωから上記マスクm_kτωを用いて時間周波数領域の分離音を計算する分離音作成ステップと、
   を含み、
   上記確率モデルは、各音源に関する特徴ベクトルx_τωの分布の重み付き和で表される混合モデルであり、
   上記確率モデルの混合重みは、時間フレームτに依存し、角周波数ωに依存しない重みであり、
   上記確率モデルのモデルパラメータは、上記混合重みと、各音源に関する特徴ベクトルx_τωの分布のパラメータである
   音源分離方法。
10. 請求項１から６のいずれかに記載のモデル推定装置としてコンピュータを機能させるためのプログラム。
11. 請求項７に記載の音源分離装置としてコンピュータを機能させるためのプログラム。

Images