(57)【要約】【課題】パーミュテーション問題を生じず、パーミュテーション解決の後処理を行わなくても音源分離ができる。【解決手段】特徴抽出部２は、K個の音源からの信号が混合された混合信号をM個のマイクロフォンで観測した観測信号から各観測信号の時間周波数成分からなる観測信号ベクトルyτωに対応する特徴ベクトルxτωを抽出する。モデル推定部３は、特徴ベクトルxτωを所定の確率モデルにあてはめ各音源の尤度の時系列が周波数ビン間で同期しているほど高い評価値を与える評価関数を用いて確率モデルのモデルパラメータを推定する。確率モデルは各音源に関する特徴ベクトルxτωの分布の重み付き和で表される混合モデルである。確率モデルの混合重みは時間フレームτに依存し角周波数ωに依存しない重みである。確率モデルのモデルパラメータは混合重みと各音源に関する特徴ベクトルxτωの分布のパラメータである。【選択図】図１
2014215385.tif 000002
【特許請求の範囲】
【請求項１】
kを音源のインデックスとし、τを時間フレームのインデックスとし、ωを角周波数とし、
K個の音源からの信号が混合された混合信号をM個のマイクロフォンで観測した観測信号から、各観測信号の時間周波数成分からなる観測信号ベクトルyτωに対応する特徴ベクトルxτωを抽出する特徴抽出部と、
上記特徴ベクトルxτωを所定の確率モデルにあてはめ、各音源の尤度の時系列が周波数ビン間で同期しているほど高い評価値を与える評価関数を用いて、上記確率モデルのモデルパラメータを推定するモデル推定部と、
を含み、
上記確率モデルは、各音源に関する特徴ベクトルxτωの分布の重み付き和で表される混合モデルであり、
上記確率モデルの混合重みは、時間フレームτに依存し、角周波数ωに依存しない重みであり、
上記確率モデルのモデルパラメータは、上記混合重みと、各音源に関する上記特徴ベクトルxτωの分布のパラメータである
モデル推定装置。
【請求項２】
請求項１に記載のモデル推定装置であって、
音源kに関する上記特徴ベクトルxτωの分布は、平均方向をakωとし、密度パラメータをκkωとするワトソン分布であり、
音源kに関する上記特徴ベクトルxτωの分布のパラメータは、上記平均方向akωと上記密度パラメータκkωである
モデル推定装置。
【請求項３】
請求項１または２に記載のモデル推定装置であって、
上記混合重みの事前分布は、音源kに依存しないハイパーパラメータφを各混合重みの指数とする上記混合重みについてのディリクレ分布である
モデル推定装置。
【請求項４】
請求項１から３のいずれかに記載のモデル推定装置であって、
上記モデル推定部は、音源kに関する上記特徴ベクトルxτωの分布と音源kの時間フレームτにおける混合重みとの積に基づいて、特徴ベクトルxτωが与えられたもとで、xτωに対応する上記観測信号ベクトルyτωが音源kに対応するクラスタに属する条件付き確率を計算する事後確率計算部と、
上記条件付き確率と上記ハイパーパラメータφとに基づいて、上記混合重みを更新する混合重み更新手段と、
上記条件付き確率と上記特徴ベクトルxτωとに基づいて、音源kに対する相関行列Rkωを計算する相関行列更新手段と、
上記相関行列Rkωの正規化された主成分ベクトルを新たな値として上記平均方向akωを更新する平均方向更新手段と、
上記相関行列Rkωの最大固有値に基づいて、上記密度パラメータκkωを更新する密度パラメータ更新手段と、
周波数ビンごとに、上記評価関数が最大になるように、上記平均方向akωと上記密度パラメータκkωとをそれぞれ音源間で並べ替えるパーミュテーション解決手段と、
を含むモデル推定装置。
【請求項５】
請求項４に記載のモデル推定装置であって、
γkτωを上記条件付き確率とし、αkτを上記混合重みとし、d(τ,ω)を上記観測信号ベクトルyτωに寄与する音源番号とし、Fを周波数ビンの数とし、・Hを・のエルミート転置とし、λkωを上記相関行列Rkωの最大固有値とし、
上記事後確率計算部は、次式により上記条件付き確率を計算し、
【数４１】
2014215385.tif 000045
上記混合重み更新手段は、次式により求めたα’kτを新たな値として上記混合重みを更新し、
【数４２】
2014215385.tif 000046
上記相関行列更新手段は、次式により求めたR’kωを新たな値として上記相関行列Rkωを更新し、
【数４３】
2014215385.tif 000047
上記密度パラメータ更新手段は、次式により求めたκ’kωを新たな値として上記密度パラメータκkωを更新する
【数４４】
2014215385.tif 000048
モデル推定装置。
【請求項６】
kを音源のインデックスとし、τを時間フレームのインデックスとし、ωを角周波数とし、
K個の音源からの信号が混合された混合信号をM個のマイクロフォンで観測した観測信号から、各観測信号の時間周波数成分からなる観測信号ベクトルyτωに対応する特徴ベクトルxτωを抽出する特徴抽出部と、
上記特徴ベクトルxτωを所定の確率モデルにあてはめ、各音源の尤度の時系列が周波数ビン間で同期しているほど高い評価値を与える評価関数を用いて、上記確率モデルのモデルパラメータを推定するモデル推定部と、
上記モデルパラメータを用いて音源kに対応するマスクmkτωを求めるマスク作成部と、
上記観測信号ベクトルyτωから上記マスクmkτωを用いて時間周波数領域の分離音を計算する分離音作成部と、
を含み、
上記確率モデルは、各音源に関する特徴ベクトルxτωの分布の重み付き和で表される混合モデルであり、
上記確率モデルの混合重みは、時間フレームτに依存し、角周波数ωに依存しない重みであり、
上記確率モデルのモデルパラメータは、上記混合重みと、各音源に関する上記特徴ベクトルxτωの分布のパラメータである
音源分離装置。
【請求項７】
kを音源のインデックスとし、τを時間フレームのインデックスとし、ωを角周波数とし、
特徴抽出部が、K個の音源からの信号が混合された混合信号をM個のマイクロフォンで観測した観測信号から、各観測信号の時間周波数成分からなる観測信号ベクトルyτωに対応する特徴ベクトルxτωを抽出する特徴抽出ステップと、
モデル推定部が、上記特徴ベクトルxτωを所定の確率モデルにあてはめ、各音源の尤度の時系列が周波数ビン間で同期しているほど高い評価値を与える評価関数を用いて、上記確率モデルのモデルパラメータを推定するモデル推定ステップと、
を含み、
上記確率モデルは、各音源に関する特徴ベクトルxτωの分布の重み付き和で表される混合モデルであり、
上記確率モデルの混合重みは、時間フレームτに依存し、角周波数ωに依存しない重みであり、
上記確率モデルのモデルパラメータは、上記混合重みと、各音源に関する上記特徴ベクトルxτωの分布のパラメータである
モデル推定方法。
【請求項８】
kを音源のインデックスとし、τを時間フレームのインデックスとし、ωを角周波数とし、
特徴抽出部が、K個の音源からの信号が混合された混合信号をM個のマイクロフォンで観測した観測信号から、各観測信号の時間周波数成分からなる観測信号ベクトルyτωに対応する特徴ベクトルxτωを抽出する特徴抽出ステップと、
モデル推定部が、上記特徴ベクトルxτωを所定の確率モデルにあてはめ、各音源の尤度の時系列が周波数ビン間で同期しているほど高い評価値を与える評価関数を用いて、上記確率モデルのモデルパラメータを推定するモデル推定ステップと、
マスク作成部が、上記モデルパラメータを用いて音源kに対応するマスクmkτωを求めるマスク作成ステップと、
分離音作成部が、上記観測信号ベクトルyτωから上記マスクmkτωを用いて時間周波数領域の分離音を計算する分離音作成ステップと、
を含み、
上記確率モデルは、各音源に関する特徴ベクトルxτωの分布の重み付き和で表される混合モデルであり、
上記確率モデルの混合重みは、時間フレームτに依存し、角周波数ωに依存しない重みであり、
上記確率モデルのモデルパラメータは、上記混合重みと、各音源に関する特徴ベクトルxτωの分布のパラメータである
音源分離方法。
【請求項９】
請求項１から５のいずれかに記載のモデル推定装置としてコンピュータを機能させるためのプログラム。
【請求項１０】
請求項６に記載の音源分離装置としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
この発明は、複数の音が混合した混合信号をもとに、混合前の原音（音源信号）を復元する音源分離技術に関し、特に、音響系に関する事前情報を用いずに音源分離を行うブラインド音源分離技術のうちクラスタリングに基づく音源分離を行う技術に関する。
【背景技術】
【０００２】
まず、ブラインド音源分離技術全体における、クラスタリングに基づくブラインド音源分離技術の位置づけを説明する。ブラインド音源分離技術には、単一のマイクロフォンを用いる方法と複数のマイクロフォンを備えるマイクロフォンアレイを用いる方法とがある。後者の方法は、音源分離の手がかりとして音源位置の違いを利用できるため、一般に前者の方法よりも高い音源分離性能を達成することができる。複数のマイクロフォンを用いた音源分離技術としては、独立成分分析による方法とクラスタリングに基づく方法とが一般的である。後者の方法は、音源の数がマイクロフォンの数よりも多い場合にも適用できるという利点がある。
【０００３】
クラスタリングに基づく音源分離技術は、混合信号の各時間周波数成分に寄与する音源信号は高々一つであるという仮定（W-Disjoint Orthogonality: WDO）に基づく（非特許文献１参照）。この仮定は、例えば音声の場合のように、各音源信号が、
・ごく少数の時間周波数成分のみが大きく、
・その他の時間周波数成分はすべて０に近い
というスパース性を満たす場合に、特に精度よく成立することが知られている。
【０００４】
ここで、記法を定義する。m番目のマイクロフォンで観測された混合信号の時間周波数変換をymτωで表す。ここで、τはフレーム番号を表し、ωは角周波数を表す。また、各マイクロフォンで観測された混合信号の時間周波数変換を並べたベクトルを式(1)で表す。
【数１】
2014215385.tif 000003
【０００５】
ここで、Mはマイクロフォンの個数であり、・T（上付き文字のT）はベクトル・の転置を表す。また、仮にk番目の音源のみが存在するとしたときに、１番目のマイクロフォンにおいて観測される信号の時間周波数変換をskτωで表す。また、skτωからyτωへの伝達関数をhkωで表す。skτωを１番目のマイクロフォンにおいて定義したため、hkωの第１成分は１に等しいことに注意する。
【０００６】
上述の記法を用いて、混合信号のモデルを記述する。まず、WDOを仮定せず、yτωにすべての音源信号が寄与する一般の場合を考えると、混合信号のモデルは式(2)で与えられる。
【数２】
2014215385.tif 000004
ここで、Kは音源数（既知と仮定する）である。
【０００７】
一方、WDOを仮定する場合、yτωに寄与する音源信号の番号をd(τ,ω)で表すと、
【数３】
2014215385.tif 000005
【０００８】
であるから、式(2)は式(4)のように単純になる。
【数４】
2014215385.tif 000006
【０００９】
WDOの仮定のもとでは、各音源skτωを推定する音源分離の問題は、d(τ,ω)を推定する問題と等価である。以下でその理由を説明する。まず、式(4)の第一要素を抜き出すと式(5)を得る。
【数５】
2014215385.tif 000007
【００１０】
式(5)と式(3)とにより、各音源信号は次のように表せる。
【数６】
2014215385.tif 000008
【００１１】
ただし、mkτωは時間周波数マスクと呼ばれ、式(8)で定義される。
【数７】
2014215385.tif 000009
【００１２】
以上より、d(τ,ω)を推定できれば、これを用いて式(8)によりマスクmkτωを計算でき、このマスクmkτωを式(7)のように混合信号y1τωに乗算することによりskτωを計算することができる。
【００１３】
d(τ,ω)の推定は、通常、各時間周波数成分に対して、音源の方向や位置などの特徴量を抽出して、この特徴量に基づいて時間周波数成分を各音源に対応するクラスにクラスタリングすることにより行われる。このクラスタリングの精度は、クラスタリングに基づく音源分離の性能に大きく影響する。
【００１４】
クラスタリングの従来技術として非特許文献２に記載されている方法が挙げられる。この方法では、混合信号の各時間周波数成分に対して、yτωをノルムが１になるように正規化した、式(9)で定義される特徴量を計算し、この特徴量に基づいて周波数ビンごとに時間周波数成分をクラスタリングする。
【数８】
2014215385.tif 000010
ここで、||・||はベクトルのユークリッドノルムを表す。
【００１５】
この方法では、ある周波数ビンにおける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に対応するマスクmkτωを求めるマスク作成部と、観測信号ベクトルyτωからマスクmkτωを用いて時間周波数領域の分離音を計算する分離音作成部と、を含む。確率モデルは、各音源に関する特徴ベクトルxτωの分布の重み付き和で表される混合モデルであり、確率モデルの混合重みは、時間フレームτに依存し、角周波数ωに依存しない重みであり、確率モデルのモデルパラメータは、混合重みと、各音源に関する特徴ベクトルxτωの分布のパラメータである。
【発明の効果】
【００２１】
この発明によれば、パーミュテーション問題を生じず、二段階の処理を必要としない音源分離技術を実現することができる。これにより、例えば、音源位置などが時間的に変化する時変の環境での音声強調のためのオンライン音源分離を容易に実現することが可能となる。
【図面の簡単な説明】
【００２２】
【図１】モデル推定装置の機能構成を例示する図。
【図２】モデル推定装置の処理フローを例示する図。
【図３】音源分離装置の機能構成を例示する図。
【図４】音源分離装置の処理フローを例示する図。
【図５】実験環境を説明するための図。
【図６】実験結果を示す図。
【図７】実験結果を示す図。
【図８】実験結果を示す図。
【図９】実験結果を示す図。
【発明を実施するための形態】
【００２３】
以下、この発明の実施の形態について詳細に説明する。なお、図面中において同じ機能を有する構成部には同じ番号を付し、重複説明を省略する。
【００２４】
［発明のポイント］
詳細は後述するが、この発明の音源分離技術のポイントは以下の３点である。
１．一般に音源信号は、「音源信号の時間周波数変換の振幅値の時系列{|skτω|}τが、周波数ビン間で類似する」という共通振幅変調の性質をもつ（例えば、非特許文献「G. J. Brown, “Computational Auditory Scene Analysis: A Representational Approach”, Ph.D. thesis, University of Sheffield, 1992.」を参照）。この発明では、この共通振幅変調の性質を、パーミュテーション問題を回避するための手掛かりとして利用できることに着目した。各音源信号に対する、周波数ビン間での振幅変調の共通性を、時変・周波数非依存の混合重みとしてモデル化することにより、パーミュテーションを引き起こさずに、クラスタリングできる。この時変・周波数非依存の混合重みは、従来用いられてきた時不変の混合重みとは一線を画す画期的なアイディアである。
２．混合重みの事前分布としてディリクレ分布を導入する。このディリクレ分布のハイパーパラメータを調整することにより、時変・周波数非依存の混合重みが極端に変動しすぎないように制御し、音源分離性能の低下を防ぐことができる。
３．パラメータ更新の反復ごとに、目的関数が最大となるように平均方向と密度パラメータを音源間で置換する。これにより、混合重みを時変・周波数非依存にすることによって生じる目的関数の局所解への収束を回避し、適切にモデルパラメータを推定することができる。
【００２５】
［第一実施形態］
この発明の第一実施形態は、複数の音源からの信号を複数個のマイクロフォンで観測し、モデルパラメータを推定するモデル推定装置である。
【００２６】
図１を参照して、第一実施形態のモデル推定装置Ａの機能構成例を説明する。モデル推定装置Ａは、周波数領域変換部１、特徴抽出部２及びモデル推定部３を含む。モデル推定部３は、事後確率計算部３１、パラメータ更新部３２及びパラメータ保持部３３を含む。パラメータ更新部３２は、混合重み更新手段３２１、相関行列更新手段３２２、平均方向更新手段３２３、密度パラメータ更新手段３２４及びパーミュテーション解決手段３２５を含む。
【００２７】
図２を参照して、モデル推定装置Ａの動作例を手続きの順に従って説明する。
【００２８】
周波数領域変換部１へM個のマイクロフォンにより観測された時間領域の混合信号~ytが入力される。時間領域の混合信号~ytは式(10)で定義される。
【数９】
2014215385.tif 000011
【００２９】
ここで、tは時間インデックスを表し、・T（上付き文字のT）はベクトル・の転置を表し、~ymtはm（1≦m≦M）番目のマイクロフォンで観測された時間領域の混合信号を表す。
【００３０】
周波数領域変換部１は入力された時間領域の混合信号~ytから短時間フーリエ変換などにより時間周波数領域の観測信号ベクトルyτωを生成し出力する（ステップＳ１）。時間周波数領域の観測信号ベクトルyτωは式(11)で定義される。
【数１０】
2014215385.tif 000012
【００３１】
ここで、τは時間フレームのインデックスを表し、ωは角周波数を表し、ymτωは混合信号~ymtの時間周波数領域での表現である。
【００３２】
特徴抽出部２は、周波数領域変換部１の出力する時間周波数領域の観測信号ベクトルyτωを入力とし、特徴ベクトルxτωを計算し出力する（ステップＳ２）。特徴ベクトルxτωの計算は、時間周波数領域の観測信号ベクトルyτωを正規化することにより行ってもよいし、時間周波数領域の観測信号ベクトルyτωを白色化した後に正規化することによって行ってもよいし、時間周波数領域の観測信号ベクトルyτωを正規化した後に白色化して再度正規化することにより行ってもよい。例えば、時間周波数領域の観測信号ベクトルyτωを正規化することにより特徴ベクトルxτωを計算する場合は式(12)により計算すればよい。
【数１１】
2014215385.tif 000013
【００３３】
また、例えば、時間周波数領域の観測信号ベクトルyτωを白色化した後に正規化する場合には、以下のように特徴ベクトルxτωを計算すればよい。まず、時間周波数領域の観測信号ベクトルyτωを用いて、式(13)により時間周波数領域の観測信号ベクトルyτωの標本相関行列Rωｙを計算する。
【数１２】
2014215385.tif 000014
【００３４】
ここで、Tはフレームの個数であり、・H（上付き文字のH）はエルミート（Hermite）転置である。
【００３５】
次に、標本相関行列Rωｙの固有値と固有ベクトルを計算する。計算した標本相関行列Rωｙの固有値を、大きい順に並べたものをσω1,σω2,…,σωMと表す。したがって、式(14)の関係が成り立つ。
【数１３】
2014215385.tif 000015
【００３６】
ここで、標本相関行列Rωｙはエルミート行列であるから、固有値σω1,σω2,…,σωMはすべて実数であることに注意する。また、固有値σω1,σω2,…,σωMに対応し、正規直交系をなす標本相関行列Rωｙの固有ベクトルをuω1,uω2,…,uωMで表す。ここで、標本相関行列Rωｙはエルミート行列であるから、このような固有ベクトルが存在することに注意する。
【００３７】
次に、行列Σωを式(15)により求め、行列Uωを式(16)により求める。
【数１４】
2014215385.tif 000016
【００３８】
次に、行列Uω及び行列Σωを用いて、時間周波数領域の観測信号ベクトルyτωを白色化したベクトルy’τωを式(17)により計算する。
【数１５】
2014215385.tif 000017
【００３９】
最後に、次式のようにベクトルy’τωをそのノルムで正規化することにより、特徴ベクトルxτωを計算する。
【数１６】
2014215385.tif 000018
【００４０】
モデル推定部３は、特徴ベクトルxτωを特徴ベクトルの分布を表す確率モデルに当てはめ、確率モデルを評価する所定の評価関数を用いて、信号抽出に適した確率モデルのモデルパラメータを計算する。
【００４１】
特徴ベクトルの分布を表す確率モデルは、例えば、以下のようにモデル化されたものである。音源の位置が固定の場合、各周波数ビンにおいて、特徴ベクトルxτωは理想的にはそれぞれの音源kごとに固有の値をとる。ただし、実際には、雑音・残響の影響、モデル化誤差などに起因する変動が存在するため、特徴ベクトルxτωは音源kごとにある値を中心に分布する。そのため、この発明では、音源kに関する特徴ベクトルxτωの分布を、例えば、以下のようにワトソン（Watson）分布でモデル化する。
【００４２】
【数１７】
2014215385.tif 000019
【００４３】
ここで、akωは音源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.（参考文献１）」を参照されたい。ここで、特徴ベクトルの分布が周波数ビンごとに定義されることに注意する。
【００４４】
一般に音源信号は、「音源信号の時間周波数変換の振幅値の時系列{|skτω|}τが、周波数ビン間で類似する」という共通振幅変調の性質をもつ（例えば、「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=1Kα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)で表す混合モデルで与えられる。
【数１８】
2014215385.tif 000020
【００４８】
ここで、Θは、式(21)に示すパラメータ集合である。
【数１９】
2014215385.tif 000021
【００４９】
ここで、{αkτ}kτは式(22)により定義される。
【数２０】
2014215385.tif 000022
【００５０】
他の同様の記法もこれにならって定義される。以降では、αkτを混合重みと呼ぶ。混合重みαkτが急激に変動して音源分離性能が低下することを防ぐため、混合重みαkτの事前分布として式(23)に示すディリクレ（Dirichlet）分布を用いる。これがこの発明の第２のポイントである。
【数２１】
2014215385.tif 000023
【００５１】
ここで、Γはガンマ関数であり、φはハイパーパラメータと呼ばれる。φの値を十分大きく定めることにより、混合重みαkτの変動を抑えることができる。φの値を微調整する必要はないが、例えば、φ=1,10,100,1000などの値を用いることができる。
【００５２】
混合重みαkτ以外のパラメータについては一様な事前分布を仮定する。したがって、p(Θ)=Πτp({αkτ}k)である。
【００５３】
モデル推定部３では、特徴ベクトルxτωを以上のようにモデル化された確率モデルに当てはめ、確率モデルを評価する所定の評価関数を用いて、事後確率及び信号抽出に適したパラメータ集合Θを求める。
【００５４】
以下、モデル推定部３の各部の処理を詳細に説明する。モデル推定部３は、図１に示すとおり、事後確率計算部３１、パラメータ更新部３２及びパラメータ保持部３３を含む。モデル推定部３での処理に先立ち、パラメータ集合Θの初期値をパラメータ保持部３３に用意しておく（ステップＳ０）。この初期値は、例えば、αkτ=1/K、κkω=20とし、akωは{xτω}τωから無作為に選ぶことにより設定することができる。
【００５５】
事後確率計算部３１は、パラメータ保持部３３に記憶されたパラメータ集合Θから事後確率γkτω、すなわち特徴ベクトルxτωが与えられたもとでd(τ,ω)=kとなる条件付き確率を式(24)により計算する（ステップＳ３１）。
【数２２】
2014215385.tif 000024
【００５６】
パラメータ更新部３２は、図１に示すとおり、混合重み更新手段３２１、相関行列更新手段３２２、平均方向更新手段３２３、密度パラメータ更新手段３２４及びパーミュテーション解決手段３２５を含み、現在のパラメータ集合Θを更新して新たなパラメータ集合Θ’を生成する（ステップＳ３２）。
【００５７】
混合重み更新手段３２１は、事後確率γkτωを用いて、式(25)を計算することにより、混合重みαkτを新しい値α’kτに更新する。
【００５８】
【数２３】
2014215385.tif 000025
【００５９】
ここで、Fは周波数ビンの個数を表す。φ=1のとき、α’kτは全周波数ビンにわたる事後確率γkτωの平均値となることがわかる。φの増加とともに、α’kτは定数1/Kに近づく。
【００６０】
相関行列更新手段３２２は、特徴ベクトルxτωと事後確率γkτωを用いて、式(26)を計算することにより、各音源kに対する相関行列Rkωを新しい値R’kωに更新する。
【数２４】
2014215385.tif 000026
【００６１】
平均方向更新手段３２３は、相関行列Rkωの正規化された主成分ベクトルとして、平均方向akωを新しい値a’kωに更新する。
【００６２】
密度パラメータ更新手段３２４は、相関行列Rkωの最大固有値λkωを用いて、密度パラメータкkωを式(27)により新しい値к’kωに更新する。
【数２５】
2014215385.tif 000027
【００６３】
パーミュテーション解決手段３２５は、式(28)〜(30)に示すように、各周波数ビンにおいて、平均方向a’kωと密度パラメータк’kωを、事後確率p(Θ’|{xτω}τω)が最大になるように音源間で置換する（ステップＳ３２５）。これがこの発明の第３のポイントである。
【００６４】
【数２６】
2014215385.tif 000028
ここで、Π:{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)により与えられる。
【数２７】
2014215385.tif 000029
【００６８】
ここで、{xτω}τωは互いに独立であると仮定し、Θに依存しない定数項を無視した。この目的関数を式(34)に示す制約条件のもとで最大化する。
【数２８】
2014215385.tif 000030
【００６９】
目的関数L(Θ)は、パーミュテーション問題がない場合に大きい値を取るため、L(Θ)の最大化によりパーミュテーション問題が回避できる。実際、式(33)の第一項から分かるように、目的関数L(Θ)が大きくなるのは、混合重みαkτが大きい値をとるk、τに対し、音源kに対する尤度（もっともらしさ）p(xτω|d(τ,ω)=k,akω,κkω)が大きい場合である。したがって、L(Θ)の最大化により、音源kに対する尤度の時系列{p(xτω|d(τ,ω)=k,akω,κkω)}τが周波数ビン間で同期する。このことと、上述の「観測信号に寄与する音源インデックスの時系列{d(τ,ω)}τは、周波数ビン間で類似する」という性質を考え合わせると、L(Θ)はパーミュテーション問題がない場合に大きい値を取ることがわかる。EMアルゴリズムで用いる評価関数（Q関数）は式(35)(36)により与えられる。
【００７０】
【数２９】
2014215385.tif 000031
【００７１】
更新後のパラメータ集合Θ'は次式により定義され、
【数３０】
2014215385.tif 000032
【００７２】
Q関数を式(34)の制約のもとで最大にするものとして導かれる。すなわち、混合重みαkτの新たな値α'kτを求める式(25)は、ラグランジュ（Lagrange）の未定乗数法によって、式(37)(38)により導出される。
【数３１】
2014215385.tif 000033
ここで、μはラグランジュの未定乗数である。
【００７３】
平均方向の算出方法は、クーラン・フィッシャー（Courant-Fischer）の定理を式(39)に適用することにより導出される。
【数３２】
2014215385.tif 000034
【００７４】
また、密度パラメータの更新式(27)については、まず∂Q/∂κkω=0より式(40)を得る。
【数３３】
2014215385.tif 000035
【００７５】
ここで、
【数３４】
2014215385.tif 000036
であり、λkωは相関行列Rkωの最大固有値である。上式は、近似的に次のように解くことができる（詳細は上述の参考文献１参照）。
【００７６】
【数３５】
2014215385.tif 000037
【００７７】
パラメータ保持部３３は、パラメータ更新部３２での更新処理により得られたパラメータ集合Θ’を記憶する（ステップＳ３３）。また、事後確率計算部３１での次回の処理の際には、記憶したパラメータ集合Θ’をパラメータ集合Θとして提供する。
【００７８】
ステップＳ３１からステップＳ３３までの処理は、事前に設定した最大反復回数max_iterに達するまで、またはパラメータ更新部３２における各パラメータの更新による変動幅が収束判定の閾値Δよりも小さくなるまで、反復して行う（ステップＳ９１）。最大反復回数max_iter及び閾値Δの具体的な値は、例えば、max_iter=100、Δ=10-10とすることができる。
【００７９】
ステップＳ９１において、モデル推定部３における処理が最大反復回数max_iterに達した場合、または各パラメータの更新による変動幅が閾値Δよりも小さくなった場合、モデル推定部３は反復終了後の事後確率γokτωを出力する。
【００８０】
［第二実施形態］
この発明の第二実施形態は、第一実施形態のモデル推定装置Ａを用いて音源分離装置として構成した実施形態である。
【００８１】
図３を参照して、第二実施形態の音源分離装置Ｂの機能構成例を説明する。音源分離装置Ｂは、第一実施形態のモデル推定装置Ａの各部に加えて、音源分離部４及び時間領域変換部５を含む。音源分離部４は、マスク作成部４１及び分離音作成部４２を含む。
【００８２】
図４を参照して、音源分離装置Ｂの動作例を手続きの順に従って説明する。ステップＳ０からステップＳ９１までの処理は第一実施形態のモデル推定装置Ａの動作例と同様であるので詳細な説明は省略する。
【００８３】
音源分離部４は、周波数領域変換部１の出力する混合音の時間周波数変換yτωと、事後確率計算部３１の出力する反復終了後の事後確率γokτωとを用いて、分離音の時間周波数変換^skτωを推定する。
【００８４】
マスク作成部４１は、反復終了後の事後確率γokτωを用いて、混合音に含まれる各音源を抽出するマスクmkτωを作成する（ステップＳ４１）。まず、マスク作成部４１は、反復終了後の事後確率γokτωを用いて、式(42)によりd(τ,ω)の推定値^d(τ,ω)を計算する。
【数３６】
2014215385.tif 000038
【００８５】
次に、マスク作成部４１は、マスクmkτωを式(43)により計算する。
【数３７】
2014215385.tif 000039
【００８６】
なお、マスク作成部４１は、マスクmkτωを式(44)により求めてもよい。
【数３８】
2014215385.tif 000040
【００８７】
分離音作成部４２は、式(45)により、マスクmkτωを混合音の時間周波数変換y1τωに乗算し、分離音の時間周波数変換^skτωを計算する。
【数３９】
2014215385.tif 000041
【００８８】
時間領域変換部５は、音源kごとに、時間周波数領域の分離信号^skτωを時間領域の分離信号~^sktに変換して出力する（ステップＳ５）。
【００８９】
［実験結果］
この発明の効果を確認するために二通りの実験を行った。
＜実験１＞
１番目の実験は、以下の二点を実証するための実験である。
・本発明の方法により、パーミュテーション解決の後処理を行わなくても、音源分離が可能であること。
・混合重みのディリクレ事前分布が、音源分離性能の低下を防ぐ効果を持つこと。
【００９０】
そのために、以下に示す３つの異なるクラスタリング方法に基づく音源分離手法により特徴ベクトルxτωを音源分離し、その性能を比較した。
１．本発明の方法１（ディリクレ事前分布あり）：時変・周波数非依存の混合重みを持つ混合ワトソンモデルによるクラスタリング方法に基づく音源分離手法である。混合重みの事前分布はディリクレ分布であると仮定し、ディリクレ分布のハイパーパラメータをφ=10,102,103として構成した。
２．本発明の方法２（ディリクレ事前分布なし）：本発明の方法１において、ディリクレ分布のハイパーパラメータをφ=1として、ディリクレ事前分布の効果をなくすように構成した。
３．比較手法：時不変・周波数非依存の混合重み0.5を持つ混合ワトソンモデルによる、周波数ビンごとのクラスタリングの後、非特許文献２に記載されているパーミュテーション解決技術を適用するように構成した。混合重みの事前分布は一様であると仮定した。
【００９１】
特徴ベクトルxτωとしては、観測信号ベクトルyτωを正規化したのち、白色化し、再び正規化したものを用いた。
【００９２】
図５を参照して実験環境を説明する。直方体の空間の中央付近に正三角形を描くように３つのマイクロフォンが配置され、それらのマイクロフォンを取り囲む円を描くように４つの音源が配置される。４つの音源は図５の紙面下方向を0°として反時計回りに70°、150°、245°、315°の位置に配置されている。３つのマイクロフォンの描く正三角形の中心から４つの音源までの距離は1.2メートルとした。３つのマイクロフォンと４つの音源は水平に設置され、その床面からの高さは1.2メートルとした。実験では、図５中に太字で描いた紙面右側の２つの音源（70°及び150°の位置）と、紙面左側の２つのマイクロフォン（m=1,2）を用いた。
【００９３】
表１にその他の条件を含めた実験条件をまとめる。
【表１】
2014215385.tif 000042
【００９４】
音源分離性能の評価尺度としては、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を示す。
【表２】
2014215385.tif 000043
【００９６】
本発明の各方法及び比較手法の実験結果により同等のSDRが得られたことがわかる。特に、本発明の方法１においてφ=103とすると、残響時間が130,200,250,300,370ミリ秒である場合には、比較手法と同じか0.1〜0.4dB大きいSDRが得られた。残響時間が440ミリ秒である場合には比較手法のSDRの方が高かったが、その差はわずか0.1dBであった。この結果から、本発明の方法により、パーミュテーション解決の後処理を行わなくても良好に音源分離ができていることがわかる。
【００９７】
また、本発明の方法１と本発明の方法２とを比較すると、前者の方が高いSDRを与えた。特に、本発明の方法１において、φ=103とすると、本発明の方法２と比べて、0.8〜1.3dB大きいSDRが得られており、比較手法と同等以上のSDRが得られた。このことから、混合重みのディリクレ事前分布が、音源分離性能の低下を防ぐ効果を持つことがわかる。
【００９８】
＜実験２＞
２番目の実験は、パラメータ更新部におけるパーミュテーション解決手段の有効性を検証するための実験である。これを検証するために、次の二つの方法によるパラメータ推定を比較した。
１．本発明の方法１：実験１で用いた本発明の方法１において、ディリクレ分布のハイパーパラメータをφ=103として構成した。
２．本発明の方法３：本発明の方法１において、パーミュテーション解決手段３２５における処理を行わない。
【００９９】
平均方向の初期値としては、1〜2kHzの範囲において、真の値を音源間で交換することで得られる、パーミュテーションが生じている初期値を用いた。ここで、真の値は、仮にそれぞれの音源を別個に観測できるとして、それを用いて計算した。他の実験条件は、実験１におけるものと同一である。
【０１００】
本発明の各方法と比較手法とによる平均方向の推定値を比較するために、次式に示す平均方向のマイクロフォン間位相差をプロットした。
【０１０１】
【数４０】
2014215385.tif 000044
ここで、[・]mは、ベクトルの第m要素を表す。
【０１０２】
実験２の実験結果を図６〜９に示す。実験結果は、位相差δkω（縦軸、単位は度（°））を周波数ω/(2π)（横軸、単位はキロヘルツ（kHz））の関数としてプロットした。各クラスタkに対するプロットを異なるマーカで示しており、クラスタk=1は「＋」でプロットし、クラスタk=2は「×」でプロットした。
【０１０３】
図６は、本発明の方法３の反復終了後の値である。図７は、本発明の方法１の反復終了後の値である。図８は、初期値である。図９は、真の値である。
【０１０４】
図６に示す本発明の方法３の反復終了後のマイクロフォン間位相差は、図８に示す初期値の場合と非常に近い値をとることがわかる。EMアルゴリズムは、局所解への収束性が保証されたアルゴリズムであることから、上述の結果は、パーミュテーションの生じている初期値が目的関数の局所解となっていることを示している。
【０１０５】
一方、図７に示す本発明の方法１の反復終了後のマイクロフォン間位相差は、図９に示す真の値の場合と非常に近い値をとることがわかる。
【０１０６】
この結果より、本発明のパーミュテーション解決手段はパーミュテーション問題を回避するのに有効であり、本発明の方法においてパーミュテーション解決手段を用いないとパーミュテーション問題の生じている局所解に陥ってしまうことがわかる。
【０１０７】
以上に説明したとおり、実験１の結果から、この発明によりパーミュテーション解決の後処理を行わなくても音源分離が可能であること、及び、混合重みのディリクレ事前分布が音源分離性能の低下を防ぐ効果を持つことが確認された。また、実験２の結果から、パラメータ更新部におけるパーミュテーション解決手段の有効性が確認された。
【０１０８】
［プログラム、記録媒体］
この発明は上述の実施形態に限定されるものではなく、この発明の趣旨を逸脱しない範囲で適宜変更が可能であることはいうまでもない。上記実施例において説明した各種の処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。
【０１０９】
また、上記実施形態で説明した各装置における各種の処理機能をコンピュータによって実現する場合、各装置が有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記各装置における各種の処理機能がコンピュータ上で実現される。
【０１１０】
この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。
【０１１１】
また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。
【０１１２】
このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。
【０１１３】
また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、本装置を構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。
【符号の説明】
【０１１４】
Ａ  モデル推定装置
Ｂ  音源分離装置
１  周波数領域変換部
２  特徴抽出部
３  モデル推定部
３１  事後確率計算部
３２  パラメータ更新部
３２１  混合重み更新手段
３２２  相関行列更新手段
３２３  平均方向更新手段
３２４  密度パラメータ更新手段
３２５  パーミュテーション解決手段
３３  パラメータ保持部
４  音源分離部
４１  マスク作成部
４２  分離音作成部
５  時間領域変換部
【図１】
2014215385.tif 000049
【図２】
2014215385.tif 000050
【図３】
2014215385.tif 000051
【図４】
2014215385.tif 000052
【図５】
2014215385.tif 000053
【図６】
2014215385.tif 000054
【図７】
2014215385.tif 000055
【図８】
2014215385.tif 000056
【図９】
2014215385.tif 000057