JP6290803B2 - モデル推定装置、目的音強調装置、モデル推定方法及びモデル推定プログラム - Google Patents

Description

本発明は、モデル推定装置、目的音強調装置、モデル推定方法及びモデル推定プログラムに関する。

従来から、目的音強調の技術として音源分離技術がある。音源分離技術は、複数のマイクロホンで取得した、複数の音源信号の混合信号を用いて、各音源信号を推定する技術である。特に、クラスタリングに基づく音源分離技術、独立成分分析に基づく音源分離技術がよく知られている。以下、従来技術として、クラスタリングに基づく音源分離技術ついて説明する。以下において、例えばＡがベクトルである場合には“ベクトルＡ”と表記し、例えばＡがスカラーである場合には単に“Ａ”と表記する。また、以下において、特に断らない限り、時間周波数領域での信号表現を用いる。時間フレームの番号をｔ∈｛１，２，・・・，Ｔ｝（Ｔは、フレーム総数）で表し、周波数binの番号をｆ∈｛１，２，・・・，Ｆ｝（Ｆは、ナイキスト周波数以下の周波数binの総数）で表す。

時間周波数領域での信号表現は、時間領域での信号表現に対し、短時間フーリエ変換などの時間周波数変換を適用することで得られる。逆に、時間領域での信号表現は、時間周波数領域での信号表現に対し、逆短時間フーリエ変換などの時間周波数変換の逆変換を適用することで得られる。

Ｎ個（Ｎは、自然数）の音源からの信号をＭ個（Ｍは、自然数）のマイクロホンで観測するとする。ｍ（１≦ｍ≦Ｍ）番目のマイクロホンで観測される混合信号をｙ^(m) _tfで表し、下記（１）式のように、Ｍ個のマイクロホンで観測される混合信号を混合信号ベクトルｙ_tfとしてまとめて表記する。

上記（１）式において、・^Ｔは、・の転置を表す。残響時間がフレーム長に比べて十分短い場合、混合信号ベクトルｙ_tfは、下記（２）式によりモデル化できる。

上記（２）式において、ｃ⁽ⁿ⁾ _tfは、ｎ番目の音源信号を表す。また、上記（２）式におけるベクトルｈ⁽ⁿ⁾ _fは、下記（３）式により定義される。なお、下記（３）式において、ｈ^(m,n) _fは、ｎ番目の音源信号からｍ番目のマイクロホンへの時不変の伝達関数を表す。

ベクトルｈ⁽ⁿ⁾ _fは、ステアリングベクトルと呼ばれ、ｎ番目の音源の位置に関する情報を含む。以下では、簡単のため、マイクロホン数がＭ＝２であり、残響や反響の影響は無視でき、各音源信号は平面波として伝搬すると仮定する。この場合、ベクトルｈ⁽ⁿ⁾ _fは、下記（４）式によりモデル化できる。なお、下記（４）式において、“ｊ”は虚数単位を表す。

ここで、上記（４）式におけるω_fは、周波数binの番号ｆに対応する角周波数を表し、ｄ^(m,n)は、ｍ番目のマイクロホンとｎ番目の音源との距離を表し、ｃは、音速を表す。ｎ番目の音源のマイクロホン間到来時間差δ⁽ⁿ⁾を、下記（５）式により定義する。

すると、ステアリングベクトルｈ⁽ⁿ⁾ _ｆにおけるマイクロホン間位相差arg(ｈ^(1,n) _f)−arg(ｈ^(2,n) _f)（arg(・)は、・の偏角（位相）を表す）と、ｎ番目のマイクロホン間到来時間差δ⁽ⁿ⁾との間には、下記（６）式に示す関係がある。

クラスタリングに基づく音源分離技術では、観測された混合信号ベクトルｙ_tfは、「各時間周波数点では単一の音源成分のみからなる」（以下、「スパース」と表記する）と仮定する（例えば、非特許文献１参照）。スパースは、残響の影響が小さく、音源信号が音声である場合に、精度よく成立することが知られている。スパースの仮定の下では、「時間周波数点（ｔ,ｆ）において混合信号ベクトルｙ_tfに含まれる」（以下、「アクティブ」と表記する）音源の番号をｄ_tfで表すと、上記（２）式は、下記（７）式のように書き換えられる。

スパース性の仮定の下では、観測信号から、下記（８）式の定義に基づき計算される特徴量ｚ_tfは、下記（９）式に示すように、アクティブなｄ_tf番目の音源のマイクロホン間到来時間差と一致する。

よって、ｚ_tfのクラスタリングにより音源分離が実現できる。クラスタリングは、例えば、混合モデルのフィッティングやk-meansクラスタリングなどのクラスタリング技術により行うことができる（例えば、非特許文献２参照）。

O. Yilmaz and S. Rickard, "Blind separation of speech mixtures via time-frequency masking." IEEE Trans. SP, vol. 52, no. 7, pp. 1830-1847, Jul. 2004. S. Araki, H. Sawada, R. Mukai, and S. Makino, "Underdetermined blind sparse source separation for arbitrarily arranged multiple sensors." Signal Processing, vol. 87, no. 8, pp. 1833-1847, Aug. 2007. 伊藤信貴，荒木章子，木下慶介，中谷智広,"音源位置情報に基づく劣決定ブラインド音源分離のためのパーミュテーションフリークラスタリング法"，電子情報通信学会論文誌, vol. J97-A, no. 4, pp. 234-246, Apr. 2014. T. Yoshioka, T. Nakatani, M. Miyoshi, and H.G. Okuno, "Blind separation and dereverberation of speech mixtures by joint optimization." IEEE Trans. ASLP, vol. 19, no. 1, pp. 69-84, Jan. 2011. N.Q.K. Duong, E. Vincent, and R. Gribonval, "Under-determined reverberant audio source separation using a full-rank spatial covariance model." IEEE Trans. ASLP, vol. 18, no. 7, pp. 1830-1840, Sep. 2010. A.P. Dempster, N.M. Laird, and D.B. Rubin, "Maximum likelihood from incomplete data via the EM algorithm." Journal of the Royal Statistical Society: Series B (Methodological), vol. 39, no. 1, pp. 1-38, 1977. H. Sawada, S. Araki, and S. Makino, "Underdetermined convolutive blind source separation via frequency bin-wise clustering and permutation alignment." IEEE Trans. ASLP, vol. 19, no. 3, pp. 516-527, Mar. 2011.

しかしながら、上記従来技術は、フレーム長と比べて残響時間が十分短いことを前提とするため、この前提が成立しない多くの実環境（例えば、会議室など）において、音源分離性能が低下する問題がある。

本願が開示する実施形態の一例は、上記に鑑みてなされたものであって、フレーム長と比べて残響時間が長い場合においても、より高精度な音源分離を実現することを目的とする。

本願の実施形態の一例は、モデル推定装置は、複数の音源が出力する音による残響の特性を示す回帰行列を含む、残響を含む混合信号のモデルのパラメータを保存する記憶部を備える。モデル推定装置は、音を複数のマイクロホンで観測した観測信号と、記憶部に保存される回帰行列とを用いた線形予測により、残響を含まない混合信号を推定する。モデル推定装置は、推定された混合信号を、各時間周波数点が属する音源毎のクラスタにクラスタリングし、記憶部に保存されるパラメータから、各クラスタと対応する事後確率を計算する。モデル推定装置は、推定された混合信号と、計算された事後確率とから、パラメータを推定し、推定したパラメータで記憶部に保存されるパラメータを更新する。モデル推定装置は、信号推定、クラスタリング及びパラメータ推定を、所定条件が満たされるまで繰り返す。

本願が開示する実施形態の一例によれば、例えば、フレーム長と比べて残響時間が長い場合においても、より高精度な音源分離を実現できる。

