JP2013167698A - 音源ごとに信号のスペクトル形状特徴量を推定する装置、方法、目的信号のスペクトル特徴量を推定する装置、方法、プログラム - Google Patents

Abstract

【課題】スペクトルの値が周波数間で相関を持つようなモデルを用いる場合でも、目的音のスペクトル推定が効率的に行えるように、音源ごとに信号のスペクトル形状特徴量を推定する技術を提供する。
【解決手段】各音源に対応する、スペクトル形状特徴量の事前確率密度関数（スペクトル形状モデル）と、スペクトル形状特徴量が与えられた下でのスペクトル特徴量の条件付き確率密度関数（スペクトル観測モデル）とを用いる。各時間周波数点において最大のエネルギーを持つ音響信号の音源を表す占有的音源番号を潜在変数に持ち、全音源のスペクトル形状特徴量が与えられたもとでの観測信号のスペクトル特徴量の条件付き確率密度関数と、音源ごとに定められるスペクトル形状特徴量の事前確率密度関数の積とで表される最適化関数を、スペクトル形状モデルとスペクトル観測モデルを用いて最大化し、音源ごとのスペクトル形状特徴量および音源占有度を推定する。
【選択図】図２

Description

本発明は、複数の音源から出た音響信号が混ざってマイクロホンで収音された観測信号から、音源ごとに信号のスペクトル形状特徴量を推定する技術ならびに目的信号のスペクトル特徴量を推定する技術に関する。

図１に、非特許文献１，２，３に開示された従来の目的信号スペクトル特徴量推定装置の機能構成例を示す。この目的信号スペクトル特徴量推定装置９００は、観測信号に含まれる個々の音源（以下、番号m=1,2,・・・で区別する）について、当該音源に関連付けられた離散状態番号モデル記憶部９０１−ｍと状態依存スペクトルモデル記憶部９０２−ｍと離散状態番号選択部９０３−ｍを具備し、さらに、特徴抽出部９０４、音源占有度推定部９０５、目的音スペクトル推定部９０６を具備している。図１では、簡単のため音源数が２の場合を例示しているが、一般に３個以上の音源も考慮することができ、この場合は、離散状態番号モデル記憶部９０１−ｍと状態依存スペクトルモデル記憶部９０２−ｍとスペクトル番号選択部９０３−ｍは、音源数だけ別々のものが用意されているものとする。

以下の説明では、音源数はN(m)と仮定して説明する。音源m=1は目的音源とし、音源m>1は背景音源とする。nは時間、k(=1,2,・・・,N(k))は周波数を表すものとし、或る時間周波数点(n,k)で値をとる変数をx_n,kのように表記し、その値を全周波数でまとめてできるベクトル（例えば、スペクトル特徴量など）を参照する場合にx_n(=[x_n,1,・・・,x_n,N(k)]^T)と表記するものとする。Tは、行列やベクトルの非共役転置を表す。

離散状態番号モデル記憶部９０１−ｍは、音源mについて、短時間フレームに対するスペクトル特徴量の離散状態を表す離散状態番号q^(m)に関する事前確率関数p(q^(m))を記憶している。離散状態番号q^(m)は、N_q ^(m)(>0)個の番号(1からN_q ^(m)の自然数値)のいずれかを取るものとする。

状態依存スペクトルモデル記憶部９０２−ｍは、音源mについて、離散状態番号q^(m)が与えられた場合の音源mのスペクトル特徴量s^(m)の条件付き確率密度関数p(s^(m)|q^(m))を記憶している。

特徴抽出部９０４は、マイクロホンで収音した時間領域信号x(t)を入力として、各短時間フレームnに対する観測信号のスペクトル特徴量x_nを抽出する。

離散状態番号選択部９０３−ｍは、観測信号のスペクトル特徴量x_n、音源mに関する離散状態番号の事前確率関数p(q^(m))とスペクトル特徴量の条件付き確率密度関数p(s^(m)|q^(m))、音源占有度の推定値M^_n ^(m)(=[M^_n,1 ^(m),・・・,M^_n,N(k) ^(m)]^T)を受け取り、観測音に最も適合する音源mの離散状態番号の推定値q^_n ^(m)を求める。

音源占有度推定部９０５は、観測音のスペクトル特徴量x_nと、すべての音源mに関するスペクトル特徴量の条件付き確率密度関数p(s^(m)|q^(m))と、すべての音源mのそれぞれに関して選択された離散状態番号の推定値の組合せ{q^_n ⁽¹⁾,・・・,q^_n ^(N(m))}を受け取り、各時間周波数点において各音源mが最も大きなエネルギーを持つ事後確率である音源占有度の推定値M^_n,k ^(m)を求める。

さらに、目的音スペクトル推定部９０６は、目的音の音源占有度の推定値M^_n ⁽¹⁾と目的音の離散状態番号の推定値q^_n ⁽¹⁾と、特徴抽出部９０４が出力する観測信号のスペクトル特徴量x_nと、目的音のスペクトル特徴量の条件付き確率密度関数p(s⁽¹⁾|q⁽¹⁾)とを入力として、目的音のスペクトルの推定値を求める。

従来の目的信号スペクトル特徴量推定装置９００は、音源占有度推定部９０５と各離散状態番号選択部９０３−ｍのそれぞれが、音源占有度と観測信号のスペクトル特徴量に最も適合するスペクトル番号を効率的に求められるようにするために、以下の二つの前提を置いていた。
〈１〉
各時間周波数点(n,k)において、観測信号のスペクトル特徴量x_n,kは、各音源mのスペクトル特徴量s_n,k ^(m)のうち、値が大きいものに一致する。すなわち、
x_n,k=max{s_n,k ⁽¹⁾,・・・,s_n,k ^(N(m))} （１）
〈２〉
各音源信号のスペクトル特徴量s_n ^(m)は、離散状態番号q(m)が与えられたもとで、周波数ごとに独立である。すなわち、スペクトル特徴量の条件付き確率密度関数p(s_n ^(m)|q^(m))は、以下のように分解できる。ただし、s_n ^(m)=[s_n,1 ^(m),・・・,s_n,N(k) ^(m)]^Tである。

非特許文献１，２，３に開示されるように、上記の二つの仮定の下で、各音源の離散状態番号と音源占有度を交互に繰り返し更新することで、離散状態番号と音源占有度を推定できる。このとき、各繰り返しにおいて、離散状態番号は音源ごとに別々（独立）に更新することができ、また、音源占有度は周波数ごとに別々（独立）に更新することができる。すなわち、それぞれの更新において、音源間の組合せ、および周波数間の組合せを考慮する必要がないため、少ない計算コストで各更新を行うことができる。

中谷智広、荒木章子、吉岡拓也、藤本雅清、"DOAクラスタリングと音声の対数スペクトルHMMに基づく音源分離," 日本音響学会2010年秋季研究発表会講演論文集、pp.577-580, 9月, 2010年. 中谷智広、荒木章子、吉岡拓也、藤本雅清、"音源スペクトルHMMと音源方向モデルの教師無し同時学習に基づく多チャンネル音源分離," 日本音響学会2011年春季研究発表会講演論文集, pp.805-808, 3月, 2011年. 中谷智広、荒木章子、マーク・デルクロア、吉岡拓也、藤本雅清、"非定常雑音に頑健な統合的音声認識アプローチ:音源方向GMMと対数スペクトルGMMに基づく統計モデルベース音声強調," 日本音響学会2011年秋季研究発表会講演論文集, 9月, 2011年.

従来の目的信号スペクトル特徴量推定装置は、上記<２>の仮定が成立しない場合は、これらの値を求めるための計算コストが膨大になるという問題があった。したがって、少ない計算コストで目的信号のスペクトル推定を行うために、上記の仮定に従うようにする必要があった。このため、各信号のスペクトルのモデルとして、対数パワースペクトルの分布を離散状態モデルで表現した混合ガウス分布（もしくは、その離散状態の時間遷移もモデル化した隠れマルコフモデル）などを用い、特に、上記<２>の仮定を満たすために、混合ガウス分布（もしくは、隠れマルコフモデル）を構成する各ガウス分布の共分散行列は対角行列のみに限定する必要があった。

