(57)【特許請求の範囲】
【請求項１】
複数の音源が出力する音による残響の特性を示す回帰行列を含む、該残響を含む混合信号のモデルのパラメータを保存する記憶部と、
前記音を複数のマイクロホンで観測した観測信号と、前記記憶部に保存される回帰行列とを用いた線形予測により、前記残響を含まない混合信号を推定する信号推定部と、
前記信号推定部により推定された混合信号から、各時間周波数点が属する前記音源に対応するクラスタ毎の事後確率を計算する事後確率計算部と、
前記観測信号と、前記信号推定部により推定された混合信号と、前記事後確率計算部により計算された事後確率と、前記記憶部に保存されるパラメータとから、前記パラメータを推定し、推定したパラメータで前記記憶部に保存されるパラメータを更新するパラメータ推定部と
を備え、
前記信号推定部、前記事後確率計算部及び前記パラメータ推定部は、所定条件が満たされるまで各処理を繰り返す
ことを特徴とするモデル推定装置。
【請求項２】
前記残響を含む混合信号のモデルは、前記残響を含む混合信号の分布を表す確率モデルであり、
前記確率モデルは、各前記クラスタに関する前記残響を含む混合信号の分布を表す確率モデルの重み付き和で表される混合モデルであり、
前記パラメータ推定部は、前記確率モデルを評価する所定の評価関数により前記パラメータを推定する
ことを特徴とする請求項１に記載のモデル推定装置。
【請求項３】
前記所定の評価関数は、前記パラメータ推定部により推定されたパラメータに対する前記残響を含む混合信号の尤度、又は、前記パラメータ推定部により推定されたパラメータの事後確率である
ことを特徴とする請求項２に記載のモデル推定装置。
【請求項４】
前記パラメータ推定部により推定されるパラメータは、各時間周波数点において前記残響を含む混合信号に含まれる前記複数の音源の分布を示す混合重み値を含み、
前記混合重み値は、前記残響を含む混合信号の周波数毎の混合重み値又は前記残響を含む混合信号の時刻毎の混合重み値である
ことを特徴とする請求項３に記載のモデル推定装置。
【請求項５】
前記パラメータ推定部は、各時間周波数点において前記残響を含む混合信号に含まれる前記複数の音源それぞれと対応する前記事後確率から、該複数の音源のうち該残響を含む混合信号に含まれる音源を推定し、推定した音源に対応するパラメータを前記推定したパラメータとする
ことを特徴とする請求項４に記載のモデル推定装置。
【請求項６】
請求項１〜５のいずれか１つに記載のモデル推定装置により推定された前記パラメータ及び前記事後確率と、時間周波数領域での各前記音源の残響を含む混合信号とから、時間周波数領域での各前記音源の残響を含まない音響信号の推定値を推定して出力する出力部
を備えることを特徴とする目的音強調装置。
【請求項７】
モデル推定装置が実行するモデル推定方法であって、
前記モデル推定装置は、複数の音源が出力する音による残響の特性を示す回帰行列を含む、該残響を含む混合信号のモデルのパラメータを保存する記憶部を備え、
前記音を複数のマイクロホンで観測した観測信号と、前記記憶部に保存される回帰行列とを用いた線形予測により、前記残響を含まない混合信号を推定する信号推定工程と、
前記信号推定工程により推定された混合信号から、各時間周波数点が属する前記音源に対応するクラスタ毎の事後確率を計算する事後確率計算工程と、
前記観測信号と、前記信号推定工程により推定された混合信号と、前記事後確率計算工程により計算された事後確率と、前記記憶部に保存されるパラメータとから、前記パラメータを推定し、推定したパラメータで前記記憶部に保存されるパラメータを更新するパラメータ推定工程と
を含み、
前記信号推定工程、前記事後確率計算工程及び前記パラメータ推定工程は、所定条件が満たされるまで繰り返される
ことを特徴とするモデル推定方法。
【請求項８】
請求項１〜５のいずれか１つに記載のモデル推定装置としてコンピュータを機能させるモデル推定プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、モデル推定装置、目的音強調装置、モデル推定方法及びモデル推定プログラムに関する。
【背景技術】
【０００２】
従来から、目的音強調の技術として音源分離技術がある。音源分離技術は、複数のマイクロホンで取得した、複数の音源信号の混合信号を用いて、各音源信号を推定する技術である。特に、クラスタリングに基づく音源分離技術、独立成分分析に基づく音源分離技術がよく知られている。以下、従来技術として、クラスタリングに基づく音源分離技術ついて説明する。以下において、例えばＡがベクトルである場合には“ベクトルＡ”と表記し、例えばＡがスカラーである場合には単に“Ａ”と表記する。また、以下において、特に断らない限り、時間周波数領域での信号表現を用いる。時間フレームの番号をｔ∈｛１，２，・・・，Ｔ｝（Ｔは、フレーム総数）で表し、周波数binの番号をｆ∈｛１，２，・・・，Ｆ｝（Ｆは、ナイキスト周波数以下の周波数binの総数）で表す。
【０００３】
時間周波数領域での信号表現は、時間領域での信号表現に対し、短時間フーリエ変換などの時間周波数変換を適用することで得られる。逆に、時間領域での信号表現は、時間周波数領域での信号表現に対し、逆短時間フーリエ変換などの時間周波数変換の逆変換を適用することで得られる。
【０００４】
Ｎ個（Ｎは、自然数）の音源からの信号をＭ個（Ｍは、自然数）のマイクロホンで観測するとする。ｍ（１≦ｍ≦Ｍ）番目のマイクロホンで観測される混合信号をｙ(m)tfで表し、下記（１）式のように、Ｍ個のマイクロホンで観測される混合信号を混合信号ベクトルｙtfとしてまとめて表記する。
【０００５】
【数１】
0006290803.tif 000002
【０００６】
上記（１）式において、・Ｔは、・の転置を表す。残響時間がフレーム長に比べて十分短い場合、混合信号ベクトルｙtfは、下記（２）式によりモデル化できる。
【０００７】
【数２】
0006290803.tif 000003
【０００８】
上記（２）式において、ｃ(n)tfは、ｎ番目の音源信号を表す。また、上記（２）式におけるベクトルｈ(n)fは、下記（３）式により定義される。なお、下記（３）式において、ｈ(m,n)fは、ｎ番目の音源信号からｍ番目のマイクロホンへの時不変の伝達関数を表す。
【０００９】
【数３】
0006290803.tif 000004
【００１０】
ベクトルｈ(n)fは、ステアリングベクトルと呼ばれ、ｎ番目の音源の位置に関する情報を含む。以下では、簡単のため、マイクロホン数がＭ＝２であり、残響や反響の影響は無視でき、各音源信号は平面波として伝搬すると仮定する。この場合、ベクトルｈ(n)fは、下記（４）式によりモデル化できる。なお、下記（４）式において、“ｊ”は虚数単位を表す。
【００１１】
【数４】
0006290803.tif 000005
【００１２】
ここで、上記（４）式におけるωfは、周波数binの番号ｆに対応する角周波数を表し、ｄ(m,n)は、ｍ番目のマイクロホンとｎ番目の音源との距離を表し、ｃは、音速を表す。ｎ番目の音源のマイクロホン間到来時間差δ(n)を、下記（５）式により定義する。
【００１３】
【数５】
0006290803.tif 000006
【００１４】
すると、ステアリングベクトルｈ(n)ｆにおけるマイクロホン間位相差arg(ｈ(1,n)f)−arg(ｈ(2,n)f)（arg(・)は、・の偏角（位相）を表す）と、ｎ番目のマイクロホン間到来時間差δ(n)との間には、下記（６）式に示す関係がある。
【００１５】
【数６】
0006290803.tif 000007
【００１６】
クラスタリングに基づく音源分離技術では、観測された混合信号ベクトルｙtfは、「各時間周波数点では単一の音源成分のみからなる」（以下、「スパース」と表記する）と仮定する（例えば、非特許文献１参照）。スパースは、残響の影響が小さく、音源信号が音声である場合に、精度よく成立することが知られている。スパースの仮定の下では、「時間周波数点（ｔ,ｆ）において混合信号ベクトルｙtfに含まれる」（以下、「アクティブ」と表記する）音源の番号をｄtfで表すと、上記（２）式は、下記（７）式のように書き換えられる。
【００１７】
【数７】
0006290803.tif 000008
【００１８】
スパース性の仮定の下では、観測信号から、下記（８）式の定義に基づき計算される特徴量ｚtfは、下記（９）式に示すように、アクティブなｄtf番目の音源のマイクロホン間到来時間差と一致する。
【００１９】
【数８】
0006290803.tif 000009
【００２０】
【数９】
0006290803.tif 000010
【００２１】
よって、ｚ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としてまとめて表記する。
【００３２】
【数１０】
0006290803.tif 000011
【００３３】
実施形態１のモデル推定装置は、残響を含む混合信号ベクトルｙtfを、混合信号ベクトルｙtfの分布を表す確率モデルに当てはめ、所定の確率モデルのパラメータを推定する。以下では、先ず、混合信号ベクトルｙtfの分布を表す確率モデルについて説明し、次に、混合信号ベクトルｙtfの分布を表す確率モデルのパラメータを推定するアルゴリズムを導出する。以下、残響を含む混合信号ベクトルのモデル化、及び、パラメータ推定アルゴリズムの導出それぞれについて、理論的背景を説明する。
【００３４】
（実施形態１の残響を含む混合信号ベクトルのモデル化）
ｎ（１≦ｎ≦Ｎ）番目の音源のみが存在し、残響および他の音源が存在しないと仮定した場合に、Ｍ個のマイクロホンで観測される予定の信号を並べたベクトル（以下、「ｎ番目の音源の残響を含まないマイクロホン像」と表記する）をベクトルｓ(n)tf∈集合ＣＭで表す。ここで、ベクトルｓ(n)tfは、複素数を要素とするＭ次元のベクトルである。残響が存在しないと仮定した場合に、Ｍ個のマイクロホンで観測される予定の混合信号を並べたベクトル（以下、「残響を含まない混合信号ベクトル」と表記する）をｘtf∈集合ＣＭで表す。残響を含まない混合信号ベクトルｘtfがスパースであると仮定すれば、混合信号ベクトルｘtfは、下記（１１）式によりモデル化できる。
【００３５】
【数１１】
0006290803.tif 000012
【００３６】
従来のクラスタリングに基づく音源分離では、残響を含む混合信号ベクトルｙtfがスパースであると仮定するのに対し、実施形態１は、残響を含まない混合信号ベクトルｘtfがスパースであると仮定する。これにより、残響下でも正確なモデル化が可能である。上記（１１）式による混合信号ベクトルｘtfのモデルに基づき、残響を含まない混合信号ベクトルｘtfの分布は、下記（１２）式の混合分布によりモデル化される。
【００３７】
【数１２】
0006290803.tif 000013
【００３８】
上記（１２）式において、ｐ(ベクトルｓ(n)tf|Θ)は、ｎ番目の音源の残響を含まないマイクロホン像のベクトルｓ(n)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.」に詳述されている。
【００４０】
【数１３】
0006290803.tif 000014
【００４１】
ここで、上記（１３）式において、ｋはタップ番号を表し、Ｋはタップ数を表し、行列Ｇkf∈集合ＣM×Mは、複素数を要素とするＭ行Ｍ列の回帰行列を表し、行列ＧHkfは、回帰行列Ｇkfのエルミート転置を表す。また、上記（１３）式において、Δは、所定の遅延を表すが、好ましくは、音源信号が自己相関を持つ時間（音声の場合、２０〜３０ｍｓ程度）に相当するように設定する。遅延Δを導入することで、推定された回帰行列Ｇkfを用いて残響除去を行う際に、音源信号の自己相関が除去されることを防ぐ。また、便宜上、ｔ＜０に対しては、混合信号ベクトルｙtf＝０（ゼロベクトル）と定義する。便宜上、上記（１３）式のモデルを確率モデルとして表すと、下記（１４）式を得る。なお、下記（１４）式において、δは、ディラックのデルタ関数である。
【００４２】
【数１４】
0006290803.tif 000015
【００４３】
上記（１２）式及び上記（１４）式の確率モデルを用いると、残響を含む混合信号ベクトルｙtfの分布を表す確率モデルを、下記（１５）式及び下記（１６）式のように導出できる。
【００４４】
【数１５】
0006290803.tif 000016
【００４５】
残響を含む混合信号ベクトルｙtfの分布を表す、上記（１６）式の確率モデルの導出においては、各音源の残響を含まないマイクロホン像のベクトルｓ(n)tfの分布を表す確率モデルｐ(ベクトルｓ(n)tf|Θ)と、アクティブな音源の番号ｄtfの確率モデルＰ(ｄtf|Θ)との具体形について、何の仮定も置いていないことに注意する。すなわち、これらの確率モデルを任意の確率分布によりモデル化しても、残響を含む混合信号ベクトルｙtfの分布を表す確率モデルは、上記（１６）式により与えられる。
【００４６】
上記（１６）式によれば、残響を含む混合信号ベクトルｙtfの分布を表す確率モデルを定めることは、アクティブな音源の番号ｄtfの確率モデルＰ(ｄtf|Θ)と、各音源の残響を含まないマイクロホン像のベクトルｓ(n)tfの分布を表す確率モデルｐ(ベクトルｓ(n)tf|Θ)とを定めることに帰着することが分かる。これらの確率モデルは、任意の確率分布を用いてモデル化できるが、以下では、実施形態１における、これらのモデル化について説明する。
【００４７】
ｎ番目の音源の残響を含まないマイクロホン像のベクトルｓ(n)tfの分布を表す確率モデルｐ(ベクトルｓ(n)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.」に詳述されている。
【００４８】
【数１６】
0006290803.tif 000017
【００４９】
ここで、上記（１７）式において、φ(n)tfは、ベクトルｓ(n)tfの時変のパワースペクトルをモデル化するパラメータであり、行列Ｂ(n)fは、ベクトルｓ(n)tfの時不変の空間共分散行列をモデル化するパラメータである。また、上記（１７）式の右辺は、下記（１８）式により表される複素ガウス分布の確率密度関数である。下記（１８）式は、確率変数がベクトルα、平均がベクトルμ、共分散行列Σである複素ガウス分布の確率密度関数を表す。下記（１８）式において、πは円周率、det(πΣ)は、行列πΣの行列式を表す。
【００５０】
【数１７】
0006290803.tif 000018
【００５１】
また、実施形態１では、アクティブな音源の番号ｄtfの確率モデルＰ(ｄtf|Θ)を、周波数依存の混合重みα(n)ｆを用いて、下記（１９）式によりモデル化する。
【００５２】
【数１８】
0006290803.tif 000019
【００５３】
実施形態１における、残響を含む混合信号ベクトルｙtfの分布を表す確率モデルの具体形は、一般の場合である上記（１６）式に、ｎ番目の音源の残響を含まないマイクロホン像のベクトルｓ(n)tfの分布を表す確率モデルｐ(ベクトルｓ(n)tf|Θ)の具体形である上記（１７）式と、アクティブな音源の番号ｄtfの確率モデルＰ(ｄtf|Θ)の具体形である上記（１９）式とを代入することで、下記（２０）式のように得られる。  
【００５４】
【数１９】
0006290803.tif 000020
【００５５】
ここで、パラメータの集合Θは、具体的には、下記（２１）式により定義される。
【００５６】
【数２０】
0006290803.tif 000021
【００５７】
（実施形態１のパラメータ推定アルゴリズムの導出）
残響を含む混合信号ベクトルｙtfの確率モデルを示す上記（１６）式に基づくと、例えば、最尤法又はＭＡＰ（Maximum  A  Posteriori）推定法に従って、パラメータの集合Θを推定することができる。
【００５８】
最尤法では、残響を含む混合信号ベクトルｙtfの尤度ｐ(Ｙ|Θ)を評価関数とし、尤度ｐ(Ｙ|Θ)を最大化することでパラメータの集合の推定値Θ＝arg  maxΘ｛ｐ(Ｙ|Θ)｝を求める。ここで、集合Ｙは、Ｙ:＝｛ベクトルｙtf｝tf:＝｛ベクトルｙtf  |∀ｔ,ｆ｝と定義する。
【００５９】
一方、ＭＡＰ推定法では、パラメータの集合Θの事後確率ｐ(Θ|Ｙ)を評価関数とし、事後確率ｐ(Θ|Ｙ)を最大化することでパラメータの集合の推定値Θ＝arg  maxΘ｛ｐ(Θ|Ｙ)｝を求める。さらに、ベイズの定理より、ｐ(Θ|Ｙ)＝｛ｐ(Ｙ|Θ)ｐ(Θ)｝／ｐ(Ｙ)であることと、ｐ(Ｙ)は定数であることに注意すると、ＭＡＰ推定法によるパラメータの集合Θの推定値は、下記（２２）式のように書きなおせる。なお、下記（２２）式において、ｐ(Θ)はパラメータの集合Θの事前確率を表す。
【００６０】
【数２１】
0006290803.tif 000022
【００６１】
残響を含む混合信号ベクトルｙtfの尤度ｐ(Ｙ|Θ)は、上記（１５）式の左辺に現れる、残響を含む混合信号ベクトルｙtfの分布を表す確率モデルを用いて、下記（２３）式で表される。
【００６２】
【数２２】
0006290803.tif 000023
【００６３】
パラメータの集合Θの事前確率ｐ(Θ)は、任意の確率モデルを用いてモデル化することができるが、例えば一様な分布を用いることができる。一様分布を用いる場合、上記（２２）式に基づく、ＭＡＰ推定法によるパラメータの集合Θの推定値は、最尤推定と一致する。もしくは、混合重みの事前分布として、下記（２４）式のようなディリクレ分布を仮定する。
【００６４】
【数２３】
0006290803.tif 000024
【００６５】
そして、混合重み以外のパラメータに対しては、一様な事前分布を仮定してもよい。この場合、パラメータの集合Θの事前分布Ｐ(Θ)は、上記（２４）式に示す混合重みの事前分布に比例する。ここで、上記（２４）式におけるψは、ハイパーパラメータと呼ばれる所定の定数である。ψは、任意の正数に設定することができるが、例えばψ＝６００とすればよい。
【００６６】
以下では、上記（２２）式に基づくＭＡＰ推定法により、パラメータの集合Θを推定するためのアルゴリズムの一例として、集合Ｄ:＝｛ｄ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.」に詳述されている。
【００６７】
ＥＭアルゴリズムとは、以下に定義するＥステップとＭステップを、収束条件が満たされるまで反復するものである。Ｅステップでは、下記（２５）式で定義されるＱ関数：Ｑ（Θ;Θ´）を計算する。
【００６８】
【数２４】
0006290803.tif 000025
【００６９】
ここで、lnＰ(Ｙ,Ｄ|Θ)は、完全データの集合｛Ｙ,Ｄ｝の対数尤度を表し、Ｐ(Ｄ|Ｙ,Θ´)は、パラメータの集合Θの現在の推定値Θ´に対する集合Ｄの事後確率を表し、<・>P(Ｄ|Ｙ,Θ´)は、Ｐ(Ｄ|Ｙ,Θ´)に関する期待値演算を表す。
【００７０】
一方、Ｍステップでは、Ｑ関数の最大化により、パラメータの集合Θを更新する。ＥＭアルゴリズムの各反復において、評価関数ｐ(Ｙ|Θ)ｐ(Θ)の単調非減少性が保証されている。Ｑ関数の具体形を計算するために、下記（２６）式〜（３０）式のように、ｐ(Ｙ,Ｄ|Θ)、ｐ(Ｄ|Ｙ,Θ´)を求める。
【００７１】
【数２５】
0006290803.tif 000026
【００７２】
ここで、上記（３０）式におけるγ(n)tfは、下記（３１）式〜（３３）式で定義する。
【００７３】
【数２６】
0006290803.tif 000027
【００７４】
ただし、簡単のため、上記（３３）式において、α(n)f、回帰行列Ｇkf、φ(n)tf、行列Ｂ(n)fの現在の推定値を、それぞれ、単にα(n)f、回帰行列Ｇkf、φ(n)tf、行列Ｂ(n)fと表記した。
【００７５】
上記（２８）式、上記（３０）式を、上記（２５）式へ代入することで、Ｑ関数の具体形が、下記（３４）式、（３５）式のように得られる。
【００７６】
【数２７】
0006290803.tif 000028
【００７７】
混合重みα(n)fの更新式は、拘束条件ΣNn=1α(n)f＝１に注意して、ラグランジュの未定乗数法を用いることで得られる。φ(n)tf及び行列Ｂ(n)fの更新式は、上記（３５）式に示すＱ関数のφ(n)tf、及び、行列Ｂ(n)fの複素共役である(行列Ｂ(n)f)＊に関する偏微分を０とおくことで得られる。
【００７８】
回帰行列Ｇkfの更新式は、上記（３５）式に示すＱ関数から、回帰行列Ｇkfのエルミート転置である行列ＧHkfに依存する項のみを抜き出すと、下記（３６）式、（３７）式のようになる。
【００７９】
【数２８】
0006290803.tif 000029
【００８０】
上記（３７）式の、行列〜ＧHfに関する偏微分を０とおいて整理すると、下記（４１）式のようになる。
【００８１】
【数２９】
0006290803.tif 000030
【００８２】
上記（４１）式の両辺に、vec作用素を作用させ、下記（ａ）式で示される、行列Ａ、行列Ｂ、行列Ｘについての、クロネッカー積に関する公式を適用すると、下記（４２）式のようになる。
【００８３】
【数３０】
0006290803.tif 000031
【００８４】
【数３１】
0006290803.tif 000032
【００８５】
ただし、上記（４２）式において、vec[a1・・・aP]及び行列〜ＧHfを、それぞれ下記（４３）式、（４４）式のように定義する。
【００８６】
【数３２】
0006290803.tif 000033
【００８７】
よって、上記（４２）式から、下記（４５）式のように、vec[行列〜ＧHf]が求まる。
【００８８】
【数３３】
0006290803.tif 000034
【００８９】
＜実施形態１の一態様＞
以下、上述の実施形態１の理論的背景に基づく、実施形態１の一態様を説明する。なお、実施形態１の一態様において、音源数Ｎは既知と仮定する。
【００９０】
（実施形態１に係るモデル推定装置の構成）
図１は、実施形態１に係るモデル推定装置の構成の一例を示す図である。実施形態１に係るモデル推定装置１０Ａは、残響除去処理部１１Ａ、クラスタリング部１２Ａを有する。残響除去処理部１１Ａは、初期化部１１Ａ−１、共分散行列更新部１１Ａ−２、回帰行列更新部１１Ａ−３、残響除去部１１Ａ−４を有する。共分散行列更新部１１Ａ−２及び回帰行列更新部１１Ａ−３及び混合重み更新部１２Ａ−２は、パラメータ推定部の一例である。残響除去部１１Ａ−４は、信号推定部の一例である。事後確率更新部１２Ａ−１は、事後確率計算部の一例である。
【００９１】
初期化部１１Ａ−１は、まず、パラメータの集合Θの初期値を計算する。この初期値は、例えば、以下のように計算することができる。まず、アクティブな音源の番号ｄtfの推定値＾ｄtfを、残響モデルを含まない従来のクラスタリングに基づく音源分離技術を用いて計算する。残響モデルを含まない従来のクラスタリングに基づく音源分離技術は、文献４「伊藤信貴，荒木章子，木下慶介，中谷智広,“音源位置情報に基づく劣決定ブラインド音源分離のためのパーミュテーションフリークラスタリング法”，電子情報通信学会論文誌,  vol.  J97-A,  no.  4,  pp.  234.246,  Apr.  2014.」に詳述されている。
【００９２】
次に、初期化部１１Ａ−１は、推定値＾ｄtfを用いて、下記（４６）式〜（４９）式により、各パラメータを初期化する。なお、下記（４６）式及び（４８）式における集合Ｃ(n)fは、Ｃ(n)f:＝｛ｔ|ｄtf＝ｎ｝で定義される行列である。また、下記（４６）式及び（４８）式における＃Ｃ(n)fは、集合Ｃ(n)fの要素数を表す。また、下記（４９）式における“tr[・]”は、行列[・]のトレースを表す。
【００９３】
【数３４】
0006290803.tif 000035
【００９４】
共分散行列更新部１１Ａ−２は、各音源ｎ（ｎ＝１,・・・,Ｎ）の残響を含まないマイクロホン像のベクトルｓ(n)tfの共分散行列φ(n)tfＢ(n)fのパラメータφ(n)tf及び行列Ｂ(n)fを、それぞれ下記（５０）式、（５１）式により更新する。
【００９５】
【数３５】
0006290803.tif 000036
【００９６】
回帰行列更新部１１Ａ−３は、回帰行列Ｇkfを、下記（５２）式、（５３）式により更新する。
【００９７】
【数３６】
0006290803.tif 000037
【００９８】
ここで、上記（５３）式の左辺に現れる行列〜Ｇf及び上記（５３）式の右辺に現れるベクトル〜ｙt-Δ-1,fは、下記（５４）式、（５５）式のように定義される。
【００９９】
【数３７】
0006290803.tif 000038
【０１００】
残響除去部１１Ａ−４は、残響を含まない混合信号ベクトルの推定値＾ｘtfを、下記（５６）式により更新する。
【０１０１】
【数３８】
0006290803.tif 000039
【０１０２】
クラスタリング部１２Ａは、事後確率更新部１２Ａ−１、混合重み更新部１２Ａ−２を有する。事後確率更新部１２Ａ−１は、時間周波数点(ｔ,ｆ)でｎ（ｎ＝１,・・・,Ｎ）番目の音源信号がアクティブである事後確率γ(n)tfを、下記（５７）式により更新する。なお、γ(n)tf:＝Ｐ(ｄtf＝ｎ|ベクトルｙtf,Θ)と定義する。
【０１０３】
【数３９】
0006290803.tif 000040
【０１０４】
混合重み更新部１２Ａ−２は、混合重みα(n)fを、下記（５８）式により更新する。
【０１０５】
【数４０】
0006290803.tif 000041
【０１０６】
なお、性能向上のため、モデル推定装置１０Ａの全処理に先立ち、残響を含む混合信号ベクトルｙtfに対し、前処理として、下記に示す白色化をおこなってもよい。
【０１０７】
【数４１】
0006290803.tif 000042
【０１０８】
なお、実施形態１は、クラスタリング部１２Ａの事後確率更新部１２Ａ−１が、上記（５７）式に基づき、時間周波数点(ｔ,ｆ)でｎ（ｎ＝１,・・・,Ｎ）番目の音源信号がアクティブである事後確率γ(n)tfを計算するとした。しかし、これに限らず、k-meansクラスタリング等の従来技法を用い、時間周波数点(ｔ,ｆ)でｎ（ｎ＝１,・・・,Ｎ）番目の音源信号がアクティブである事後確率γ(n)tfを計算するとしてもよい。
【０１０９】
（実施形態１に係るモデル推定装置の処理）
図２は、実施形態１に係るモデル推定装置の処理手順の一例を示すフローチャートである。以下に述べるモデル推定装置１０Ａの処理は、所定の収束判定条件が満たされるまで反復される。所定の収束条件は、例えば、「所定の反復回数に達している、又は、事後確率更新部１２Ａ−１、混合重み更新部１２Ａ−２の各更新部のうち１つ以上の更新部による更新前後のパラメータ値の差分が所定の閾値未満である」などとすればよい。
【０１１０】
先ず、ステップＳ１１では、初期化部１１Ａ−１は、パラメータの集合Θの初期値を、上記（４６）式〜（４９）式に基づき計算し、モデル推定装置１０Ａの主記憶装置に保存する。次に、ステップＳ１２では、残響除去部１１Ａ−４は、モデル推定装置１０Ａの主記憶装置に現在保存されている回帰行列Ｇkfに基づき、上記（５６）式により、残響を含まない混合信号ベクトルの推定値＾ｘtfを更新する（“残響除去”処理）。
【０１１１】
次に、ステップＳ１３では、事後確率更新部１２Ａ−１は、時間周波数点(ｔ,ｆ)でｎ（ｎ＝１,・・・,Ｎ）番目の音源信号がアクティブである事後確率γ(n)tfを、上記（５７）式により計算し、モデル推定装置１０Ａの主記憶装置に保存する。また、ステップＳ１３では、混合重み更新部１２Ａ−２は、混合重みα(n)fを、上記（５８）式により計算し、モデル推定装置１０Ａの主記憶装置に保存する（以上、“クラスタリング”処理）。
【０１１２】
次に、モデル推定装置１０Ａは、収束判定条件が満たされているか否かを判定する（ステップＳ１４）。モデル推定装置１０Ａは、収束判定条件が満たされている場合（ステップＳ１４Ｙｅｓ）、処理を終了する。モデル推定装置１０Ａは、収束判定条件が満たされていない場合（ステップＳ１４Ｎｏ）、ステップＳ１５へ処理を移す。
【０１１３】
ステップＳ１５では、共分散行列更新部１１Ａ−２は、各音源ｎ（ｎ＝１,・・・,Ｎ）の残響を含まないマイクロホン像のベクトルｓ(n)tfの共分散行列φ(n)tfＢ(n)fのパラメータφ(n)tf及び行列Ｂ(n)fを、それぞれ上記（５０）式、（５１）式により計算し、モデル推定装置１０Ａの主記憶装置に更新保存する。また、ステップＳ１５では、回帰行列更新部１１Ａ−３は、共分散行列更新部１１Ａ−２により計算されたパラメータφ(n)tf及び行列Ｂ(n)fに基づき、回帰行列Ｇkfを、上記（５２）式、（５３）式により計算し、モデル推定装置１０Ａの主記憶装置に更新保存する。
【０１１４】
また、ステップＳ１５では、事後確率更新部１２Ａ−１は、モデル推定装置１０Ａの主記憶装置に現在保存されているパラメータの集合Θ、及び、最後に実行したステップＳ１２による残響を含まない混合信号ベクトルの推定値＾ｘtfに基づき、上記（５７）式により、事後確率γ(n)tfを計算し、モデル推定装置１０Ａの主記憶装置に更新保存する。また、ステップＳ１５では、混合重み更新部１２Ａ−２は、事後確率更新部１２Ａ−１により計算された事後確率γ(n)tfに基づき、上記（５８）式により、混合重みα(n)fを更新し、モデル推定装置１０Ａの主記憶装置に更新保存する。以上のステップＳ１５の処理が終了すると、モデル推定装置１０Ａは、ステップＳ１２へ処理を移す。
【０１１５】
［実施形態２］
以下、実施形態２について、実施形態２の理論的背景を説明後、実施形態２の一態様を説明する。
【０１１６】
＜実施形態２の理論的背景＞
実施形態１のように、上記（１９）式に示す周波数依存の混合重みを用いる場合、評価関数である事後確率には、パーミュテーション（置換）の不定性がある。すなわち、｛１,・・・,Ｎ｝上の置換Πfにより、パラメータの集合Θのα(n)f、φ(n)tf、行列Ｂ(n)tfの順序を、下記（６２）式のように入れ替えた場合を考える。
【０１１７】
【数４２】
0006290803.tif 000043
【０１１８】
このとき、下記（６３）式が成り立つ。
【０１１９】
【数４３】
0006290803.tif 000044
【０１２０】
すなわち、事後確率を最大化するだけでは、推定されたΘにおける番号ｎは、周波数毎に、異なる音源に対応してしまうというパーミュテーション問題がある。よって、推定されたΘをそのまま用いては、適切に目的音強調を行うことはできない。従って、実施形態１に基づいて目的音強調装置を構成する際には、番号ｎが周波数によらず同一の音源に対応するように置換Πfを決定する、パーミュテーション解決の処理が別途必要となる。
【０１２１】
これに対し、実施形態２のモデル推定装置は、時間依存の混合重みを用いる。これにより、上記文献４に開示されている通り、事後確率の最大化により、パーミュテーション問題を生じずにモデル推定が可能である。
【０１２２】
以下、実施形態２の理論的背景を、実施形態１との差異に重点を置きながら説明する。
【０１２３】
（実施形態２の残響を含む混合信号ベクトルのモデル化）
実施形態２では、アクティブな音源の番号ｄtfの確率モデルＰ(ｄtf|Θ)を、周波数依存の混合重みではなく、時間依存の混合重みα(n)tを用いて、下記（６４）式でモデル化する。
【０１２４】
【数４４】
0006290803.tif 000045
【０１２５】
従って、実施形態２における残響を含む混合信号ベクトルｙtfの分布を表す確率モデル（上記（１６）式参照）の具体形は、下記（６５）式のように得られる。  
【０１２６】
【数４５】
0006290803.tif 000046
【０１２７】
パラメータの集合Θは、具体的には、下記（６６）式で表される。
【０１２８】
【数４６】
0006290803.tif 000047
【０１２９】
（実施形態２のパラメータ推定アルゴリズムの導出）
ＥＭアルゴリズムにより、事後確率を最大化する点は、実施形態２は、実施形態１と同様である。しかし、実施形態２は、ＥＭアルゴリズムの各反復において、Ｅステップ、Ｍステップの処理に加えて、Ｐ（Permutation）ステップの処理を行う。Ｐステップでは、各周波数binの番号ｆにて、目的関数である事後確率が最大となるように、共分散行列φ(n)tfＢ(n)fを音源間で置換することにより、パーミュテーションを解決する。すなわち、Πfを｛１,・・・,Ｎ｝上の置換として、下記（６７）式〜（６９）式の処理を行う。
【０１３０】
【数４７】
0006290803.tif 000048
【０１３１】
なお、Ｅステップ及びＭステップにおける更新式の導出は、実施形態１と同様であるので、説明を省略する。
【０１３２】
＜実施形態２の一態様＞
以下、上述の実施形態２の理論的背景に基づく、実施形態２の一態様を説明する。なお、実施形態２の一態様において、音源数Ｎは既知と仮定する。しかし、実施形態２は、真の音源数Ｎ0が既知でなくても、その上限は分かっていると仮定し、仮定する音源数Ｎを、真の音源数Ｎ0の上限より大きく設定することで、音源数が既知である場合と同様に実施可能である。
【０１３３】
（実施形態２に係るモデル推定装置の構成）
図３は、実施形態２に係るモデル推定装置の構成の一例を示す図である。実施形態２に係るモデル推定装置１０Ｂは、残響除去処理部１１Ｂ、クラスタリング部１２Ｂを有する。残響除去処理部１１Ｂは、初期化部１１Ｂ−１、共分散行列更新部１１Ｂ−２、回帰行列更新部１１Ｂ−３、残響除去部１１Ｂ−４を有する。共分散行列更新部１１Ｂ−２及び回帰行列更新部１１Ｂ−３及び混合重み更新部１２Ｂ−２は、パラメータ推定部の一例である。残響除去部１１Ｂ−４は、信号推定部の一例である。事後確率更新部１２Ｂ−１は、事後確率計算部の一例である。
【０１３４】
初期化部１１Ｂ−１は、まず、パラメータの集合Θの初期値を計算する。この初期値は、例えば、以下のように計算することができる。まず、アクティブな音源の番号ｄtfの推定値＾ｄtfを、実施形態１と同様に、残響モデルを含まない従来のクラスタリングに基づく音源分離技術を用いて計算する。次に、初期化部１１Ｂ−１は、推定値＾ｄtfを用いて、上記（４７）式〜（４９）式、及び、下記（７０）式により、各パラメータを初期化する。なお、下記（７０）式における集合〜Ｃ(n)tは、Ｃ(n)t:＝｛ｆ|ｄtf＝ｎ｝で定義される行列である。また、下記（７０）式における＃Ｃ(n)tは、集合Ｃ(n)tの要素数を表す。
【０１３５】
【数４８】
0006290803.tif 000049
【０１３６】
共分散行列更新部１１Ｂ−２、回帰行列更新部１１Ｂ−３、残響除去部１１Ｂ−４は、実施形態１の共分散行列更新部１１Ａ−２、回帰行列更新部１１Ａ−３、残響除去部１１Ａ−４とそれぞれ同様である。
【０１３７】
クラスタリング部１２Ｂは、事後確率更新部１２Ｂ−１、混合重み更新部１２Ｂ−２、パーミュテーション解決部１２Ｂ−３を有する。事後確率更新部１２Ｂ−１は、時間周波数点(ｔ,ｆ)でｎ（ｎ＝１,・・・,Ｎ）番目の音源信号がアクティブである事後確率γ(n)tfを、下記（７１）式により更新する。なお、γ(n)tf:＝Ｐ(ｄtf＝ｎ|ベクトルｙtf,Θ)と定義する。
【０１３８】
【数４９】
0006290803.tif 000050
【０１３９】
混合重み更新部１２Ｂ−２は、混合重みα(n)tを、下記（７２）式により更新する。
【０１４０】
【数５０】
0006290803.tif 000051
【０１４１】
パーミュテーション解決部１２Ｂ−３は、各周波数binの番号ｆにて、目的関数である事後確率が最大となるように、共分散行列φ(n)tfＢ(n)fを音源間で置換することにより、パーミュテーションを解決する。すなわち、Πfを｛１,・・・,Ｎ｝上の置換として、下記（７３）式〜（７５）式により、共分散行列φ(n)tfＢ(n)fを置換する。
【０１４２】
【数５１】
0006290803.tif 000052
【０１４３】
なお、性能向上のため、モデル推定装置１０Ｂの全処理に先立ち、残響を含む混合信号ベクトルｙtfに対し、前処理として、上記（５９）式〜（６１）式に示す白色化をおこなってもよい。
【０１４４】
なお、実施形態２は、クラスタリング部１２Ｂの事後確率更新部１２Ｂ−１が、上記（７１）式に基づき、時間周波数点(ｔ,ｆ)でｎ（ｎ＝１,・・・,Ｎ）番目の音源信号がアクティブである事後確率γ(n)tfを計算するとした。しかし、これに限らず、k-meansクラスタリング等の従来技法を用い、時間周波数点(ｔ,ｆ)でｎ（ｎ＝１,・・・,Ｎ）番目の音源信号がアクティブである事後確率γ(n)tfを計算するとしてもよい。
【０１４５】
［実施形態３］
実施形態３は、実施形態２のモデル推定装置１０Ｂを用いて、上記文献４に記載の音源数推定技術により、音源数も推定する構成にしたものである。実施形態３は、真の音源数Ｎ0は分からないがその上限は分かっていると仮定し、仮定する音源数Ｎを、真の音源数Ｎ0の上限より大きく設定する。
【０１４６】
（実施形態３に係るモデル推定装置の構成）
図４は、実施形態３に係るモデル推定装置の構成の一例を示す図である。実施形態３に係るモデル推定装置１０Ｃは、実施形態２に係るモデル推定装置１０Ｂと比較して、音源数推定部１３をさらに有する。
【０１４７】
音源数推定部１３は、クラスタリング部１２Ｂによって計算されたｎ番目の音源がアクティブである事後確率γ(n)tfを用いて、番号ｎ＝１,・・・,Ｎのうち、真の音源に対応する番号ｎ(１),・・・,  ｎ(Ｎ0)を判定し、真の音源に対応する番号のパラメータのみを出力する。具体的には、音源数推定部１３は、ｎ番目の音源がアクティブである事後確率γ(n)tfを用いて、ｎ番目の音源がアクティブである事後確率の総和を、例えば下記（７６）式により算出する。
【０１４８】
【数５２】
0006290803.tif 000053
【０１４９】
そして、音源数推定部１３は、各ｎ番目の音源がアクティブである事後確率の総和ρ(n)を２つにクラスタリングし、総和の大きい方のクラスタに属するρ(n)の番号ｎ＝ｎ(１),・・・,  ｎ(＾Ｎ0)を求め、真の音源に対応する番号とみなす。例えば、音源数推定部１３は、ρ(n)に対して、クラスタ数２のk-meansクラスタリングを適用してクラスタリングする。
【０１５０】
最後に、音源数推定部１３は、真の音源に対応するｎ＝ｎ(１),・・・,  ｎ(＾Ｎ0)に対応する、下記（７７）式に示すパラメータのみを出力する。なお、下記（７７）式において、ｌ＝１,・・・,＾Ｎ0である。
【０１５１】
【数５３】
0006290803.tif 000054
【０１５２】
（実施形態３に係るモデル推定装置の処理）
図５は、実施形態３に係るモデル推定装置の処理手順の一例を示すフローチャートである。以下に述べるモデル推定装置１０Ｃの処理は、実施形態１又は２と同様の所定の収束判定条件が満たされるまで反復される。
【０１５３】
先ず、ステップＳ２１では、初期化部１１Ｂ−１は、パラメータの集合Θの初期値を、上記（４７）式〜（４９）式、及び、（７０）式に基づき計算し、モデル推定装置１０Ｃの主記憶装置に保存する。次に、ステップＳ２２では、残響除去部１１Ｂ−４は、モデル推定装置１０Ｃの主記憶装置に現在保存されている回帰行列Ｇkfに基づき、上記（５６）式により、残響を含まない混合信号ベクトルの推定値＾ｘtfを更新する（“残響除去”処理）。
【０１５４】
次に、ステップＳ２３では、事後確率更新部１２Ｂ−１は、時間周波数点(ｔ,ｆ)でｎ（ｎ＝１,・・・,Ｎ）番目の音源信号がアクティブである事後確率γ(n)tfを、上記（７１）式により計算し、モデル推定装置１０Ｃの主記憶装置に保存する。また、ステップＳ２３では、混合重み更新部１２Ｂ−２は、混合重みα(n)tを、上記（７２）式により計算し、モデル推定装置１０Ｃの主記憶装置に保存する（以上、“クラスタリング”処理）。また、ステップＳ２３では、パーミュテーション解決部１２Ｂ−３は、Πfを｛１,・・・,Ｎ｝上の置換として、上記（７３）式〜（７５）式により、共分散行列φ(n)tfＢ(n)fを置換する。
【０１５５】
次に、モデル推定装置１０Ｃは、収束判定条件が満たされているか否かを判定する（ステップＳ２４）。モデル推定装置１０Ｃは、収束判定条件が満たされている場合（ステップＳ２４Ｙｅｓ）、ステップＳ２６へ処理を移す。モデル推定装置１０Ｃは、収束判定条件が満たされていない場合（ステップＳ２４Ｎｏ）、ステップＳ２５へ処理を移す。
【０１５６】
ステップＳ２５の処理は、図２に示す実施形態１のステップＳ１５の処理と同様である。ステップＳ２６では、音源数推定部１３は、ｎ番目の音源がアクティブである事後確率γ(n)tfを用いて、真の音源数を推定し、推定結果を出力する。
【０１５７】
［実施形態４］
実施形態４に係る目的音強調装置は、実施形態１〜３に係るモデル推定装置１０Ａ〜１０Ｃのいずれかを有する目的音強調装置１００である。
【０１５８】
（実施形態４に係る目的音強調装置の構成）
図６は、実施形態４に係る目的音強調装置の構成の一例を示す図である。実施形態４に係る目的音強調装置１００は、周波数領域変換部２０、モデル推定装置１０Ａ（あるいは１０Ｂ又は１０Ｃ）、強調音計算部３０、時間領域変換部４０を有する。
【０１５９】
周波数領域変換部２０は、時間領域での残響を含む混合信号ベクトル〜ｙτを、短時間フーリエ変換などの時間周波数変換により、時間周波数領域での残響を含む混合信号ベクトルｙtfに変換する。ここで、混合信号ベクトル〜ｙτは、下記（７８）式により定義される。
【０１６０】
【数５４】
0006290803.tif 000055
【０１６１】
ただし、上記（７８）式において、〜ｙ(m)τは、時間領域でのｍ（ｍ＝１,・・・,Ｍ）番目のマイクロホンで観測された残響を含む混合信号であり、τはサンプル番号を表す。モデル推定装置１０Ａ（あるいは１０Ｂ又は１０Ｃ）は、パラメータの集合Θと、各音源ｎがアクティブである事後確率γ(n)tfを計算する。
【０１６２】
強調音計算部３０は、周波数領域変換部２０から出力された時間周波数領域での残響を含む混合信号ベクトルｙtfと、モデル推定装置１０Ａ（あるいは１０Ｂ又は１０Ｃ）から出力されたパラメータの集合Θと各音源ｎがアクティブである事後確率γ(n)tfとを用いて、時間周波数領域での各音源の残響を含まないマイクロホン像の推定値＾ｓ(n)tfを、下記（７９）式及び（８０）式により計算し、出力する。
【０１６３】
【数５５】
0006290803.tif 000056
【０１６４】
なお、目的音強調装置１００において実施形態１のモデル推定装置１０Ａを用いる場合は、上記（７９）式及び（８０）式の処理に先立って、γ(n)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.」に記載の方法により行うことができる。
【０１６５】
時間領域変換部４０は、強調音計算部３０から出力された時間周波数領域での各音源の残響を含まないマイクロホン像の推定値のベクトル＾ｓ(n)tfに、逆短時間フーリエ変換などの時間周波数変換の逆変換を適用して、時間領域での各音源の残響を含まないマイクロホン像の推定値のベクトル〜＾ｓ(n)τを計算する。ここで、ベクトル〜＾ｓ(n)τは、下記（８１）式により定義される。ただし、〜＾ｓ(m,n)τは、ベクトル＾ｓ(n)τの第ｍ要素＾ｓ(m,n)tfの逆短時間フーリエ変換である。
【０１６６】
【数５６】
0006290803.tif 000057
【０１６７】
なお、強調音計算部３０において、残響除去と音源分離を同時に実現する例を示したが、残響のみを除去するために、時間周波数領域での残響を含まない混合信号の推定値＾ｘtfに、逆短時間フーリエ変換などの時間周波数変換の逆変換を適用して、時間領域での残響を含まない混合信号の推定値のベクトル〜＾ｘτを得る構成としてもよい。ここで、ベクトル〜＾ｘτは、下記（８２）式で定義される。ただし、〜＾ｓ(m,n)τは、ベクトル＾ｓ(n)τの第ｍ要素＾ｓ(m,n)tfの逆短時間フーリエ変換である。
【０１６８】
【数５７】
0006290803.tif 000058
【０１６９】
（実施形態４に係る目的音強調装置の処理）
図７は、実施形態４に係る目的音強調装置の処理手順の一例を示すフローチャートである。実施形態４に係る目的音強調装置１００において、先ず、ステップＳ３１では、周波数領域変換部２０は、各マイクロホンで観測された信号をそれぞれ時間周波数領域の信号に変換する。次に、ステップＳ３２では、モデル推定装置１０Ａ（あるいは１０Ｂ又は１０Ｃ）は、モデル推定を行う。次に、ステップＳ３３では、強調音計算部３０は、強調音を計算により推定する。次に、ステップＳ３４では、時間領域変換部４０は、強調音計算部３０により推定された強調音を周波数領域から時間領域に変換する。
【実施例】
【０１７０】
以下、実施形態４を例に取り、開示の実施形態の実施例及びその効果について説明する。図８及び図９は、実施形態４の効果の一例を説明する図である。実施形態４に係る目的音強調装置１００（以下「提案法」）と、従来の残響モデルを含まないクラスタリングベースの音源分離手法（例えば、文献４に記載の手法、以下「従来法」）の性能を比較する実験をおこなった。ただし、実施形態４に係る目的音強調装置１００のモデル推定装置としては、実施形態２に係るモデル推定装置１０Ｂを用いた。
【０１７１】
マイクロホンで観測される残響を含む混合信号は、残響を含まない音声波形に、実験室で計測したインパルス応答（例えば、上述の文献５参照）を畳み込むことにより生成した。図８は、インパルス応答を計測した際のマイクロホンと音源の位置を示す。なお、提案法及び従来法の両方において、パラメータΘの推定に先立って、残響を含む混合信号ベクトルｙtfに対し、上記（５９）式〜（６１）式に示す白色化をおこなった。また、音源数Ｎは既知とした。また、他の実験条件は、下記（表１）に示すとおりとした。なお、図８に示す実験室は、４．４５ｍ×３．５５ｍ×（高さ）２．５０ｍの空間であった。また、図８に示すSource1及び2とMicrophone1及び2の、実験室の床面に対する高さは、１．２ｍとした。
【０１７２】
【表１】
0006290803.tif 000059
【０１７３】
提案法及び従来法の性能は、下記（８３）式で定義されるＳＩＲ（Signal-to-Interference  Ratio）により評価した。
【０１７４】
【数５８】
0006290803.tif 000060
【０１７５】
ここで、〜＾ｓ(1,n,ν)τは、〜＾ｓ(1,n)に含まれるν番目の音源成分を表す。Τ:＝８ｋHz×８ｓ＝６４０００は、サンプリング点の総数を表し、Σν≠ｎは、ｎ以外のνの値に対する総和を表す。
【０１７６】
ここで、〜＾ｓ(1,n,ν)τの求め方を説明する。観測された残響を含む混合信号ベクトルｙtfは、ν番目の音源の残響を含むマイクロホン像のベクトルｘ(ν)tfを用いて、下記（８４）式のように分解できる。
【０１７７】
【数５９】
0006290803.tif 000061
【０１７８】
従って、ｎ番目の音源の残響を含まないマイクロホン像の推定値のベクトル＾ｓ(n)tfは、下記（８５）式及び（８６）式のように分解できる。
【０１７９】
【数６０】
0006290803.tif 000062
【０１８０】
ここで、上記（８６）式において、＾ｓ(n,ν)tfは、＾ｓ(n)tfに含まれるν番目の音源成分を表す。よって、下記（８７）式により、＾ｓ(n,ν)tfを求め、＾ｓ(n,ν)tfを逆短時間フーリエ変換して〜＾ｓ(n,ν)τを求め、〜＾ｓ(n,ν)τの第１要素として〜＾ｓ(1,n,ν)τが求まる。
【０１８１】
【数６１】
0006290803.tif 000063
【０１８２】
図９に、各残響時間に対し、音声波形の組み合わせを変えて８回の試行を行った際のＳＩＲの平均値をプロットしたグラフを示す。残響時間が最も小さい条件（残響時間１３０ｍｓ程度）では、提案法と従来法は同等の性能を示した。しかし、図９に示すように、残響時間が大きくなるにつれて、従来法に対する提案法の性能改善量が増加する傾向があった。特に、残響時間が３７０ｍｓ程度の場合に、性能改善量は、試行中、最大の約４ｄＢとなった。
【０１８３】
以上から、実施形態１〜４は、独立成分分析に基づく音源分離技術と比較して、音源数が未知の場合でも適用できる等の利点があるクラスタリングに基づく音源分離技術において、線形予測に基づく残響除去とクラスタリングに基づく音源分離を交互に反復する。実施形態１〜４は、線形予測に基づく残響除去により推定された残響を含まない混合信号に対して、クラスタリングに基づく音源分離を適用することで、音源分離の性能を向上させることができる。さらに、実施形態１〜４は、改善された音源分離結果を用いることで、残響除去の性能を改善することができる。よって、実施形態１〜４は、残響除去と上記音源分離の反復により、残響時間がフレーム長に比べて長い場合でも、より高精度な音源分離を実現することができる。
【０１８４】
（モデル推定装置及び目的音強調装置の装置構成について）
図１、図３、図４に示すモデル推定装置１０Ａ〜１０Ｃ及び図６に示す目的音強調装置１００の各構成要素は機能概念的なものであり、必ずしも物理的に図示のように構成されていることを要さない。すなわち、モデル推定装置１０Ａ〜１０Ｃ及び目的音強調装置１００の機能の分散及び統合の具体的形態は図示のものに限られず、全部又は一部を、各種の負荷や使用状況等に応じて、任意の単位で機能的又は物理的に分散又は統合して構成することができる。
【０１８５】
また、モデル推定装置１０Ａ〜１０Ｃ及び目的音強調装置１００において行われる各処理は、全部又は任意の一部が、ＣＰＵ（Central  Processing  Unit）等の処理装置及び処理装置により解析実行されるプログラムにて実現されてもよい。また、モデル推定装置１０Ａ〜１０Ｃ及び目的音強調装置１００において行われる各処理は、ワイヤードロジックによるハードウェアとして実現されてもよい。
【０１８６】
また、実施形態において説明した各処理のうち、自動的に行われるものとして説明した処理の全部又は一部を手動的に行うこともできる。もしくは、実施形態において説明した各処理のうち、手動的に行われるものとして説明した処理の全部又は一部を公知の方法で自動的に行うこともできる。この他、上述及び図示の処理手順、制御手順、具体的名称、各種のデータやパラメータを含む情報については、特記する場合を除いて適宜変更することができる。
【０１８７】
（プログラムについて）
図１０は、プログラムが実行されることにより、モデル推定装置及び目的音強調装置が実現されるコンピュータの一例を示す図である。コンピュータ１０００は、例えば、メモリ１０１０、ＣＰＵ１０２０を有する。また、コンピュータ１０００は、ハードディスクドライブインタフェース１０３０、ディスクドライブインタフェース１０４０、シリアルポートインタフェース１０５０、ビデオアダプタ１０６０、ネットワークインタフェース１０７０を有する。コンピュータ１０００において、これらの各部はバス１０８０によって接続される。
【０１８８】
メモリ１０１０は、ＲＯＭ（Read  Only  Memory）１０１１及びＲＡＭ（Random  Access  Memory）１０１２を含む。ＲＯＭ１０１１は、例えば、ＢＩＯＳ（Basic  Input  Output  System）等のブートプログラムを記憶する。ハードディスクドライブインタフェース１０３０は、ハードディスクドライブ１０３１に接続される。ディスクドライブインタフェース１０４０は、ディスクドライブ１０４１に接続される。例えば磁気ディスクや光ディスク等の着脱可能な記憶媒体が、ディスクドライブ１０４１に挿入される。シリアルポートインタフェース１０５０は、例えばマウス１０５１、キーボード１０５２に接続される。ビデオアダプタ１０６０は、例えばディスプレイ１０６１に接続される。
【０１８９】
ハードディスクドライブ１０３１は、例えば、ＯＳ１０９１、アプリケーションプログラム１０９２、プログラムモジュール１０９３、プログラムデータ１０９４を記憶する。すなわち、モデル推定装置１０Ａ〜１０Ｃ及び目的音強調装置１００の各処理を規定するプログラムは、コンピュータ１０００によって実行される指令が記述されたプログラムモジュール１０９３として、例えばハードディスクドライブ１０３１に記憶される。例えば、モデル推定装置１０Ａ〜１０Ｃ及び目的音強調装置１００における機能構成と同様の情報処理を実行するためのプログラムモジュール１０９３が、ハードディスクドライブ１０３１に記憶される。
【０１９０】
また、実施形態の処理で用いられる設定データは、プログラムデータ１０９４として、例えばメモリ１０１０やハードディスクドライブ１０３１に記憶される。そして、ＣＰＵ１０２０が、メモリ１０１０やハードディスクドライブ１０３１に記憶されたプログラムモジュール１０９３やプログラムデータ１０９４を必要に応じてＲＡＭ１０１２に読み出して実行する。
【０１９１】
なお、プログラムモジュール１０９３やプログラムデータ１０９４は、ハードディスクドライブ１０３１に記憶される場合に限らず、例えば着脱可能な記憶媒体に記憶され、ディスクドライブ１０４１等を介してＣＰＵ１０２０によって読み出されてもよい。あるいは、プログラムモジュール１０９３やプログラムデータ１０９４は、ネットワーク（ＬＡＮ（Local  Area  Network）、ＷＡＮ（Wide  Area  Network）等）を介して接続された他のコンピュータに記憶されてもよい。そして、プログラムモジュール１０９３やプログラムデータ１０９４は、ネットワークインタフェース１０７０を介してＣＰＵ１０２０によって読み出されてもよい。
【０１９２】
上記実施形態及びその他の実施形態は、本願が開示する技術に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。
【符号の説明】
【０１９３】
１０Ａ、１０Ｂ、１０Ｃ  モデル推定装置
１１Ａ、１１Ｂ          残響除去処理部
１１Ａ−１、１１Ｂ−１  初期化部
１１Ａ−２、１１Ｂ−２  共分散行列更新部
１１Ａ−３、１１Ｂ−３  回帰行列更新部
１１Ａ−４、１１Ｂ−４  残響除去部
１２Ａ、１２Ｂ          クラスタリング部
１２Ａ−１、１２Ｂ−１  事後確率更新部
１２Ａ−２、１２Ｂ−２  混合重み更新部
１２Ｂ−３              パーミュテーション解決部
１３                    音源数推定部
２０                    周波数領域変換部
３０                    強調音計算部
４０                    時間領域変換部
１００                  目的音強調装置
１０００  コンピュータ
１０１０  メモリ
１０２０  ＣＰＵ
【図１】
0006290803.tif 000064
【図２】
0006290803.tif 000065
【図３】
0006290803.tif 000066
【図４】
0006290803.tif 000067
【図５】
0006290803.tif 000068
【図６】
0006290803.tif 000069
【図７】
0006290803.tif 000070
【図８】
0006290803.tif 000071
【図９】
0006290803.tif 000072
【図１０】
0006290803.tif 000073