図１は、実施形態１に係るモデル推定装置の構成の一例を示す図である。 図２は、実施形態１に係るモデル推定装置の処理手順の一例を示すフローチャートである。 図３は、実施形態２に係るモデル推定装置の構成の一例を示す図である。 図４は、実施形態３に係るモデル推定装置の構成の一例を示す図である。 図５は、実施形態３に係るモデル推定装置の処理手順の一例を示すフローチャートである。 図６は、実施形態４に係る目的音強調装置の構成の一例を示す図である。 図７は、実施形態４に係る目的音強調装置の処理手順の一例を示すフローチャートである。 図８は、実施形態４の効果の一例を説明する図である。 図９は、実施形態４の効果の一例を説明する図である。 図１０は、プログラムが実行されることにより、モデル推定装置及び目的音強調装置が実現されるコンピュータの一例を示す図である。

［実施形態］
以下、本願が開示するモデル推定装置、目的音強調装置、モデル推定方法及びモデル推定プログラムの実施形態を説明する。なお、以下の実施形態は、一例を示すに過ぎず、本願が開示する技術を限定するものではない。また、以下に示す各実施形態は、矛盾しない範囲で適宜組合せてもよい。

なお、以下の実施形態では、例えばＡがベクトルである場合には“ベクトルＡ”と表記し、例えばＡがスカラーである場合には単に“Ａ”と表記する。また、例えばＡが集合である場合には、“集合Ａ”と表記するものとする。また、例えばベクトルＡの関数ｆは、ｆ（ベクトルＡ）と表記するものとする。また、ベクトル又はスカラーであるＡに対し、“^〜Ａ”と記載する場合は「“Ａ”の直上に“〜”が記された記号」と同等であるとする。また、ベクトル又はスカラーであるＡに対し、“＾Ａ”と記載する場合は「“Ａ”の直上に“＾”が記された記号」と同等であるとする。また、ベクトル又はスカラーであるＡに対し、“〜＾Ａ”と記載する場合は「“Ａ”の直上に“＾”が記され、さらにその直上に“〜”が付された記号」と同等であるとする。また、ベクトル又はスカラーであるＡに対し、Ａ^TはＡの転置を表す。また、行列Ａに対し、行列Ａ^−１は行列Ａの逆行列を表し、detＡは行列Ａの行列式を表し、trＡは行列Ａの対角和（トレース）を表す。また、行列Ａに対し、行列Ａ^Hは、行列Ａのエルミート転置を表し、行列Ａ^＊は、行列Ａの複素共役を表す。また、集合Ａに対し、＃Ａは集合Ａの要素数を表す。また、exp（・）は、指数関数であり、ln(・)は対数関数である。

［実施形態１］
以下、実施形態１について、実施形態１の理論的背景を説明後、実施形態１の一態様を説明する。

＜実施形態１の理論的背景＞
実施形態１は、残響下で、Ｎ個（Ｎは、自然数）の音源からの信号をＭ個（Ｍは、自然数）のマイクロホンで観測するとする。ｍ（１≦ｍ≦Ｍ）番目のマイクロホンで観測された残響を含む混合信号をｙ^(m) _tfで表し、下記（１０）式のように、Ｍ個のマイクロホンで観測される混合信号を混合信号ベクトルｙ_tfとしてまとめて表記する。

実施形態１のモデル推定装置は、残響を含む混合信号ベクトルｙ_tfを、混合信号ベクトルｙ_tfの分布を表す確率モデルに当てはめ、所定の確率モデルのパラメータを推定する。以下では、先ず、混合信号ベクトルｙ_tfの分布を表す確率モデルについて説明し、次に、混合信号ベクトルｙ_tfの分布を表す確率モデルのパラメータを推定するアルゴリズムを導出する。以下、残響を含む混合信号ベクトルのモデル化、及び、パラメータ推定アルゴリズムの導出それぞれについて、理論的背景を説明する。

（実施形態１の残響を含む混合信号ベクトルのモデル化）
ｎ（１≦ｎ≦Ｎ）番目の音源のみが存在し、残響および他の音源が存在しないと仮定した場合に、Ｍ個のマイクロホンで観測される予定の信号を並べたベクトル（以下、「ｎ番目の音源の残響を含まないマイクロホン像」と表記する）をベクトルｓ⁽ⁿ⁾ _tf∈集合Ｃ^Ｍで表す。ここで、ベクトルｓ⁽ⁿ⁾ _tfは、複素数を要素とするＭ次元のベクトルである。残響が存在しないと仮定した場合に、Ｍ個のマイクロホンで観測される予定の混合信号を並べたベクトル（以下、「残響を含まない混合信号ベクトル」と表記する）をｘ_tf∈集合Ｃ^Ｍで表す。残響を含まない混合信号ベクトルｘ_tfがスパースであると仮定すれば、混合信号ベクトルｘ_tfは、下記（１１）式によりモデル化できる。

従来のクラスタリングに基づく音源分離では、残響を含む混合信号ベクトルｙ_tfがスパースであると仮定するのに対し、実施形態１は、残響を含まない混合信号ベクトルｘ_tfがスパースであると仮定する。これにより、残響下でも正確なモデル化が可能である。上記（１１）式による混合信号ベクトルｘ_tfのモデルに基づき、残響を含まない混合信号ベクトルｘ_tfの分布は、下記（１２）式の混合分布によりモデル化される。

上記（１２）式において、ｐ(ベクトルｓ⁽ⁿ⁾ _tf|Θ)は、ｎ番目の音源の残響を含まないマイクロホン像のベクトルｓ⁽ⁿ⁾ _tfの分布を表す確率モデルを表す。また、上記（１２）式において、Ｐ(ｄ_tf|Θ)は、混合重みと呼ばれ、アクティブな音源の番号ｄ_tfの確率モデルを表す。また、上記（１２）式において、Θは、確率モデルのパラメータの集合を表す。集合Θの定義は、後述する。

一方、残響を含む混合信号ベクトルｙ_tfは、残響を含まない混合信号ベクトルｘ_tfにより駆動されたマルチチャネル自己回帰過程により、下記（１３）式のようにモデル化できる。混合信号ベクトルｙ_tfのモデル化については、文献１「T. Yoshioka, T. Nakatani, M. Miyoshi, and H.G. Okuno, “Blind separation and dereverberation of speech mixtures by joint optimization.” IEEE Trans. ASLP, vol. 19, no. 1, pp. 69.84, Jan. 2011.」に詳述されている。

ここで、上記（１３）式において、ｋはタップ番号を表し、Ｋはタップ数を表し、行列Ｇ_kf∈集合Ｃ^M×Mは、複素数を要素とするＭ行Ｍ列の回帰行列を表し、行列Ｇ^H _kfは、回帰行列Ｇ_kfのエルミート転置を表す。また、上記（１３）式において、Δは、所定の遅延を表すが、好ましくは、音源信号が自己相関を持つ時間（音声の場合、２０〜３０ｍｓ程度）に相当するように設定する。遅延Δを導入することで、推定された回帰行列Ｇ_kfを用いて残響除去を行う際に、音源信号の自己相関が除去されることを防ぐ。また、便宜上、ｔ＜０に対しては、混合信号ベクトルｙ_tf＝０（ゼロベクトル）と定義する。便宜上、上記（１３）式のモデルを確率モデルとして表すと、下記（１４）式を得る。なお、下記（１４）式において、δは、ディラックのデルタ関数である。

上記（１２）式及び上記（１４）式の確率モデルを用いると、残響を含む混合信号ベクトルｙ_tfの分布を表す確率モデルを、下記（１５）式及び下記（１６）式のように導出できる。

残響を含む混合信号ベクトルｙ_tfの分布を表す、上記（１６）式の確率モデルの導出においては、各音源の残響を含まないマイクロホン像のベクトルｓ⁽ⁿ⁾ _tfの分布を表す確率モデルｐ(ベクトルｓ⁽ⁿ⁾ _tf|Θ)と、アクティブな音源の番号ｄ_tfの確率モデルＰ(ｄ_tf|Θ)との具体形について、何の仮定も置いていないことに注意する。すなわち、これらの確率モデルを任意の確率分布によりモデル化しても、残響を含む混合信号ベクトルｙ_tfの分布を表す確率モデルは、上記（１６）式により与えられる。