しかし、この限定は、スペクトルモデルが表現できる信号の分布を限定してしまい、信号が本来持つ分布を精度よく表現できるとは限らない。この結果、限定された分布で表現できるスペクトルの精度以上に、目的音および背景音のスペクトルの推定精度を高くできないという問題があった。例えば、より具体的には、音声認識でしばしば用いられるメル周波数ケプストラム係数（以下、MFCCという）に関する混合ガウス分布（もしくは、MFCCに関する隠れマルコフモデル）は、現在、音声信号のスペクトル形状のモデルとして最良なものの一つとして考えられているが、従来の方法では、これを信号のスペクトルモデルとして用いることはできなかった。この理由は、このモデルでは、各離散状態に対応する各ガウス分布がMFCCの分布を規定しており、かつ、MFCCの各係数は周波数間にまたがるスペクトル形状全体に関連するものであるため、各離散状態に対応する信号のスペクトルの分布は、一般的には、周波数間で独立になりえないからである。

このような課題に鑑みて、本発明は、従来技術における制限的なモデルと異なり、制限を持たないより一般的なモデル（スペクトルの値が周波数間で相関を持つような高精度なモデル）を用いる場合でも、目的信号のスペクトル推定が効率的に行えるように、複数の音源から出た音響信号が混ざってマイクロホンで収音された観測信号から、音源ごとに信号のスペクトル形状特徴量を推定する技術ならびに目的信号のスペクトル特徴量を推定する技術を提供することを目的とする。

本発明のスペクトル形状特徴量推定は、複数の音源それぞれからの音響信号が混ざって収音された観測信号から、各音響信号のスペクトル特徴量に対応するスペクトル形状特徴量を推定するものであって、各音源に対応する、スペクトル形状特徴量の事前確率密度関数（スペクトル形状モデル）と、スペクトル形状特徴量が与えられたもとでのスペクトル特徴量の条件付き確率密度関数（スペクトル観測モデル）とを用いる。具体的には、各時間周波数点において最大のエネルギーを持つ音響信号の音源を表す占有的音源番号を潜在変数に持つ最適化関数を、スペクトル形状モデルとスペクトル観測モデルを用いて最大化するとともに、音源ごとのスペクトル形状特徴量、および、全ての音源のスペクトル形状特徴量が与えられたもとで各音源が占有的音源番号で表される音源である事後確率（音源占有度）、を推定する。最適化関数は、全ての音源のスペクトル形状特徴量が与えられたもとでの観測信号のスペクトル特徴量の条件付き確率密度関数と、音源ごとに定められるスペクトル形状特徴量の事前確率密度関数の積とで表される。

本発明の目的信号スペクトル特徴量推定技術は、本発明のスペクトル形状特徴量推定技術によって推定された、音源ごとのスペクトル形状特徴量および音源占有度のうち、目的信号の音源に対応するスペクトル形状特徴量と、目的信号の音源の音源占有度と、目的信号の音源に対応するスペクトル観測モデルと、観測信号のスペクトル特徴量とから、目的信号のスペクトル特徴量を推定する。

本発明に拠ると、詳細は実施形態の説明に譲るが、スペクトル形状モデルとして、スペクトルの値が周波数間で相関を持つような高精度なモデルを導入することができる。そして、低い計算コストで、従来例よりも高精度に目的信号のスペクトルが推定できる。

従来例の目的信号スペクトル特徴量推定装置の機能構成例を示す図。 目的信号スペクトル特徴量推定装置の機能構成例を示す図。 目的信号スペクトル特徴量推定装置の機能構成例を示す図。

《本技術の概略》
まず本技術の概略を説明してから、本技術の詳細を説明する。
図２にスペクトル形状特徴量推定装置／目的信号スペクトル特徴量推定装置の構成例を示す。この目的信号スペクトル特徴量推定装置１００は、観測信号に含まれる複数の音源のうち、m番目の音源ごとに、当該音源に関連付けられたスペクトル形状モデル記憶部１０１−ｍとスペクトル観測モデル記憶部１０２−ｍとスペクトル形状推定部１０３−ｍを具備し、さらに、特徴抽出部１０４、音源占有度更新部１０５、および目的音スペクトル推定部１０６を具備している。スペクトル形状特徴量推定装置１００ｐは、目的音スペクトル推定部１０６を具備しない点で、目的信号スペクトル特徴量推定装置１００と異なる。図２では、簡単のため音源数が２の場合を例示しているが、３個以上の音源を考慮する場合は、スペクトル形状モデル記憶部１０１−ｍとスペクトル観測モデル記憶部１０２−ｍとスペクトル形状推定部１０３−ｍは、音源数だけ別々のものが用意されているものとする。

実施形態におけるスペクトル形状特徴量推定装置は、それ単体で独立に存在するよりは、得られた最尤スペクトル形状特徴量などを用いて目的信号のスペクトル特徴量を推定する装置（実施形態における目的信号スペクトル特徴量推定装置）を構成する構成要素として存在することが実用的な場合がある。さらに云えば、スペクトル形状特徴量推定装置は、目的信号スペクトル特徴量推定装置とは容易に分離可能に目的信号スペクトル特徴量推定装置を構成する構成要素ではなく、目的信号スペクトル特徴量推定装置自体を或る機能に着眼して片面的に評価したものと云うこともできる。要するに、スペクトル形状特徴量推定装置は、目的信号スペクトル特徴量推定装置そのものであることが凡そ実用的と言うことができる場合がある。
ただし、スペクトル形状特徴量推定装置が、単体独立の構成要素として存在すること、目的信号スペクトル特徴量推定装置とは容易に分離可能に目的信号スペクトル特徴量推定装置を構成する構成要素であることを排除する趣旨ではない。例えば各音源のスペクトル形状特徴量を推定すること自体を目的とするならば、スペクトル形状特徴量推定装置を単体独立の構成要素として実現することに何らの妨げは無い。
ここでは、目的信号スペクトル特徴量推定装置／スペクトル形状特徴量推定装置は、例えば専用のハードウェアで構成された専用機やパーソナルコンピュータのような汎用機といったコンピュータで実現されるとし、スペクトル形状特徴量推定装置は、目的信号スペクトル特徴量推定装置を構成する構成要素として説明する。
目的信号スペクトル特徴量推定装置／スペクトル形状特徴量推定装置を単体単独の構成要素として、これをコンピュータで実現する場合のハードウェア構成例は後述する。

以下の説明では、音源数はN(m)と仮定して説明する。音源m=1は目的音源とし、音源m>1は背景音源とする。nは時間、k(=1,2,・・・,N(k))は周波数を表すものとし、或る時間周波数点(n,k)で値をとる変数をx_n,kのように表記し、その値を全周波数でまとめてできるベクトル（例えば、スペクトル特徴量など）を参照する場合にx_n(=[x_n,1,・・・,x_n,N(k)]^T)と表記するものとする。Tは、行列やベクトルの非共役転置を表す。

スペクトル形状モデル記憶部１０１−ｍは、音源mについて、短時間フレームに対するスペクトル形状特徴量c^(m)の事前確率密度関数p(c^(m))を記憶している。この事前確率密度関数を、以下、スペクトル形状モデルと呼ぶ。

スペクトル観測モデル記憶部１０２−ｍは、音源mについて、スペクトル形状特徴量c^(m)が得られたもとでのスペクトル特徴量s^(m)の条件付き確率密度関数p(s^(m)|c^(m))を記憶している。与えられたスペクトル形状特徴量c^(m)に対応してスペクトル特徴量s^(m)をユニークには決定できないため、確率的な揺らぎが含まれると仮定している。この揺らぎをモデル化したものがスペクトル特徴量s^(m)の条件付き確率密度関数p(s^(m)|c^(m))である。以下、この関数をスペクトル観測モデルと呼ぶ。

特徴抽出部１０４は、マイクロホンで収音した時間領域信号x(t)を入力として、各短時間フレームnにおける観測信号のスペクトル特徴量x_nを抽出する。

スペクトル形状推定部１０３−ｍは、観測信号のスペクトル特徴量x_n、音源mの音源占有度M^_n ^(m)、音源mに関するスペクトル形状モデルp(c^(m))、音源mのスペクトル観測モデルp(s^(m)|c^(m))を受け取り、音源mのスペクトル形状特徴量c^_n ^(m)を推定する。

