(57)【要約】  目的音源と背景雑音が混在する状況において、異なる位置で収録された観測信号に基づいて観測特徴量ベクトルを計算し、目的音源に対応するマスクと、背景雑音に対応するマスクを推定し、観測信号および目的音源に対応するマスクに基づいて背景雑音を含んだ目的音源の空間相関行列を計算し、観測信号および背景雑音に対応するマスクに基づいて背景雑音の空間相関行列を計算し、さらにそれぞれの空間相関行列を所定の係数で重み付けした行列に基づいて目的音源の空間相関行列を推定する。
2017094862.tif 000002
【特許請求の範囲】
【請求項１】
Ｎ個の目的音源（ただし、Ｎは１以上の整数）に対応するＮ個の第１の音響信号と、背景雑音に対応する第２の音響信号とが混在する状況において、それぞれ異なる位置で収録されたＭ個（ただし、Ｍは２以上の整数）の観測信号に基づいて計算される観測特徴量ベクトルに基づいて、時間周波数点ごとの、前記第１の音響信号が前記観測信号の特徴量に含まれる割合である第１のマスクと、時間周波数点ごとの、前記第２の音響信号が前記観測信号の特徴量に含まれる割合である第２のマスクと、を推定し、前記第１のマスクおよび前記第２のマスクに基づいて前記目的音源の空間相関行列を推定する空間相関行列推定装置であって、
前記観測信号および前記第１のマスクに基づいて計算された第１の特徴量行列を第１の係数で重み付けした第１の空間相関行列と、前記観測信号および前記第２のマスクに基づいて計算された第２の特徴量行列を第２の係数で重み付けした第２の空間相関行列と、に基づいて前記目的音源の空間相関行列を推定する雑音除去部を有することを特徴とする空間相関行列推定装置。
【請求項２】
前記雑音除去部は、背景雑音の空間相関行列が時間的に変化しないとの条件下において、前記目的音源の空間相関行列の推定値に含まれる背景雑音由来の成分が０となるように、前記第１の係数及び前記第２の係数を計算することを特徴とする請求項１に記載の空間相関行列推定装置。
【請求項３】
前記雑音除去部は、前記第１の係数と前記第２の係数との比率が、前記第１のマスクの時間平均値の逆数と前記第２のマスクの時間平均値の逆数との比率に等しくなるように、前記第１の係数及び前記第２の係数を計算することを特徴とする請求項１または２に記載の空間相関行列推定装置。
【請求項４】
前記第１の空間相関行列は、Ｎ＝１である場合、前記観測特徴量ベクトルに基づいて計算された観測特徴量行列の周波数ごとの時間平均であることを特徴とする請求項１から３のいずれか１項に記載の空間相関行列推定装置。
【請求項５】
前記観測信号に短時間信号分析を適用し、時間周波数点ごとの信号特徴量を抽出し、前記信号特徴量を成分とするＭ次元縦ベクトルである観測特徴量ベクトルを時間周波数点ごとに計算する時間周波数分析部と、
前記観測特徴量ベクトルに基づいて、時間周波数点ごとに、前記観測特徴量ベクトルと該観測特徴量ベクトルのエルミート転置とを乗じて観測特徴量行列を計算する観測特徴量行列計算部と、
前記目的音源のそれぞれについて、時間周波数点ごとに前記観測特徴量行列と前記第１のマスクを乗じて得られる行列の周波数ごとの時間平均を第１の特徴量行列として計算し、前記第１の特徴量行列に前記第１の係数を乗じることで前記第１の空間相関行列を推定する雑音下目的音空間相関行列推定部と、
前記背景雑音について、時間周波数点ごとに前記観測特徴量行列と前記第２のマスクを乗じて得られる行列の周波数ごとの時間平均を第２の特徴量行列として計算し、前記第２の特徴量行列に前記第２の係数を乗じることで前記第２の空間相関行列を推定する雑音空間相関行列推定部と、
をさらに有し、
前記雑音除去部は、前記第１の空間相関行列から前記第２の空間相関行列を減じることで前記目的音源の空間相関行列を推定し、
前記第１の係数と前記第２の係数との比率は、前記第１のマスクの時間平均値の逆数と前記第２のマスクの時間平均値の逆数との比率に等しいことを特徴とする請求項１に記載の空間相関行列推定装置。
【請求項６】
前記観測特徴量ベクトルの確率分布を、共分散行列が時刻ごとに異なる値を取るスカラーパラメータと時不変のパラメータとを要素にもつ正定値エルミート行列の積で表される、平均０のＭ次元複素ガウス分布であるＮ＋１個の要素分布からなる混合分布で周波数ごとにモデル化し、前記混合分布が前記観測特徴量ベクトルの分布に近くなるように前記混合分布のパラメータを推定することで得られる前記要素分布のそれぞれの事後確率を、前記第１のマスクおよび前記第２のマスクとするマスク推定部をさらに有することを特徴とする請求項１から５のいずれか１項に記載の空間相関行列推定装置。
【請求項７】
前記マスク推定部は、前記要素分布のうち、前記時不変のパラメータを要素に持つ正定値エルミート行列の固有値の分布の形状が最も平坦である要素分布の事後確率を前記第２のマスクとすることを特徴とする請求項６に記載の空間相関行列推定装置。
【請求項８】
Ｎ個の目的音源（ただし、Ｎは１以上の整数）に対応するＮ個の第１の音響信号と、背景雑音に対応する第２の音響信号とが混在する状況において、それぞれ異なる位置で収録されたＭ個（ただし、Ｍは２以上の整数）の観測信号に基づいて計算される観測特徴量ベクトルに基づいて、時間周波数点ごとの、前記第１の音響信号が前記観測信号の特徴量に含まれる割合である第１のマスクと、時間周波数点ごとの、前記第２の音響信号が前記観測信号の特徴量に含まれる割合である第２のマスクと、を推定し、前記第１のマスクおよび前記第２のマスクに基づいて前記目的音源の空間相関行列を推定する空間相関行列推定方法であって、
前記観測信号および前記第１のマスクに基づいて計算された第１の特徴量行列を第１の係数で重み付けした第１の空間相関行列と、前記観測信号および前記第２のマスクに基づいて計算された第２の特徴量行列を第２の係数で重み付けした第２の空間相関行列と、に基づいて前記目的音源の空間相関行列を推定する雑音除去工程を含んだことを特徴とする空間相関行列推定方法。
【請求項９】
前記雑音除去工程は、背景雑音の空間相関行列が時間的に変化しないとの条件下において、前記目的音源の空間相関行列の推定値に含まれる背景雑音由来の成分が０となるように、前記第１の係数及び前記第２の係数を計算することを特徴とする請求項８に記載の空間相関行列推定方法。
【請求項１０】
前記雑音除去工程は、前記第１の係数と前記第２の係数との比率が、前記第１のマスクの時間平均値の逆数と前記第２のマスクの時間平均値の逆数との比率に等しくなるように、前記第１の係数及び前記第２の係数を計算することを特徴とする請求項８または９に記載の空間相関行列推定方法。
【請求項１１】
前記観測信号に短時間信号分析を適用し、時間周波数点ごとの信号特徴量を抽出し、前記信号特徴量を成分とするＭ次元縦ベクトルである観測特徴量ベクトルを時間周波数点ごとに計算する時間周波数分析工程と、
前記観測特徴量ベクトルに基づいて、時間周波数点ごとに、前記観測特徴量ベクトルと該観測特徴量ベクトルのエルミート転置とを乗じて観測特徴量行列を計算する観測特徴量行列計算工程と、
前記目的音源のそれぞれについて、時間周波数点ごとに前記観測特徴量行列と前記第１のマスクを乗じて得られる行列の周波数ごとの時間平均を第１の特徴量行列として計算し、前記第１の特徴量行列に前記第１の係数を乗じることで前記第１の空間相関行列を推定する雑音下目的音空間相関行列推定工程と、
前記背景雑音について、時間周波数点ごとに前記観測特徴量行列と前記第２のマスクを乗じて得られる行列の周波数ごとの時間平均を第２の特徴量行列として計算し、前記第２の特徴量行列に前記第２の係数を乗じることで前記第２の空間相関行列を推定する雑音空間相関行列推定工程と、
をさらに含み、
前記雑音除去工程は、前記第１の空間相関行列から前記第２の空間相関行列を減じることで前記目的音源の空間相関行列を推定し、
前記第１の係数と前記第２の係数との比率は、前記第１のマスクの時間平均値の逆数と前記第２のマスクの時間平均値の逆数との比率に等しいことを特徴とする請求項８に記載の空間相関行列推定方法。
【請求項１２】
コンピュータに、請求項８から１１のいずれか１項に記載の空間相関行列推定方法を実行させるための空間相関行列推定プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、空間相関行列推定装置、空間相関行列推定方法および空間相関行列推定プログラムに関する。
【背景技術】
【０００２】
従来、目的音源から出た音響信号と背景雑音による音響信号とが混在する状況において、複数のマイクロホンで収音された観測信号から、各目的音源のみが観測信号に含まれている場合の空間相関行列を推定する方法が提案されている。また、空間相関行列を推定する際には、各音響信号が観測された音響信号に含まれる割合であるマスクが用いられる場合がある。【０００３】
なお、空間相関行列とは、マイクロホン間の信号の自己相関、および相互相関を表す行列であり、例えば目的音源の位置を推定することや、観測信号から目的音源のみを取り出すビームフォーマを設計することに用いられる。【０００４】
ここで、図６を用いて、従来の空間相関行列推定装置について説明する。図６は、従来の空間相関行列推定装置の構成を示す図である。図６に示すように、まず、時間周波数分析部１０ａは、観測信号から抽出した時間周波数点ごとの観測特徴量ベクトルを計算する。次に、マスク推定部２０ａは、観測特徴量ベクトルを基に目的音源および背景雑音に対応するマスクを推定する。また、観測特徴量行列計算部３０ａは、観測特徴量ベクトルと当該観測特徴量ベクトルのエルミート転置とを乗じて観測特徴量行列を計算する。【０００５】
そして、目的音特徴量行列時間平均計算部４０ａは、観測特徴量行列に目的音源に対応するマスクを乗じて得られた行列の時間平均である平均目的音特徴量行列を計算する。また、雑音特徴量行列時間平均計算部５０ａは、観測特徴量行列に背景雑音に対応するマスクを乗じて得られた行列の時間平均である平均雑音特徴量行列を計算する。最後に、目的音特徴量雑音除去部６０ａは、平均目的音特徴量行列から平均雑音特徴量行列を減じることで目的音源の空間相関行列を推定する。
【先行技術文献】
【非特許文献】
【０００６】
【非特許文献１】Mehrez  Souden,  Shoko  Araki,  Keisuke  Kinoshita,  Tomohiro  Nakatani,  Hiroshi  Sawada、“A  multichannel  MMSE-based  framework  for  speech  source  separation  and  noise  reduction,”  IEEE  Trans.  Audio,  Speech,  and  Language  Processing,  vol.  21,  no.  9,pp.  1913-1928,  2013.
【非特許文献２】Ozgur  Yilmaz,  and  Scott  Rickard,  “Blind  separation  of  speech  mixture  via  time-frequency  masking,”  IEEE  Trans.  Signal  Processing,  vol.  52,  no.  7,  pp.  1830-1847,  2004.
【非特許文献３】Dang  Hai  Tran  Vu  and  Reinhold  Haeb-Umbach,  “Blind  speech  separation  employing  directional  statistics  in  an  expectation  maximization  framework,”  Proc.IEEE  Int.  Conf.  Acoust.,  Speech,  Signal  Process.  (ICASSP-2010),  pp.  241-244,  2010.
【非特許文献４】Tomohiro  Nakatani,  Shoko  Araki,  Takuya  Yoshioka,  Marc  Delcroix,  and  Masakiyo  Fujimoto,  “Dominance  based  integration  of  spatial  and  spectral  features  for  speech  enhancement,”  IEEE  Transactions  on  Audio,  Speech,  and  Language  Processing,  vol.  21,  no.  12,  pp.  2516-2531,  Dec.  2013.
【発明の概要】
【発明が解決しようとする課題】
【０００７】
しかしながら、従来の空間相関行列の推定方法には、背景雑音の影響を観測信号から正確に取り除くことができないため、目的音源の空間相関行列を精度良く推定できない場合があるという問題があった。【０００８】
例えば、従来の空間相関行列の推定方法では、平均目的音特徴量行列から平均雑音特徴量行列を減じた結果を目的音源の空間相関行列として推定しているが、これは経験的に得られた方法であり、平均目的音特徴量行列に含まれる雑音の影響の量が平均雑音特徴量行列に一致するとは限らないため、雑音の影響がキャンセルされる保証はない。そのため、従来の空間相関行列の推定方法では、目的音源の空間相関行列を精度良く推定できない場合がある。
【課題を解決するための手段】
【０００９】
本発明の空間相関行列推定装置は、Ｎ個の目的音源（ただし、Ｎは１以上の整数）に対応するＮ個の第１の音響信号と、背景雑音に対応する第２の音響信号とが混在する状況において、それぞれ異なる位置で収録されたＭ個（ただし、Ｍは２以上の整数）の観測信号に基づいて計算される観測特徴量ベクトルに基づいて、時間周波数点ごとの、前記第１の音響信号が前記観測信号の特徴量に含まれる割合である第１のマスクと、時間周波数点ごとの、前記第２の音響信号が前記観測信号の特徴量に含まれる割合である第２のマスクと、を推定し、前記第１のマスクおよび前記第２のマスクに基づいて前記目的音源の空間相関行列を推定する空間相関行列推定装置であって、前記観測信号および前記第１のマスクに基づいて計算された第１の特徴量行列を第１の係数で重み付けした第１の空間相関行列と、前記観測信号および前記第２のマスクに基づいて計算された第２の特徴量行列を第２の係数で重み付けした第２の空間相関行列と、に基づいて前記目的音源の空間相関行列を推定する雑音除去部を有することを特徴とする。【００１０】
また、本発明の空間相関行列推定方法は、Ｎ個の目的音源（ただし、Ｎは１以上の整数）に対応するＮ個の第１の音響信号と、背景雑音に対応する第２の音響信号とが混在する状況において、それぞれ異なる位置で収録されたＭ個（ただし、Ｍは２以上の整数）の観測信号に基づいて計算される観測特徴量ベクトルに基づいて、時間周波数点ごとの、前記第１の音響信号が前記観測信号の特徴量に含まれる割合である第１のマスクと、時間周波数点ごとの、前記第２の音響信号が前記観測信号の特徴量に含まれる割合である第２のマスクと、を推定し、前記第１のマスクおよび前記第２のマスクに基づいて前記目的音源の空間相関行列を推定する空間相関行列推定方法であって、前記観測信号および前記第１のマスクに基づいて計算された第１の特徴量行列を第１の係数で重み付けした第１の空間相関行列と、前記観測信号および前記第２のマスクに基づいて計算された第２の特徴量行列を第２の係数で重み付けした第２の空間相関行列と、に基づいて前記目的音源の空間相関行列を推定する雑音除去工程を含んだことを特徴とする。
【発明の効果】
【００１１】
本発明によれば、背景雑音の影響を観測信号から正確に取り除き、目的音源の空間相関行列を精度良く推定できる。
【図面の簡単な説明】
【００１２】
【図１】図１は、第１の実施形態に係る空間相関行列推定装置の構成の一例を示す図である。
【図２】図２は、第１の実施形態に係る空間相関行列推定装置のマスク推定部の構成の一例を示す図である。
【図３】図３は、第１の実施形態に係る空間相関行列推定装置の処理の一例を示す図である。
【図４】図４は、第１の実施形態に係る空間相関行列推定装置のマスク推定処理の一例を示す図である。
【図５】図５は、プログラムが実行されることにより空間相関行列推定装置が実現されるコンピュータの一例を示す図である。
【図６】図６は、従来の空間相関行列推定装置の構成を示す図である。
【発明を実施するための形態】
【００１３】
以下に、本願に係る空間相関行列推定装置、空間相関行列推定方法および空間相関行列推定プログラムの実施形態を図面に基づいて詳細に説明する。なお、この実施形態により本発明が限定されるものではない。
【００１４】
［第１の実施形態］
まず、第１の実施形態に係る空間相関行列推定装置の構成、処理の流れおよび効果を説明する。なお、第１の実施形態においては、Ｎ個の目的音源（ただし、Ｎは１以上の整数）に対応するＮ個の第１の音響信号と、背景雑音に対応する第２の音響信号とが混在する状況において、それぞれ異なる位置で収録されたＭ個（ただし、Ｍは２以上の整数）の観測信号が空間相関行列推定装置に入力されるものとする。
【００１５】
［第１の実施形態の構成］
図１を用いて、第１の実施形態の構成について説明する。図１は、第１の実施形態に係る空間相関行列推定装置の構成の一例を示す図である。図１に示すように、空間相関行列推定装置１は、時間周波数分析部１０、マスク推定部２０、観測特徴量行列計算部３０、雑音下目的音空間相関行列推定部４０、雑音空間相関行列推定部５０および目的音空間相関行列雑音除去部６０を有する。
【００１６】
まず、空間相関行列推定装置１の各部の概要について説明する。時間周波数分析部１０は、入力された観測特徴量に基づいて観測特徴量ベクトルを計算する。具体的には、時間周波数分析部１０は、各観測信号ｙ（ｍ）（τ）に短時間信号分析を適用し、時間周波数点ごとの信号特徴量を抽出し、信号特徴量を成分とするＭ次元縦ベクトルである観測特徴量ベクトルｘ（ｔ、ｆ）を時間周波数点ごとに計算する。
【００１７】
また、マスク推定部２０は、時間周波数点ごとの、第１の音響信号が観測信号の特徴量に含まれる割合である第１のマスクφｎ（ｔ，ｆ）と、時間周波数点ごとの、第２の音響信号が観測信号の特徴量に含まれる割合である第２のマスクφｖ（ｔ，ｆ）と、を推定する。そして、観測特徴量行列計算部３０は、観測特徴量ベクトルに基づいて、時間周波数点ごとに、観測特徴量ベクトルと該観測特徴量ベクトルのエルミート転置とを乗じて観測特徴量行列Ｒｘｘ（ｔ，ｆ）を計算する。
【００１８】
雑音下目的音空間相関行列推定部４０は、観測信号および第１のマスクに基づいて計算された第１の特徴量行列を第１の係数で重み付けした第１の空間相関行列を計算する。具体的には、雑音下目的音空間相関行列推定部４０は、目的音源のそれぞれについて、時間周波数点ごとに観測特徴量行列と第１のマスクを乗じて得られる行列の周波数ごとの時間平均を第１の特徴量行列Ｒ´ｎ＋ｖ（ｔ，ｆ）として計算し、第１の特徴量行列に第１の係数αを乗じた結果を第１の空間相関行列Ｒｎ＋ｖ（ｔ，ｆ）とする。
【００１９】
雑音空間相関行列推定部５０は、観測信号および第２のマスクに基づいて計算された第２の特徴量行列を第２の係数で重み付けした第２の空間相関行列を計算する。具体的には、背景雑音について、時間周波数点ごとに観測特徴量行列と第２のマスクを乗じて得られる行列の周波数ごとの時間平均を第２の特徴量行列Ｒ´ｖ（ｔ，ｆ）として計算し、第２の特徴量行列に第２の係数βを乗じた結果を第２の空間相関行列Ｒｖ（ｔ，ｆ）とする。
【００２０】
雑音除去部として機能する目的音空間相関行列雑音除去部６０は、第１の空間相関行列および第２の空間相関行列に基づいて目的音源の空間相関行列を推定する。具体的には、目的音空間相関行列雑音除去部６０は、第１の空間相関行列から第２の空間相関行列を減じた結果を目的音源の空間相関行列Ｒｎ（ｔ，ｆ）とする。なお、第１の係数と第２の係数との比率は、例えば、第１のマスクの時間平均値の逆数と第２のマスクの時間平均値の逆数との比率に等しい。
【００２１】
次に、空間相関行列推定装置１の各部の詳細について説明する。目的音源はスパース性を有し、各時間周波数点において高々１つの目的音源だけが存在すると仮定される。また、背景雑音はすべての時間周波数点に存在すると仮定される。これより、時間周波数分析部１０が、入力された観測特徴量から短時間フーリエ変換等の短時間信号分析を用いて計算される観測特徴量ベクトルは、式（１）および式（２）のどちらかに一致することになる。
【００２２】
【数１】
2017094862.tif 000003
【００２３】
【数２】
2017094862.tif 000004
【００２４】
なお、式（１）および式（２）のｔとｆは、それぞれ時間と周波数の番号であり、ｔは１〜Ｔの整数、ｆは０〜Ｆの整数をとることとする。ここで、式（１）は当該時間周波数点において目的音源のうちｎ番目の音源だけが存在する場合、式（２）は目的音源が一つも存在しない場合を表し、ｓｎ（ｔ，ｆ）とｖ（ｔ，ｆ）は、観測特徴量ベクトルを目的音源ｎの成分と背景雑音に対応する成分の和に分解したものである。
【００２５】
マスク推定部２０は、既知のマスク推定技術を用いて、マスクを推定する。マスク推定部２０がｎ番目の目的音源に関して推定したマスクをφｎ（ｔ，ｆ）、背景雑音に関して推定したマスクをφｖ（ｔ，ｆ）とする。以下、添え字ｎはどの目的音源に対応するかを示す番号として、添え字ｖは雑音に対応することを表す記号とする。
【００２６】
雑音下目的音空間相関行列推定部４０は、式（３）によりｎ番目の目的音源に対応する第１の特徴量行列、すなわち平均目的音特徴量行列Ｒ´ｎ＋ｖ（ｆ）を計算する。
【００２７】
【数３】
2017094862.tif 000005
【００２８】
また、雑音空間相関行列推定部５０は、式（４）により背景雑音に対応する第２の特徴量行列、すなわち平均雑音特徴量行列Ｒ´ｖ（ｆ）を計算する。
【００２９】
【数４】
2017094862.tif 000006
【００３０】
ただし、観測特徴量行列Ｒｘｘ（ｔ，ｆ）は、式（５）のように表される。なお、式（５）のＨは行列のエルミート転置を表す。
【００３１】
【数５】
2017094862.tif 000007
【００３２】
式（１）および式（２）に示したように、背景雑音は、すべての時間周波数点に含まれているので、Ｒ´ｎ＋ｖ（ｆ）にも、雑音の影響は必ず含まれてしまう。Ｒ´ｎ＋ｖ（ｆ）の添え字ｎ＋ｖは、Ｒ´ｎ＋ｖ（ｆ）に目的音源ｎと雑音の両方の影響が含まれていることを表している。
【００３３】
ここで、式（１）に対応する時間周波数点だけを集めて空間相関行列を求めることができれば、その空間相関行列は目的音源ｎと背景雑音の影響だけを受けたものになる。一方、背景雑音の空間相関行列は、式（２）に対応する時間周波数点だけを集めて空間相関行列を求めることで得ることができる。
【００３４】
そこで、従来の空間相関行列推定方法においては、式（６）のように、それぞれ求めた空間相関行列の差を求めることで、目的音源の空間相関行列を求めていた。
【００３５】
【数６】
2017094862.tif 000008
【００３６】
一方、本発明の第１の実施形態においては、さらにこれらの空間相関行列に重み付けをしたうえで差を求める。ここで、各目的音源および背景雑音が無相関であると仮定すると、Ｒｘｘ（ｔ，ｆ）は式（７）により表される。
【００３７】
【数７】
2017094862.tif 000009
【００３８】
式（７）で、背景雑音由来の成分はｖ（ｔ，ｆ）ｖＨ（ｔ，ｆ）であることと、式（３）および式（４）を考慮すると、式（６）に残存する背景雑音由来の成分は式（８）により表される。
【００３９】
【数８】
2017094862.tif 000010
【００４０】
これより、式（８）の値が０になる場合に、目的音源の空間相関行列の推定値に残存する背景雑音の影響が０になるといえる。そこで、目的音空間相関行列雑音除去部６０は、式（９）のように、第１の係数αによって重み付けされた第１の空間相関行列、すなわち平均目的音特徴量行列Ｒ´ｎ＋ｖ（ｆ）、および第２の係数βによって重み付けされた第２の空間相関行列、すなわち平均雑音特徴量行列Ｒ´ｖ（ｔ，ｆ）を用いて目的音源の空間相関行列を計算する。
【００４１】
【数９】
2017094862.tif 000011
【００４２】
なお、Ｒ´ｎ＋ｖ（ｆ）を第１の係数αで重み付けしたＲｎ＋ｖ（ｆ）は雑音下目的音空間相関行列推定部４０によって計算され、Ｒ´ｖ（ｆ）を第２の係数βで重み付けしたＲｖ（ｆ）は雑音空間相関行列推定部５０によって計算される。
【００４３】
このとき、式（８）の目的音源の空間相関行列推定値に残存する背景雑音由来の成分は式（１０）により表される。
【００４４】
【数１０】
2017094862.tif 000012
【００４５】
式（１０）の値が０に一致するための必要十分条件は、式（１１）が成立することである。
【００４６】
【数１１】
2017094862.tif 000013
【００４７】
式（１１）のΣｔφｎ（ｔ，ｆ）ｖ（ｔ，ｆ）ｖＨ（ｔ，ｆ）／Σｔφｎ（ｔ，ｆ）とΣｔφｖ（ｔ，ｆ）ｖ（ｔ，ｆ）ｖＨ（ｔ，ｆ）／Σｔφｖ（ｔ，ｆ）は、雑音特徴量行列ｖ（ｔ，ｆ）ｖＨ（ｔ，ｆ）の重み付き時間平均を、異なる重みを用いて計算したものである。いま、背景雑音の空間相関行列は、時間的に大きく変化しないと仮定すると、これら２つの重み付き時間平均値は近似的に一致するといえる。その結果、式（１１）はさらに式（１２）のように書き換えられる。
【００４８】
【数１２】
2017094862.tif 000014
【００４９】
そして、式（１２）および式（９）より、式（１３）が得られる。
【００５０】
【数１３】
2017094862.tif 000015
【００５１】
式（１３）で、Ｔ／Σｔφｎ（ｔ，ｆ）は、目的音源ｎに対応するマスクの時間平均の逆数、Ｔ／Σｔφｖ（ｔ，ｆ）は、背景雑音に対応するマスクの時間平均の逆数であり、ｃはスカラー定数とする。ｃは、どの時間区間で目的音源の空間相関行列を求めるかで決まる定数で、全時間区間の場合はｃ＝Σｔφｎ（ｔ，ｆ）／Ｔとし、主に目的音源ｎが存在する時間区間で求める場合はｃ＝１とすればよい。
【００５２】
ｃ＝Σｔφｎ（ｔ，ｆ）／Ｔとした場合は、式（９）中で、α＝１とした場合に相当し、式（６）に対し、Ｒ´ｖ（ｆ）のゲインだけを変えて、目的音源に関する空間相関行列のゲインを変えずに雑音の影響を取り除く場合に相当する。
【００５３】
式（１３）を、さらに、式（３）、式（４）とともに整理すると、式（１４）〜（１６）が得られる。
【００５４】
【数１４】
2017094862.tif 000016
【００５５】
【数１５】
2017094862.tif 000017
【００５６】
【数１６】
2017094862.tif 000018
【００５７】
例えば、ｃ＝１の場合、式（１６）は、式（１７）のように表される。このように、背景雑音の空間相関行列は時間的に大きく変化しないとの仮定のもと、適切な係数を乗じたうえで差を求めることで、ｎ番目の目的音源に関して背景雑音の影響を正確に取り除いた空間相関行列を推定することができる。
【００５８】
【数１７】
2017094862.tif 000019
【００５９】
式（１４）は、雑音下目的音空間相関行列推定部４０が雑音下目的音空間相関行列Ｒｎ＋ｖ（ｆ）を推定する処理に相当する。また、式（１５）は、雑音空間相関行列推定部５０が雑音空間相関行列Ｒｖ（ｆ）を推定する処理に相当する。また、式（１６）は、目的音空間相関行列雑音除去部６０が目的音の空間相関行列Ｒｎ（ｆ）を推定する処理に相当する。
【００６０】
また、音源数Ｎ＝１であるときに、ｃを式（１８）のように定めた場合には、目的音源の空間相関行列は、式（１９）〜（２１）によって計算されてもよい。
【００６１】
【数１８】
2017094862.tif 000020
【００６２】
【数１９】
2017094862.tif 000021
【００６３】
【数２０】
2017094862.tif 000022
【００６４】
【数２１】
2017094862.tif 000023
【００６５】
式（１９）〜（２１）では目的音源のマスクφｎ（ｔ，ｆ）が使用されていないことから、目的音源のマスクを推定することなく目的音源の空間相関行列を推定することが可能であるといえる。この場合、式（１９）に示すように、雑音下目的音空間相関行列は、Ｎ＝１である場合、観測特徴量行列の周波数ごとの時間平均である。
【００６６】
マスク推定部２０は、観測特徴量ベクトルの確率分布を、共分散行列が時刻ごとに異なる値を取るスカラーパラメータと時不変のパラメータを要素にもつ正定値エルミート行列との積で表される、平均０のＭ次元複素ガウス分布であるＮ＋１個の要素分布からなる混合分布で周波数ごとにモデル化する。そして、マスク推定部２０は、混合分布が観測特徴量ベクトルの分布に近くなるように混合分布のパラメータを推定することで得られる要素分布のそれぞれの事後確率を、第１のマスクおよび第２のマスクとする。
【００６７】
これにより、マスク推定部２０は、観測特徴量ベクトルの分布の形状が、超球面上の円状では正確に近似できない場合でも、分布の形状を正確に近似し、正確なマスク推定を行う。
【００６８】
目的音源ｎが存在する時間周波数点の観測特徴量ベクトルの確率密度関数に対応する要素分布をｐｎ（ｘ（ｔ，ｆ）；Θ）、雑音のみが存在する時間周波数点の観測特徴量ベクトルの確率密度関数に対応する要素分布をｐｖ（ｘ（ｔ，ｆ）；Θ）とすると、マスク推定部２０は、それぞれの要素分布を、式（２２）および式（２３）のようにモデル化する。
【００６９】
【数２２】
2017094862.tif 000024
【００７０】
【数２３】
2017094862.tif 000025
【００７１】
ここで、Ｎｃ（ｘ；μ，Σ）は、平均ベクトルμ，共分散行列ΣのＭ次元複素ガウス分布である。式（２２）および式（２３）の要素分布の式中で、ｒｎ（ｔ，ｆ）、ｒｖ（ｔ，ｆ）は、各音響信号の大きさに対応するスカラーパラメータであり、時間周波数点ごとに異なる値を取ることができるように設定されている。
【００７２】
一方、Ｂｎ（ｆ）とＢｖ（ｆ）は、音響信号が空間的にどの方向から到来するかを表現する行列であり、時不変のパラメータを要素に持つ行列として規定されている。Ｂｎ（ｆ）とＢｖ（ｆ）は、要素分布の形状を決定するパラメータであり、上記のモデルでは、特段の制約を設けていない。このため、各要素分布は、Ｍ次元複素ガウス分布が表しうるあらゆる形状を持つことができ、超球面上の円状の分布に限定されない。
【００７３】
また、Θ＝｛ｒｎ（ｔ，ｆ），ｒｖ（ｔ，ｆ），Ｂｎ（ｆ），Ｂｖ（ｆ），λｎ（ｆ），λｖ（ｆ）｝は、上記の複素ガウス分布を要素分布として構成される混合分布のモデルパラメータの集合を表す。λｎ（ｆ）、λｖ（ｆ）は、それぞれ目的音源ｎの存在する時間周波数点に対応する要素分布の混合比、背景雑音のみが存在する時間周波数点に対応する要素分布の混合比であり、Σｎλｎ（ｆ）＋λｖ（ｆ）＝１、１＞λｎ（ｆ）＞０、１＞λｖ（ｆ）＞０を満たす。また、上記の要素分布からなる混合分布は式（２４）のように表される。
【００７４】
【数２４】
2017094862.tif 000026
【００７５】
マスク推定部２０は、上記混合モデルを用いて、すべての時間周波数点における観測特徴量ベクトルをモデル化し、上記の混合分布が観測特徴量ベクトルの確率分布に近くなるように各モデルパラメータを推定する。
【００７６】
マスク推定部２０は、モデルパラメータが推定されたのちに、目的音源ｎ，背景雑音のそれぞれに対応するマスクを、各要素分布の事後確率分布として式（２５）または式（２６）によって推定する。
【００７７】
【数２５】
2017094862.tif 000027
【００７８】
【数２６】
2017094862.tif 000028
【００７９】
各要素分布は、Ｍ次元複素ガウス分布の範囲であらゆる形状を持つことができるため、各要素分布は観測特徴量ベクトルの分布の形状が、超球面上の円では正確に近似できない場合でも、その形状を正確に近似することが可能になる。
【００８０】
ところで、一般に各目的音源ｎに対応する音響信号は、マイクロホン位置から見て音源のある方向（音源方向）から主に到来するという性質を持つ。このため、目的音源ｎに対応する要素分布の正定値エルミート行列は、音源方向に対応する部分空間に最大の固有値を持ち、それ以外の部分空間の固有値は比較的小さな値を持つという性質を持つ。
【００８１】
一方、背景雑音は、通常、あらゆる方向から音が到来するため、背景雑音に対応する要素分布の正定値エルミート行列は、全ての方向に対応する部分空間に行列の成分が分散する。このため、固有値が特定の部分空間に偏るようなことは生じにくい。
【００８２】
そこで、マスク推定部２０はさらに、要素分布のうち、時不変のパラメータを要素に持つ正定値エルミート行列の固有値の分布の形状が最も平坦である要素分布の事後確率を背景雑音に対応する第２のマスクとする。これにより、マスク推定部２０は、推定したマスクのうち、どれが背景雑音に対応するものであるかを自動的に推定することができる。
【００８３】
（実施例１）
第１の実施形態について、具体例を用いて説明する。まず、Ｎ＝１の場合、空間相関行列推定装置１は、例えば背景雑音下でＭ＝２以上のマイクで収録された、１人の人が話している声について、雑音の影響を除いた空間相関行列を推定する。また、Ｎ＞１の場合、空間相関行列推定装置１は、例えばＭ＞１個のマイクロホンで収録された、Ｎ人による会話について、雑音の影響を除いた空間相関行列を話者ごとに推定する。
【００８４】
ここで、マイクロホンｍで収録された観測信号は、ｙ（ｍ）（τ）と書くことにする。ｙ（ｍ）（τ）は、各音源信号ｎに由来する音響信号ｚｎ（ｍ）（τ）と背景雑音に由来する音響信号ｕ（ｍ）（τ）の和で構成されていることから、式（２７）のようにモデル化される。
【００８５】
【数２７】
2017094862.tif 000029
【００８６】
時間周波数分析部１０は、すべてのマイクロホンで収録された上記観測信号を受け取り、各観測信号ｙ（ｍ）（τ）ごとに短時間信号分析を適用して時間周波数ごとの信号特徴量ｘ（ｍ）（ｔ，ｆ）を求める。短時間信号分析としては、短時間離散フーリエ変換や短時間離散コサイン変換等の様々な方法を用いることができる。
【００８７】
時間周波数分析部１０は、さらに、各時間周波数で得られた信号特徴量ｘ（ｍ）（ｔ，ｆ）をすべてのマイクロホンに関してまとめたベクトルとして、式（２８）に示すような、観測特徴量ベクトルｘ（ｔ，ｆ）を構成する。
【００８８】
【数２８】
2017094862.tif 000030
【００８９】
次に、観測特徴量行列計算部３０は、観測特徴量ベクトルｘ（ｔ，ｆ）を受け取り、時間周波数点ごとに、観測特徴量行列Ｒｘｘ（ｔ，ｆ）を式（２９）によって求める。
【００９０】
【数２９】
2017094862.tif 000031
【００９１】
また、マスク推定部２０は、観測特徴量ベクトルｘ（ｔ，ｆ）を受け取り、時間周波数点ごとに、各目的音源と背景雑音がどのような割合で混ざっているかをマスクの値として推定する。なお、式（３０）に示すように、時間周波数点で、全目的音源と背景雑音に関するマスクの総和は１になると仮定する。
【００９２】
【数３０】
2017094862.tif 000032
【００９３】
雑音下目的音空間相関行列推定部４０は、各目的音源に関するマスクの推定値φｎ（ｔ，ｆ）と観測特徴量行列Ｒｘｘ（ｔ，ｆ）を受け取り、各目的音源ｎに対し、周波数ｆごとに、雑音下目的音空間相関行列Ｒｎ＋ｖ（ｆ）を式（３１）のように求める。
【００９４】
【数３１】
2017094862.tif 000033
【００９５】
雑音空間相関行列推定部５０は、背景雑音に関するマスクの推定値φｖ（ｔ，ｆ）と観測特徴量行列Ｒｘｘ（ｔ，ｆ）を受け取り、周波数ｆごとに、雑音空間相関行列Ｒｖ（ｆ）を式（３２）のように求める。
【００９６】
【数３２】
2017094862.tif 000034
【００９７】
目的音空間相関行列雑音除去部６０は、雑音下目的音空間相関行列の推定値Ｒｎ＋ｖ（ｆ）と雑音空間相関行列の指定値Ｒｖ（ｆ）を受け取り、各目的音源ｎに対し、周波数ｆごとに、目的音の空間相関行列Ｒｎ（ｆ）を式（３３）により求める。
【００９８】
【数３３】
2017094862.tif 000035
【００９９】
求めた空間相関行列は様々な用途に利用できる。例えば、目的音源ｎの空間相関行列の最大固有値に対応する固有ベクトルは、目的音源ｎからマイクロホンまでの空間伝達特性を表すステアリングベクトルと一致する。さらに、このようにして推定されたステアリングベクトルｈｎ（ｆ）と、式（３４）に示す観測信号自身の空間相関行列Ｒｘ（ｆ）とから、最小分散無歪応答（ＭＶＤＲ：Minimum  Variance  Distortionless  Response）フィルタｗｎ（ｆ）を式（３５）のように求めることができる。
【０１００】
【数３４】
2017094862.tif 000036
【０１０１】
【数３５】
2017094862.tif 000037
【０１０２】
このＭＶＤＲフィルタを観測特徴量ベクトルｘ（ｔ，ｆ）に適用することで、目的音源ｎ以外の音源や背景雑音の成分を抑圧し、式（３６）に示すように、目的音源ｎに対応する信号特徴量の推定値ｓｎ（ｔ，ｆ）を得ることができる。
【０１０３】
【数３６】
2017094862.tif 000038
【０１０４】
また、目的音源ｎの空間相関行列Ｒｎ（ｆ）と観測信号の空間相関行列Ｒｘ（ｆ）が求められているとき、多チャンネルウィナフィルタＷｎ（ｆ）を式（３７）のように構成することができる。
【０１０５】
【数３７】
2017094862.tif 000039
【０１０６】
この多チャンネルウィナフィルタＷｎ（ｆ）を観測特徴量ベクトルｘ（ｔ，ｆ）に適用することで、目的音源ｎ以外の音源や背景雑音の成分を抑圧し、式（３８）に示すように、目的音源ｎに対応する特徴量ベクトルの推定値ｓｎ（ｔ，ｆ）を得ることができる。
【０１０７】
【数３８】
2017094862.tif 000040
【０１０８】
（実施例２）
次に、マスク推定部２０の具体例について図２を用いて説明する。図２は、第１の実施形態に係る空間相関行列推定装置のマスク推定部の構成の一例を示す図である。マスク推定部２０は、観測特徴量ベクトルの確率分布を、混合複素ガウス分布を用いてモデル化することで、マスクを推定する。
【０１０９】
まず、マスク推定部２０は、各周波数ｆにおける観測信号ｘ（ｔ，ｆ）の生成分布について、混合複素ガウス分布を用いて式（３９）のようにモデル化を行う。
【０１１０】
【数３９】
2017094862.tif 000041
【０１１１】
ここでΘ＝｛λｎ（ｆ），λｖ（ｆ），ｒｎ（ｔ，ｆ），ｒｖ（ｔ，ｆ），Ｂｎ（ｆ），Ｂｖ（ｆ）｝は、混合複素ガウス分布のパラメータ集合である。λｎ（ｆ）とλｖ（ｆ）は、ｎ番目の音源と背景雑音にそれぞれ対応する複素ガウス分布の混合重みを表すパラメータであり、式（４０）を満たす。ｒｎ（ｔ，ｆ）とｒｖ（ｔ，ｆ）は、それぞれｎ番目の音源と背景雑音の、時間周波数点（ｔ，ｆ）におけるパワーの期待値を表すスカラーパラメータである。
【０１１２】
【数４０】
2017094862.tif 000042
【０１１３】
Ｂｎ（ｆ）とＢｖ（ｆ）はそれぞれパワーで正規化されたｎ番目の音源と背景雑音の時不変な空間相関行列である。ここでＢｎ（ｆ）とＢｖ（ｆ）は観測特徴量ベクトルの分布を決定するパラメータとなるが、このパラメータをフルランクの行列として求めることで、超球面上の円状では正確に近似できない場合においても、観測特徴量ベクトルの分布をより正確に近似することができる。
【０１１４】
事後確率推定部２０１は、式（３９）の確率分布に基づき、観測信号ｘ（ｔ，ｆ）がそれぞれの要素分布から生起された確率を求めることで、マスクの推定を行う。まず、パラメータ初期化部２０３は、各パラメータの初期値を設定し、設定した初期値をパラメータ保持部２０４に保持しておく。パラメータ初期化部２０３は、例えば乱数によりパラメータの初期値を決定する。
【０１１５】
次に、事後確率推定部２０１は、入力データ（観測信号）と現在の分布パラメータを用いて、それぞれの要素分布に関する事後確率を式（４１）および式（４２）のように計算する。ここで計算された事後確率が各周波数点のマスクに相当する。
【０１１６】
【数４１】
2017094862.tif 000043
【０１１７】
【数４２】
2017094862.tif 000044
【０１１８】
次に、パラメータ更新部２０２は、ＥＭアルゴリズムに基づいて分布パラメータを更新する。このとき、パラメータ更新部２０２は、最尤推定のためのコスト関数を式（４３）のように設定する。
【０１１９】
【数４３】
2017094862.tif 000045
【０１２０】
また、パラメータ更新部２０２は、事後確率推定部２０１で推定した事後確率を用いて、Ｑ関数を式（４４）のように設定する。
【０１２１】
【数４４】
2017094862.tif 000046
【０１２２】
ここで、Θｔは、ｔ回目の反復更新で得られたパラメータを示す。また、φｎ（ｔ，ｆ）とφｖ（ｔ，ｆ）は、式（３６）および式（３７）で与えられる。パラメータ更新部２０２は、式（４５）に示す条件下で、式（４４）のＱ関数をそれぞれのパラメータで偏微分したものを０と置くことで、式（４６）〜式（４８）に示すパラメータ更新則を導く。
【０１２３】
【数４５】
2017094862.tif 000047
【０１２４】
【数４６】
2017094862.tif 000048
【０１２５】
【数４７】
2017094862.tif 000049
【０１２６】
【数４８】
2017094862.tif 000050
【０１２７】
これにより、パラメータ更新部２０２は、分布パラメータΘを更新する。なお、Θに対して適切な事前分布を設定することで、既知の方法を用い、より精度良いマスク推定を実現することもできる。
【０１２８】
また、パラメータ更新部２０２は、分布パラメータの更新をオンラインで実施してもよい。この場合、パラメータ更新部２０２は、式（４７）で与えられる更新則を、時刻ｔ´において、１つ前の時刻ｔ´−１における推定値Ｂｎ（ｔ´−１，ｆ）を用いて式（４９）のように表す。
【０１２９】
【数４９】
2017094862.tif 000051
【０１３０】
また、パラメータ更新部２０２は、式（４８）で与えられる更新則を、同様に式（５０）のように表す。
【０１３１】
【数５０】
2017094862.tif 000052
【０１３２】
次に、パラメータ更新部２０２は、更新則を用いて更新した新たなパラメータを、パラメータ保持部２０４にコピーする。そして、マスク推定部２０は、事後確率推定部２０１、パラメータ更新部２０２、パラメータ保持部２０４の処理が決められた回数（例えば３０回）実行されるまで、もしくは計算結果が収束するまで反復する。
【０１３３】
（実施例３）
実施例３では、実施例２のマスク推定方法で発生するパーミュテーション問題の解決方法について説明する。実施例２において、マスク推定部２０は、周波数ｆごとに、マスクφｎ（ｔ，ｆ）とφｖ（ｔ，ｆ）を求めていた。しかし、各周波数で推定されたマスクにおいて、雑音に対応するマスクが目的音源のマスクと入れ替わっていたり、異なる周波数間で、同じ目的音源に対応するマスクが、異なる目的音源番号に対応付けられたりすることが起きる。
【０１３４】
このため、目的音源ごとに空間相関行列を正しく推定するためには、マスク推定部２０は、背景雑音に対応するマスクがどれであるかを正しく定め、また、異なる周波数間において同じ目的音源を同じ音源番号に対応付ける必要がある。ここでは、この問題をパーミュテーション問題と呼ぶ。
【０１３５】
パーミュテーション問題を解決するためには、マスク推定部２０は、以下の（１）および（２）の操作を行うことが必要である。
（１）各周波数において、どのマスクが背景雑音に対応するかを定める。
（２）異なる周波数間で、同じ目的音源に対応するマスクが、同じ音源番号に関連付けられるようにする。
【０１３６】
まず、（１）の操作について説明する。いま、実施例２の方法に従い、各周波数ｆにおいて、Ｎ個のＢｎ（ｆ）と１個のＢｖ（ｆ）が求められているとする。以下、説明を簡単にするため、Ｂ０（ｆ）＝Ｂｖ（ｆ）と表記する。ここで、マスク推定部２０は、Ｎ＋１個のＢｎ（ｆ）（Ｎ≧ｎ≧０）のうち、どのＢｎ（ｆ）が背景雑音に対応するかを以下の（１−１）〜（１−３）により決定する。
【０１３７】
（１−１）
各ｎに対し、Ｂｎ（ｆ）のＭ個の固有値を求め、それらを値の大きいものから順に並べてできるベクトルγｎ（ｆ）を式（５１）のように構成する。
【０１３８】
【数５１】
2017094862.tif 000053
【０１３９】
（１−２）
γｎ（ｆ）の分布の平坦さの度合いを評価する関数Ｅ（・）を用意し、その値が最も大きいｎに対応する番号ｎｖを、式（５２）により求める。
【０１４０】
【数５２】
2017094862.tif 000054
【０１４１】
（１−３）
ｎｖに対応するマスクを、背景雑音に対応するマスクとして定める。Ｅ（・）の定め方としては、例えば、式（５３）に示すような、ベクトルの要素を足して１になるように正規化したγｎ（ｆ）のエントロピーを求める関数として、式（５４）のように定めることが可能である。
【０１４２】
【数５３】
2017094862.tif 000055
【０１４３】
【数５４】
2017094862.tif 000056
【０１４４】
ここで、Ｈ（・）は、要素を足して１になるベクトルｕ＝［ｕ１，ｕ２，…，ｕＭ］のエントロピーを求める関数で、式（５５）のように定義される。
【０１４５】
【数５５】
2017094862.tif 000057
【０１４６】
次に、（２）の操作について説明する。まず、マスク推定部２０は、推定されたＮ個のマスクについて、全ての周波数において、同一の目的音源ｎに対応するマスクφｎ（ｔ，ｆ）が、同一の目的音源の番号ｎに関連付ける必要がある。具体的な手段としては、下記の（２−１）〜（２−４）が考えられる。
【０１４７】
（２−１）
会話への参加人数Ｎが既知であるとし、マスク推定部２０は、実施例２の方法で推定されたマスクのうち、背景雑音のマスクを除いたＮ個のマスクをφｎ（ｔ，ｆ）（ｎ＝１，…，Ｎ）とする。
【０１４８】
ここで、マスクは、当該目的信号が、各時間周波数点にどの程度含まれているかの割合を表すものであるため、ある１つの音源のマスクの時系列は全ての周波数で同期する傾向がある。この性質を用いて、マスク推定部２０は、得られたマスクの全てのｎおよびｆでの時系列φｎ（ｔ，ｆ）（ｔ＝１，・・・，Ｔ）をＮ個のクラスタにクラスタリングすることで、パーミュテーション問題を解決する。クラスタリングには、例えばｋ−ｍｅａｎｓ法を用いたり、参考文献１（H.  Sawada,  S.  Araki,  S.  Makino,  “Underdetermined  Convolutive  Blind  Source  Separation  via  Frequency  Bin-Wise  Clustering  and  Permutation  Alignment,”  IEEE  Trans.  Audio,  Speech,  and  Language  Processing,  vol.19,  no.3,  pp.516-527,  March  2011.）に記載された方法を用いることができる。
【０１４９】
（２−２）
マスク推定部２０は、式（４１）および式（４２）によるマスクの推定において、Ｂｎ（ｆ）を、話者の位置ごとにあらかじめ学習した空間相関行列Ｂｎｔｒａｉｎｅｄ（ｆ）に固定する。Ｂｎｔｒａｉｎｅｄ（ｆ）は、例えば、学習データとして話者の位置ごとの観測信号をあらかじめ用意し、その学習データにて実施例２の方法でマスクを推定し、式（４７）の結果として得られたＢｎ（ｆ）である。
【０１５０】
この手段は、椅子の位置がほぼ固定されている会議室等での会話に有効であり、各席に対応する話者を目的音源ｎとして、それに対応するマスクφｎ（ｔ，ｆ）を推定できる。
【０１５１】
（２−３）
手段（２−３）において、マスク推定部２０は、手段（２−２）において、Ｂｎ（ｆ）の初期値をＢｎｔｒａｉｎｅｄ（ｆ）とし、実施例２の方法でマスクを推定する。手段（２−２）は、椅子の位置がほぼ固定されているが、椅子にキャスターがついていること等により話者の位置が会話中に少しずつ変動するような場合に有効である。
【０１５２】
（２−４）
手段（２−４）において、マスク推定部２０は、Ｂｎｔｒａｉｎｅｄ（ｆ）をＢｎ（ｆ）の事前情報として用いながらマスクの推定を行う。具体的に、マスク推定部２０は、式（４７）の推定を、η（０〜１までの間の実数）を重みとして、式（５６）により行う。
【０１５３】
【数５６】
2017094862.tif 000058
【０１５４】
手段（２−３）は、手段（２−２）と同様に、椅子の位置がほぼ固定されているが、椅子にキャスターがついていること等により話者の位置が会話中に少しずつ変動するような場合に有効である。
【０１５５】
（実施例４）
実施例４として、空間相関行列推定装置１によって得られた目的音源の空間相関行列を用いて、方向推定を行なう場合について説明する。まず、実施例１と同様の手順で、音源ｎに関するステアリングベクトルが式（５７）のように得られているとする。
【０１５６】
【数５７】
2017094862.tif 000059
【０１５７】
次に、参考文献２（S.  Araki,  H.  Sawada,  R.  Mukai  and  S.  Makino,“DOA  estimation  for  multiple  sparse  sourceswith  normalized  observation  vector  clustering,”,  ICASSP2006,  Vol.  5,  pp.33-36,  2006.）のように、Ｍ個のマイク配置が既知でマイクｍの３次元座標をｄｍとし、マイクアレイから見た音源ｎの方位角をθｎ、仰角をψｎとすると、ｑｎ＝［ｃｏｓ（θｎ）ｃｏｓ（ψｎ），ｃｏｓ（θｎ）ｓｉｎ（ψｎ），ｓｉｎ（ψｎ）］Ｔは、式（５８）により計算できる。
【０１５８】
【数５８】
2017094862.tif 000060
【０１５９】
ここで、ｃは音速、ｆバーは周波数インデックスｆに対応する周波数（Ｈｚ）、ξｎ（ｆ）＝［ａｒｇ（ｈｎ１／ｈｎＪ），…，ａｒｇ（ｈｎＭ／ｈｎＪ）］Ｔ、Ｄ＝［ｄ１-ｄＪ，…，ｄＭ-ｄＪ］Ｔ、Ｊは基準マイクのインデックス（１〜Ｍから任意に選択）であり、＋は一般化逆行列を示す。
【０１６０】
そして、式（５８）によって得られた到来方向ｑｎ（ｆ）について、空間的エリアシングが発生しない周波数範囲のｑｎ（ｆ）の平均値をもって、音源ｎの到来方向ｑｎとする。またｑではなく、方位角、仰角等の平均値を計算しても良い。
【０１６１】
［第１の実施形態の処理］
図３を用いて、第１の実施形態の空間相関行列推定装置１の処理について説明する。図３は、第１の実施形態に係る空間相関行列推定装置の処理の一例を示す図である。まず、図３に示すように、時間周波数分析部１０は、観測信号を取得し（ステップＳ１０）、短時間フーリエ変換等の短時間信号分析を用いて時間周波数点ごとの信号特徴量を計算し（ステップＳ１１）、観測特徴量ベクトルを構成する（ステップＳ１２）。
【０１６２】
次に、観測特徴量行列計算部３０は、観測特徴量ベクトルを基に、時間周波数点ごとの観測特徴量行列を計算する（ステップＳ１３）。そして、マスク推定部２０は、観測特徴量ベクトルを基に、マスクを推定する（ステップＳ１４）。
【０１６３】
雑音下目的音空間相関行列推定部４０は、観測特徴量行列に目的音に対応するマスクを適用し、所定の係数で重み付けすることで、雑音下目的音空間相関行列を推定する（ステップＳ１５）。また、雑音空間相関行列推定部５０は、観測特徴量行列に背景雑音に対応するマスクを適用し、所定の係数で重み付けすることで、雑音空間相関行列を推定する（ステップＳ１６）。
【０１６４】
このとき、雑音下目的音空間相関行列の推定に用いられる係数と、雑音空間相関行列の推定に用いられる係数との比は、例えば、目的音に対応するマスクの時間平均の逆数と、背景雑音に対応するマスクの時間平均の逆数との比に等しい。
【０１６５】
最後に、目的音空間相関行列雑音除去部６０は、例えば雑音下目的音空間相関行列から雑音空間相関行列を減じることにより、目的音の空間相関行列を推定する（ステップＳ１７）。
【０１６６】
また、図３のステップＳ１４のマスク推定処理の例を、図４を用いて説明する。図４は、第１の実施形態に係る空間相関行列推定装置のマスク推定処理の一例を示す図である。まず、マスク推定部２０は、観測信号の生成分布を、混合複素ガウス分布を用いてモデル化する（ステップＳ１４１）。
【０１６７】
パラメータ初期化部２０３は、乱数等によりモデルのパラメータの初期値を設定する（ステップＳ１４２）。次に、事後確率推定部２０１は、観測信号とパラメータを用いて各要素分布に関する事後確率を計算する（ステップＳ１４３）。ここで、事後確率の計算が３０回行われていない場合（ステップＳ１４４、Ｎｏ）は、パラメータ更新部２０２は、計算した事後確率を用いてパラメータを更新する（ステップＳ１４５）。さらに、マスク推定部２０は、ステップＳ１４３に戻り処理を繰り返す。
【０１６８】
そして、事後確率の計算が３０回行われた場合は（ステップＳ１４４、Ｙｅｓ）、パラメータ更新部２０２は、最後のパラメータ更新を行う。最後に、マスク推定部２０は、計算した事後確率をマスクとして推定する（ステップＳ１４６）。
【０１６９】
［第１の実施形態の効果］
本発明の効果を確認するために、従来の方法および第１の実施形態を用いた確認実験について説明する。
【０１７０】
（確認実験１）
確認実験１では、バスの中、カフェ等の背景雑音の存在する環境下において、１人の話者（Ｎ＝１）がタブレットに向かって文章を読み上げている状況で、タブレットに装着されたＭ＝６個のマイクで信号を収録した。このとき、収録した信号に対して、各方法を用いて音声認識を行った場合の音声認識精度は下記の通りであった。下記の結果より、第１の実施形態を適用することで、音声認識精度が向上することが確認できた。
（１）そのまま音声認識をした場合：８７．１１（％）
（２）Ｗａｔｓｏｎ分布でマスク推定をした後、ＭＶＤＲを適応した場合（従来の方法）：８９．４０（％）
（３）第１の実施形態を適用し、オフラインでマスク推定した後、ＭＶＤＲを適応した場合（実施例１、オフライン）：９１．５４（％）
（４）第１の実施形態を適用し、事前学習したパラメータを初期値として、オンラインでマスク推定した後、ＭＶＤＲを適応した場合（実施例１、オンライン）：９１．８０（％）
【０１７１】
（確認実験２）
確認実験２では、通常の会議室において、４人の話者（Ｎ＝４）が直径１．２ｍの円卓を囲んで自由に会話している状況で、円卓中央のＭ＝８個のマイクで信号を収録した。このとき、収録した信号に対して、各方法を用いて音声認識を行った場合の音声認識精度は下記の通りであった。下記の結果より、第１の実施形態を適用することで、音声認識精度が向上することが確認できた。
（１）そのまま音声認識をした場合：２０．９（％）
（２）第１の実施形態を適用し、オフラインでマスク推定した後、ＭＶＤＲを適応した場合（実施例１、オフライン）：５４．０（％）
（３）第１の実施形態を適用し、オンラインでマスク推定した後、ＭＶＤＲを適応した場合（実施例１、オンライン）：５２．０（％）
【０１７２】
時間周波数分析部１０は、入力された観測特徴量に基づいて観測特徴量ベクトルを計算する。また、マスク推定部２０は、時間周波数点ごとの、第１の音響信号が観測信号の特徴量に含まれる割合である第１のマスクと、時間周波数点ごとの、第２の音響信号が観測信号の特徴量に含まれる割合である第２のマスクと、を推定する。そして、観測特徴量行列計算部３０は、観測特徴量ベクトルに基づいて、時間周波数点ごとに、観測特徴量ベクトルと該観測特徴量ベクトルのエルミート転置とを乗じて観測特徴量行列を計算する。
【０１７３】
雑音下目的音空間相関行列推定部４０は、観測信号および第１のマスクに基づいて計算された第１の特徴量行列を第１の係数で重み付けした第１の空間相関行列を計算する。また、雑音空間相関行列推定部５０は、観測信号および第２のマスクに基づいて計算された第２の特徴量行列を第２の係数で重み付けした第２の空間相関行列を計算する。そして、目的音空間相関行列雑音除去部６０は、第１の空間相関行列および第２の空間相関行列に基づいて目的音源の空間相関行列を推定する。
【０１７４】
このように、第１の実施形態によれば、第１の係数および第２の係数による適切な重み付けが行われているため、第１の特徴量行列および第２の特徴量行列をそのまま用いる場合と比較して、背景雑音の影響を観測信号から正確に取り除き、目的音源の空間相関行列を精度良く推定できる。
【０１７５】
また、第１の係数と第２の係数との比率は、例えば、第１のマスクの時間平均値の逆数と第２のマスクの時間平均値の逆数との比率に等しいこととしてもよい。これにより、推定される目的音源の空間相関行列に、背景雑音の空間相関行列が時間的に大きく変化しないことが盛り込まれ、推定精度が向上する。
【０１７６】
また、マスク推定部２０は、観測特徴量ベクトルの確率分布を、共分散行列が時刻ごとに異なる値を取るスカラーパラメータと時不変のパラメータとを要素にもつ正定値エルミート行列の積で表される、平均０のＭ次元複素ガウス分布であるＮ＋１個の要素分布からなる混合分布で周波数ごとにモデル化する。
【０１７７】
そして、マスク推定部２０は、混合分布が観測特徴量ベクトルの分布に近くなるように混合分布のパラメータを推定することで得られる要素分布のそれぞれの事後確率を、第１のマスクおよび第２のマスクとする。これにより、観測特徴量ベクトルの分布の形状が、超球面上の円状では正確に近似できない場合でも、正確にマスクを推定することができる。
【０１７８】
マスク推定部２０はさらに、要素分布のうち、時不変のパラメータを要素に持つ正定値エルミート行列の固有値の分布の形状が最も平坦である要素分布の事後確率を背景雑音に対応する第２のマスクとする。これにより、マスク推定部が推定したマスクのうち、どれが背景雑音に対応するものであるかを自動的に推定することができる。
【０１７９】
［システム構成等］
また、図示した各装置の各構成要素は機能概念的なものであり、必ずしも物理的に図示の如く構成されていることを要しない。すなわち、各装置の分散・統合の具体的形態は図示のものに限られず、その全部または一部を、各種の負荷や使用状況等に応じて、任意の単位で機能的または物理的に分散・統合して構成することができる。さらに、各装置にて行なわれる各処理機能は、その全部または任意の一部が、ＣＰＵ（Central  Processing  Unit）および当該ＣＰＵにて解析実行されるプログラムにて実現され、あるいは、ワイヤードロジックによるハードウェアとして実現され得る。
【０１８０】
また、本実施形態において説明した各処理のうち、自動的におこなわれるものとして説明した処理の全部または一部を手動的におこなうこともでき、あるいは、手動的におこなわれるものとして説明した処理の全部または一部を公知の方法で自動的におこなうこともできる。この他、上記文書中や図面中で示した処理手順、制御手順、具体的名称、各種のデータやパラメータを含む情報については、特記する場合を除いて任意に変更することができる。
【０１８１】
［プログラム］
一実施形態として、空間相関行列推定装置は、パッケージソフトウェアやオンラインソフトウェアとして上記の空間相関行列推定を実行する空間相関行列推定プログラムを所望のコンピュータにインストールさせることによって実装できる。例えば、上記の空間相関行列推定プログラムを情報処理装置に実行させることにより、情報処理装置を空間相関行列推定装置として機能させることができる。ここで言う情報処理装置には、デスクトップ型またはノート型のパーソナルコンピュータが含まれる。また、その他にも、情報処理装置にはスマートフォン、携帯電話機やＰＨＳ（Personal  Handyphone  System）等の移動体通信端末、さらには、ＰＤＡ（Personal  Digital  Assistant）等のスレート端末等がその範疇に含まれる。
【０１８２】
また、空間相関行列推定装置は、ユーザが使用する端末装置をクライアントとし、当該クライアントに上記の空間相関行列推定に関するサービスを提供するサーバ装置として実装することもできる。例えば、空間相関行列推定装置は、観測信号を入力とし、目的音源の空間相関行列を出力とする空間相関行列推定サービスを提供するサーバ装置として実装される。この場合、空間相関行列推定装置は、Ｗｅｂサーバとして実装することとしてもよいし、アウトソーシングによって上記の空間相関行列推定に関するサービスを提供するクラウドとして実装することとしてもかまわない。
【０１８３】
図５は、プログラムが実行されることにより空間相関行列推定装置が実現されるコンピュータの一例を示す図である。コンピュータ１０００は、例えば、メモリ１０１０、ＣＰＵ１０２０を有する。また、コンピュータ１０００は、ハードディスクドライブインタフェース１０３０、ディスクドライブインタフェース１０４０、シリアルポートインタフェース１０５０、ビデオアダプタ１０６０、ネットワークインタフェース１０７０を有する。これらの各部は、バス１０８０によって接続される。
【０１８４】
メモリ１０１０は、ＲＯＭ（Read  Only  Memory）１０１１およびＲＡＭ（Random  Access  Memory）１０１２を含む。ＲＯＭ１０１１は、例えば、ＢＩＯＳ（Basic  Input  Output  System）等のブートプログラムを記憶する。ハードディスクドライブインタフェース１０３０は、ハードディスクドライブ１０９０に接続される。ディスクドライブインタフェース１０４０は、ディスクドライブ１１００に接続される。例えば磁気ディスクや光ディスク等の着脱可能な記憶媒体が、ディスクドライブ１１００に挿入される。シリアルポートインタフェース１０５０は、例えばマウス１１１０、キーボード１１２０に接続される。ビデオアダプタ１０６０は、例えばディスプレイ１１３０に接続される。
【０１８５】
ハードディスクドライブ１０９０は、例えば、ＯＳ１０９１、アプリケーションプログラム１０９２、プログラムモジュール１０９３、プログラムデータ１０９４を記憶する。すなわち、空間相関行列推定装置１の各処理を規定するプログラムは、コンピュータにより実行可能なコードが記述されたプログラムモジュール１０９３として実装される。プログラムモジュール１０９３は、例えばハードディスクドライブ１０９０に記憶される。例えば、空間相関行列推定装置１における機能構成と同様の処理を実行するためのプログラムモジュール１０９３が、ハードディスクドライブ１０９０に記憶される。なお、ハードディスクドライブ１０９０は、ＳＳＤ（Solid  State  Drive）により代替されてもよい。
【０１８６】
また、上述した実施形態の処理で用いられる設定データは、プログラムデータ１０９４として、例えばメモリ１０１０やハードディスクドライブ１０９０に記憶される。そして、ＣＰＵ１０２０が、メモリ１０１０やハードディスクドライブ１０９０に記憶されたプログラムモジュール１０９３やプログラムデータ１０９４を必要に応じてＲＡＭ１０１２に読み出して実行する。
【０１８７】
なお、プログラムモジュール１０９３やプログラムデータ１０９４は、ハードディスクドライブ１０９０に記憶される場合に限らず、例えば着脱可能な記憶媒体に記憶され、ディスクドライブ１１００等を介してＣＰＵ１０２０によって読み出されてもよい。あるいは、プログラムモジュール１０９３およびプログラムデータ１０９４は、ネットワーク（ＬＡＮ（Local  Area  Network）、ＷＡＮ（Wide  Area  Network）等）を介して接続された他のコンピュータに記憶されてもよい。そして、プログラムモジュール１０９３およびプログラムデータ１０９４は、他のコンピュータから、ネットワークインタフェース１０７０を介してＣＰＵ１０２０によって読み出されてもよい。
【符号の説明】
【０１８８】
１  空間相関行列推定装置
１０  時間周波数分析部
２０  マスク推定部
３０  観測特徴量行列計算部
４０  雑音下目的音空間相関行列推定部
５０  雑音空間相関行列推定部
６０  目的音空間相関行列雑音除去部
２０１  事後確率推定部
２０２  パラメータ更新部
２０３  パラメータ初期化部
２０４  パラメータ保持部
【図１】
2017094862.tif 000061
【図２】
2017094862.tif 000062
【図３】
2017094862.tif 000063
【図４】
2017094862.tif 000064
【図５】
2017094862.tif 000065
【図６】
2017094862.tif 000066
【国際調査報告】
2017094862.tif 000067
2017094862.tif 000068