上記（１６）式によれば、残響を含む混合信号ベクトルｙ_tfの分布を表す確率モデルを定めることは、アクティブな音源の番号ｄ_tfの確率モデルＰ(ｄ_tf|Θ)と、各音源の残響を含まないマイクロホン像のベクトルｓ⁽ⁿ⁾ _tfの分布を表す確率モデルｐ(ベクトルｓ⁽ⁿ⁾ _tf|Θ)とを定めることに帰着することが分かる。これらの確率モデルは、任意の確率分布を用いてモデル化できるが、以下では、実施形態１における、これらのモデル化について説明する。

ｎ番目の音源の残響を含まないマイクロホン像のベクトルｓ⁽ⁿ⁾ _tfの分布を表す確率モデルｐ(ベクトルｓ⁽ⁿ⁾ _tf|Θ)は、例えば、下記（１７）式の時変ガウス分布でモデル化できる。このモデル化については、文献２「N.Q.K. Duong, E. Vincent, and R. Gribonval, “Under-determined reverberant audio source separation using a full-rank spatial covariance model.” IEEE Trans. ASLP, vol. 18, no. 7, pp. 1830.1840, Sep. 2010.」に詳述されている。

ここで、上記（１７）式において、φ⁽ⁿ⁾ _tfは、ベクトルｓ⁽ⁿ⁾ _tfの時変のパワースペクトルをモデル化するパラメータであり、行列Ｂ⁽ⁿ⁾ _fは、ベクトルｓ⁽ⁿ⁾ _tfの時不変の空間共分散行列をモデル化するパラメータである。また、上記（１７）式の右辺は、下記（１８）式により表される複素ガウス分布の確率密度関数である。下記（１８）式は、確率変数がベクトルα、平均がベクトルμ、共分散行列Σである複素ガウス分布の確率密度関数を表す。下記（１８）式において、πは円周率、det(πΣ)は、行列πΣの行列式を表す。

また、実施形態１では、アクティブな音源の番号ｄ_tfの確率モデルＰ(ｄ_tf|Θ)を、周波数依存の混合重みα⁽ⁿ⁾ _ｆを用いて、下記（１９）式によりモデル化する。

実施形態１における、残響を含む混合信号ベクトルｙ_tfの分布を表す確率モデルの具体形は、一般の場合である上記（１６）式に、ｎ番目の音源の残響を含まないマイクロホン像のベクトルｓ⁽ⁿ⁾ _tfの分布を表す確率モデルｐ(ベクトルｓ⁽ⁿ⁾ _tf|Θ)の具体形である上記（１７）式と、アクティブな音源の番号ｄ_tfの確率モデルＰ(ｄ_tf|Θ)の具体形である上記（１９）式とを代入することで、下記（２０）式のように得られる。

ここで、パラメータの集合Θは、具体的には、下記（２１）式により定義される。

（実施形態１のパラメータ推定アルゴリズムの導出）
残響を含む混合信号ベクトルｙ_tfの確率モデルを示す上記（１６）式に基づくと、例えば、最尤法又はＭＡＰ（Maximum A Posteriori）推定法に従って、パラメータの集合Θを推定することができる。

最尤法では、残響を含む混合信号ベクトルｙ_tfの尤度ｐ(Ｙ|Θ)を評価関数とし、尤度ｐ(Ｙ|Θ)を最大化することでパラメータの集合の推定値Θ＝arg max_Θ｛ｐ(Ｙ|Θ)｝を求める。ここで、集合Ｙは、Ｙ:＝｛ベクトルｙ_tf｝_tf:＝｛ベクトルｙ_tf |∀ｔ,ｆ｝と定義する。

一方、ＭＡＰ推定法では、パラメータの集合Θの事後確率ｐ(Θ|Ｙ)を評価関数とし、事後確率ｐ(Θ|Ｙ)を最大化することでパラメータの集合の推定値Θ＝arg max_Θ｛ｐ(Θ|Ｙ)｝を求める。さらに、ベイズの定理より、ｐ(Θ|Ｙ)＝｛ｐ(Ｙ|Θ)ｐ(Θ)｝／ｐ(Ｙ)であることと、ｐ(Ｙ)は定数であることに注意すると、ＭＡＰ推定法によるパラメータの集合Θの推定値は、下記（２２）式のように書きなおせる。なお、下記（２２）式において、ｐ(Θ)はパラメータの集合Θの事前確率を表す。

残響を含む混合信号ベクトルｙ_tfの尤度ｐ(Ｙ|Θ)は、上記（１５）式の左辺に現れる、残響を含む混合信号ベクトルｙ_tfの分布を表す確率モデルを用いて、下記（２３）式で表される。

パラメータの集合Θの事前確率ｐ(Θ)は、任意の確率モデルを用いてモデル化することができるが、例えば一様な分布を用いることができる。一様分布を用いる場合、上記（２２）式に基づく、ＭＡＰ推定法によるパラメータの集合Θの推定値は、最尤推定と一致する。もしくは、混合重みの事前分布として、下記（２４）式のようなディリクレ分布を仮定する。

そして、混合重み以外のパラメータに対しては、一様な事前分布を仮定してもよい。この場合、パラメータの集合Θの事前分布Ｐ(Θ)は、上記（２４）式に示す混合重みの事前分布に比例する。ここで、上記（２４）式におけるψは、ハイパーパラメータと呼ばれる所定の定数である。ψは、任意の正数に設定することができるが、例えばψ＝６００とすればよい。

以下では、上記（２２）式に基づくＭＡＰ推定法により、パラメータの集合Θを推定するためのアルゴリズムの一例として、集合Ｄ:＝｛ｄ_tf｝を隠れ変数とみなしたＥＭ（Expectation-Maximization）アルゴリズムを導出する。なお、ＥＭについては、文献３「A.P. Dempster, N.M. Laird, and D.B. Rubin, “Maximum likelihood from incomplete data via the EM algorithm.” Journal of the Royal Statistical Society: Series B (Methodological), vol. 39, no. 1, pp. 1.38, 1977.」に詳述されている。

ＥＭアルゴリズムとは、以下に定義するＥステップとＭステップを、収束条件が満たされるまで反復するものである。Ｅステップでは、下記（２５）式で定義されるＱ関数：Ｑ（Θ;Θ´）を計算する。

ここで、lnＰ(Ｙ,Ｄ|Θ)は、完全データの集合｛Ｙ,Ｄ｝の対数尤度を表し、Ｐ(Ｄ|Ｙ,Θ´)は、パラメータの集合Θの現在の推定値Θ´に対する集合Ｄの事後確率を表し、<・>_{P(Ｄ|Ｙ,Θ´)}は、Ｐ(Ｄ|Ｙ,Θ´)に関する期待値演算を表す。

一方、Ｍステップでは、Ｑ関数の最大化により、パラメータの集合Θを更新する。ＥＭアルゴリズムの各反復において、評価関数ｐ(Ｙ|Θ)ｐ(Θ)の単調非減少性が保証されている。Ｑ関数の具体形を計算するために、下記（２６）式〜（３０）式のように、ｐ(Ｙ,Ｄ|Θ)、ｐ(Ｄ|Ｙ,Θ´)を求める。

ここで、上記（３０）式におけるγ⁽ⁿ⁾ _tfは、下記（３１）式〜（３３）式で定義する。

ただし、簡単のため、上記（３３）式において、α⁽ⁿ⁾ _f、回帰行列Ｇ_kf、φ⁽ⁿ⁾ _tf、行列Ｂ⁽ⁿ⁾ _fの現在の推定値を、それぞれ、単にα⁽ⁿ⁾ _f、回帰行列Ｇ_kf、φ⁽ⁿ⁾ _tf、行列Ｂ⁽ⁿ⁾ _fと表記した。

上記（２８）式、上記（３０）式を、上記（２５）式へ代入することで、Ｑ関数の具体形が、下記（３４）式、（３５）式のように得られる。