音源占有度更新部１０５は、特徴抽出部１０４から観測信号のスペクトル特徴量x_nを受け取り、すべてのスペクトル形状推定部１０３−ｍのそれぞれから、各音源mのスペクトル形状特徴量c^_n ^(m)を受け取り、すべてのスペクトル観測モデルから各音源mのスペクトル観測モデルp(s^(m)|c^(m))を受け取り、各音源mの音源占有度の推定値M^_n ^(m)を推定する。

目的音スペクトル推定部１０６は、特徴抽出部１０４から観測信号のスペクトル特徴量x_nを受け取り、スペクトル形状推定部１０３−１から目的信号のスペクトル形状特徴量の推定値c^_n ⁽¹⁾を受け取り、音源占有度更新部１０５から目的信号の音源占有度の推定値M^_n ⁽¹⁾を受け取り、スペクトル観測モデル記憶部１０２−１から目的信号のスペクトル特徴量の条件付き確率密度関数p(s⁽¹⁾|c⁽¹⁾)を受け取り、目的信号のスペクトルの推定値を求める。

さらに、本技術では、各音源mのスペクトル形状特徴量を効率的に求めるために、以下の前提を置く。
<１>
各時間周波数点(n,k)において、観測信号のスペクトル特徴量x_n,kは、音源mのスペクトル特徴量s_n,k ^(m)のうち、値が大きいものに一致する。すなわち、式（３）が成立する。以下、この仮定をLogMaxモデルと呼ぶ。
x_n,k=max{s_n,k ⁽¹⁾,・・・,s_n,k ^(N(m))} （３）
〈２〉
各音源信号のスペクトル特徴量s_n ^(m)は、各音源のスペクトル形状特徴量c_n ^(m)が与えられたもとで、周波数ごとに独立である。すなわち、スペクトル特徴量の条件付き確率密度関数p(s_n ^(m)|c_n ^(m))は、以下のように分解できる。ただし、s_n ^(m)=[s_n,1 ^(m),・・・,s_n,N(k) ^(m)]^Tである。

この仮定のうち、仮定<１>は従来例と同じであり、仮定<２>が従来例と異なる。

上記の構成において、従来例との違いで最も重要な点は、従来例では、各音源のスペクトル特徴量のモデルを、離散状態モデルと状態依存スペクトルモデルの二つの組合せで表現していたのに対し、本発明では、スペクトル形状モデルとスペクトル観測モデルの二つの組合せで表現している点である。

従来例では、スペクトルの形状に関する特徴量はスペクトル観測過程と混在して状態依存スペクトルモデルで表現されており、そのモデルにおいて、スペクトルの値が周波数間で独立であると仮定する必要があった。

他方、本発明では、周波数間で独立であると仮定する必要があるのはスペクトル観測モデルのみであり、スペクトル形状モデルについては何も規定されない。その結果、スペクトルの値が周波数間で相関を持つようなスペクトル形状モデルを導入することができる。例えば、MFCCの混合ガウスモデルなど、高精度なスペクトル形状モデルを利用することができる。

また、上記の二つの仮定の下で、各音源のスペクトル形状特徴量と音源占有度を交互に繰り返し更新することで、スペクトル形状特徴量と音源占有度を推定できる。この時、各繰り返しにおいて、スペクトル形状特徴量は音源ごとに別々に更新することができ、また、音源占有度の更新は周波数ごとに別々に更新することができる。すなわち、それぞれの更新において、音源間で値の組合せ、および周波数間で値の組合せを考慮する必要がないため、少ない計算コストで各更新を行うことができる。

《本技術の詳細》
c_n=H(s_n)を、信号のスペクトル特徴量s_nからスペクトル形状特徴量c_nを抽出する特徴量変換関数とする。ここで、nは短時間フレームの番号を表し、スペクトル特徴量s_nは各周波数に対応する連続スカラー値s_n,kを要素にもつベクトルs_n=[s_n,1,s_n,2,・・・,s_n,N(k)]^Tと定義され、スペクトル形状特徴量c_nは形状パラメータの各次元に対応する連続スカラー値c_n,hを要素に持つベクトルc_n=[c_n,1,c_n,2,・・・,c_n,N(h)]^Tと定義され、H(・)は、スペクトルの形状に関する特徴量を抽出する任意の関数とする。例えば、s_nを音声認識などでしばしば用いられる対数メルフィルタバンクの出力、c_nをそのMFCCとすると、H(s_n)はs_nの離散コサイン変換D(s_n)に対応する。すなわち、H(s_n)=D(s_n)となる。また、s_nを対数パワースペクトルとし、c_nを対応するMFCCとすると、H(s_n)は、s_nに対数関数の逆変換、すなわち指数関数(exp(・)と表記)を適用したのち、メルフィルタバンク処理(mfb(・)と表記)を施し、再度、対数変換(log(・))を適用することに対応する。すなわち、H(s_n)=log(mfb(exp(s_n)))である。

以下の議論では、基本的に、複数の短時間フレームを観測信号から特徴抽出部１０４が抽出することを前提とするが、下記の説明（例えば実施形態１，２）においては、議論を簡単にするため、一つの短時間フレームnに閉じた処理として説明する。複数の短時間フレームに対しては、各短時間フレームに対して同じ処理を個別に適用するものとする。その他の実施形態では、必ずしもその限りではない。

[最適化関数の定義]
いま、s_n ^(m)を音源mの短時間フレームnにおけるスペクトル特徴量とする。また、すべての音源mに関するスペクトル特徴量およびスペクトル形状特徴量をひとまとめにして、以下のように記述することにする。ここで、N(m)は音源数を表す。
S_n={s_n ⁽¹⁾,・・・,s_n ^(N(m))} （５）
S_n,k={s_n,k ⁽¹⁾,・・・,s_n,k ^(N(m))} （６）
C_n={c_n ⁽¹⁾,・・・,c_n ^(N(m))} ただし、c_n ^(m)=H(s_n ^(m)) （７）

本技術では、目的音スペクトル推定のためのステップとして、以下で定義されるように、各音源mのスペクトル形状特徴量を事後確率最大化(MAP)推定で求める。

このMAP推定は、後述するように、主としてスペクトル形状推定部１０３−ｍと音源占有度推定部１０５による繰り返し処理により実現される。

MAP推定を行うために、本技術では、音源mのスペクトル形状モデル（すなわち、スペクトル形状特徴量の事前確率密度関数）p(c_n ^(m))、および、スペクトル観測モデル（すなわち、スペクトル形状特徴量が与えられたもとのでスペクトル特徴量の条件付き確率密度関数）p(s_n ^(m)|c_n ^(m))が事前に与えられているものとする（もしくは、後述するように目的音スペクトルを推定する対象となる観測信号から推定できるものとする）。

スペクトル形状モデルp(c_n ^(m))に関して、特段の前提を置いていないので、c_n ^(m)の確率密度関数を表すものであれば、どんなものでもモデルとして利用することができる。代表的なものとして、c_n ^(m)に関する正規分布、ガウス混合モデル、隠れマルコフモデルなどがあげられる。

スペクトル観測モデルp(s_n ^(m)|c_n ^(m))は、前述の特徴量変換H(s)の逆過程として定義される。しかし、一般には、c_n=H(s_n)は多対１の変換となるため、その逆変換はユニークには定められない。したがって、その定め方には任意性がある。ここでは、一例をあげる。いま、s=G(c)を、c=H(s)を満たすsの集合の中の代表点に変換する関数とし、H(s)の疑似逆変換と呼ぶ。さらに、e=s-G(H(s))を逆変換誤差と呼ぶ。そして、eを期待値0、共分散行列Ξの正規分布に従うと仮定すると、スペクトル観測モデルは、以下のように定義できる。ここで、Nd(x;μ,Σ)は、平均μ、共分散行列Σの正規分布の確率密度関数を表す。
p(s_n ^(m)|c_n ^(m))=Nd(s_n ^(m);G(s_n ^(m)),Ξ) （９）

さらに、このスペクトル観測モデルp(s_n ^(m)|c_n ^(m))は、式（４）に基づき、各周波数で定義される条件付き確率密度関数の積に分解できることが仮定される。これは、上記のモデルの場合、共分散行列は対角行列Ξ=diag{ξ_k}で規定されることを意味する。s=G(c)となるs=[s₁,s₂,・・・,s_N(k)]のk番目の要素を与える関数をs_k=G_k(c)と表すとすると、式（９）は、以下のように分解可能である。

