(57)【要約】  マスク推定装置（１０）の特徴抽出部（１１）は、複数の音響信号を異なる位置で観測した複数の観測信号から、観測信号の時間周波数成分を時間周波数点毎にまとめた特徴量ベクトルを抽出する。マスク更新部（１２）は、特徴量ベクトルと、特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化する際の各要素分布の混合重みと、各要素分布の形状を制御可能なモデルパラメータである形状パラメータとから、各時間周波数点に対して各要素分布が寄与する割合を示すマスクを推定する。混合重み更新部（１３）は、更新されたマスクから、混合重みを更新する。パラメータ更新部（１４）は、特徴量ベクトルと、マスクとを用いて、形状パラメータを更新する。
2017141542.tif 000002
【特許請求の範囲】
【請求項１】
複数の音響信号が混在する状況下で、異なる位置で取得した複数の観測信号から、前記観測信号の時間周波数成分を時間周波数点毎にまとめた特徴量ベクトルを抽出する特徴抽出部と、
前記特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化する際の各要素分布の混合重みと、各要素分布の形状を制御可能なモデルパラメータである形状パラメータとを記憶する記憶部と、
前記特徴量ベクトルと、前記混合重みと、前記形状パラメータとを用いて、各時間周波数点に対して各要素分布が寄与する割合を要素分布毎かつ時間周波数点毎に示すマスクを推定して更新するマスク更新部と、
前記マスク更新部により更新されたマスクから、前記混合重みを更新する混合重み更新部と、
前記特徴量ベクトルと、前記マスクとを用いて、前記形状パラメータを更新するパラメータ更新部と
を備えたことを特徴とするマスク推定装置。
【請求項２】
前記要素分布は、回転対称ではない分布形状も表現することができることを特徴とする請求項１に記載のマスク推定装置。
【請求項３】
前記要素分布の確率密度関数は、ｚを独立変数、ａを単位ベクトル、・Ｈをエルミート転置として、ａとｚの内積ａＨｚの絶対値｜ａＨｚ｜の関数ではない確率密度関数も表現することができることを特徴とする請求項１又は２に記載のマスク推定装置。
【請求項４】
前記要素分布の確率密度関数は、ｚを独立変数、Ａを行列、・Ｈをエルミート転置として、二次形式ｚＨＡｚの関数であり、
前記形状パラメータは、前記行列Ａを定めるパラメータである
ことを特徴とする請求項１〜３のいずれか１項に記載のマスク推定装置。
【請求項５】
前記形状パラメータは、各要素分布の分布形状を表す楕円の各軸方向の分布広がりの情報を含む
ことを特徴とする請求項１〜４のいずれか１項に記載のマスク推定装置。
【請求項６】
前記要素分布は、複素ビンガム分布であり、
前記形状パラメータは、複素ビンガム分布のパラメータ行列である
ことを特徴とする請求項１〜５のいずれか１項に記載のマスク推定装置。
【請求項７】
前記要素分布は、複素角度ガウス分布であり、
前記形状パラメータは、複素角度ガウス分布のパラメータ行列である
ことを特徴とする請求項１〜５のいずれか１項に記載のマスク推定装置。
【請求項８】
マスク推定装置が実行するマスク推定方法であって、
前記マスク推定装置は、複数の音響信号が混在する状況下で、異なる位置で取得した複数の観測信号から、前記観測信号の時間周波数成分を時間周波数点毎にまとめることにより抽出された特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化する際の各要素分布の混合重みと、各要素分布の形状を制御可能なモデルパラメータである形状パラメータとを記憶する記憶部を備え、
前記特徴量ベクトルを抽出する特徴抽出ステップと、
前記特徴量ベクトルと、前記混合重みと、前記形状パラメータとを用いて、各時間周波数点に対して各要素分布が寄与する割合を要素分布毎かつ時間周波数点毎に示すマスクを推定して更新するマスク更新ステップと、
前記マスク更新ステップにより更新されたマスクから、前記混合重みを更新する混合重み更新ステップと、
前記特徴量ベクトルと、前記マスクとを用いて、前記形状パラメータを更新するパラメータ更新ステップと
を含んだことを特徴とするマスク推定方法。
【請求項９】
前記要素分布は、回転対称ではない分布形状も表現することができることを特徴とする請求項８に記載のマスク推定方法。
【請求項１０】
前記要素分布の確率密度関数は、ｚを独立変数、ａを単位ベクトル、・Ｈをエルミート転置として、ａとｚの内積ａＨｚの絶対値｜ａＨｚ｜の関数ではない確率密度関数も表現することができることを特徴とする請求項８又は９に記載のマスク推定方法。
【請求項１１】
前記要素分布の確率密度関数は、ｚを独立変数、Ａを行列、・Ｈをエルミート転置として、二次形式ｚＨＡｚの関数であり、
前記形状パラメータは、前記行列Ａを定めるパラメータである
ことを特徴とする請求項８〜１０のいずれか１項に記載のマスク推定方法。
【請求項１２】
前記形状パラメータは、各要素分布の分布形状を表す楕円の各軸方向の分布広がりの情報を含む
ことを特徴とする請求項８〜１１のいずれか１項に記載のマスク推定方法。
【請求項１３】
請求項１〜７のいずれか１つに記載のマスク推定装置としてコンピュータを機能させるためのマスク推定プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、マスク推定装置、マスク推定方法及びマスク推定プログラムに関する。
【背景技術】
【０００２】
複数のマイクロホンにより取得された観測信号から目的音を推定する際に用いられる、各時間周波数点における各音響信号の寄与度を表すマスクを推定する技術がある。マスクは、観測信号の雑音除去や音源分離等に用いられる。【０００３】
例えば、各音響信号の到来方向が互いに異なることに着目し、各時間周波数点において各音響信号音がどの方向から到来するかの情報を表す特徴量ベクトルをクラスタリングすることにより、マスクを推定する技術がある。
【先行技術文献】
【非特許文献】
【０００４】
【非特許文献１】M.  Souden,  S.  Araki,  K.  Kinoshita,  T.  Nakatani,  and  H.  Sawada,  “A  Multichannel  MMSE-Based  Framework  for  Speech  Source  Separation  and  Noise  Reduction,”  IEEE  TRANSACTIONS  ON  AUDIO,  SPEECH  AND  LANGUAGE  PROCESSING,  VOL.  21,  NO.  9,  SEPTEMBER  2013,  pp.  1913  -  1928.
【発明の概要】
【発明が解決しようとする課題】
【０００５】
しかしながら、上記技術では、特徴量ベクトルをモデル化する混合分布の位置パラメータ及び広がりパラメータを更新することにより、混合分布の要素分布の位置及び広がりを調整する。そのため、要素分布は、ある軸に関して回転対称な分布形状を表すに過ぎない。特徴量ベクトルの分布形状は、マイクロホン位置や音響空間の伝達特性等の様々な影響を受け、必ずしも回転対称になるとは限らず、例えば楕円状の分布形状となる場合がある。そのため、上記技術では、混合分布を特徴量ベクトルの分布に十分近づけることができず、マスク推定の精度が必ずしも高くないという問題がある。【０００６】
本願が開示する実施形態の一例は、上記に鑑みてなされたものであって、例えば、マスク推定の精度をより向上させるマスク推定装置、マスク推定方法及びマスク推定プログラムを提供することを目的とする。
【課題を解決するための手段】
【０００７】
本願の実施形態の一例において、例えば、マスク推定装置は、複数の音響信号が混在する状況下で、異なる位置で取得した複数の観測信号から、観測信号の時間周波数成分を時間周波数点毎にまとめた特徴量ベクトルを抽出する。マスク推定装置は、特徴量ベクトルと、特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化する際の各要素分布の混合重みと、各要素分布の形状を制御可能なモデルパラメータである形状パラメータとを用いて、各時間周波数点に対して各要素分布が寄与する割合を示すマスクを更新する。マスク推定装置は、更新されたマスクから、混合重みを更新する。マスク推定装置は、特徴量ベクトルと、マスクとから、形状パラメータを更新する。
【発明の効果】
【０００８】
本願が開示する実施形態の一例によれば、例えば、特徴量ベクトルの分布が回転対称でない場合でも高精度なマスク推定が可能になる。
【図面の簡単な説明】
【０００９】
【図１】図１は、実施形態１（実施形態２）に係るマスク推定装置の一例を示す図である。
【図２】図２は、実施形態１（実施形態２）に係るマスク推定処理の一例を示すフローチャートである。
【図３】図３は、実施形態３に係る目的音推定システムの一例を示す図である。
【図４】図４は、実施形態３に係る目的音推定処理の一例を示すフローチャートである。
【図５】図５は、プログラムが実行されることにより、実施形態に係るマスク推定装置及び目的音推定装置を含む目的音推定システムが実現されるコンピュータの一例を示す図である。
【図６】図６は、背景技術に係るマスク推定装置の一例を示す図である。
【発明を実施するための形態】
【００１０】
以下、本願が開示する技術に関する実施形態の一例について、図面を参照して説明する。なお、以下の実施形態により、本願の開示技術が限定されるものではない。また、以下の実施形態は、矛盾しない範囲で適宜組合せてもよい。
【００１１】
なお、以下の記載において、“Ａ”に対して“＾Ａ”と表記する場合は、「Ａの直上に＾が付された記号」と同等であるとする。また、“Ａ”がベクトルである場合には「ベクトルＡ」、“Ａ”が行列である場合には「行列Ａ」、“Ａ”が集合である場合には「集合Ａ」と表記する。
【００１２】
［背景技術について］
本願の実施形態の開示に先立ち、本願が開示する技術の基礎となる背景技術について説明する。背景技術は、文献１「M.  Souden,  S.  Araki,  K.  Kinoshita,  T.  Nakatani,  and  H.  Sawada,  “A  Multichannel  MMSE-Based  Framework  for  Speech  Source  Separation  and  Noise  Reduction,”  IEEE  TRANSACTIONS  ON  AUDIO,  SPEECH  AND  LANGUAGE  PROCESSING,  VOL.  21,  NO.  9,  SEPTEMBER  2013,  pp.  1913  -  1928.」に開示されている。
【００１３】
図６は、背景技術に係るマスク推定装置の一例を示す図である。背景技術に係るマスク推定装置１０Ａは、記憶部２０Ａが接続される。マスク推定装置１０Ａは、特徴抽出部１１Ａ、マスク更新部１２Ａ、混合重み更新部１３Ａ、位置パラメータ更新部１４Ａ、広がりパラメータ更新部１５Ａを有する。
【００１４】
記憶部２０Ａは、後述するように特徴抽出部１１Ａにより抽出された特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化したときの、各要素分布の混合重み、各要素分布の位置を表す位置パラメータと、各要素分布の広がりを表す広がりパラメータを記憶する。
【００１５】
特徴抽出部１１Ａは、Ｎ個（Ｎ＞１）の音響信号が混在する状況において、異なる位置で収録されたＭ個（Ｍ＞１）の観測信号を受け取り、時間周波数点毎にすべての観測信号の時間周波数成分をまとめたＭ次元縦ベクトルを構成し、Ｍ次元縦ベクトルを用いて所定のノルムを持つＭ次元の特徴量ベクトルを抽出する。
【００１６】
マスク更新部１２Ａは、特徴抽出部１１Ａにより抽出された特徴量ベクトルを特徴抽出部１１Ａから受け取り、記憶部２０Ａに保存されている混合重みと、位置パラメータと、広がりパラメータを記憶部２０Ａから受け取り、各時間周波数点に各要素分布がどの程度寄与しているかの割合を表すマスクを更新する。
【００１７】
混合重み更新部１３Ａは、マスク更新部１２Ａにより更新されたマスクを受け取り、混合重みを更新する。位置パラメータ更新部１４Ａは、特徴抽出部１１Ａにより抽出された特徴量ベクトルと、マスク更新部１２Ａにより更新されたマスクを受け取り、位置パラメータを更新する。広がりパラメータ更新部１５Ａは、特徴抽出部１１Ａにより抽出された特徴量ベクトルと、マスク更新部１２Ａにより更新されたマスクを受け取り、広がりパラメータを更新する。
【００１８】
上記文献１では、音響信号の数Ｎを既知と仮定すると共に、Ｎ−１個が目的音であり、１個が背景雑音であると仮定する。Ｎ＝２の場合は、例えば、背景雑音の存在下で一人の音声をＭ個のマイクロホンで収録する場合に対応する。また、Ｎ＞２の場合は、例えば、背景雑音の存在下でＮ−１人による会話をＭ個のマイクロホンで収録する場合に対応する。
【００１９】
いま、Ｍ個の観測信号のそれぞれに短時間フーリエ変換等の時間周波数分析を適用して得られる各観測信号の時間周波数成分を、時間周波数点毎にまとめたＭ次元縦ベクトル（以下、観測信号ベクトルと呼ぶ）をｙ(t,f)で表す。ｔ（ｔ＝１，・・・Ｔ。ｔは正整数）は時刻、ｆ（ｆ＝１，・・・Ｆ。ｆは正整数）は周波数ビン番号である。
【００２０】
上記文献１では、目的音はスパース性を有し、各時間周波数点において高々一つの目的音のみが存在すると仮定され、また、背景雑音はすべての時間周波数点において存在すると仮定されている。このとき、観測信号ベクトルｙ(t,f)は、各時間周波数点において次のＮ個の状態のうちのいずれかを取る。
【００２１】
状態Ｓn：目的音のうちｎ番目の目的音のみが存在する状態（ｎ＝１，・・・，Ｎ−１。ｎは正整数）
状態ＳN：目的音が一つも存在しない状態
【００２２】
このため、観測信号ベクトルｙ(t,f)は、下記（１−１）式又は（１−２）式のいずれかでモデル化できるとされている。
【００２３】
【数１】
2017141542.tif 000003
【００２４】
上記（１−１）式は、当該時間周波数点においてｎ番目の目的音のみが存在する場合を表し、上記（１−２）式は、当該時間周波数点において目的音が一つも存在しない場合を表す。上記（１−１）式及び（１−２）式におけるベクトルｓn(t,f)はｎ番目の目的音に対応する成分を表し、ベクトルｖ(t,f)は背景雑音に対応する成分を表す。
【００２５】
マスク推定装置１０Ａの目的は、観測信号ベクトルｙ(t,f)が与えられたときに、各時間周波数点(t,f)に各状態Ｓn（ｎ＝１，・・・，Ｎ）がどの程度寄与しているかの割合を表すマスクγn(t,f)を推定することである。マスクγn(t,f)は、γ1(t,f)＋・・・＋γN(t,f)＝１，０≦γn(t,f)≦１を満たす。マスクγn(t,f)は、例えば、当該時間周波数点での状態がＳnの場合にγn(t,f)＝１、そうでない場合にγn(t,f)＝０と定義することができる。あるいは、より一般的には、マスクγn(t,f)を当該時間周波数点における状態Ｓnの事後確率と定義することもできる。以下の説明では、マスクγn(t,f)は、当該時間周波数点における状態Ｓnの事後確率と定義する。
【００２６】
マスクγn(t,f)が推定できると、各目的音ｓn(t,f)の推定等、さまざまな応用に用いることができる。例えば、マスクγn(t,f)を用いて当該目的音が鳴っている時間周波数成分を集めてくることにより、当該目的音の性質を表す統計量を推定することができる。この統計量を用いて設計したフィルタにより各目的音ｓn(t,f)を推定することができる。
【００２７】
上記文献１では、各音響信号の到来方向は互いに異なることに着目し、到来方向をクラスタリングすることにより、マスクγn(t,f)を推定する。特徴抽出部１１Ａは、各時間周波数点で音がどの方向から到来するかの情報を表す特徴量ベクトルを抽出する。観測信号ベクトルｙ(t,f)の大きさ（ノルム）は音響信号自体に依存するが、観測信号ベクトルｙ(t,f)の方向はその音源位置のみで決まると仮定する。この仮定に基づき、音源位置を表す特徴量ベクトルとして、一定のノルムを持つＭ次元の特徴量ベクトルｚ(t,f)を観測信号ベクトルｙ(t,f)から抽出する。特徴量ベクトルｚ(t,f)は、例えば下記（２）式により抽出できる。
【００２８】
【数２】
2017141542.tif 000004
【００２９】
ここで、||・||は２ノルムを表す。当該時間周波数点がどの状態Ｓnを取るかによって、音の到来方向が異なるため、特徴量ベクトルｚ(t,f)は、状態Ｓn毎に異なる確率分布を持つ。
【００３０】
そこで、特徴量ベクトルｚ(t,f)に対してクラスタリング技術を適用することにより、各時間周波数点に各状態がどの程度寄与しているかの割合を表すマスクを推定することができる。これが、上記文献１における基本的な考え方である。
【００３１】
上記文献１では、周波数毎に、特徴量ベクトルz(t,f)の確率分布を、下記（３）式の混合分布でモデル化する。
【００３２】
【数３】
2017141542.tif 000005
【００３３】
上記（３）式の混合分布が、特徴量ベクトルｚ(t,f)の分布を近似するようにモデルパラメータ（集合）Θを推定することにより、クラスタリングを行う。上記（３）式の混合分布の要素分布は、下記（４）式で表される。
【００３４】
【数４】
2017141542.tif 000006
【００３５】
上記（４）式の要素分布は、状態が既知であるという条件下での特徴量ベクトルｚ(t,f)の条件付き確率分布をモデル化し、α1(f)＋・・・＋αK(f)＝１，０≦αk(f)≦１を満たすαk(f)は、ｋ番目の要素分布の混合重みである。また、||ａk(f)||＝１を満たす位置パラメータ（ベクトル）ａk(f)は、ｋ番目の要素分布の位置を表す。また、広がりパラメータσ2k(f)は、ｋ番目の要素分布の広がりを表す。なお、Ｋは混合数であり、モデルパラメータ（集合）Θは上記（３）式の混合分布のモデルパラメータの集合｛αk(f),ａk(f),σ2k(f)｝を表し、・Hは行列のエルミート転置を表す。上記（３）式における要素分布ｐW(ｚ(t,f);ａk(f),σ2k(f))を山に例えると、上記（３）式の混合分布のモデルパラメータである混合重みαk(f)は山の高さ、位置パラメータａk(f)は山の位置、広がりパラメータσ2k(f)は山の裾野の広さに該当する。上記文献１では、混合数Ｋは、音響信号の数Ｎ（既知と仮定）と等しいとする。
【００３６】
モデルパラメータΘが求まれば、ベイズの定理より、下記（５）式により、特徴量ベクトルｚ(t,f)が観測された条件下で当該時間周波数点がｋ番目の要素分布に対応する事後確率（すなわちマスク）が求まる。
【００３７】
【数５】
2017141542.tif 000007
【００３８】
ここで、モデルパラメータΘをいかに推定するか問題となる。上記文献１では、モデルパラメータΘを用いて、上記（５）式によりマスクγk(t,f)を更新するステップと、マスクγk(t,f)を用いて、モデルパラメータΘを更新するステップの２つのステップを交互に反復することにより、モデルパラメータΘ及びマスクγk(t,f)を推定する。モデルパラメータΘを更新するステップでは、マスクγk(t,f)を用いて、下記（６−１）式〜（６−３）式により、モデルパラメータΘを更新する。
【００３９】
【数６】
2017141542.tif 000008
【００４０】
ここで、行列Ｒk(f)は、下記（７）式により計算される。
【００４１】
【数７】
2017141542.tif 000009
【００４２】
上記の反復する処理は、期待値最大化法による対数尤度の最適化として、理論的に導出することができる。なお、上記（５）式はマスク更新部１２Ａがマスクγk(t,f)を更新する処理に相当し、上記（６−１）式は混合重み更新部１３Ａが混合重みαk(f)を更新する処理に相当し、上記（６−２）式は位置パラメータ更新部１４Ａが位置パラメータａk(f)を更新する処理に相当し、上記（６−３）式は広がりパラメータ更新部１５Ａが広がりパラメータσ2k(f)を更新する処理に相当する。
【００４３】
ここで、上記文献１の方法では、上記（３）式の混合分布に基づき上記（５）式によりマスクを推定するため、上記（３）式の混合分布が特徴量ベクトルｚ(t,f)の分布をどの程度正確に近似できるかが、マスクの推定精度に大きく影響する。文献１の方法では、上記（３）式の混合分布を特徴量ベクトルｚ(t,f)の分布に近づけるために、上記（３）式の混合分布の位置パラメータａk(f)と広がりパラメータσ2k(f)を更新することにより、上記（４）式の要素分布の位置と広がりを調整する。そのため、上記（４）式の要素分布は、ある軸に関して回転対称な分布形状しか表せないという制限がある。
【００４４】
特徴量ベクトルｚ(t,f)の分布形状は、マイクロホン配置や部屋の音響伝達特性等のさまざまな影響を受け、必ずしも回転対称になるとは限らない。そのため、上記文献１の方法では、必ずしも上記（３）式の混合分布が特徴量ベクトルｚ(t,f)の分布を十分近似することができず、上記（５）式によるマスク推定の精度が高くないという問題がある。
【００４５】
そこで、開示の実施形態の第１の態様のマスク推定装置は、Ｎ個（Ｎ＞１）の音響信号が混在する状況で収録したＭ個の観測信号（Ｍ＞１）に基づくＭ次元の特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化したときの、各要素分布の混合重みと、各要素分布の形状を制御可能なモデルパラメータである形状パラメータをもとにマスクを更新する。なお、形状とは、図形の属性のうち、位置及び広がりを示す属性を除いた属性とする。形状には、例えば、楕円の長軸及び短軸に相当する情報等がある。
【００４６】
第１の態様のマスク推定装置によれば、形状パラメータの更新により要素分布の形状を調整することで、上記文献１の方法と比べ、混合分布を用いて特徴量ベクトルの分布をより正確に近似することができ、マスクをより正確に推定することができる。
【００４７】
また、開示の実施形態の第２の態様のマスク推定装置は、第１の態様のマスク推定装置において、特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化したときの各要素分布を複素ビンガム分布とし、複素ビンガム分布のパラメータ行列を形状パラメータとする。このとき、特徴量ベクトルｚ(t,f)の確率分布ｐ(ｚ(t,f);Θ)は、複素ビンガム分布を要素分布とする、下記（８）式に示す混合分布（以下、混合複素ビンガム分布と呼ぶ）によりモデル化される。
【００４８】
【数８】
2017141542.tif 000010
【００４９】
なお、上記（８）式におけるｐB(ｚ;Ｂ)は、下記（９）式により定義される複素ビンガム分布である。なお、exp・は、指数関数を表す。
【００５０】
【数９】
2017141542.tif 000011
【００５１】
上記（９）式における行列Ｂは、複素ビンガム分布のパラメータ行列であり、ｃ(Ｂ)は規格化定数である。パラメータΘは、上記（８）式の混合複素ビンガム分布のモデルパラメータの集合｛αk(f),Ｂk(f)｝を表す。
【００５２】
なお、上記（８）式の混合分布のモデルパラメータである混合重みαk(f)は要素分布ｋの高さを表し、行列Ｂk(f)は要素分布ｋの位置、広がりに加えて分布形状（分布形状を表す楕円の各軸方向の分布広がり等）も制御可能な形状パラメータである。行列Ｂk(f)の第一固有ベクトルは要素分布ｋの位置を表し、行列Ｂk(f)の第一固有値と第二固有値の差の絶対値は要素分布ｋの広がりの小ささを表し、行列Ｂk(f)の第一固有値と第ｍ固有値（３≦ｍ≦Ｍ）の差の絶対値は、要素分布ｋの分布形状を表す楕円の第（ｍ−１）軸方向の分布広がりの小ささを表す。
【００５３】
いま、上記（９）式の複素ビンガム分布は、パラメータ行列Ｂがσ2とａの２つのパラメータのみを用いてＢ＝(１/σ2)・ａａHと表される特殊な場合（分布形状を表す楕円の第一軸方向の分布広がり、第二軸方向の分布広がり、・・・・・・、第（Ｍ−１）軸方向の分布広がりがすべて等しい場合）を考えると、下記（１０）式のようになる。
【００５４】
【数１０】
2017141542.tif 000012
【００５５】
上記（１０）式は、上記（４）式と同じ形をしており、ベクトルａに平行な軸に関して回転対称となる。これに対し、上記（９）式の複素ビンガム分布は、パラメータ行列Ｂに特段の制約は設けられておらず、パラメータ行列Ｂを更新することにより、上記（９）式の複素ビンガム分布の分布形状を調整して回転対称でない分布形状も表現することができる。
【００５６】
これにより、上記（９）式の複素ビンガム分布を用いれば、特徴量ベクトルｚ(t,f)の分布形状が回転対称でない場合でも、上記（８）式の混合分布で特徴量ベクトルｚ(t,f)の分布を十分近似することができる。よって、上記文献１に記載の方法よりも高精度にマスクを推定することが可能になる。
【００５７】
上記（８）式によるモデル化の下、期待値最大化法に基づき、モデルパラメータΘを用いて、マスクγk(t,f)を更新するステップと、マスクγk(t,f)を用いて、モデルパラメータΘを更新するステップの２つのステップを交互に反復することにより、マスクγk(t,f)とモデルパラメータΘを推定するアルゴリズム（上述）を導出することができる。
【００５８】
また、開示の実施形態の第３の態様のマスク推定装置は、第１の態様のマスク推定装置において、特徴量ベクトルの確率分布を複数の要素分布からなる混合分布でモデル化したときの各要素分布をComplex  Angular  Central  Gaussian  (cACG)  distribution（以下、複素角度ガウス分布と呼ぶ）とし、複素角度ガウス分布のパラメータ行列を形状パラメータとする。
【００５９】
このとき、特徴量ベクトルｚ(t,f)の確率分布ｐ(ｚ(t,f);Θ)は、複素角度ガウス分布を要素分布とする、下記（１１）式に示す混合分布（以下、混合複素角度ガウス分布と呼ぶ）によりモデル化される。
【００６０】
【数１１】
2017141542.tif 000013
【００６１】
ここで、上記（１１）式におけるｐA(ｚ;Σ)は、パラメータ行列を行列Σとする、下記（１２）式により定義される複素角度ガウス分布である。
【００６２】
【数１２】
2017141542.tif 000014
【００６３】
行列Σk(f)は、要素分布ｋの位置、広がりに加えて分布形状（分布形状を表す楕円の各軸方向の分布広がり等）も制御可能な形状パラメータである。行列Σk(f)の第一固有ベクトルは要素分布ｋの位置を表し、行列Σk(f)の第一固有値を第二固有値で割った比は要素分布ｋの広がりの小ささを表し、行列Σk(f)の第一固有値を第ｍ固有値（３≦ｍ≦Ｍ）で割った比は、要素分布ｋの分布形状を表す楕円の第（ｍ−１）軸方向の分布広がりの小ささを表す。モデルパラメータΘは、上記（１１）式における混合複素角度ガウス分布のモデルパラメータの集合｛αk(f);Σk(f)｝を表す。
【００６４】
上記の第２の態様のマスク推定装置における上記（９）式の複素ビンガム分布と同様に、上記（１２）式の複素角度ガウス分布においてもパラメータ行列Σに特段の制約は設けられておらず、回転対称でない分布形状を表現することができる。そのため、上記（１２）式の複素角度ガウス分布を用いれば、特徴量ベクトルｚ(t,f)の分布形状が回転対称でない場合でも分布形状をより正確にモデル化でき、上記文献１の方法よりもより高精度にクラスタリングを行うことが可能になる。
【００６５】
［実施形態１］
実施形態１では、要素分布として複素ビンガム分布を用いる。実施形態１では、上記文献１と同様に、音響信号の数Ｎを既知と仮定し、Ｎ−１個の目的音と、１個の背景雑音が混ざったものをＭ個のマイクロホンで収録するものとする。
【００６６】
（実施形態１に係るマスク推定装置）
図１は、実施形態１に係るマスク推定装置の一例を示す図である。実施形態１に係るマスク推定装置１０は、記憶部２０が接続されている。マスク推定装置１０は、特徴抽出部１１、マスク更新部１２、混合重み更新部１３、パラメータ更新部１４を有する。τを時間を表す番号（サンプル番号）として、マイクロホンｍで収録された時間領域の観測信号をｙm(τ)と書くことにする。ｙm(τ)は、ｎ番目の目的音源信号ｓmn(τ)(ｎ＝１〜Ｎ−１。ｎは正整数)と、背景雑音ｖm(τ)の和として、下記（１３）式によりモデル化できる。
【００６７】
【数１３】
2017141542.tif 000015
【００６８】
特徴抽出部１１は、複数のマイクロホンで収録された観測信号を受け取り、各観測信号ｙm(τ)毎に時間周波数分析を適用して各観測信号の時間周波数成分ｙm(t,f)を求める（ｍはマイクロホンの番号を表し、１〜Ｍの整数とする）。時間周波数分析としては、短時間フーリエ変換や短時間離散コサイン変換などさまざまな方法を用いることができる。
【００６９】
次に、特徴抽出部１１は、時間周波数点毎にすべての観測信号の時間周波数成分をまとめた、下記（１４）式により表されるＭ次元縦ベクトル（観測信号べクトルと呼ぶ）ｙ(t,f)を構成する。
【００７０】
【数１４】
2017141542.tif 000016
【００７１】
次に、特徴抽出部１１は、観測信号ベクトルｙ(t,f)を用いて所定のノルムを持つＭ次元の特徴量ベクトルｚ(t,f)を抽出する。特徴量ベクトルｚ(t,f)としては、文献２「H.  Sawada,  S.  Araki,  and  S.  Makino,  “Underdetermined  convolutive  blind  sourceseparation  via  frequency  bin-wise  clustering  and  permutation  alignment,”  IEEE  Trans-actions  on  Audio,  Speech  and  Language  Processing(ASLP),  vol.  19,  no.  3,  pp.  516  -  527,  Mar.  2011.」、文献３「D.H.  Tran  Vu  and  R.  Haeb-Umbach,  “Blind  speech  separation  em-ploying  directional  statistics  in  an  expectation  maximization  framework,”  Proceedings  of  the  International  Conference  on  Acoustics,  Speech  and  Signal  Processing  (ICASSP),  pp.  241  -  244,  Mar.  2010.」、文献４「S.  Araki,  H.  Sawada,  R.  Mukai,  and  S.  Makino,  “Underdeter-mined  blind  sparse  source  separation  for  arbitrarily  arranged  multiple  sensors,”  Signal  Processing,  vol.  87,  no.  8,  pp.  1833  -  1847,  Aug.  2007.」に記載の特徴量ベクトル等、さまざまなものを用いることができる。例えば、上記文献２に倣うと、特徴量ベクトルｚ(t,f)は、下記（１５）式により抽出することができる。
【００７２】
【数１５】
2017141542.tif 000017
【００７３】
次に、マスク更新部１２は、特徴量ベクトルｚ(t,f)を受け取り、各時間周波数点(t,f)にｋ番目の要素分布がどの程度寄与しているかの割合を表すマスクγk(t,f)を計算し出力する。
【００７４】
以下、マスク更新部１２における処理の流れを説明する。まず、マスク更新部１２は、記憶部２０にあらかじめ設定し保持しておいたモデルパラメータΘ:=｛αk(f),Ｂk(f)｝の初期値を記憶部２０から読み出す。モデルパラメータΘの初期値の設定は、乱数を用いる等、さまざまな方法により行うことができる。
【００７５】
次に、マスク更新部１２は、特徴量ベクトルｚ(t,f)と、モデルパラメータΘの現在の推定値を受け取って、ｋ番目の要素分布に対応するマスクを、例えば、事後確率γk(t,f)として、下記（１６）式により更新する。
【００７６】
【数１６】
2017141542.tif 000018
【００７７】
ここで、上記（１６）式における行列Ｂk(f)の固有値をｅk1(f),・・・,  ｅkM(f)とすると、規格化定数ｃ（Ｂk(f)）は、下記（１７）式のようになる。また、混合数Ｋは、Ｋ＝Ｎとする。
【００７８】
【数１７】
2017141542.tif 000019
【００７９】
次に、混合重み更新部１３は、事後確率γk(t,f)を受け取って、下記（１８）式により、混合重みを更新する。
【００８０】
【数１８】
2017141542.tif 000020
【００８１】
次に、パラメータ更新部１４は、特徴量ベクトルｚ(t,f)と、事後確率γk(t,f)を受け取って、正定値エルミート行列Ｒk(f)を、下記（１９）式により計算する。
【００８２】
【数１９】
2017141542.tif 000021
【００８３】
次に、パラメータ更新部１４は、下記（２０）式のように、行列Ｒk(f)の固有値分解を行う。
【００８４】
【数２０】
2017141542.tif 000022
【００８５】
上記（２０）式の固有値分解によれば、行列Ｒk(f)の固有ベクトルにより構成されるユニタリ行列Ｕk(f)と、行列Ｒk(f)の固有値により構成される対角行列Ｄk(f)が求まる。ここで、行列Ｄk(f)のｍ番目の対角成分（固有値）ｄkm(f)は、ｄk1(f)≦・・・≦ｄkM(f)のように、昇順に並べてあるものとする。
【００８６】
次に、パラメータ更新部１４は、下記（２１）式により、パラメータ行列Ｂk(f)を更新する。
【００８７】
【数２１】
2017141542.tif 000023
【００８８】
ここで、上記（２１）式における行列Ｅk(f)は、ｍ番目の対角成分がｅkm(f)である対角行列である。ｅkm(f)は、下記（２２）式で与えられる。
【００８９】
【数２２】
2017141542.tif 000024
【００９０】
マスク更新部１２、混合重み更新部１３、パラメータ更新部１４での処理は、終了条件が満たされるまで反復して行われる。終了条件としては、例えば、「所定の回数（例えば２０回）反復が行われた」、「更新前後での（後述する）対数尤度関数の増加量が所定の閾値以下であった」等、さまざまな条件を用いることができる。
【００９１】
なお、記憶部２０は、混合重み更新部１３により更新された混合重み、パラメータ更新部１４により更新された形状パラメータを保存し、次回、マスク更新部１２、パラメータ更新部１４での処理の際に、保存した混合重み、形状パラメータを提供する。
【００９２】
（マスク推定装置１０における処理の理論的背景）
以下、マスク推定装置１０における処理の理論的背景について説明する。マスク推定装置１０における処理は、下記（２３）式の対数尤度関数Ｌ(Θ)をモデルパラメータΘに関して最大化する最適化問題を解くことにより導かれるものである。
【００９３】
【数２３】
2017141542.tif 000025
【００９４】
上記（２３）式の対数尤度関数Ｌ(Θ)は、期待値最大化法に基づいて、以下に述べるＥステップとＭステップを交互に反復することで最適化することができる。
【００９５】
||ｚ||＝１より、任意の実数ｈに対してｐB(ｚ;Ｂ+ｈＩ)＝ｐB(ｚ;Ｂ)（ただし、Ｉは単位行列）が成り立つことから、以下では、一般性を失わず、行列Ｂk(f)の最大固有値ｅkM(f)を０とする。これは、複素ビンガム分布は、各固有値の値そのものではなく、各固有値の間の差分により分布形状が定まるものであることを意味する。このことから、複素ビンガム分布は、最大固有値を指定値に固定する条件を追加したとしても、固有値間の値の差が保たれれば不変であるため、最大固有値ｅkM(f)＝０と固定することができる。
【００９６】
Ｅステップでは、モデルパラメータΘの現在の推定値を用いて、ｋ番目の要素分布の事後確率γk(t,f)を、下記（２４−２）式により更新する。ただし、規格化定数ｃ(・)は、上記（１７）式により定義される。
【００９７】
【数２４】
2017141542.tif 000026
【００９８】
Ｍステップでは、Ｅステップで更新した事後確率γk(t,f)を用いて、下記（２５−１）式及び（２５−２）式のように定義されるＱ関数を最大化することにより、モデルパラメータΘを更新する。なお、行列Ｒk(f)は、上記（１９）式により定義され、trは行列のトレースを表す。
【００９９】
【数２５】
2017141542.tif 000027
【０１００】
すなわち、ラグランジュの未定乗数法により、制約条件α1(f)＋・・・＋αK(f)＝１の下で、関数Ｑ(Θ)を最大化するαk(f)が、上記（１８）式のように得られる。以下では、要素分布である複素ビンガム分布の形状パラメータＢk(f)の更新式を導出する。行列Ｒk(f)の固有値分解を下記（２６−１）式のように定義し、行列Ｂk(f)の固有値分解を下記（２６−２）式のように定義する。
【０１０１】
【数２６】
2017141542.tif 000028
【０１０２】
ここで、上記（２６−１）式における行列Ｕk(f)は行列Ｒk(f)の固有ベクトルからなるユニタリ行列であり、上記（２６−２）式における行列Ｖk(f)は行列Ｂk(f)の固有ベクトルからなるユニタリ行列である。また、上記（２６−１）式における行列Ｄk(f)は行列Ｒk(f)の固有値からなる対角行列であり、上記（２６−２）式における行列Ｅk(f)は行列Ｂk(f)の固有値からなる対角行列である。
【０１０３】
そして、行列Ｒk(f)の固有値をｄk1(f)≦・・・≦ｄkM(f)、行列Ｂ(f)の固有値をｅk1(f)≦・・・≦ｅkM(f)＝０とするとき、行列Ｄk(f)及び行列Ｅk(f)は、それぞれ下記（２７−１）式及び（２７−２）式のように表される。なお、行列Ｂ(f)の最大固有値ｅkM(f)を０とおいた。
【０１０４】
なお、ある音源に着目して空間相関行列Ｒk(f)と複素ビンガム分布のパラメータ行列Ｂ(f)を固有値分解してマイクロホンの個数分の固有値を求める際、ｍ番目に大きい第ｍ固有値（ｍ＝１〜Ｍ）は、次の意味を持つ。すなわち、第１固有ベクトルは、着目音源の分布のピークの位置を表す。また、第ｍ固有ベクトル（ｍ＝２〜Ｍ）は、着目音源の分布を表す楕円の第（ｍ−１）軸の方向（行列Ｒk(f)の第ｍ固有値が軸の長さに相当）を表す。
【０１０５】
【数２７】
2017141542.tif 000029
【０１０６】
ここで、文献５「K.V.  Mardia  and  P.E.  Jupp,  “Directional  Statistics,”  John  Wiley  &  Sons,  WestSussex,  2000.」に開示されている結果より、行列Ｂk(f)が関数Ｑ(Θ)を最大化するとき、Ｖk(f)＝Ｕk(f)となる。したがって、上記（２５−２）式において、下記（２８）式のようになることを用いて、ｅkm(f)（ｍ＝１,・・・,Ｍ−１）に関する関数Ｑ(Θ)の偏微分を０と置くと、下記（２９）式を得る。
【０１０７】
【数２８】
2017141542.tif 000030
【０１０８】
【数２９】
2017141542.tif 000031
【０１０９】
上記（２９）式をｅkm(f)について解くために、行列Ｒ(f)の固有値ｄk1(f),・・・,ｄkM(f)はすべて異なると仮定する。この仮定は、現実のデータに対してもほとんど常に成り立つ。このとき、上記（１７）式により、規格化定数ｃ（Ｂk(f)）が与えられる。
【０１１０】
いま、各音響信号が特定の方向から到来するため複素ビンガム分布は十分に集中していると仮定し、ｅkm(f)＝κk(f)ｅkm0(f)とおいてκk(f)→∞の極限を取ると、下記（３０−１）式〜（３０−３）式のように、規格化定数ｃ（Ｂk(f)）が近似できる。
【０１１１】
【数３０】
2017141542.tif 000032
【０１１２】
すなわち、複素ビンガム分布では、「音源の位置が大きく変化せずに音を発生するなら、音源を表現する特徴量ベクトルの分布形状が、比較的急峻なピークを持つこととなる（分布が十分に集中することになる）」ことに着目すると、上記（３０−２）式でのκk(f)は非常に大きな値を有するとの仮定を導入することができ、この結果、規格化定数ｃ（Ｂk(f)）を（３０−３）式の形式で表現できるようになる。したがって、上記（２９）式をｅkm(f)について解くことにより、上記（２２）式が得られる。
【０１１３】
（実施形態１に係るマスク推定処理）
図２は、実施形態１に係るマスク推定処理の一例を示すフローチャートである。図２に示すように、先ず、特徴抽出部１１は、複数のマイクロホンで収録された観測信号から、Ｍ次元の特徴量ベクトルｚ(t,f)を抽出する（ステップＳ１１）。
【０１１４】
次に、マスク更新部１２は、特徴量ベクトルｚ(t,f)、混合重み、及び形状パラメータからマスクγk(t,f)を計算し更新する（ステップＳ１２）。次に、混合重み更新部１３は、事後確率γk(t,f)から混合重みを更新する（ステップＳ１３）。次に、パラメータ更新部１４は、特徴量ベクトルｚ(t,f)と、事後確率γk(t,f)から、パラメータ行列を更新する（ステップＳ１４）。
【０１１５】
次に、マスク更新部１２は、終了条件充足か否かを判定する（ステップＳ１５）。マスク更新部１２は、終了条件充足の場合（ステップＳ１５：Ｙｅｓ）、マスク推定処理を終了し、終了条件充足でない場合（ステップＳ１５：Ｎｏ）、ステップＳ１２へ処理を移す。
【０１１６】
［実施形態２］
実施形態２では、要素分布として複素角度ガウス分布を用いる。実施形態１では、上記文献１と同様に、音響信号の数Ｎを既知と仮定し、Ｎ−１個の目的音と、１個の背景雑音が混ざったものをＭ個のマイクロホンで収録するものとする。
【０１１７】
図１を参照して、実施形態２に係るマスク推定装置１０−２の一例を説明する。マスク推定装置１０−２は、記憶部２０−２が接続されている。マスク推定装置１０−２は、特徴抽出部１１、マスク更新部１２−２、混合重み更新部１３−２、パラメータ更新部１４−２を有する。
【０１１８】
特徴抽出部１１は、実施形態１と同様である。マスク更新部１２−２は、特徴量ベクトルｚ(t,f)を受け取り、各時間周波数点(t,f)に、ｋ番目の要素分布がどの程度寄与しているかの割合を表すマスクγk(t,f)を計算し出力する。
【０１１９】
以下、マスク更新部１２−２における処理の流れを説明する。まず、マスク更新部１２−２は、記憶部２０−２にあらかじめ設定し保持しておいたモデルパラメータΘ:=｛αk(f),Σk(f)｝の初期値を記憶部２０−２から読み出す。行列Σk(f)は、複素角度ガウス分布のパラメータ行列である。
【０１２０】
次に、マスク更新部１２−２は、特徴量ベクトルｚ(t,f)とモデルパラメータΘの現在の推定値を受け取って、ｋ番目の要素分布に対応するマスクを、例えば、事後確率γk(t,f)として、下記（３１）式により更新する。なお、混合数ＫはＫ＝Ｎとする。
【０１２１】
【数３１】
2017141542.tif 000033
【０１２２】
次に、混合重み更新部１３−２は、事後確率γk(t,f)を受け取って、下記（３２）式により、混合重みを更新する。
【０１２３】
【数３２】
2017141542.tif 000034
【０１２４】
次に、パラメータ更新部１４−２は、特徴量ベクトルｚ(t,f)と、パラメータ行列Σk(t,f)と、事後確率γk(t,f)を受け取って、パラメータ行列Σk(t,f)を、下記（３３）式により更新する。
【０１２５】
【数３３】
2017141542.tif 000035
【０１２６】
マスク更新部１２−２、混合重み更新部１３−２、パラメータ更新部１４−２での処理は、実施形態１と同様に、終了条件が満たされるまで反復して行われる。
【０１２７】
なお、記憶部２０−２は、混合重み更新部１３−２により更新された混合重み、パラメータ更新部１４−２により更新された形状パラメータを保持し、次回、マスク更新部１２−２、パラメータ更新部１４−２での処理の際に、保存したマスク、形状パラメータを提供する。
【０１２８】
（マスク推定装置１０−２における処理の理論的背景）
以下、マスク推定装置１０−２における処理の理論的背景について説明する。マスク推定装置１０−２における処理は、下記（３４）式の対数尤度関数Ｌ(Θ)をモデルパラメータΘに関して最大化する最適化問題を解くことにより導かれるものである。
【０１２９】
【数３４】
2017141542.tif 000036
【０１３０】
上記（３４）式の対数尤度関数Ｌ(Θ)は、期待値最大化法に基づいて、以下に述べるＥステップとＭステップを交互に反復することで最適化することができる。
【０１３１】
Ｅステップでは、モデルパラメータΘの現在の推定値を用いて、ｋ番目の要素分布の事後確率γk(t,f)を、下記（３５−２）式により更新する。
【０１３２】
【数３５】
2017141542.tif 000037
【０１３３】
Ｍステップでは、Ｅステップで更新した事後確率γk(t,f)を用いて、下記（３６−１）式及び（３６−２）式のように定義されるＱ関数を最大化することにより、モデルパラメータΘを更新する。なお、ln・は、対数関数を表す。
【０１３４】
【数３６】
2017141542.tif 000038
【０１３５】
すなわち、ラグランジュの未定乗数法により、制約条件α1(f)＋・・・＋αK(f)＝１の下で、関数Ｑ(Θ)を最大化するαk(f)が、上記（３２）式のように得られる。また、関数Ｑ(Θ)のΣk(f)に関する偏微分は、下記（３７）式のようになる。
【０１３６】
【数３７】
2017141542.tif 000039
【０１３７】
そして、上記（３７）式で表される関数Ｑ(Θ)のΣk(f)に関する偏微分を０とおくと、上記（３３）式が得られる。
【０１３８】
なお、実施形態２に係るマスク推定処理は、実施形態１に係るマスク推定処理と同様に、図２に示すマスク推定処理の一例を示すフローチャートの処理順序に沿って実行される。
【０１３９】
［実施形態１及び実施形態２の変形例］
実施形態１及び実施形態２では、Ｎ個の音響信号のうちＮ−１個が目的音であり、１個が背景雑音であるとする。しかし、これに限らず、実施形態１及び実施形態２の変形例では、Ｎ個の音響信号がすべて目的音である場合にも、同一の処理により各音響信号（目的音）に対応するマスクを推定することができる。
【０１４０】
また、実施形態１及び実施形態２では、パラメータΘは確定的であるとみなし、最尤法によりパラメータΘを推定する。しかし、これに限らず、実施形態１及び実施形態２の変形例、例えば、パラメータΘを確率変数とみなして事前分布を与え、事後確率最大化によりパラメータΘを推定してもよい。
【０１４１】
［実施形態３］
実施形態３は、実施形態１のマスク推定装置１０又は実施形態２のマスク推定装置１０−２のいずれかを用いて目的音を推定する目的音推定システムを開示する。図３は、実施形態３に係る目的音推定システムの一例を示す図である。
【０１４２】
実施形態３に係る目的音推定システム１００は、実施形態１のマスク推定装置１０又は実施形態２のマスク推定装置１０−２のいずれか１つと、目的音推定装置３０とを有する。目的音推定装置３０は、観測信号ベクトルｙ(t,f)と、マスク推定装置１０又はマスク推定装置１０−２からのマスクγk(t,f)を入力とし、ｎ番目の目的音の推定値＾ｓn(t,f)（ｎ＝１〜Ｎ−１）を出力する。
【０１４３】
目的音推定装置３０は、行列計算部３１、ウィーナーフィルタ構成部３２、目的音推定部３３を有する。
【０１４４】
実施形態３では、周波数ｆごとにマスクを求めるため、同一番号ｎの音響信号に対応するクラスタが、異なる周波数ｆでは異なるクラスタ番号を持つということが起こる場合がある。この問題を、パーミュテーション問題という。このパーミュテーション問題のため、実施形態１又は実施形態２において求めたマスクγk(t,f)を用いて各目的音を推定するために、同一番号ｎの音響信号に対応するクラスタが、周波数ｆによらず同じクラスタ番号を持つように対応を取る。これをパーミュテーション解決という。パーミュテーション解決は、上記文献２「H.  Sawada,  S.  Araki,  and  S.  Makino,  “Underdetermined  convolutive  blind  source  separation  via  frequency  bin-wise  clustering  and  permutation  alignment,”  IEEE  Transactions  on  Audio,  Speech  and  Language  Processing(ASLP),  vol.  19,  no.  3,  pp.  516  -  527,  Mar.  2011.」に記載の方法等、さまざまな方法により行うことができる。
【０１４５】
いま、ｎ番目の音響信号に対応するパーミュテーション解決後のマスクを、改めてγn(t,f)と書くことにする。また、Ｎ個のマスクγn(t,f)のうち、いずれが背景雑音に対応するかが分かっていると仮定する。例えば、マスクを用いてマスキングにより作成した分離音を聞き比べることにより、人為的に背景雑音に対応するクラスタを同定することができる。
【０１４６】
以下では、一般性を失わずに、背景雑音に対応するクラスタがＮ番目のクラスタであると仮定する。行列計算部３１は、各目的音ｎ（ｎ＝１，・・・，Ｎ−１）に対して、下記（３８）式により雑音下目的音共分散行列Φn+v(f)を計算する。
【０１４７】
【数３８】
2017141542.tif 000040
【０１４８】
次に、行列計算部３１は、下記（３９）式により雑音共分散行列Φv(f)を計算する。
【０１４９】
【数３９】
2017141542.tif 000041
【０１５０】
次に、行列計算部３１は、各目的音ｎ（ｎ＝１，・・・，Ｎ−１）に対して、下記（４０）式により目的音共分散行列Φn(f)を求める。
【０１５１】
【数４０】
2017141542.tif 000042
【０１５２】
次に、行列計算部３１は、下記（４１）式により、観測共分散行列Φy(f)を求める。
【０１５３】
【数４１】
2017141542.tif 000043
【０１５４】
次に、ウィーナーフィルタ構成部３２は、各目的音ｎ（ｎ＝１，・・・，Ｎ−１）に対して、マルチチャネル・ウィーナー・フィルタＷn(f)を、下記（４２）式により求める。
【０１５５】
【数４２】
2017141542.tif 000044
【０１５６】
そして、目的音推定部３３は、下記（４３）式のように、マルチチャネル・ウィーナー・フィルタＷn(f)を観測信号ベクトルｙ(t,f)に適用することで、背景雑音や、ｎ番目以外の目的音の成分を抑圧し、目的音ｎの成分の推定値＾ｓn(t,f)を得ることができる。
【０１５７】
【数４３】
2017141542.tif 000045
【０１５８】
（実施形態３に係る目的音推定処理）
図４は、実施形態３に係る目的音推定処理の一例を示すフローチャートである。先ず、目的音推定装置３０の行列計算部３１は、マスク推定装置１０（１０−２）からマスク情報を取得する（ステップＳ２１）。
【０１５９】
次に、行列計算部３１は、雑音下目的音共分散行列Φn+v(f)を計算する（ステップＳ２２）。次に、行列計算部３１は、雑音共分散行列Φv(f)を計算する（ステップＳ２３）。次に、行列計算部３１は、目的音共分散行列Φn(f)を計算する（ステップＳ２４）。次に、行列計算部３１は、観測共分散行列Φy(f)を計算する（ステップＳ２５）。
【０１６０】
次に、ウィーナーフィルタ構成部３２は、マルチチャネル・ウィーナー・フィルタＷn(f)を構成する（ステップＳ２６）。次に、目的音推定部３３は、ステップＳ２６で構成したマルチチャネル・ウィーナー・フィルタＷn(f)  を観測信号ベクトルｙ(t,f)に適用することで、目的音ｎの成分の推定値＾ｓn(t,f)を得て、出力する（ステップＳ２７）。
【０１６１】
（性能実証結果）
実施形態３の性能実証として、実施形態１のマスク推定装置１０を含む目的音推定システム１００を用い、バスの中、カフェなどの背景雑音が存在する環境において、１人の話者がタブレットに向かって文章を読み上げている状況（実施形態１においてＮ＝２の場合）で、タブレットに装着されたＭ＝６個のマイクロホンで収録した信号に対して音声認識を行った。音声認識性能（単語誤り率）は、次のとおりであった。
【０１６２】
すなわち、マスク推定を行わず音声認識をした場合の単語誤り率は、14.29(%)であった。また、混合複素ワトソン分布でマスク推定をした後、マルチチャネル・ウィーナー・フィルタを適用した音声認識の場合の単語誤り率は、9.51(%)であった。実施形態１のマスク推定装置１０を含む目的音推定システム１００において、混合複素ビンガム分布でマスク推定した後、マルチチャネル・ウィーナー・フィルタを適用した場合の単語誤り率は、8.53(%)であった。以上から、実施形態３は、従来技術と比較して、音声認識性能が向上することが分かる。
【０１６３】
上記の実施形態のマスク推定装置１０（１０−２）及び目的音推定装置３０を含む目的音推定システム１００において行われる各処理は、全部又は任意の一部が、ＣＰＵ（Central  Processing  Unit）等の処理装置及び処理装置により解析実行されるプログラムにて実現されてもよい。また、マスク推定装置１０（１０−２）及び目的音推定装置３０を含む目的音推定システム１００において行われる各処理は、ワイヤードロジックによるハードウェアとして実現されてもよい。
【０１６４】
また、実施形態において説明した各処理のうち、自動的に行われるものとして説明した処理の全部又は一部を手動的に行うこともできる。もしくは、実施形態において説明した各処理のうち、手動的に行われるものとして説明した処理の全部又は一部を公知の方法で自動的に行うこともできる。また、実施形態においてフローチャートを参照して説明した各処理は、最終的な実行結果が同一である限りにおいて、適宜実行順序の入れ替えや並列実行を行うことができる。この他、上述及び図示の処理手順、制御手順、具体的名称、各種のデータやパラメータを含む情報については、特記する場合を除いて適宜変更することができる。
【０１６５】
（プログラムについて）
図５は、プログラムが実行されることにより、実施形態に係るマスク推定装置及び目的音推定装置を含む目的音推定システムが実現されるコンピュータの一例を示す図である。コンピュータ１０００は、例えば、メモリ１０１０、ＣＰＵ１０２０を有する。また、コンピュータ１０００は、ハードディスクドライブインタフェース１０３０、ディスクドライブインタフェース１０４０、シリアルポートインタフェース１０５０、ビデオアダプタ１０６０、ネットワークインタフェース１０７０を有する。コンピュータ１０００において、これらの各部はバス１０８０によって接続される。
【０１６６】
メモリ１０１０は、ＲＯＭ１０１１及びＲＡＭ１０１２を含む。ＲＯＭ１０１１は、例えば、ＢＩＯＳ等のブートプログラムを記憶する。ハードディスクドライブインタフェース１０３０は、ハードディスクドライブ１０３１に接続される。ディスクドライブインタフェース１０４０は、ディスクドライブ１０４１に接続される。例えば磁気ディスクや光ディスク等の着脱可能な記憶媒体が、ディスクドライブ１０４１に挿入される。シリアルポートインタフェース１０５０は、例えばマウス１０５１、キーボード１０５２に接続される。ビデオアダプタ１０６０は、例えばディスプレイ１０６１に接続される。
【０１６７】
ハードディスクドライブ１０３１は、例えば、ＯＳ１０９１、アプリケーションプログラム１０９２、プログラムモジュール１０９３、プログラムデータ１０９４を記憶する。すなわちマスク推定装置１０（１０−２）及び目的音推定装置３０の各処理を規定するプログラムは、コンピュータ１０００によって実行される指令が記述されたプログラムモジュール１０９３として、例えばハードディスクドライブ１０３１に記憶される。例えば、マスク推定装置１０（１０−２）及び目的音推定装置３０における機能構成と同様の情報処理を実行するためのプログラムモジュール１０９３が、ハードディスクドライブ１０３１に記憶される。
【０１６８】
また、実施形態の処理で用いられる設定データは、プログラムデータ１０９４として、例えばメモリ１０１０やハードディスクドライブ１０３１に記憶される。そして、ＣＰＵ１０２０が、メモリ１０１０やハードディスクドライブ１０３１に記憶されたプログラムモジュール１０９３やプログラムデータ１０９４を必要に応じてＲＡＭ１０１２に読み出して実行する。
【０１６９】
なお、プログラムモジュール１０９３やプログラムデータ１０９４は、ハードディスクドライブ１０３１に記憶される場合に限らず、例えば着脱可能な記憶媒体に記憶され、ディスクドライブ１０４１等を介してＣＰＵ１０２０によって読み出されてもよい。あるいは、プログラムモジュール１０９３やプログラムデータ１０９４は、ネットワーク（ＬＡＮ（Local  Area  Network）、ＷＡＮ（Wide  Area  Network）等）を介して接続された他のコンピュータに記憶されてもよい。そして、プログラムモジュール１０９３やプログラムデータ１０９４は、ネットワークインタフェース１０７０を介してＣＰＵ１０２０によって読み出されてもよい。
【０１７０】
実施形態は、本願が開示する技術に含まれると同様に、特許請求の範囲に記載された発明とその均等の範囲に含まれるものである。
【符号の説明】
【０１７１】
１０、１０−２、１０Ａ  マスク推定装置
１１、１１Ａ  特徴抽出部
１２、１２−２、１２Ａ  マスク更新部
１３、１３−２、１３Ａ  混合重み更新部
１４、１４−２  パラメータ更新部
１４Ａ  位置パラメータ更新部
１５Ａ  広がりパラメータ更新部
２０、２０−２、２０Ａ  記憶部
３０  目的音推定装置
３１  行列計算部
３２  ウィーナーフィルタ構成部
３３  目的音推定部
１０００  コンピュータ
１０１０  メモリ
１０２０  ＣＰＵ
【図１】
2017141542.tif 000046
【図２】
2017141542.tif 000047
【図３】
2017141542.tif 000048
【図４】
2017141542.tif 000049
【図５】
2017141542.tif 000050
【図６】
2017141542.tif 000051
【国際調査報告】
2017141542.tif 000052
2017141542.tif 000053
2017141542.tif 000054
2017141542.tif 000055