混合重みα⁽ⁿ⁾ _fの更新式は、拘束条件Σ^N _n=1α⁽ⁿ⁾ _f＝１に注意して、ラグランジュの未定乗数法を用いることで得られる。φ⁽ⁿ⁾ _tf及び行列Ｂ⁽ⁿ⁾ _fの更新式は、上記（３５）式に示すＱ関数のφ⁽ⁿ⁾ _tf、及び、行列Ｂ⁽ⁿ⁾ _fの複素共役である(行列Ｂ⁽ⁿ⁾ _f)^＊に関する偏微分を０とおくことで得られる。

回帰行列Ｇ_kfの更新式は、上記（３５）式に示すＱ関数から、回帰行列Ｇ_kfのエルミート転置である行列Ｇ^H _kfに依存する項のみを抜き出すと、下記（３６）式、（３７）式のようになる。

上記（３７）式の、行列〜Ｇ^H _fに関する偏微分を０とおいて整理すると、下記（４１）式のようになる。

上記（４１）式の両辺に、vec作用素を作用させ、下記（ａ）式で示される、行列Ａ、行列Ｂ、行列Ｘについての、クロネッカー積に関する公式を適用すると、下記（４２）式のようになる。

ただし、上記（４２）式において、vec[a₁・・・a_P]及び行列〜Ｇ^H _fを、それぞれ下記（４３）式、（４４）式のように定義する。

よって、上記（４２）式から、下記（４５）式のように、vec[行列〜Ｇ^H _f]が求まる。

＜実施形態１の一態様＞
以下、上述の実施形態１の理論的背景に基づく、実施形態１の一態様を説明する。なお、実施形態１の一態様において、音源数Ｎは既知と仮定する。

（実施形態１に係るモデル推定装置の構成）
図１は、実施形態１に係るモデル推定装置の構成の一例を示す図である。実施形態１に係るモデル推定装置１０Ａは、残響除去処理部１１Ａ、クラスタリング部１２Ａを有する。残響除去処理部１１Ａは、初期化部１１Ａ−１、共分散行列更新部１１Ａ−２、回帰行列更新部１１Ａ−３、残響除去部１１Ａ−４を有する。共分散行列更新部１１Ａ−２及び回帰行列更新部１１Ａ−３及び混合重み更新部１２Ａ−２は、パラメータ推定部の一例である。残響除去部１１Ａ−４は、信号推定部の一例である。事後確率更新部１２Ａ−１は、事後確率計算部の一例である。

初期化部１１Ａ−１は、まず、パラメータの集合Θの初期値を計算する。この初期値は、例えば、以下のように計算することができる。まず、アクティブな音源の番号ｄ_tfの推定値＾ｄ_tfを、残響モデルを含まない従来のクラスタリングに基づく音源分離技術を用いて計算する。残響モデルを含まない従来のクラスタリングに基づく音源分離技術は、文献４「伊藤信貴，荒木章子，木下慶介，中谷智広,“音源位置情報に基づく劣決定ブラインド音源分離のためのパーミュテーションフリークラスタリング法”，電子情報通信学会論文誌, vol. J97-A, no. 4, pp. 234.246, Apr. 2014.」に詳述されている。

次に、初期化部１１Ａ−１は、推定値＾ｄ_tfを用いて、下記（４６）式〜（４９）式により、各パラメータを初期化する。なお、下記（４６）式及び（４８）式における集合Ｃ⁽ⁿ⁾ _fは、Ｃ⁽ⁿ⁾ _f:＝｛ｔ|ｄ_tf＝ｎ｝で定義される行列である。また、下記（４６）式及び（４８）式における＃Ｃ⁽ⁿ⁾ _fは、集合Ｃ⁽ⁿ⁾ _fの要素数を表す。また、下記（４９）式における“tr[・]”は、行列[・]のトレースを表す。

共分散行列更新部１１Ａ−２は、各音源ｎ（ｎ＝１,・・・,Ｎ）の残響を含まないマイクロホン像のベクトルｓ⁽ⁿ⁾ _tfの共分散行列φ⁽ⁿ⁾ _tfＢ⁽ⁿ⁾ _fのパラメータφ⁽ⁿ⁾ _tf及び行列Ｂ⁽ⁿ⁾ _fを、それぞれ下記（５０）式、（５１）式により更新する。

回帰行列更新部１１Ａ−３は、回帰行列Ｇ_kfを、下記（５２）式、（５３）式により更新する。

ここで、上記（５３）式の左辺に現れる行列〜Ｇ_f及び上記（５３）式の右辺に現れるベクトル〜ｙ_t-Δ-1,fは、下記（５４）式、（５５）式のように定義される。

残響除去部１１Ａ−４は、残響を含まない混合信号ベクトルの推定値＾ｘ_tfを、下記（５６）式により更新する。

クラスタリング部１２Ａは、事後確率更新部１２Ａ−１、混合重み更新部１２Ａ−２を有する。事後確率更新部１２Ａ−１は、時間周波数点(ｔ,ｆ)でｎ（ｎ＝１,・・・,Ｎ）番目の音源信号がアクティブである事後確率γ⁽ⁿ⁾ _tfを、下記（５７）式により更新する。なお、γ⁽ⁿ⁾ _tf:＝Ｐ(ｄ_tf＝ｎ|ベクトルｙ_tf,Θ)と定義する。

混合重み更新部１２Ａ−２は、混合重みα⁽ⁿ⁾ _fを、下記（５８）式により更新する。

なお、性能向上のため、モデル推定装置１０Ａの全処理に先立ち、残響を含む混合信号ベクトルｙ_tfに対し、前処理として、下記に示す白色化をおこなってもよい。

なお、実施形態１は、クラスタリング部１２Ａの事後確率更新部１２Ａ−１が、上記（５７）式に基づき、時間周波数点(ｔ,ｆ)でｎ（ｎ＝１,・・・,Ｎ）番目の音源信号がアクティブである事後確率γ⁽ⁿ⁾ _tfを計算するとした。しかし、これに限らず、k-meansクラスタリング等の従来技法を用い、時間周波数点(ｔ,ｆ)でｎ（ｎ＝１,・・・,Ｎ）番目の音源信号がアクティブである事後確率γ⁽ⁿ⁾ _tfを計算するとしてもよい。

（実施形態１に係るモデル推定装置の処理）
図２は、実施形態１に係るモデル推定装置の処理手順の一例を示すフローチャートである。以下に述べるモデル推定装置１０Ａの処理は、所定の収束判定条件が満たされるまで反復される。所定の収束条件は、例えば、「所定の反復回数に達している、又は、事後確率更新部１２Ａ−１、混合重み更新部１２Ａ−２の各更新部のうち１つ以上の更新部による更新前後のパラメータ値の差分が所定の閾値未満である」などとすればよい。

先ず、ステップＳ１１では、初期化部１１Ａ−１は、パラメータの集合Θの初期値を、上記（４６）式〜（４９）式に基づき計算し、モデル推定装置１０Ａの主記憶装置に保存する。次に、ステップＳ１２では、残響除去部１１Ａ−４は、モデル推定装置１０Ａの主記憶装置に現在保存されている回帰行列Ｇ_kfに基づき、上記（５６）式により、残響を含まない混合信号ベクトルの推定値＾ｘ_tfを更新する（“残響除去”処理）。

次に、ステップＳ１３では、事後確率更新部１２Ａ−１は、時間周波数点(ｔ,ｆ)でｎ（ｎ＝１,・・・,Ｎ）番目の音源信号がアクティブである事後確率γ⁽ⁿ⁾ _tfを、上記（５７）式により計算し、モデル推定装置１０Ａの主記憶装置に保存する。また、ステップＳ１３では、混合重み更新部１２Ａ−２は、混合重みα⁽ⁿ⁾ _fを、上記（５８）式により計算し、モデル推定装置１０Ａの主記憶装置に保存する（以上、“クラスタリング”処理）。

次に、モデル推定装置１０Ａは、収束判定条件が満たされているか否かを判定する（ステップＳ１４）。モデル推定装置１０Ａは、収束判定条件が満たされている場合（ステップＳ１４Ｙｅｓ）、処理を終了する。モデル推定装置１０Ａは、収束判定条件が満たされていない場合（ステップＳ１４Ｎｏ）、ステップＳ１５へ処理を移す。