なお、上記では、eの期待値がE{e}=0となることを仮定した（ここで、E{・}は確率変数の期待値をとる関数である）。他方、この値が0でない場合も、G'(c)=G(c)-E{e}としたG'(c)を疑似逆変換とすることで、逆変換誤差の期待値を0としたモデル化が可能である。さらに、上記では、スペクトル観測モデル（すなわち、逆変換誤差）が正規分布に従うと仮定したが、正規分布に限らず一般の分布を用いても同様の定義が可能である。

次に、本発明で導入した式（３）の仮定を確率的に扱うために、従来例と同様に、以下の関係式を導入する。

ここで、d_n=[d_n,1,・・・,d_n,N(k)]は、各時間周波数において、最もエネルギーの大きな音源の番号を示す占有的音源番号を表す。例えば、d_n,k=1は、時間周波数点(n,k)において、m=1番の音源（すなわち、目的音）が最も占有的な音源であることを示す。δ(・)は、ディラックのデルタ関数である。式（１１）は、観測信号のスペクトル特徴量x_n,kと占有的音源番号で示された音源のスペクトル特徴量s_n,k ^(m)［m=d_n,k］は一致することを意味し、式（１２）は、占有的音源番号が最もスペクトル特徴量の値が大きい音源の番号に一致することを意味している。これは、式（３）と等価であることは明らかである。すると、x_n,k,d_n,kおよびC_nの間の関係は、式（４）,式（１１）および式（１２）に基づき以下のように導出できる。

式（１３）は、全音源のスペクトル形状特徴量が与えられたもとでの、各周波数における、観測信号のスペクトル特徴量と占有的音源番号の条件付き同時確率密度関数であって、以下の二種類の確率関数の積によって定められる関数を表している。
１．占有的音源番号に一致する音源のスペクトル観測モデルにおいて、当該音源のスペクトル特徴量が観測信号のスペクトル特徴量と同一の値をとると規定された場合の確率関数（右辺の第一要素）
２．占有的音源番号に一致する音源以外の音源のスペクトル観測モデルにおいて、当該音源のスペクトル特徴量が観測信号のスペクトル特徴量の値以下の値をとると規定された場合の確率関数（右辺の第二要素中の各積分項）

最終的に、式（１３）を用いて、式（８）中のMAP関数p(x_n,C_n)は、以下のように定義される。

ここで、d_n,kは隠れ変数として扱われている。式（１４）の右辺第一項は、スペクトル観測モデルとLogMaxモデルに対応し、式（１４）の右辺第二項は、各音源のスペクトル形状モデルに対応する。本技術では、式（１４）を最大にするスペクトル形状特徴量を全音源に関して求めることで、スペクトル形状モデル、スペクトル観測モデル、LogMaxモデルのすべてを考慮した推定を実現する。

より詳しくは、式（１４）は、観測信号のスペクトル特徴量x_nと全音源のスペクトル形状特徴量C_nの同時確率密度関数であり、右辺のように二つの要素の積からなる。式（１４）の右辺の一つ目の要素は、全音源のスペクトル形状特徴量が与えられたもとでの観測信号のスペクトル特徴量の条件付き確率密度関数（すなわち、p(x_n|C_n))である。以下、これを、観測信号スペクトル観測モデルと呼ぶ。式（１４）の右辺の二つ目の要素は、音源ごとに定められるスペクトル形状特徴量の事前確率密度関数であるスペクトル形状モデルの積（すなわち、Π_l p(c_n ^(l))）である。
観測信号スペクトル観測モデルは、式（１４）の右辺にあるように以下の特徴を持つ。
１．観測信号スペクトル観測モデルは、各周波数に対応する条件付き確率密度関数の積に分解可能（すなわち、p(x_n|C_n)=Π_k p(x_n,k|C_n))。
２．各周波数に対応する条件付き確率密度関数は、当該周波数において、どの音源が最も占有的な音源であるかを示す占有的音源番号を潜在変数として持つ。（すなわち、p(x_n,k|C_n)=Σ_d p(x_n,k,d_n,k|C_n)）
３．各周波数に対応する条件付き確率密度関数は、音源ごとに定められるスペクトル観測モデルとLogMaxモデルに基づき、式（１３）のように定められる。

本技術は、式（１４）を最大化する各音源のスペクトル形状特徴量を推定するとともに、各音源が占有的音源番号に一致する事後確率関数である各音源の音源占有度を推定する。さらに、目的信号について推定されたスペクトル形状特徴量および音源占有度と、目的信号のスペクトル観測モデルに基づき、目的信号のスペクトルを推定する。

[最適化のアルゴリズム]
式（１４）を最大化するアルゴリズムとしては、共役勾配法、準ニュートン法などの一般の非線形最適化アルゴリズムを適用することができる。これらの方法は、各非線形最適化アルゴリズムに基づき、一般的な方法で導出できるので、ここでは説明を省略する。

他方、式（１４）は隠れ変数を含む関数であるので、期待値最大化(EM)アルゴリズムを用いて効率的に最大化することもできる。以下では、この最適化アルゴリズムについて詳しく述べる。EMアルゴリズムで用いられる補助関数は、以下のように定義される。

いま、音源mの音源占有度の推定値をM^_n,k ^(m)=E{p(d_n,k=m|x_n,k,C_n)|C^_n}と表記すると、式（１７）はさらに以下のように展開できる。

EMアルゴリズムに従うと、E-stepにおいて式（１８）を計算し、M-stepにおいて式（１９）を最大化するC_nを求め、E-stepとM-stepを繰り返すことによって、前記の最適化関数を最大化するスペクトル形状特徴量C_nを求めることができる。

ここで注意すべきは、E-stepにおいて、式（１８）は、各周波数で独立に計算できることである。このため、音源占有度の更新において、周波数間の値の組合せを考慮する必要がない。さらに、M-stepにおいて、式（１９）は、式（２０）のように、音源ごとの最適化関数に分解して最大化することができる。したがって、式（２０）に従い、音源ごとに独立にスペクトル形状特徴量を更新することができる。このため、音源間の値の組合せを考慮する必要がない。その結果、EMアルゴリズムの繰り返し推定に基づき、音源mの音源占有度とスペクトル形状特徴量を、効率的に更新することができる。

[目的信号のスペクトル特徴量の推定]
目的音スペクトル推定部１０６は、観測信号のスペクトル特徴量x_n,kと、目的信号の音源占有度の推定値M^_n,k ⁽¹⁾と、スペクトル形状特徴量の推定値c^_n ⁽¹⁾と、スペクトル観測モデルp(s_n ⁽¹⁾|c_n ⁽¹⁾)と、を入力として、目的信号のスペクトルの推定値s^_n,k ⁽¹⁾を最小自乗誤差推定により求める。推定の方法は次式によって行う。

[実施形態１]
上記の説明に基づき、本技術では、以下の手順で目的音スペクトルの推定が行える。
ステップ１．各短時間フレームnに対して、特徴抽出部１０４が、観測信号のスペクトル特徴量x_nを抽出する。
ステップ２．各音源ｍに対応するスペクトル形状推定部１０３−ｍが、スペクトル形状特徴量の推定値c^_n ⁽¹⁾を初期化する。例えば、観測信号のスペクトル特徴量x_nと特徴量変換関数H(・)を用いて、c^_n ⁽¹⁾=H(x_n)とする。
ステップ３．(a)(b)を収束するまで繰り返す。
(a)音源占有度更新部１０５は、各周波数ごとに独立に、式（１８）を計算することで、各音源mの音源占有度の推定値M^_n,k ^(m)を更新する(E-step)。
(b)各音源mについて、スペクトル形状推定部１０３−ｍが、式（２０）を最大化するc_n ^(m)を求めることで、スペクトル形状特徴量の推定値c^_n ^(m)を更新する(M-step)。このとき、式（２０）は、一般には、非線形関数となるため、その最大化は、共役勾配法、準ニュートン法、ニュートン法などの一般的な非線形最適化法により実現される。
ステップ４．目的音スペクトル推定部１０６は、式（２１）により、目的信号のスペクトル特徴量の推定値s^_n,k ⁽¹⁾を求める。

なお、一例として、上記３．(b)において、ニュートン法に従ってスペクトル形状特徴量を更新する場合の更新式は、以下のようになる。

[実施形態２]
実施形態１より、さらに具体的な実施形態について説明する。まず、観測信号のスペクトル徴量としてメルフィルタバンクの出力を用い、スペクトル形状特徴量としてMFCCを用い、スペクトル形状モデルとしてMFCCの混合ガウス分布を用いるとする。

各音源mのスペクトル形状特徴量c_n ^(m)は、混合ガウス分布を用いて、以下のようにモデル化されているものとする。ここで、p(c_n ^(m)|i)とα_iは、それぞれ、混合番号iに対するガウス分布とその混合比を表す。各ガウス分布の平均μ_i ^(m)と共分散行列Σ_i ^(m)は、事前に、音響信号のデータベースなどを用いて求められているものとする。このモデルにおいて、混合番号iは隠れ変数として扱われる。

実施形態２において、実施形態１との最も大きな違いは、EMアルゴリズムの隠れ変数として新たに各音源mのスペクトル形状モデルにおける混合ガウス分布の混合番号iを含んでいることである。これにともない、EMアルゴリズムで用いる補助関数は、以下のように修正される。

さらに、Z^_n,i ^(m)=E{p(i|c_n ^(m))|c^_n ^(m)}と表記すると、式（２０）は以下のように修正される。

この結果、以下の手順で目的音スペクトルの推定が行えることになる。
ステップ１．各短時間フレームnに対して、特徴抽出部１０４が、観測信号に関するメルフィルタバンク出力x_nを抽出する。
ステップ２．各音源mに対応するスペクトル形状推定部１０３−ｍが、MFCCの推定値c^_n ^(m)を初期化する。例えば、H(x_n)を離散コサイン変換とし、c^_n ^(m)=H(x_n)とする。
ステップ３．(a)(b)(c)を収束するまで繰り返す。
(a)音源占有度更新部１０５は、各周波数ごとに独立に、式（１８）を計算することで、各音源mの音源占有度の推定値M^_n,k ^(m)を更新する(E-step1)。
(b)各音源mについて、スペクトル形状推定部１０３−ｍが、式（２８）に基づき、Z^_n,i ^(m)を更新する(E-step2)。
(c)各音源mについて、スペクトル形状推定部１０３−ｍが、式（２７）を最大化するMFCC c_n ^(m)を求めることで、MFCCの推定値c^_n ^(m)を更新する。このとき、式（２７）は、一般には、非線形関数となるため、その最大化は、共役勾配法、準ニュートン法、ニュートン法などの一般的な非線形最適化法により実現される(M-step)。
ステップ４．目的音スペクトル推定部１０６は、式（２１）により、目的信号のフィルタバンク出力（＝スペクトル特徴量）の推定値s^_n,k ^(m)を求める。

次に、実施形態２において、さらに、各音源mに対するスペクトル観測モデルp(s_n,k ^(m)|c_n ^(m))について、より具体的なモデルを用いた場合について説明する。

この実施形態２では、疑似逆行列G(c)を、以下によって定義される線形回帰モデルでモデル化するものとする。ただし、Aは行列(N(k)×N(m))、bはベクトル(N(k)×1)を表す。
G(c)=Ac+b (29)

Aとbの値は、事前に音響信号のデータベースにより学習されるか、観測信号を用いて学習されるものとする。すなわち、いま学習用のデータベース（もしくは、観測信号）から、複数のスペクトル特徴量x_nと対応するスペクトル形状特徴量c_n=H(x_n)の組合せが与えられているときに、Aとbは、以下のように定められるものとする。

また、逆変換誤差eの分散ξ_kの値は、式（３０）で得られたAとbを用いて得られる平均自乗回帰誤差E{|x_n-(Ac_n+b)|²}の、各周波数kに対応する値として定められるものとする。上記の定義に基づき、スペクトル観測モデルp(s_n,k ^(m)|c_n ^(m))を、式（１０）により定めることができる。なお、Aとbの推定に用いることができるデータに基づき、各音源ごとに個別のスペクトル観測モデルを用意してもよいし、すべての音源で共通のスペクトル観測モデルを用意してもよい。

この場合のアルゴリズムは、実施形態２のアルゴリズムと比べて、疑似逆変換G(c)の定義が変わるのみである。

G(c)の定義に関する、もうひとつの別の例として、Aを離散コサイン変換行列に対するムーアペンローズ疑似逆変換行列、b=0と定めることも考えられる。これは、離散コサイン変換の逆変換のうち、変換後の特徴量のノルムが最小になるものを選択したことになる。このように、本技術では、様々な基準に基づく逆変換を、疑似逆変換G(c)として採用することができる。

[実施形態２の変形例]
実施形態２において、観測信号や各音源のスペクトル特徴量として対数パワースペクトルを用いる場合について説明する。これは、上記のアルゴリズムにおいて、メルフィルタバンク出力に関する処理を対数パワースペクトルに関する処理に、単純に置き換えるだけで実現できる。上記のアルゴリズムに基づき、変更点をまとめると以下の通りである。
＝１＝
ステップ１．において、特徴抽出部１０４は、メルフィルタバンク出力のかわりに対数パワースペクトルを観測信号のスペクトル特徴量x_nとして抽出する。
＝２＝
ステップ２．の初期化においては、特徴量変換関数として、一例として挙げた、対数パワースペクトルからMFCCに変換する関数H(x_n)=log(mfb(exp(x_n)))を用いる。
＝３＝
ステップ３．(a), ３．(c)において用いる、各音源mに対するスペクトル観測モデルは、MFCCから対数パワースペクトルを求める逆変換G(c)に基づき定められるとする。
＝４＝
ステップ４．においては、目的信号のスペクトル推定値として、対数パワースペクトルを推定する。このため、観測信号のスペクトル特徴量x_nとして、観測信号の対数パワースペクトを用い、ステップ３．(a), ステップ３．(c)と同様に、スペクトル観測モデルは、MFCCから対数パワースペクトルを求める逆変換G(c)に基づき定められるものを用いる。

上記のアルゴリズムでは、最終的に得られるものが目的信号の対数パワースペクトルの推定値s^_n,k ^(m)である。したがって、ここから目的信号のパワースペクトルの推定値も、exp(s^_n,k ^(m))として求めることができる。さらに観測信号の位相情報を用いて、離散逆フーリエ変換とオーバラップ加算合成技術などにより、目的信号の波形の推定値をも得ることができる。

[実施形態３]
実施形態２において、スペクトル形状モデルとしてMFCCに関する混合ガウス分布を用いる代わりに、MFCCに関する隠れマルコフモデルを用いる場合を構成する。本実施形態では、簡単のためすべての観測信号が与えられてからすべての推定処理を行うバッチ処理を前提として説明する。ただし、隠れマルコフモデルを逐次処理で動作させることは周知の技術を用いれば可能であり、それらの技術を用いて本実施形態を逐次的に動作させることも可能である。

まず、{c_n ^(m)}をすべての短時間フレームnに関するc_n ^(m)の集合を表すとする。各音源mのスペクトル形状モデルp({c_n ^(m)})は、隠れマルコフモデルを用いると、以下のようにモデル化できる。ただし、i₀は初期状態、i_n for n>0は、短時間フレームnにおける状態、π_i0は初期状態がi₀となる確率、π_i,jは、状態iから状態jへの遷移確率、p(c_n ^(m)|i)は状態iにおけるスペクトル形状特徴量の出力確率密度関数、Σ_Iは、すべての状態系列の組合せに関する総和を表す。

隠れマルコフモデルでは、観測特徴量の系列が与えられたもとで、各短時間フレームにおいて各状態iをとる事後確率関数Z_n,i=p(i_n=i|{c_n ^(m)})は、forward-backwardアルゴリズムやViterbiアルゴリズムなどを用いて効率的に求められることが知られている。

スペクトル形状モデルに隠れマルコフモデルを用いた場合、EMアルゴリズムで用いる補助関数は、以下のように修正される。

さらに、Z^_n,i ^(m)=E{p(i_n|{c_n ^(m)})|{c^_n ^(m)}}と表記すると、式（１９）、式（２０）は以下のように修正される。