ステップＳ１５では、共分散行列更新部１１Ａ−２は、各音源ｎ（ｎ＝１,・・・,Ｎ）の残響を含まないマイクロホン像のベクトルｓ⁽ⁿ⁾ _tfの共分散行列φ⁽ⁿ⁾ _tfＢ⁽ⁿ⁾ _fのパラメータφ⁽ⁿ⁾ _tf及び行列Ｂ⁽ⁿ⁾ _fを、それぞれ上記（５０）式、（５１）式により計算し、モデル推定装置１０Ａの主記憶装置に更新保存する。また、ステップＳ１５では、回帰行列更新部１１Ａ−３は、共分散行列更新部１１Ａ−２により計算されたパラメータφ⁽ⁿ⁾ _tf及び行列Ｂ⁽ⁿ⁾ _fに基づき、回帰行列Ｇ_kfを、上記（５２）式、（５３）式により計算し、モデル推定装置１０Ａの主記憶装置に更新保存する。

また、ステップＳ１５では、事後確率更新部１２Ａ−１は、モデル推定装置１０Ａの主記憶装置に現在保存されているパラメータの集合Θ、及び、最後に実行したステップＳ１２による残響を含まない混合信号ベクトルの推定値＾ｘ_tfに基づき、上記（５７）式により、事後確率γ⁽ⁿ⁾ _tfを計算し、モデル推定装置１０Ａの主記憶装置に更新保存する。また、ステップＳ１５では、混合重み更新部１２Ａ−２は、事後確率更新部１２Ａ−１により計算された事後確率γ⁽ⁿ⁾ _tfに基づき、上記（５８）式により、混合重みα⁽ⁿ⁾ _fを更新し、モデル推定装置１０Ａの主記憶装置に更新保存する。以上のステップＳ１５の処理が終了すると、モデル推定装置１０Ａは、ステップＳ１２へ処理を移す。

［実施形態２］
以下、実施形態２について、実施形態２の理論的背景を説明後、実施形態２の一態様を説明する。

＜実施形態２の理論的背景＞
実施形態１のように、上記（１９）式に示す周波数依存の混合重みを用いる場合、評価関数である事後確率には、パーミュテーション（置換）の不定性がある。すなわち、｛１,・・・,Ｎ｝上の置換Π_fにより、パラメータの集合Θのα⁽ⁿ⁾ _f、φ⁽ⁿ⁾ _tf、行列Ｂ⁽ⁿ⁾ _tfの順序を、下記（６２）式のように入れ替えた場合を考える。

このとき、下記（６３）式が成り立つ。

すなわち、事後確率を最大化するだけでは、推定されたΘにおける番号ｎは、周波数毎に、異なる音源に対応してしまうというパーミュテーション問題がある。よって、推定されたΘをそのまま用いては、適切に目的音強調を行うことはできない。従って、実施形態１に基づいて目的音強調装置を構成する際には、番号ｎが周波数によらず同一の音源に対応するように置換Π_fを決定する、パーミュテーション解決の処理が別途必要となる。

これに対し、実施形態２のモデル推定装置は、時間依存の混合重みを用いる。これにより、上記文献４に開示されている通り、事後確率の最大化により、パーミュテーション問題を生じずにモデル推定が可能である。

以下、実施形態２の理論的背景を、実施形態１との差異に重点を置きながら説明する。

（実施形態２の残響を含む混合信号ベクトルのモデル化）
実施形態２では、アクティブな音源の番号ｄ_tfの確率モデルＰ(ｄ_tf|Θ)を、周波数依存の混合重みではなく、時間依存の混合重みα⁽ⁿ⁾ _tを用いて、下記（６４）式でモデル化する。

従って、実施形態２における残響を含む混合信号ベクトルｙ_tfの分布を表す確率モデル（上記（１６）式参照）の具体形は、下記（６５）式のように得られる。

パラメータの集合Θは、具体的には、下記（６６）式で表される。

（実施形態２のパラメータ推定アルゴリズムの導出）
ＥＭアルゴリズムにより、事後確率を最大化する点は、実施形態２は、実施形態１と同様である。しかし、実施形態２は、ＥＭアルゴリズムの各反復において、Ｅステップ、Ｍステップの処理に加えて、Ｐ（Permutation）ステップの処理を行う。Ｐステップでは、各周波数binの番号ｆにて、目的関数である事後確率が最大となるように、共分散行列φ⁽ⁿ⁾ _tfＢ⁽ⁿ⁾ _fを音源間で置換することにより、パーミュテーションを解決する。すなわち、Π_fを｛１,・・・,Ｎ｝上の置換として、下記（６７）式〜（６９）式の処理を行う。

なお、Ｅステップ及びＭステップにおける更新式の導出は、実施形態１と同様であるので、説明を省略する。

＜実施形態２の一態様＞
以下、上述の実施形態２の理論的背景に基づく、実施形態２の一態様を説明する。なお、実施形態２の一態様において、音源数Ｎは既知と仮定する。しかし、実施形態２は、真の音源数Ｎ₀が既知でなくても、その上限は分かっていると仮定し、仮定する音源数Ｎを、真の音源数Ｎ₀の上限より大きく設定することで、音源数が既知である場合と同様に実施可能である。

（実施形態２に係るモデル推定装置の構成）
図３は、実施形態２に係るモデル推定装置の構成の一例を示す図である。実施形態２に係るモデル推定装置１０Ｂは、残響除去処理部１１Ｂ、クラスタリング部１２Ｂを有する。残響除去処理部１１Ｂは、初期化部１１Ｂ−１、共分散行列更新部１１Ｂ−２、回帰行列更新部１１Ｂ−３、残響除去部１１Ｂ−４を有する。共分散行列更新部１１Ｂ−２及び回帰行列更新部１１Ｂ−３及び混合重み更新部１２Ｂ−２は、パラメータ推定部の一例である。残響除去部１１Ｂ−４は、信号推定部の一例である。事後確率更新部１２Ｂ−１は、事後確率計算部の一例である。

初期化部１１Ｂ−１は、まず、パラメータの集合Θの初期値を計算する。この初期値は、例えば、以下のように計算することができる。まず、アクティブな音源の番号ｄ_tfの推定値＾ｄ_tfを、実施形態１と同様に、残響モデルを含まない従来のクラスタリングに基づく音源分離技術を用いて計算する。次に、初期化部１１Ｂ−１は、推定値＾ｄ_tfを用いて、上記（４７）式〜（４９）式、及び、下記（７０）式により、各パラメータを初期化する。なお、下記（７０）式における集合〜Ｃ⁽ⁿ⁾ _tは、Ｃ⁽ⁿ⁾ _t:＝｛ｆ|ｄ_tf＝ｎ｝で定義される行列である。また、下記（７０）式における＃Ｃ⁽ⁿ⁾ _tは、集合Ｃ⁽ⁿ⁾ _tの要素数を表す。

共分散行列更新部１１Ｂ−２、回帰行列更新部１１Ｂ−３、残響除去部１１Ｂ−４は、実施形態１の共分散行列更新部１１Ａ−２、回帰行列更新部１１Ａ−３、残響除去部１１Ａ−４とそれぞれ同様である。

クラスタリング部１２Ｂは、事後確率更新部１２Ｂ−１、混合重み更新部１２Ｂ−２、パーミュテーション解決部１２Ｂ−３を有する。事後確率更新部１２Ｂ−１は、時間周波数点(ｔ,ｆ)でｎ（ｎ＝１,・・・,Ｎ）番目の音源信号がアクティブである事後確率γ⁽ⁿ⁾ _tfを、下記（７１）式により更新する。なお、γ⁽ⁿ⁾ _tf:＝Ｐ(ｄ_tf＝ｎ|ベクトルｙ_tf,Θ)と定義する。

混合重み更新部１２Ｂ−２は、混合重みα⁽ⁿ⁾ _tを、下記（７２）式により更新する。