したがって、実施形態２の補助関数との主な違いは、式（３５）のZ^_n,i ^(m)の計算方法のみとなる。Z^_n,i ^(m)は、M-stepにおいて更新されたスペクトル形状特徴量の推定値{c^_n ^(m)}を出力値として持つ隠れマルコフモデルが短時間フレームnで状態iをとる事後確率p(i_n|{c^_n ^(m)})と同一であるので、前記のように、forward-backwardアルゴリズムやViterbiアルゴリズムなどを用いて効率的に求めることができる。

最終的に、実施形態２とのアルゴリズムの違いは、ステップ３．(b)の処理が以下のように修正されるのみである。
ステップ３．(b)各音源mごとに、スペクトル形状推定部１０３−ｍが、forward-backwardアルゴリズムなどを用いて、Z^_n,i ^(m)=p(i_n|{c^_n ^(m)})を求める。

[実施形態４]
これまでの実施形態では、すべての音源に関して、スペクトル形状モデルが事前に与えられていると仮定しているが、実際の環境では、必ずしもそのような状況は期待できない。例えば、目的信号のスペクトル形状モデルは事前に学習しておくことは可能だが、背景雑音のスペクトル形状モデルは事前に学習しておくことができない場合などがある。この場合、観測信号から得られる情報のみで、事前学習されていない音源のスペクトル形状モデルを学習する必要がある。そのような状況に対応するための方法として、実施形態４を説明する。

実施形態４では、実施形態２や実施形態３で、スペクトル形状特徴量c_n ^(m)のみをEMアルゴリズムにより求めるべきパラメータとして推定していたのに対し、スペクトル形状モデルである混合ガウス分布のパラメータ(混合比、各ガウス分布の平均と共分散)や隠れマルコフモデルのパラメータ（初期状態確率、状態遷移確率、出力確率分布の平均と共分散）なども求めるべきパラメータに含めて推定する。まず、スペクトル形状特徴量を求める問題を、以下の最適化関数を最大化する問題として定義する。ここで、θは、求めるべき全パラメータを含んだ集合であり、スペクトル形状特徴量C_n以外に、スペクトル形状モデルの未知パラメータも含まれているとする。

ここでは、スペクトル形状モデルのパラメータが、一つ以上の音源に関して、事前に与えられていない場合のアルゴリズムについて説明する。簡単のため、実施形態２や実施形態３のアルゴリズムとの違いについてのみ説明する。実施形態２に対する実施形態３の修正の場合と同様に、この場合も、アルゴリズムの修正は、ステップ３．(b)についてのみになる。具体的には、実施形態４では、スペクトル形状モデルが事前に与えられていない音源mについて、ステップ３．(b)の処理が以下の二つの処理で構成される。
ステップ３．(b)-1:初期化処理、もしくはM-stepにおいて更新されたc^_n ^(m)を用いて、まず、スペクトル形状モデルのパラメータを推定する。すなわち、スペクトル形状特徴量としてc^_n ^(m)が生成される尤度が最大になるように、スペクトル形状モデルのパラメータを推定する。パラメータの推定には、スペクトル形状モデルが採用している確率分布に適したアルゴリズムを用いればよい。例えば、混合ガウス分布や隠れマルコフモデルの場合は、EMアルゴリズムを用いて効率的に、パラメータ推定が行えることが知られている。
ステップ３．(b)-2:スペクトル形状モデルのパラメータ推定後に、実施形態２や実施形態３と同様の手続きにより、Z^_n,k ^(m)を求める。

以上の処理により、スペクトル形状モデルが事前に与えられていない場合でも、ステップ３．(b)の処理を遂行できることになる。また、上記以外の目的信号スペクトル推定の処理においても、上記で求めたスペクトル形状モデルを利用することができるので、特に、アルゴリズムの変更は生じない。

なお、上記のスペクトル形状モデルのパラメータ推定は、スペクトル形状特徴量の推定値c^_n ^(m)が与えられたもとで、この推定値をスペクトル形状モデルが出力する尤度が大きくなるように実施される。したがって、実施形態４の構成に利用できるスペクトル形状モデルは上記で挙げたものに限定されず、与えられた出力系列に基づきパラメータ推定が可能な任意のスペクトル形状モデルが含まれる。

[実施形態５]
本技術は、非特許文献１−３に示された従来例と同じく、音源占有度をEMアルゴリズムのE-stepで更新しながら繰り返し推定する構成をとる。したがって、観測信号が複数のマイクロホンから同時に収録されている場合には、従来例と同じ方法に基づき、観測信号から各音源の位置情報に基づく特徴量を抽出しつつ、E-stepで、音源位置の情報を考慮して音源占有度の更新を行うような構成が可能である。本実施形態では、この構成について説明する。

以下、図３を参照して、本実施形態の説明を行う。図中では、複数のマイクロホン信号を区別するために、x⁽¹⁾(t)のように右肩にマイクロホン番号を記入している。簡単のため、マイクロホン数が２の場合を例示しているが、従来例のようにマイクロホンの数は２以上であればいくつでもよい。特徴抽出部１０４は、これまでの実施形態と同様に、どれか一つの観測信号から観測信号のスペクトル特徴量x_nを抽出する。ただし、複数のマイクロホンを用いて、何らかの処理（例えば、他の音声強調技術を用いて目的音を強調するなど）を加えた結果得られるスペクトル特徴量を抽出するような構成も可能である。各音源に対するスペクトル形状モデル、スペクトル観測モデル、スペクトル形状推定部、および目的音スペクトル推定部は、これまでに実施形態として説明したもののいずれかと同じ、もしくは、同様の方法で構成されるものを用いるとする。本実施形態では、これらの機能部についての説明は、特に必要な場合を除き省略する。

他方、本実施形態では、新たな機能部として、スペクトル形状特徴量推定装置１００ｐ（つまり目的信号スペクトル特徴量推定装置１００）は、音源位置特徴抽出部１１１と、各音源mに関連付けられた音源位置状態推定部１１２−ｍを具備している。図中では、音源位置状態推定部１１２−ｍの数が２の場合を例示しているが、実際には、音源数に応じて音源位置状態推定部１１２−ｍを用意することを前提としている。これまでの実施形態と同様、m=1が目的音に対応する音源位置状態推定部１１２−ｍである。

音源位置特徴抽出部１１１は、複数のマイクロホンによる観測信号x^(m)(t)から音源位置特徴量a_n=[a_n,1,a_n,2,・・・,a_n,N(k)]^Tを抽出する。ここで、a_n,kは、時間周波数点(n,k)における、音源位置特徴量である。非特許文献１−３などと同様に、音源位置特徴量としては、音源位置に依存して異なる値を取る傾向を持つものを採用すればよく、例えば、マイクロホン間位相差、マイクロホン間時間差、マイクロホン間強度差など様々な選択肢がある。本明細書では、一例として、以下の正規化複素スペクトルを音源位置特徴量として採用する場合について説明する。ただし、X_n,k=[X_n,k ⁽¹⁾,X_n,k ⁽²⁾,・・・,X_n,k ^(N(m))]^Tとし、X_n,k ^(m)を短時間フレームnに対応するx^(m)(t)の短時間フーリエ変換のk番目の周波数要素である。また、|・|はベクトルのノルムを表す。
a_n,k=X_n,k/|X_n,k| (37)

次に、各音源位置状態推定部１１２−ｍは、音源位置特徴量a_n,kと音源位置パラメータの推定値φ^^(m)とに基づき、以下の値を求める。p(a_n,k;φ^(m))は、周波数kにおける、音源mの音源位置特徴量に関する確率密度関数を表し、以下では、w^_n,k ^(m)を音源mの音源位置状態の推定値と呼ぶ。
w^_n,k ^(m)=p(a_n,k;φ^(m)) (38)

音源位置特徴量に関する確率密度関数としては、音源位置特徴量に依存して、様々なものが知られている。例えば、上記で示した正規化複素スペクトルを音源位置特徴量とする場合の確率密度関数は、参考文献１などに記載されている。
（参考文献１）Hiroshi Sawada, Shoko Araki, and Shoji 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, 2011.

音源位置パラメータの推定値φ^^(m)は、観測信号と同じ収録条件のもと複数のマイクロホンで収音した音響信号を学習データとして、事前に学習するか、もしくは、観測信号から直接学習することができる。これらの学習方法は、音源位置特徴量に依存して、様々なものが知られている。例えば、上記で示した正規化複素スペクトルを音源位置特徴量とする場合の学習方法は、参考文献１に詳しい。また、本技術の目的信号スペクトル推定において、音源占有度の推定値とスペクトル形状特徴量の推定値を繰り返し更新するのと並行して、各繰り返しにおいて、各音源位置状態推定部１１２−ｍが音源占有度更新部１０５から音源占有度の推定値を受け取り、音源位置パラメータの推定値を更新することも可能である。目的信号スペクトル推定において、音源位置パラメータの推定値を繰り返し更新する方法については、非特許文献２に記載されている。

さらに、音源占有度更新部１０５は、観測信号のスペクトル特徴量x_nを受け取り、各音源mのそれぞれに関するスペクトル形状特徴量の推定値c^_n ^(m)とスペクトル観測モデルp(s_n ^(m)|c_n ^(m))を受け取るのに加えて、各音源mの音源位置状態の推定値w^_n,k ^(m)を受け取り、各音源の音源占有度の推定値M^_n,k ^(m)を更新する。この更新は、非特許文献２と同様、以下のように実現できる。

実施形態２でスペクトル特徴量として対数パワースペクトルを用いる場合の変形として、本実施形態のアルゴリズムをまとめると以下のようになる。
ステップ１．(a)(b)の処理によって特徴抽出する。
(a)各短時間フレームnに対して、特徴抽出部１０４が、観測信号に関する対数パワースペクトルx_nを抽出する。
(b)音源位置特徴抽出部１１１が観測信号に関する音源位置特徴量を式（３７）により抽出する。
ステップ２．(a)(b)の処理によって初期化を行う。
(a)各音源mに対応するスペクトル形状推定部１０３−ｍが、MFCCの推定値c^_n ^(m)を初期化する。例えば、H(x_n)=log(mfb(exp(x_n)))とし、c^_n ^(m)=H(x_n)とする。
(b)各音源mごとに、音源位置状態推定部１１２−ｍが、式（３８）により、音源位置状態の推定値w^_n,k ^(m)を求める。
ステップ３．(a)-(c)を収束するまで繰り返す。
(a)音源占有度更新部１０５は、各周波数ごとに独立に、式（３９）を計算することで、各音源mの音源占有度の推定値M^_n,k ^(m)を更新する(E-step1)。
(b)各音源mについて、スペクトル形状推定部１０３−ｍが、式（２８）に基づき、Z^_n,i ^(m)を更新する(E-step2)。
(c)各音源mについて、スペクトル形状推定部１０３−ｍが、式（２７）を最大化するMFCC c_n ^(m)を求めることで、MFCCの推定値c^_n ^(m)を更新する。このとき、式（２７）は、一般には、非線形関数となるため、その最大化は、共役勾配法、準ニュートン法、ニュートン法などの一般的な非線形最適化法により実現される(M-step)。
ステップ４．目的音スペクトル推定部１０６は、式（２１）により、対数パワースペクトル（＝スペクトル特徴量）の推定値s^_n,k ^(m)を求める。

[確認実験]
目的信号スペクトル推定技術を評価する目的で確認実験を行った。実験条件を説明する。残響のある部屋で、二本のマイクロホンを用いて、マイクロホンの正面にいる話者の音声が様々な周囲の背景音と同時に収録された音を、観測信号として用いた。この観測信号を用いて、非特許文献３（従来例）および実施形態５（本技術）に示された目的信号スペクトル推定法の比較実験を行った。従来例と本技術はともに、対数パワースペクトルを観測信号のスペクトル特徴量とし、本技術では、MFCCをスペクトル形状特徴量とし、スペクトル形状モデルとしてガウス混合モデルを採用した。また、従来例と本技術はともに、非特許文献３と同じ音源位置特徴量のモデルを用いて、そのパラメータを事前学習により用意した。

観測信号および目的信号のスペクトルを推定した信号に対して、自動音声認識を適用した結果を示す。観測信号をそのまま音声認識した場合の単語正解率は、69.4%であったのに対し、従来例と本技術で推定した目的信号のスペクトルに対し音声認識を適用した場合の単語正解率は、それぞれ、83.7%と87.2%であった。従来例でも、大幅な音声認識率の改善が得られているが、本技術により、さらに大幅な改善が得られた。

以上の結果より、本技術は、MFCCに関する混合ガウス分布をスペクトル形状モデルとして用いることで、従来例に比べて、大幅に目的信号のスペクトル推定精度を改善できることが確認された。

＜スペクトル形状特徴量推定装置／目的信号スペクトル特徴量推定装置のハードウェア構成例＞
上述の実施形態に関わるスペクトル形状特徴量推定装置／目的信号スペクトル特徴量推定装置（以下、単に推定装置という）は、キーボードなどが接続可能な入力部、液晶ディスプレイなどが接続可能な出力部、ＣＰＵ（Central Processing Unit）〔キャッシュメモリなどを備えていてもよい。〕、メモリであるＲＡＭ（Random Access Memory）やＲＯＭ（Read Only Memory）と、ハードディスクである外部記憶装置、並びにこれらの入力部、出力部、ＣＰＵ、ＲＡＭ、ＲＯＭ、外部記憶装置間のデータのやり取りが可能なように接続するバスなどを備えている。また必要に応じて、推定装置に、ＣＤ−ＲＯＭなどの記憶媒体を読み書きできる装置（ドライブ）などを設けるとしてもよい。このようなハードウェア資源を備えた物理的実体としては、汎用コンピュータなどがある。

推定装置の外部記憶装置には、スペクトル形状特徴量などを推定するためのプログラム並びにこのプログラムの処理において必要となるデータなどが記憶されている〔外部記憶装置に限らず、例えばプログラムを読み出し専用記憶装置であるＲＯＭに記憶させておくなどでもよい。〕。また、これらのプログラムの処理によって得られるデータなどは、ＲＡＭや外部記憶装置などに適宜に記憶される。以下、データやその格納領域のアドレスなどを記憶する記憶装置を単に「記憶部」と呼ぶことにする。

推定装置の記憶部には、観測信号からスペクトル特徴量を抽出するためのプログラム、最適化関数に基づいて最も尤もらしいスペクトル形状特徴量を推定するためのプログラム（具体的には、スペクトル形状特徴量を更新するためのプログラムと、音源占有度を更新するためのプログラム）、目的信号のスペクトル特徴量を推定するためのプログラム、必要に応じて更に、音源位置特徴量を推定するためのプログラム、音源位置状態を推定するためのプログラムが記憶されている。

推定装置では、記憶部に記憶された各プログラムとこの各プログラムの処理に必要なデータが必要に応じてＲＡＭに読み込まれて、ＣＰＵで解釈実行・処理される。この結果、ＣＰＵが所定の機能（特徴抽出部、スペクトル形状推定部、音源占有度更新部、目的音スペクトル推定部、音源位置特徴抽出部、音源位置状態推定部）を実現することでスペクトル推定が実現される。なお、目的音スペクトル推定部は、スペクトル形状特徴量推定装置の必須の構成要素ではない。

＜補記＞
本発明は上述の実施形態に限定されるものではなく、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。また、上記実施形態において説明した処理は、記載の順に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されるとしてもよい。

また、上記実施形態において説明したハードウェアエンティティ（推定装置）における処理機能をコンピュータによって実現する場合、ハードウェアエンティティが有すべき機能の処理内容はプログラムによって記述される。そして、このプログラムをコンピュータで実行することにより、上記ハードウェアエンティティにおける処理機能がコンピュータ上で実現される。

この処理内容を記述したプログラムは、コンピュータで読み取り可能な記録媒体に記録しておくことができる。コンピュータで読み取り可能な記録媒体としては、例えば、磁気記録装置、光ディスク、光磁気記録媒体、半導体メモリ等どのようなものでもよい。具体的には、例えば、磁気記録装置として、ハードディスク装置、フレキシブルディスク、磁気テープ等を、光ディスクとして、ＤＶＤ（Digital Versatile Disc）、ＤＶＤ−ＲＡＭ（Random Access Memory）、ＣＤ−ＲＯＭ（Compact Disc Read Only Memory）、ＣＤ−Ｒ（Recordable）／ＲＷ（ReWritable）等を、光磁気記録媒体として、ＭＯ（Magneto-Optical disc）等を、半導体メモリとしてＥＥＰ−ＲＯＭ（Electronically Erasable and Programmable-Read Only Memory）等を用いることができる。