パーミュテーション解決部１２Ｂ−３は、各周波数binの番号ｆにて、目的関数である事後確率が最大となるように、共分散行列φ⁽ⁿ⁾ _tfＢ⁽ⁿ⁾ _fを音源間で置換することにより、パーミュテーションを解決する。すなわち、Π_fを｛１,・・・,Ｎ｝上の置換として、下記（７３）式〜（７５）式により、共分散行列φ⁽ⁿ⁾ _tfＢ⁽ⁿ⁾ _fを置換する。

なお、性能向上のため、モデル推定装置１０Ｂの全処理に先立ち、残響を含む混合信号ベクトルｙ_tfに対し、前処理として、上記（５９）式〜（６１）式に示す白色化をおこなってもよい。

なお、実施形態２は、クラスタリング部１２Ｂの事後確率更新部１２Ｂ−１が、上記（７１）式に基づき、時間周波数点(ｔ,ｆ)でｎ（ｎ＝１,・・・,Ｎ）番目の音源信号がアクティブである事後確率γ⁽ⁿ⁾ _tfを計算するとした。しかし、これに限らず、k-meansクラスタリング等の従来技法を用い、時間周波数点(ｔ,ｆ)でｎ（ｎ＝１,・・・,Ｎ）番目の音源信号がアクティブである事後確率γ⁽ⁿ⁾ _tfを計算するとしてもよい。

［実施形態３］
実施形態３は、実施形態２のモデル推定装置１０Ｂを用いて、上記文献４に記載の音源数推定技術により、音源数も推定する構成にしたものである。実施形態３は、真の音源数Ｎ₀は分からないがその上限は分かっていると仮定し、仮定する音源数Ｎを、真の音源数Ｎ₀の上限より大きく設定する。

（実施形態３に係るモデル推定装置の構成）
図４は、実施形態３に係るモデル推定装置の構成の一例を示す図である。実施形態３に係るモデル推定装置１０Ｃは、実施形態２に係るモデル推定装置１０Ｂと比較して、音源数推定部１３をさらに有する。

音源数推定部１３は、クラスタリング部１２Ｂによって計算されたｎ番目の音源がアクティブである事後確率γ⁽ⁿ⁾ _tfを用いて、番号ｎ＝１,・・・,Ｎのうち、真の音源に対応する番号ｎ(１),・・・, ｎ(Ｎ₀)を判定し、真の音源に対応する番号のパラメータのみを出力する。具体的には、音源数推定部１３は、ｎ番目の音源がアクティブである事後確率γ⁽ⁿ⁾ _tfを用いて、ｎ番目の音源がアクティブである事後確率の総和を、例えば下記（７６）式により算出する。

そして、音源数推定部１３は、各ｎ番目の音源がアクティブである事後確率の総和ρ⁽ⁿ⁾を２つにクラスタリングし、総和の大きい方のクラスタに属するρ⁽ⁿ⁾の番号ｎ＝ｎ(１),・・・, ｎ(＾Ｎ₀)を求め、真の音源に対応する番号とみなす。例えば、音源数推定部１３は、ρ⁽ⁿ⁾に対して、クラスタ数２のk-meansクラスタリングを適用してクラスタリングする。

最後に、音源数推定部１３は、真の音源に対応するｎ＝ｎ(１),・・・, ｎ(＾Ｎ₀)に対応する、下記（７７）式に示すパラメータのみを出力する。なお、下記（７７）式において、ｌ＝１,・・・,＾Ｎ₀である。

（実施形態３に係るモデル推定装置の処理）
図５は、実施形態３に係るモデル推定装置の処理手順の一例を示すフローチャートである。以下に述べるモデル推定装置１０Ｃの処理は、実施形態１又は２と同様の所定の収束判定条件が満たされるまで反復される。

先ず、ステップＳ２１では、初期化部１１Ｂ−１は、パラメータの集合Θの初期値を、上記（４７）式〜（４９）式、及び、（７０）式に基づき計算し、モデル推定装置１０Ｃの主記憶装置に保存する。次に、ステップＳ２２では、残響除去部１１Ｂ−４は、モデル推定装置１０Ｃの主記憶装置に現在保存されている回帰行列Ｇ_kfに基づき、上記（５６）式により、残響を含まない混合信号ベクトルの推定値＾ｘ_tfを更新する（“残響除去”処理）。

次に、ステップＳ２３では、事後確率更新部１２Ｂ−１は、時間周波数点(ｔ,ｆ)でｎ（ｎ＝１,・・・,Ｎ）番目の音源信号がアクティブである事後確率γ⁽ⁿ⁾ _tfを、上記（７１）式により計算し、モデル推定装置１０Ｃの主記憶装置に保存する。また、ステップＳ２３では、混合重み更新部１２Ｂ−２は、混合重みα⁽ⁿ⁾ _tを、上記（７２）式により計算し、モデル推定装置１０Ｃの主記憶装置に保存する（以上、“クラスタリング”処理）。また、ステップＳ２３では、パーミュテーション解決部１２Ｂ−３は、Π_fを｛１,・・・,Ｎ｝上の置換として、上記（７３）式〜（７５）式により、共分散行列φ⁽ⁿ⁾ _tfＢ⁽ⁿ⁾ _fを置換する。

次に、モデル推定装置１０Ｃは、収束判定条件が満たされているか否かを判定する（ステップＳ２４）。モデル推定装置１０Ｃは、収束判定条件が満たされている場合（ステップＳ２４Ｙｅｓ）、ステップＳ２６へ処理を移す。モデル推定装置１０Ｃは、収束判定条件が満たされていない場合（ステップＳ２４Ｎｏ）、ステップＳ２５へ処理を移す。

ステップＳ２５の処理は、図２に示す実施形態１のステップＳ１５の処理と同様である。ステップＳ２６では、音源数推定部１３は、ｎ番目の音源がアクティブである事後確率γ⁽ⁿ⁾ _tfを用いて、真の音源数を推定し、推定結果を出力する。

［実施形態４］
実施形態４に係る目的音強調装置は、実施形態１〜３に係るモデル推定装置１０Ａ〜１０Ｃのいずれかを有する目的音強調装置１００である。

（実施形態４に係る目的音強調装置の構成）
図６は、実施形態４に係る目的音強調装置の構成の一例を示す図である。実施形態４に係る目的音強調装置１００は、周波数領域変換部２０、モデル推定装置１０Ａ（あるいは１０Ｂ又は１０Ｃ）、強調音計算部３０、時間領域変換部４０を有する。

周波数領域変換部２０は、時間領域での残響を含む混合信号ベクトル〜ｙ_τを、短時間フーリエ変換などの時間周波数変換により、時間周波数領域での残響を含む混合信号ベクトルｙ_tfに変換する。ここで、混合信号ベクトル〜ｙ_τは、下記（７８）式により定義される。

ただし、上記（７８）式において、〜ｙ^(m) _τは、時間領域でのｍ（ｍ＝１,・・・,Ｍ）番目のマイクロホンで観測された残響を含む混合信号であり、τはサンプル番号を表す。モデル推定装置１０Ａ（あるいは１０Ｂ又は１０Ｃ）は、パラメータの集合Θと、各音源ｎがアクティブである事後確率γ⁽ⁿ⁾ _tfを計算する。

強調音計算部３０は、周波数領域変換部２０から出力された時間周波数領域での残響を含む混合信号ベクトルｙ_tfと、モデル推定装置１０Ａ（あるいは１０Ｂ又は１０Ｃ）から出力されたパラメータの集合Θと各音源ｎがアクティブである事後確率γ⁽ⁿ⁾ _tfとを用いて、時間周波数領域での各音源の残響を含まないマイクロホン像の推定値＾ｓ⁽ⁿ⁾ _tfを、下記（７９）式及び（８０）式により計算し、出力する。

なお、目的音強調装置１００において実施形態１のモデル推定装置１０Ａを用いる場合は、上記（７９）式及び（８０）式の処理に先立って、γ⁽ⁿ⁾ _tfの番号ｎが周波数によらず同一の音源に対応するように、パーミュテーション解決を行う必要がある。このパーミュテーション解決は、例えば文献５「H. Sawada, S. Araki, and S. Makino, “Underdetermined convolutive blind source separation via frequency bin-wise clustering and permutation alignment.” IEEE Trans. ASLP, vol. 19, no. 3, pp. 516.527, Mar. 2011.」に記載の方法により行うことができる。

時間領域変換部４０は、強調音計算部３０から出力された時間周波数領域での各音源の残響を含まないマイクロホン像の推定値のベクトル＾ｓ⁽ⁿ⁾ _tfに、逆短時間フーリエ変換などの時間周波数変換の逆変換を適用して、時間領域での各音源の残響を含まないマイクロホン像の推定値のベクトル〜＾ｓ⁽ⁿ⁾ _τを計算する。ここで、ベクトル〜＾ｓ⁽ⁿ⁾ _τは、下記（８１）式により定義される。ただし、〜＾ｓ^(m,n) _τは、ベクトル＾ｓ⁽ⁿ⁾ _τの第ｍ要素＾ｓ^(m,n) _tfの逆短時間フーリエ変換である。

なお、強調音計算部３０において、残響除去と音源分離を同時に実現する例を示したが、残響のみを除去するために、時間周波数領域での残響を含まない混合信号の推定値＾ｘ_tfに、逆短時間フーリエ変換などの時間周波数変換の逆変換を適用して、時間領域での残響を含まない混合信号の推定値のベクトル〜＾ｘ_τを得る構成としてもよい。ここで、ベクトル〜＾ｘ_τは、下記（８２）式で定義される。ただし、〜＾ｓ^(m,n) _τは、ベクトル＾ｓ⁽ⁿ⁾ _τの第ｍ要素＾ｓ^(m,n) _tfの逆短時間フーリエ変換である。

（実施形態４に係る目的音強調装置の処理）
図７は、実施形態４に係る目的音強調装置の処理手順の一例を示すフローチャートである。実施形態４に係る目的音強調装置１００において、先ず、ステップＳ３１では、周波数領域変換部２０は、各マイクロホンで観測された信号をそれぞれ時間周波数領域の信号に変換する。次に、ステップＳ３２では、モデル推定装置１０Ａ（あるいは１０Ｂ又は１０Ｃ）は、モデル推定を行う。次に、ステップＳ３３では、強調音計算部３０は、強調音を計算により推定する。次に、ステップＳ３４では、時間領域変換部４０は、強調音計算部３０により推定された強調音を周波数領域から時間領域に変換する。

以下、実施形態４を例に取り、開示の実施形態の実施例及びその効果について説明する。図８及び図９は、実施形態４の効果の一例を説明する図である。実施形態４に係る目的音強調装置１００（以下「提案法」）と、従来の残響モデルを含まないクラスタリングベースの音源分離手法（例えば、文献４に記載の手法、以下「従来法」）の性能を比較する実験をおこなった。ただし、実施形態４に係る目的音強調装置１００のモデル推定装置としては、実施形態２に係るモデル推定装置１０Ｂを用いた。

マイクロホンで観測される残響を含む混合信号は、残響を含まない音声波形に、実験室で計測したインパルス応答（例えば、上述の文献５参照）を畳み込むことにより生成した。図８は、インパルス応答を計測した際のマイクロホンと音源の位置を示す。なお、提案法及び従来法の両方において、パラメータΘの推定に先立って、残響を含む混合信号ベクトルｙ_tfに対し、上記（５９）式〜（６１）式に示す白色化をおこなった。また、音源数Ｎは既知とした。また、他の実験条件は、下記（表１）に示すとおりとした。なお、図８に示す実験室は、４．４５ｍ×３．５５ｍ×（高さ）２．５０ｍの空間であった。また、図８に示すSource1及び2とMicrophone1及び2の、実験室の床面に対する高さは、１．２ｍとした。

提案法及び従来法の性能は、下記（８３）式で定義されるＳＩＲ（Signal-to-Interference Ratio）により評価した。

ここで、〜＾ｓ^(1,n,ν) _τは、〜＾ｓ^(1,n)に含まれるν番目の音源成分を表す。Τ:＝８ｋHz×８ｓ＝６４０００は、サンプリング点の総数を表し、Σ_ν≠ｎは、ｎ以外のνの値に対する総和を表す。

ここで、〜＾ｓ^(1,n,ν) _τの求め方を説明する。観測された残響を含む混合信号ベクトルｙ_tfは、ν番目の音源の残響を含むマイクロホン像のベクトルｘ^(ν) _tfを用いて、下記（８４）式のように分解できる。

従って、ｎ番目の音源の残響を含まないマイクロホン像の推定値のベクトル＾ｓ⁽ⁿ⁾ _tfは、下記（８５）式及び（８６）式のように分解できる。

ここで、上記（８６）式において、＾ｓ^(n,ν) _tfは、＾ｓ⁽ⁿ⁾ _tfに含まれるν番目の音源成分を表す。よって、下記（８７）式により、＾ｓ^(n,ν) _tfを求め、＾ｓ^(n,ν) _tfを逆短時間フーリエ変換して〜＾ｓ^(n,ν) _τを求め、〜＾ｓ^(n,ν) _τの第１要素として〜＾ｓ^(1,n,ν) _τが求まる。

図９に、各残響時間に対し、音声波形の組み合わせを変えて８回の試行を行った際のＳＩＲの平均値をプロットしたグラフを示す。残響時間が最も小さい条件（残響時間１３０ｍｓ程度）では、提案法と従来法は同等の性能を示した。しかし、図９に示すように、残響時間が大きくなるにつれて、従来法に対する提案法の性能改善量が増加する傾向があった。特に、残響時間が３７０ｍｓ程度の場合に、性能改善量は、試行中、最大の約４ｄＢとなった。

以上から、実施形態１〜４は、独立成分分析に基づく音源分離技術と比較して、音源数が未知の場合でも適用できる等の利点があるクラスタリングに基づく音源分離技術において、線形予測に基づく残響除去とクラスタリングに基づく音源分離を交互に反復する。実施形態１〜４は、線形予測に基づく残響除去により推定された残響を含まない混合信号に対して、クラスタリングに基づく音源分離を適用することで、音源分離の性能を向上させることができる。さらに、実施形態１〜４は、改善された音源分離結果を用いることで、残響除去の性能を改善することができる。よって、実施形態１〜４は、残響除去と上記音源分離の反復により、残響時間がフレーム長に比べて長い場合でも、より高精度な音源分離を実現することができる。

（モデル推定装置及び目的音強調装置の装置構成について）
図１、図３、図４に示すモデル推定装置１０Ａ〜１０Ｃ及び図６に示す目的音強調装置１００の各構成要素は機能概念的なものであり、必ずしも物理的に図示のように構成されていることを要さない。すなわち、モデル推定装置１０Ａ〜１０Ｃ及び目的音強調装置１００の機能の分散及び統合の具体的形態は図示のものに限られず、全部又は一部を、各種の負荷や使用状況等に応じて、任意の単位で機能的又は物理的に分散又は統合して構成することができる。

また、モデル推定装置１０Ａ〜１０Ｃ及び目的音強調装置１００において行われる各処理は、全部又は任意の一部が、ＣＰＵ（Central Processing Unit）等の処理装置及び処理装置により解析実行されるプログラムにて実現されてもよい。また、モデル推定装置１０Ａ〜１０Ｃ及び目的音強調装置１００において行われる各処理は、ワイヤードロジックによるハードウェアとして実現されてもよい。

また、実施形態において説明した各処理のうち、自動的に行われるものとして説明した処理の全部又は一部を手動的に行うこともできる。もしくは、実施形態において説明した各処理のうち、手動的に行われるものとして説明した処理の全部又は一部を公知の方法で自動的に行うこともできる。この他、上述及び図示の処理手順、制御手順、具体的名称、各種のデータやパラメータを含む情報については、特記する場合を除いて適宜変更することができる。