また、このプログラムの流通は、例えば、そのプログラムを記録したＤＶＤ、ＣＤ−ＲＯＭ等の可搬型記録媒体を販売、譲渡、貸与等することによって行う。さらに、このプログラムをサーバコンピュータの記憶装置に格納しておき、ネットワークを介して、サーバコンピュータから他のコンピュータにそのプログラムを転送することにより、このプログラムを流通させる構成としてもよい。

このようなプログラムを実行するコンピュータは、例えば、まず、可搬型記録媒体に記録されたプログラムもしくはサーバコンピュータから転送されたプログラムを、一旦、自己の記憶装置に格納する。そして、処理の実行時、このコンピュータは、自己の記録媒体に格納されたプログラムを読み取り、読み取ったプログラムに従った処理を実行する。また、このプログラムの別の実行形態として、コンピュータが可搬型記録媒体から直接プログラムを読み取り、そのプログラムに従った処理を実行することとしてもよく、さらに、このコンピュータにサーバコンピュータからプログラムが転送されるたびに、逐次、受け取ったプログラムに従った処理を実行することとしてもよい。また、サーバコンピュータから、このコンピュータへのプログラムの転送は行わず、その実行指示と結果取得のみによって処理機能を実現する、いわゆるＡＳＰ（Application Service Provider）型のサービスによって、上述の処理を実行する構成としてもよい。なお、本形態におけるプログラムには、電子計算機による処理の用に供する情報であってプログラムに準ずるもの（コンピュータに対する直接の指令ではないがコンピュータの処理を規定する性質を有するデータ等）を含むものとする。

また、この形態では、コンピュータ上で所定のプログラムを実行させることにより、ハードウェアエンティティを構成することとしたが、これらの処理内容の少なくとも一部をハードウェア的に実現することとしてもよい。

Claims (11)

1. 複数の音源それぞれからの音響信号が混ざって収音された観測信号から、各音響信号のスペクトル特徴量に対応するスペクトル形状特徴量を推定するスペクトル形状特徴量推定装置であって、
   上記各音源に対応する、スペクトル形状特徴量の事前確率密度関数（スペクトル形状モデル）と、スペクトル形状特徴量が与えられたもとでのスペクトル特徴量の条件付き確率密度関数（スペクトル観測モデル）とを記憶している記憶部と、
   各時間周波数点において最大のエネルギーを持つ音響信号の音源を表す占有的音源番号を潜在変数に持つ最適化関数を、上記スペクトル形状モデルと上記スペクトル観測モデルを用いて最大化し、音源ごとのスペクトル形状特徴量、および、全ての上記音源のスペクトル形状特徴量が与えられたもとで上記各音源が上記占有的音源番号で表される音源である事後確率（音源占有度）、を推定する推定手段とを含み、
   上記最適化関数は、全ての上記音源のスペクトル形状特徴量が与えられたもとでの観測信号のスペクトル特徴量の条件付き確率密度関数と、上記音源ごとに定められるスペクトル形状特徴量の事前確率密度関数の積とで表される
   スペクトル形状特徴量推定装置。
2. 請求項１に記載のスペクトル形状特徴量推定装置であって、
   上記推定手段は、
   音源ごとのスペクトル形状特徴量に基づいて、上記音源占有度を更新する音源占有度更新部と、
   上記音源占有度に基づいて、上記音源ごとのスペクトル形状特徴量を更新するスペクトル形状推定部とを含む
   ことを特徴とするスペクトル形状特徴量推定装置。
3. 請求項２に記載のスペクトル形状特徴量推定装置であって、
   上記スペクトル形状特徴量はメル周波数ケプストラム係数であり、
   上記スペクトル形状モデルはメル周波数ケプストラム係数の混合ガウス分布である
   ことを特徴とするスペクトル形状特徴量推定装置。
4. 請求項２または請求項３に記載のスペクトル形状特徴量推定装置であって、
   複数のマイクロホンによって得られた上記観測信号から上記各音源の位置に関する特徴量である音源位置特徴量を抽出する音源位置特徴抽出部と、
   上記音源ごとに、上記音源位置特徴量に関する確率密度関数で表される音源位置状態の推定値を求める音源位置状態推定部とを含み、
   上記音源占有度更新部は、上記各音源の上記音源位置状態の推定値をも用いて上記音源占有度を更新する
   ことを特徴とするスペクトル形状特徴量推定装置。
5. 請求項１から請求項４のいずれかに記載のスペクトル形状特徴量推定装置によって推定された、音源ごとのスペクトル形状特徴量および音源占有度のうち、目的信号の音源に対応するスペクトル形状特徴量と、目的信号の音源の音源占有度と、目的信号の音源に対応するスペクトル観測モデルと、観測信号のスペクトル特徴量とから、目的信号のスペクトル特徴量を推定する目的音スペクトル推定部を含む
   目的信号スペクトル特徴量推定装置。
6. 複数の音源それぞれからの音響信号が混ざって収音された観測信号から、各音響信号のスペクトル特徴量に対応するスペクトル形状特徴量を推定するスペクトル形状特徴量推定方法であって、
   記憶部には、上記各音源に対応する、スペクトル形状特徴量の事前確率密度関数（スペクトル形状モデル）と、スペクトル形状特徴量が与えられたもとでのスペクトル特徴量の条件付き確率密度関数（スペクトル観測モデル）とが記憶されており、
   各時間周波数点において最大のエネルギーを持つ音響信号の音源を表す占有的音源番号を潜在変数に持つ最適化関数を、上記スペクトル形状モデルと上記スペクトル観測モデルを用いて最大化し、音源ごとのスペクトル形状特徴量、および、全ての上記音源のスペクトル形状特徴量が与えられたもとで上記各音源が上記占有的音源番号で表される音源である事後確率（音源占有度）、を推定する推定ステップを含み、
   上記最適化関数は、全ての上記音源のスペクトル形状特徴量が与えられたもとでの観測信号のスペクトル特徴量の条件付き確率密度関数と、上記音源ごとに定められるスペクトル形状特徴量の事前確率密度関数の積とで表される
   スペクトル形状特徴量推定方法。
7. 請求項６に記載のスペクトル形状特徴量推定方法であって、
   上記推定ステップは、
   音源ごとのスペクトル形状特徴量に基づいて、上記音源占有度を更新する音源占有度更新ステップと、
   上記音源占有度に基づいて、上記音源ごとのスペクトル形状特徴量を更新するスペクトル形状推定ステップとを含む
   ことを特徴とするスペクトル形状特徴量推定方法。
8. 請求項７に記載のスペクトル形状特徴量推定方法であって、
   上記スペクトル形状特徴量はメル周波数ケプストラム係数であり、
   上記スペクトル形状モデルはメル周波数ケプストラム係数の混合ガウス分布である
   ことを特徴とするスペクトル形状特徴量推定方法。
9. 請求項７または請求項８に記載のスペクトル形状特徴量推定方法であって、
   複数のマイクロホンによって得られた上記観測信号から上記各音源の位置に関する特徴量である音源位置特徴量を抽出する音源位置特徴抽出ステップと、
   上記音源ごとに、上記音源位置特徴量に関する確率密度関数で表される音源位置状態の推定値を求める音源位置状態推定ステップとを含み、
   上記音源占有度更新ステップは、上記各音源の上記音源位置状態の推定値をも用いて上記音源占有度を更新する
   ことを特徴とするスペクトル形状特徴量推定方法。
10. 請求項６から請求項９のいずれかに記載のスペクトル形状特徴量推定方法によって推定された、音源ごとのスペクトル形状特徴量および音源占有度のうち、目的信号の音源に対応するスペクトル形状特徴量と、目的信号の音源の音源占有度と、目的信号の音源に対応するスペクトル観測モデルと、観測信号のスペクトル特徴量とから、目的信号のスペクトル特徴量を推定する目的音スペクトル推定ステップを含む
    目的信号スペクトル特徴量推定方法。
11. コンピュータを、請求項１から請求項４のいずれかに記載のスペクトル形状特徴量推定装置として、あるいは、請求項５に記載の目的信号スペクトル特徴量推定装置として、機能させるためのプログラム。

Images