（プログラムについて）
図１０は、プログラムが実行されることにより、モデル推定装置及び目的音強調装置が実現されるコンピュータの一例を示す図である。コンピュータ１０００は、例えば、メモリ１０１０、ＣＰＵ１０２０を有する。また、コンピュータ１０００は、ハードディスクドライブインタフェース１０３０、ディスクドライブインタフェース１０４０、シリアルポートインタフェース１０５０、ビデオアダプタ１０６０、ネットワークインタフェース１０７０を有する。コンピュータ１０００において、これらの各部はバス１０８０によって接続される。

メモリ１０１０は、ＲＯＭ（Read Only Memory）１０１１及びＲＡＭ（Random Access Memory）１０１２を含む。ＲＯＭ１０１１は、例えば、ＢＩＯＳ（Basic Input Output System）等のブートプログラムを記憶する。ハードディスクドライブインタフェース１０３０は、ハードディスクドライブ１０３１に接続される。ディスクドライブインタフェース１０４０は、ディスクドライブ１０４１に接続される。例えば磁気ディスクや光ディスク等の着脱可能な記憶媒体が、ディスクドライブ１０４１に挿入される。シリアルポートインタフェース１０５０は、例えばマウス１０５１、キーボード１０５２に接続される。ビデオアダプタ１０６０は、例えばディスプレイ１０６１に接続される。

ハードディスクドライブ１０３１は、例えば、ＯＳ１０９１、アプリケーションプログラム１０９２、プログラムモジュール１０９３、プログラムデータ１０９４を記憶する。すなわち、モデル推定装置１０Ａ〜１０Ｃ及び目的音強調装置１００の各処理を規定するプログラムは、コンピュータ１０００によって実行される指令が記述されたプログラムモジュール１０９３として、例えばハードディスクドライブ１０３１に記憶される。例えば、モデル推定装置１０Ａ〜１０Ｃ及び目的音強調装置１００における機能構成と同様の情報処理を実行するためのプログラムモジュール１０９３が、ハードディスクドライブ１０３１に記憶される。

また、実施形態の処理で用いられる設定データは、プログラムデータ１０９４として、例えばメモリ１０１０やハードディスクドライブ１０３１に記憶される。そして、ＣＰＵ１０２０が、メモリ１０１０やハードディスクドライブ１０３１に記憶されたプログラムモジュール１０９３やプログラムデータ１０９４を必要に応じてＲＡＭ１０１２に読み出して実行する。

なお、プログラムモジュール１０９３やプログラムデータ１０９４は、ハードディスクドライブ１０３１に記憶される場合に限らず、例えば着脱可能な記憶媒体に記憶され、ディスクドライブ１０４１等を介してＣＰＵ１０２０によって読み出されてもよい。あるいは、プログラムモジュール１０９３やプログラムデータ１０９４は、ネットワーク（ＬＡＮ（Local Area Network）、ＷＡＮ（Wide Area Network）等）を介して接続された他のコンピュータに記憶されてもよい。そして、プログラムモジュール１０９３やプログラムデータ１０９４は、ネットワークインタフェース１０７０を介してＣＰＵ１０２０によって読み出されてもよい。

上記実施形態及びその他の実施形態は、本願が開示する技術に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。

１０Ａ、１０Ｂ、１０Ｃ モデル推定装置
１１Ａ、１１Ｂ 残響除去処理部
１１Ａ−１、１１Ｂ−１ 初期化部
１１Ａ−２、１１Ｂ−２ 共分散行列更新部
１１Ａ−３、１１Ｂ−３ 回帰行列更新部
１１Ａ−４、１１Ｂ−４ 残響除去部
１２Ａ、１２Ｂ クラスタリング部
１２Ａ−１、１２Ｂ−１ 事後確率更新部
１２Ａ−２、１２Ｂ−２ 混合重み更新部
１２Ｂ−３ パーミュテーション解決部
１３ 音源数推定部
２０ 周波数領域変換部
３０ 強調音計算部
４０ 時間領域変換部
１００ 目的音強調装置
１０００ コンピュータ
１０１０ メモリ
１０２０ ＣＰＵ

Claims (8)

1. 複数の音源が出力する音による残響の特性を示す回帰行列を含む、該残響を含む混合信号のモデルのパラメータを保存する記憶部と、
   前記音を複数のマイクロホンで観測した観測信号と、前記記憶部に保存される回帰行列とを用いた線形予測により、前記残響を含まない混合信号を推定する信号推定部と、
   前記信号推定部により推定された混合信号から、各時間周波数点が属する前記音源に対応するクラスタ毎の事後確率を計算する事後確率計算部と、
   前記観測信号と、前記信号推定部により推定された混合信号と、前記事後確率計算部により計算された事後確率と、前記記憶部に保存されるパラメータとから、前記パラメータを推定し、推定したパラメータで前記記憶部に保存されるパラメータを更新するパラメータ推定部と
   を備え、
   前記信号推定部、前記事後確率計算部及び前記パラメータ推定部は、所定条件が満たされるまで各処理を繰り返す
   ことを特徴とするモデル推定装置。
2. 前記残響を含む混合信号のモデルは、前記残響を含む混合信号の分布を表す確率モデルであり、
   前記確率モデルは、各前記クラスタに関する前記残響を含む混合信号の分布を表す確率モデルの重み付き和で表される混合モデルであり、
   前記パラメータ推定部は、前記確率モデルを評価する所定の評価関数により前記パラメータを推定する
   ことを特徴とする請求項１に記載のモデル推定装置。
3. 前記所定の評価関数は、前記パラメータ推定部により推定されたパラメータに対する前記残響を含む混合信号の尤度、又は、前記パラメータ推定部により推定されたパラメータの事後確率である
   ことを特徴とする請求項２に記載のモデル推定装置。
4. 前記パラメータ推定部により推定されるパラメータは、各時間周波数点において前記残響を含む混合信号に含まれる前記複数の音源の分布を示す混合重み値を含み、
   前記混合重み値は、前記残響を含む混合信号の周波数毎の混合重み値又は前記残響を含む混合信号の時刻毎の混合重み値である
   ことを特徴とする請求項３に記載のモデル推定装置。
5. 前記パラメータ推定部は、各時間周波数点において前記残響を含む混合信号に含まれる前記複数の音源それぞれと対応する前記事後確率から、該複数の音源のうち該残響を含む混合信号に含まれる音源を推定し、推定した音源に対応するパラメータを前記推定したパラメータとする
   ことを特徴とする請求項４に記載のモデル推定装置。
6. 請求項１〜５のいずれか１つに記載のモデル推定装置により推定された前記パラメータ及び前記事後確率と、時間周波数領域での各前記音源の残響を含む混合信号とから、時間周波数領域での各前記音源の残響を含まない音響信号の推定値を推定して出力する出力部
   を備えることを特徴とする目的音強調装置。
7. モデル推定装置が実行するモデル推定方法であって、
   前記モデル推定装置は、複数の音源が出力する音による残響の特性を示す回帰行列を含む、該残響を含む混合信号のモデルのパラメータを保存する記憶部を備え、
   前記音を複数のマイクロホンで観測した観測信号と、前記記憶部に保存される回帰行列とを用いた線形予測により、前記残響を含まない混合信号を推定する信号推定工程と、
   前記信号推定工程により推定された混合信号から、各時間周波数点が属する前記音源に対応するクラスタ毎の事後確率を計算する事後確率計算工程と、
   前記観測信号と、前記信号推定工程により推定された混合信号と、前記事後確率計算工程により計算された事後確率と、前記記憶部に保存されるパラメータとから、前記パラメータを推定し、推定したパラメータで前記記憶部に保存されるパラメータを更新するパラメータ推定工程と
   を含み、
   前記信号推定工程、前記事後確率計算工程及び前記パラメータ推定工程は、所定条件が満たされるまで繰り返される
   ことを特徴とするモデル推定方法。
8. 請求項１〜５のいずれか１つに記載のモデル推定装置としてコンピュータを機能させるモデル推定プログラム。

Images