(57)【特許請求の範囲】
【請求項１】
各短時間フレームｎにおける周波数信号ｘｎのスペクトル値σｎを推定するスペクトル推定装置であって、
前記周波数信号ｘｎの対数スペクトルρｎの状態を表す状態パラメータθｎの事前確率密度関数ｐ（θｎ；Θθ）に関するモデルパラメータであるスペクトル状態モデルΘθと、前記状態パラメータθｎが既知の条件下での前記対数スペクトルρｎの条件付き確率密度関数ｐ（ρｎ｜θｎ；Θρ）に関するモデルパラメータである状態依存スペクトルモデルΘρとを記憶する記憶部と、
前記対数スペクトルρｎの推定値ρ＾ｎ、前記スペクトル状態モデルΘθ及び前記状態依存スペクトルモデルΘρを用いて、対数尤度重みｗθｎを推定するスペクトル状態推定部と、
前記周波数信号ｘｎ、前記対数尤度重みｗθｎ及び前記状態依存スペクトルモデルΘρを用いて、目的関数を最大化する対数スペクトルを前記推定値ρ＾ｎとして求める事後確率最大化スペクトル推定部とを含み、
収束条件を満たすまで、前記スペクトル状態推定部及び事後確率最大化スペクトル推定部における処理を繰り返す、
スペクトル推定装置。
【請求項２】
請求項１記載のスペクトル推定装置であって、
前記事後確率最大化スペクトル推定部が、一つのスカラー変数ｚとそのスカラー変数に関する指数関数ｅｘｐ（ｚ）と一つのスカラー定数ａとの和によって規定される非線形方程式に関して、各短時間フレームｎにおける周波数ｋ毎の周波数信号ｘｎ，ｋと前記対数尤度重みｗθｎと前記状態依存スペクトルモデルΘρに依存して前記スカラー定数ａを定めるとともに、前記非線形方程式が０に一致する前記スカラー変数ｚの値を求め、その求めた前記スカラー変数ｚと前記周波数信号ｘｎ，ｋと前記対数尤度重みｗθｎと前記状態依存スペクトルモデルΘρとに基づき、前記推定値ρ＾ｎを更新する、
スペクトル推定装置。
【請求項３】
請求項１または請求項２記載のスペクトル推定装置であって、
前記対数スペクトルρｎが取りうる値を規定する事前確率密度関数ｐ（ρｎ；Θθ，Θρ）が混合ガウス分布に従い、前記状態パラメータθｎは短時間フレームｎにおいてＮθ個の有限状態の何れかの状態をとり、前記スペクトル状態モデルΘθは全ての状態ｉに関する混合比βｉからなり、前記状態依存スペクトルモデルΘρは全ての状態ｉに関する全ての周波数ｋに関する平均μｉｋと分散ξｉｋとからなるものとし、
前記スペクトル状態推定部は、前記対数尤度重みｗｉを
【数３５】
0005689844.tif 000036
または、
【数３６】
0005689844.tif 000037
として推定する、
スペクトル推定装置。
【請求項４】
請求項１または請求項２記載のスペクトル推定装置であって、
前記状態パラメータθｎを前記周波数信号に対応するメル周波数ケプストラム係数ｃｎとし、前記状態パラメータの事前確率密度関数ｐ（ｃｎ；Θθ）として前記メル周波数ケプストラム係数ｃｎの混合ガウス分布を用い、ｊをガウス分布の番号とし、前記スペクトル状態モデルΘθを全てのｊに関する混合比γｊと平均μｊと共分散行列Σｊとの集合とし、前記メル周波数ケプストラム係数ｃｎから前記推定値ρ＾ｎへの擬似逆変換を規定する行列Ａ及びベクトルｂと、逆変換誤差ｅがガウス分布に従うと仮定したときの共分散行列Ξとを前記状態依存スペクトルモデルΘρとし、
前記スペクトル状態推定部は、
前記メル周波数ケプストラム係数である状態パラメータの推定値ｃ＾ｎと前記混合比γｊと前記平均μｊと前記共分散行列Σｊとを用いて、期待値Ｅ｛ｊ｜ｃ＾ｎ｝を
【数３７】
0005689844.tif 000038
として算出する期待値算出部と、
前記推定値ρ＾ｎと前記期待値Ｅ｛ｊ｜ｃ＾ｎ｝と前記行列Ａと前記ベクトルｂと前記共分散行列Ξと前記平均μｊと前記共分散行列Σｊとを用いて、前記状態パラメータである前記メル周波数ケプストラム係数の推定値ｃ＾ｎを
【数３８】
0005689844.tif 000039
として更新する状態パラメータ更新部と、
δをディラックデルタ関数とし、前記メル周波数ケプストラム係数である状態パラメータの推定値ｃ＾ｎを用いて、前記メル周波数ケプストラム係数である前記状態パラメータの推定値ｃ＾ｎに対する対数尤度重みｗｃｎを
【数３９】
0005689844.tif 000040
として算出する対数尤度重み算出部とを含み、
収束条件を満たすまで、前記期待値算出部及び前記状態パラメータ更新部における処理を繰り返す、
スペクトル推定装置。
【請求項５】
各短時間フレームｎにおける周波数信号ｘｎのスペクトル値σｎを推定するスペクトル推定方法であって、
前記周波数信号ｘｎの対数スペクトルρｎの状態を表す状態パラメータθｎの事前確率密度関数ｐ（θｎ；Θθ）に関するモデルパラメータであるスペクトル状態モデルΘθと、前記状態パラメータθｎが既知の条件下での前記対数スペクトルρｎの条件付き確率密度関数ｐ（ρｎ｜θｎ；Θρ）に関するモデルパラメータである状態依存スペクトルモデルΘρとを記憶しておき、
前記対数スペクトルρｎの推定値ρ＾ｎ、前記スペクトル状態モデルΘθ及び前記状態依存スペクトルモデルΘρを用いて、対数尤度重みｗθｎを推定するスペクトル状態推定ステップと、
前記周波数信号ｘｎ、前記対数尤度重みｗθｎ及び前記状態依存スペクトルモデルΘρを用いて、目的関数を最大化する対数スペクトルを前記推定値ρ^ｎとして求める事後確率最大化スペクトル推定ステップとを含み、
収束条件を満たすまで、前記スペクトル状態推定ステップ及び事後確率最大化スペクトル推定ステップにおける処理を繰り返す、
スペクトル推定方法。
【請求項６】
請求項５記載のスペクトル推定方法であって、
前記事後確率最大化スペクトル推定ステップが、一つのスカラー変数ｚとそのスカラー変数に関する指数関数ｅｘｐ（ｚ）と一つのスカラー定数ａとの和によって規定される非線形方程式に関して、各短時間フレームｎにおける周波数ｋ毎の周波数信号ｘｎ，ｋと前記対数尤度重みｗθｎと前記状態依存スペクトルモデルΘρに依存して前記スカラー定数ａを定めるとともに、前記非線形方程式が０に一致する前記スカラー変数ｚの値を求め、その求めた前記スカラー変数ｚと前記周波数信号ｘｎ，ｋと前記対数尤度重みｗθｎと前記状態依存スペクトルモデルΘρとに基づき、前記推定値ρ＾ｎを更新する、
スペクトル推定方法。
【請求項７】
請求項５または請求項６記載のスペクトル推定方法であって、
前記対数スペクトルρｎが取りうる値を規定する事前確率密度関数ｐ（ρｎ；Θθ，Θρ）が混合ガウス分布に従い、前記状態パラメータθｎは短時間フレームｎにおいてＮθ個の有限状態の何れかの状態をとり、前記スペクトル状態モデルΘθは全ての状態ｉに関する混合比βｉからなり、前記状態依存スペクトルモデルΘρは全ての状態ｉに関する全ての周波数ｋに関する平均μｉｋと分散ξｉｋとからなるものとし、
前記スペクトル状態推定ステップにおいて、前記対数尤度重みｗｉを
【数４０】
0005689844.tif 000041
または、
【数４１】
0005689844.tif 000042
として推定する、
スペクトル推定方法。
【請求項８】
請求項５または請求項６記載のスペクトル推定方法であって、
前記状態パラメータθｎを前記周波数信号に対応するメル周波数ケプストラム係数ｃｎとし、前記状態パラメータの事前確率密度関数ｐ（ｃｎ；Θθ）として前記メル周波数ケプストラム係数ｃｎの混合ガウス分布を用い、ｊをガウス分布の番号とし、前記スペクトル状態モデルΘθを全てのｊに関する混合比γｊと平均μｊと共分散行列Σｊとの集合とし、前記メル周波数ケプストラム係数ｃｎから前記推定値ρ＾ｎへの擬似逆変換を規定する行列Ａ及びベクトルｂと、逆変換誤差ｅがガウス分布に従うと仮定したときの共分散行列Ξとを前記状態依存スペクトルモデルΘρとし、
前記スペクトル状態推定ステップは、
前記メル周波数ケプストラム係数である状態パラメータの推定値ｃ＾ｎと前記混合比γｊと前記平均μｊと前記共分散行列Σｊとを用いて、期待値Ｅ｛ｊ｜ｃ＾ｎ｝を
【数４２】
0005689844.tif 000043
として算出する期待値算出ステップと、
前記推定値ρ＾ｎと前記期待値Ｅ｛ｊ｜ｃ＾ｎ｝と前記行列Ａと前記ベクトルｂと前記共分散行列Ξと前記平均μｊと前記共分散行列Σｊとを用いて、前記状態パラメータである前記メル周波数ケプストラム係数の推定値ｃ＾ｎを
【数４３】
0005689844.tif 000044
として更新する状態パラメータ更新ステップと、
δをディラックデルタ関数とし、前記メル周波数ケプストラム係数である状態パラメータの推定値ｃ＾ｎを用いて、前記メル周波数ケプストラム係数である前記状態パラメータの推定値ｃ＾ｎに対する対数尤度重みｗｃｎを
【数４４】
0005689844.tif 000045
として算出する対数尤度重み算出ステップとを含み、
収束条件を満たすまで、前記期待値算出ステップ及び前記状態パラメータ更新ステップにおける処理を繰り返す、
スペクトル推定方法。
【請求項９】
請求項１から請求項４の何れかに記載のスペクトル推定装置としてコンピュータを機能させるためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は一次元時系列信号を周波数分割した信号から、信号のスペクトルを推定するスペクトル推定技術に関する。
【背景技術】
【０００２】
以下の説明に用いる図面では、同じ機能を持つ構成部や同じ処理を行うステップには同一の符号を記し、重複説明を省略する。以下の説明において、テキスト中で使用する記号、「＾」等は、本来直前の文字の真上に記載されるべきものであるが、テキスト記法の制限により、当該文字の直後に記載する。式中においてはこれらの記号は本来の位置に記述している。また、ベクトルや行列の各要素単位で行われる処理は、特に断りが無い限り、そのベクトルやその行列の全ての要素に対して適用されるものとする。
【０００３】
ｎを短時間フレームの番号、ｋ（＝１〜Ｎｋ）を、観測信号を周波数分割する際の周波数の番号とし、短時間フレームにおける各周波数分割した信号をｘｎ，ｋと表す。さらに、ｘｎ，ｋを全ての周波数についてひとまとめにしてできるベクトルをｘｎ＝［ｘｎ，１，ｘｎ，２，…，ｘｎ，Ｎｋ］Ｔと表記し（ただし、下付添え字ＮｋはＮｋを表す）、以下では、短時間フレームｎの周波数信号と呼ぶ。Ｎｋは周波数分割数を表す。Ｔは、ベクトルや行列の非共役転置を表す。
【０００４】
図１は、非特許文献１などに開示されている従来のスペクトル推定装置９の機能ブロック図を示す。スペクトル推定装置９は、各短時間フレームｎにおいて、周波数信号ｘｎを受け取り、最尤法に基づき、周波数信号ｘｎのスペクトルσｎ＝［σｎ，１，σｎ，２，…，σｎ，Ｎｋ］Ｔを推定する。より具体的には、非特許文献１では、信号の周波数分割に短時間フーリエ変換を用いており、残響除去された信号の短時間フーリエ変換の推定値がｘｎと与えられているときに、ｘｎ，ｋが平均０、分散σｎ，ｋの複素正規分布に従うとの仮定の下で、最尤スペクトル推定部９１において、分散σｎ，ｋ（＝スペクトルの値）を最尤法により求める。つまり、ｘｎ，ｋの条件付き確率密度関数ｐ（ｘｎ，ｋ｜σｎ，ｋ）は、以下の式でモデル化される。
【０００５】
【数１】
0005689844.tif 000002
【０００６】
そして、対数尤度関数Ｌ（σｎ）＝Σｋｌｏｇｐ（ｘｎ，ｋ｜σｎ，ｋ）を最大にする値として、以下のようにスペクトルσｎ＝［σｎ，１，σｎ，２，…，σｎ，Ｎｋ］Ｔを推定する。
【０００７】
【数２】
0005689844.tif 000003
【０００８】
なお、推定値と推定すべき変数を区別するために、推定値には＾をつけて、σ＾ｎ等と表記することにする。
【０００９】
一方、非特許文献２等に詳述されているように、（１）式に加えて、分散σｎ，ｋのとりうる値を規定する事前確率密度関数ｐ（σｎ；Θ）を導入し、σｎの値を、周波数信号ｘｎが与えられた下での事後確率最大化（Maximum a posteriori、以下「ＭＡＰ」ともいう）推定により求める方法が説明されている。ここで、Θは、事前確率密度関数のモデルパラメータである。この場合のスペクトル推定装置８の機能ブロック図を図２に示す。ＭＡＰ推定は、以下のように定義される。
【００１０】
【数３】
0005689844.tif 000004
【００１１】
事後確率最大化スペクトル推定部８１は、スペクトル事前分布記憶部８２から事前確率密度関数のモデルパラメータΘを取り出し、（４）式により、σｎを求める。このように、σｎの事前確率密度関数ｐ（σｎ；Θ）を考慮することで、σｎがとりうる値の傾向をある程度制限できることになる。事前確率密度関数ｐ（σｎ；Θ）として、ガウス分布の分散に関する自然共役分布である逆ガンマ分布等を用いると、効率的な計算が可能なことが知られている。
【先行技術文献】
【非特許文献】
【００１２】
【非特許文献１】中谷智広、吉岡拓也、木下慶介、三好正人、Biing-Hwang Juang、“短時間フーリエ変換表現を用いた最尤推定に基づく音声信号の残響除去”、日本音響学会春季研究発表会、２００８年３月、pp.733-736
【非特許文献２】C. M. Bishop著、元田浩、栗田多喜夫、樋口知之、松本裕治訳、「パターン認識と機械学習上- ベイズ理論による統計的予測」、シュプリンガー・ジャパン、2007年、pp.95-100
【発明の概要】
【発明が解決しようとする課題】
【００１３】
非特許文献１では、周波数信号として短時間フーリエ変換の推定値を用いているが、一般に、推定値には必ず推定誤差が含まれる。また、周波数信号としてマイクロホンで収音した観測信号を用いる場合なども、一般に、観測信号には必ず何らかの雑音が含まれる。その結果、これらの周波数信号に基づき、従来の方法でスペクトル推定を行うと、誤差や雑音の影響で必ずしも精度よく推定が行えないという問題がある。特に、非特許文献１のように、短時間フーリエ変換の推定とスペクトルの推定を相互に依存させながら交互に繰り返すような場合、繰り返しにより誤差の影響が拡大して推定値が劣化する場合もある。
【００１４】
一方、非特許文献２にあるように、分散の事前確率密度関数ｐ（σｎ；Θ）を導入し、分散の値をＭＡＰ推定によって求めるようにすることで、分散がとりうる値を制限し、ある程度、誤差の影響を弱めることができる。しかし、ＭＡＰ推定において効率的に最適化が行えるのは、自然共役分布のようなごく一部の事前確率密度関数ｐ（σｎ；Θ）を用いる場合に限られ、σｎの分布を精度よくあらわすものを必ずしも利用できないという問題がある。特に、自動音声認識システムの音響モデルとして利用される確率分布と類似性の高い対数スペクトルに関する混合ガウス分布等（ガウス分布、混合ガウス分布、ガウス分布を出力確率分布に持つ隠れマルコフモデル等を含む）は、精度よく音声信号のスペクトルの分布を表すと考えられているが、（４）式の事前確率密度関数ｐ（σｎ；Θ）として用いた場合に、効率的に最適化を行う方法は知られていない。
【００１５】
この発明はこの課題に鑑みてなされたものであり、周波数信号が誤差を含む場合でも、対数スペクトルに関する混合ガウス分布等をスペクトルの事前確率密度関数として用いて、高精度かつ効率的にスペクトル推定が行える技術を提供することを目的とする。
【課題を解決するための手段】
【００１６】
上記の課題を解決するために、本発明の第一の態様によれば、スペクトル推定装置は、各短時間フレームｎにおける周波数信号ｘｎのスペクトル値σｎを推定する。スペクトル推定装置は、記憶部、スペクトル状態推定部及び事後確率最大化スペクトル推定部を含む。記憶部は、周波数信号ｘｎの対数スペクトルρｎの状態を表す状態パラメータθｎの事前確率密度関数ｐ（θｎ；Θθ）に関するモデルパラメータであるスペクトル状態モデルΘθと、状態パラメータθｎが既知の条件下での対数スペクトルρｎの条件付き確率密度関数ｐ（ρｎ｜θｎ；Θρ）に関するモデルパラメータである状態依存スペクトルモデルΘρとを記憶する。スペクトル状態推定部は、対数スペクトルの推定値ρ＾ｎ、スペクトル状態モデルΘθ及び状態依存スペクトルモデルΘρを用いて、対数尤度重みｗθｎを推定する。事後確率最大化スペクトル推定部は、周波数信号ｘｎ、対数尤度重みｗθｎ及び状態依存スペクトルモデルΘρを用いて、目的関数を最大化する対数スペクトルρｎを推定する。収束条件を満たすまで、スペクトル状態推定部及び事後確率最大化スペクトル推定部における処理を繰り返す。
【００１７】
上記の課題を解決するために、本発明の第二の態様によれば、スペクトル推定方法は、各短時間フレームｎにおける周波数信号ｘｎのスペクトル値σｎを推定する。スペクトル推定方法は、スペクトル状態推定ステップ及び事後確率最大化スペクトル推定ステップを含む。周波数信号ｘｎの対数スペクトルρｎの状態を表す状態パラメータθｎの事前確率密度関数ｐ（θｎ；Θθ）に関するモデルパラメータであるスペクトル状態モデルΘθと、状態パラメータθｎが既知の条件下での対数スペクトルρｎの条件付き確率密度関数ｐ（ρｎ｜θｎ；Θρ）に関するモデルパラメータである状態依存スペクトルモデルΘρとを記憶しておく。スペクトル状態推定ステップは、対数スペクトルの推定値ρ＾ｎ、スペクトル状態モデルΘθ及び状態依存スペクトルモデルΘρを用いて、対数尤度重みｗθｎを推定する。事後確率最大化スペクトル推定ステップは、周波数信号ｘｎ、対数尤度重みｗθｎ及び状態依存スペクトルモデルΘρを用いて、目的関数を最大化する対数スペクトルρｎを推定する。収束条件を満たすまで、スペクトル状態推定ステップ及び事後確率最大化スペクトル推定ステップにおける処理を繰り返す。
【発明の効果】
【００１８】
本発明によれば、スペクトルの分布を高精度に表現可能な対数スペクトルに関する潜在変数依存型ガウス分布をスペクトルの事前確率密度関数として用いた場合でも、効率的にスペクトルの値を推定できる。その結果、周波数信号が誤差を含むような場合でも、効率的かつ高精度に、そのスペクトルの推定が可能になるという効果を奏する。
【図面の簡単な説明】
【００１９】
【図１】従来のスペクトル推定装置の機能ブロック図。
【図２】従来のスペクトル推定装置の機能ブロック図。
【図３】第一実施形態のスペクトル推定装置の機能ブロック図。
【図４】第一実施形態のスペクトル推定装置の処理フローを示す図。
【図５】第一実施形態の事後確率最大化スペクトル推定部の機能ブロック図。
【図６】第一実施形態のスペクトル状態推定部と事後確率最大化スペクトル推定部の処理フローを示す図。
【図７】第一実施形態の変形例のスペクトル状態推定部の機能ブロック図。
【図８】第一実施形態の変形例のスペクトル状態推定部の処理フローを示す図。
【図９】第二実施形態のスペクトル状態推定部の機能ブロック図。
【図１０】第二実施形態のスペクトル状態推定部の処理フローを示す図。
【図１１】従来技術（最尤法）、第一実施形態、第二実施形態を用いた場合の比較結果を示す図。
【発明を実施するための形態】
【００２０】
以下、本発明の実施形態について説明する。
【００２１】
＜第一実施形態＞
図３はスペクトル推定装置１０の機能ブロック図を、図４はその処理フローを示す。スペクトル推定装置１０は、スペクトル状態モデル記憶部１０１、状態依存スペクトルモデル記憶部１０２、スペクトル状態推定部１０４及び事後確率最大化スペクトル推定部１０６を含む。
【００２２】
スペクトル推定装置１０は、各短時間フレームｎにおいて、周波数信号ｘｎを受け取り、そのスペクトルの推定値σ＾ｎを出力する。
【００２３】
まず、周波数信号ｘｎの対数スペクトルをρｎ＝［ρｎ，１，ρｎ，２，…，ρｎ，Ｎｋ］Ｔと表すことにする。ただし、ρｎ，ｋ＝ｌｏｇσｎ，ｋである。
【００２４】
スペクトル状態モデル記憶部１０１は、周波数信号ｘｎの対数スペクトルρｎの状態を表す状態パラメータθｎの事前確率密度関数ｐ（θｎ；Θθ）に関するモデルパラメータを記憶している。以下、このモデルパラメータをスペクトル状態モデルΘθと呼ぶ。
【００２５】
状態依存スペクトルモデル記憶部１０２は、状態パラメータθｎが既知の条件下での対数スペクトルρｎの条件付き確率密度関数ｐ（ρｎ｜θｎ；Θρ）に関するモデルパラメータを記憶している。以下、このモデルパラメータを状態依存スペクトルモデルΘρと呼ぶ。
【００２６】
スペクトル状態推定部１０４は、後述する事後確率最大化スペクトル推定部１０６が推定した対数スペクトルの推定値ρ＾ｎを受け取るとともに、スペクトル状態モデル記憶部１０１と状態依存スペクトルモデル記憶部１０２のそれぞれからスペクトル状態モデルΘθと状態依存スペクトルモデルΘρを受け取り、対数尤度重みｗθｎ（ただし、下付添え字θｎはθｎを表す）を推定し（ｓ１）、出力する。
【００２７】
事後確率最大化スペクトル推定部１０６は、周波数信号ｘｎと、対数尤度重みｗθｎと、状態依存スペクトルモデルΘρを受け取り、後述する目的関数を最大化する対数スペクトルの推定値ρ＾ｎ＝［ρ＾ｎ，１，ρ＾ｎ，２，…，ρ＾ｎ，Ｎｋ］Ｔを推定し（ｓ２）、出力する。また収束条件を満たすまで（ｓ３）、スペクトル状態推定部１０４及び事後確率最大化スペクトル推定部１０６における処理（ｓ１及びｓ２）を繰り返す。収束条件としては、例えば、（１）繰り返し回数が所定の回数を超えることや、（２）一つ前の繰り返し時に得られた対数スペクトルの推定値と現在の繰り返し時に得られた対数スペクトルの推定値との差分が閾値以下であること等が挙げられる。収束条件を満たした場合は、満たした時点の対数スペクトルの推定値ρ＾ｎからスペクトルの推定値σ＾ｎ＝［σ＾ｎ，１，σ＾ｎ，２，…，σ＾ｎ，Ｎｋ］Ｔを求め、出力する。ただし、σ＾ｎ，ｋ＝ｅｘｐ（ρ＾ｎ，ｋ）である。推定値ρ＾ｎ，ｋが得られれば推定値σ＾ｎ，ｋも与えられるので、以下では推定値ρ＾ｎ，ｋについての推定方法のみについて記述する。
【００２８】
＜第一実施形態のポイント＞
スペクトル推定装置１０では、対数スペクトルρｎがとりうる値を規定する事前確率密度関数ｐ（ρｎ；Θθ，Θρ）を導入し、対数スペクトルρｎの値を、周波数信号ｘｎが与えられた下での事後確率最大化（ＭＡＰ）推定により求める。すなわち、以下のように求める。
【００２９】
【数４】
0005689844.tif 000005
【００３０】
これにより、（１）式で定義される周波数信号ｘｎの条件付き確率密度関数ｐ（ｘｎ，ｋ｜σｎ，ｋ）に加えて、対数スペクトルρｎの事前確率密度関数ｐ（ρｎ；Θθ，Θρ）をも考慮しながら対数スペクトルρｎが推定されることになる。そのため、周波数信号ｘｎに含まれる誤差の影響を比較的受けにくいスペクトル推定が可能になる。なお、（５）式のｐ（ｘｎ｜ρｎ）は、従来の最尤法と同様に、ｐ（ｘｎ｜ρｎ）＝Πｋｐ（ｘｎ，ｋ｜ρｎ，ｋ）のように分解でき、（１）式とσｎ，ｋ＝ｅｘｐ（ρｎ，ｋ）の関係式に基づき、以下のように定義されているものとする。
【００３１】
【数５】
0005689844.tif 000006
【００３２】
さらに、第一実施形態のスペクトル推定装置１０では、高精度で効率的な推定を実現するために、以下の３つの仮定を導入する。
【００３３】
仮定（１）：周波数信号ｘｎの対数スペクトルρｎの事前確率密度関数ｐ（ρｎ；Θθ，Θρ）は、状態パラメータθｎを潜在変数として持つ以下の式でモデル化されている。
【００３４】
【数６】
0005689844.tif 000007
【００３５】
なお、上式では、状態パラメータθｎは離散値を取るものと仮定し、その周辺化のために全状態の総和をとっている。一方、本発明は、状態パラメータθｎが連続値を取る場合も含む。その場合、状態パラメータθｎの周辺化は、以下のように、状態パラメータθｎがとりうる値の全範囲にわたる積分として定義される。
【００３６】
【数７】
0005689844.tif 000008
【００３７】
本実施形態では、状態パラメータθｎは離散値を取るものとして説明する。なお、連続値を取る場合については、状態パラメータθｎに関する総和の部分を、適宜、状態パラメータθｎがとりうる値の全範囲にわたる積分として読み替えるだけでよいので、個別の説明は省略する。
【００３８】
仮定（２）：状態パラメータθｎが与えられた下での対数スペクトルρｎの条件付き確率密度関数ｐ（ρｎ｜θｎ；Θρ）は、多変量ガウス分布に従う。以下、仮定（１）及び仮定（２）に従う分布を潜在変数依存型ガウス分布と呼ぶ。
【００３９】
仮定（３）：さらに、条件付き確率密度関数ｐ（ρｎ｜θｎ；Θρ）は、各周波数ｋの対数スペクトルρｎ，ｋに関する条件付き確率密度関数ｐ（ρｎ，ｋ｜θｎ；Θρ）の積に分解できる。
【００４０】
【数８】
0005689844.tif 000009
【００４１】
なお、仮定（３）を満たすとき、条件付き確率密度関数ｐ（ρｎ｜θｎ；Θρ）は周波数分解可能であるという。なお、仮定（２）により、上式の右辺はさらに以下のように書き換えられる。
【００４２】
【数９】
0005689844.tif 000010
【００４３】
ここで、Ｎ（ｘ；μ，ξ）は、平均μ、分散ξの一次元ガウス分布の確率密度関数を表す。例えば、（８）式において、状態パラメータθｎが単一の状態しかとらないとすると事前確率密度関数ｐ（ρｎ；Θθ，Θρ）はガウス分布に一致する。状態パラメータθｎが有限個の状態のどれか一つを取ると仮定すると事前確率密度関数ｐ（ρｎ；Θθ，Θρ）は混合ガウス分布に一致する。さらに、隣り合う短時間フレーム間での状態パラメータθｎの遷移が、ある状態遷移確率に従うと仮定すると、対数スペクトルρｎに関する隠れマルコフモデルになる。上記の仮定および以下では、簡単のため、対数スペクトルρｎの事前確率密度関数ｐ（ρｎ；Θθ，Θρ）は、短時間フレームｎ毎に独立な分布として説明する。なお、本実施形態において、短時間フレーム間の状態遷移過程を導入する方法は、隠れマルコフモデルに関する既知の技術に基づき自明であるので、その説明を省略する。
【００４４】
（５）式の解は、状態パラメータθｎを隠れ変数とした期待値最大化（Expextation Maximization：以下「ＥＭ」とする）アルゴリズム（及び、その関連最適化手法）で求めることができる。このとき、補助関数Ｑ（ρｎ｜ρ＾ｎ）は、以下のように定義される。
【００４５】
【数１０】
0005689844.tif 000011
【００４６】
ここで、対数スペクトルρｎが既知の下で、周波数信号ｘｎは、状態パラメータθｎと独立であると仮定すると、上記右辺に含まれる完全データの確率密度関数ｐ（ｘｎ，ρｎ，θｎ；Θθ，Θρ）は、以下のように展開できる。
【００４７】
【数１１】
0005689844.tif 000012
【００４８】
したがって、ρｎと無関係の項を省略して（１１）式をさらに展開し、以下を得る。
【００４９】
【数１２】
0005689844.tif 000013
【００５０】
ただし、
【００５１】
【数１３】
0005689844.tif 000014
【００５２】
【数１４】
0005689844.tif 000015
【００５３】
したがって、ＥＭアルゴリズムでは、収束するまで、以下の二つの処理を交互に繰り返すことで、ＭＡＰ推定は実現される。
１．Ｅ−ｓｔｅｐ：スペクトル状態推定部１０４が、対数尤度重みｗθｎを（１８）式に従い更新する（ｓ１）。
２．Ｍ−ｓｔｅｐ：事後確率最大化スペクトル推定部１０６が、（１５）式を最大化するρｎ，ｋを対数スペクトルの推定値ρ＾ｎ，ｋとして更新する（ｓ２）。
【００５４】
なお、（１８）式は、状態パラメータθｎが連続値を取る場合は、状態パラメータθｎに関する連続関数になる。上記の繰り返しのうち、最も計算コストを増大させる可能性があるのは、補助関数Ｑ（ρｎ｜ρ＾ｎ）の値を最大化する対数スペクトルρｎを求めるＭ−ｓｔｅｐである。これに対し、本実施形態では、上記の仮定（１）〜（３）により、すなわち、対数スペクトルρｎの事前確率密度関数ｐ（ρｎ；Θθ，Θρ）が潜在変数依存型ガウス分布に従い、その条件付き確率密度関数ｐ（ρｎ｜θｎ；Θρ）が周波数分解可能であるとき、計算コストを抑えた処理が可能になる。より具体的には、以下の二つのポイントにより、計算コストを抑えることができる。
【００５５】
ポイント（１）：（１５）式は、各時間周波数に閉じて、対数スペクトルρｎ，ｋに関するスカラー１変数関数になっている。すなわち、ρｎの更新は、各時間周波数ｎ，ｋにおける対数スペクトルρｎ，ｋの更新に分解できる。
【００５６】
ポイント（２）：さらに、各時間周波数ｎ，ｋにおいて最大化をするべき関数であるＱｋ（ρｎ，ｋ｜ρ＾ｎ，ｋ）をρｎ，ｋで微分して得られる関数は、以下のような単純な形式をしている。
f(z)=exp(z)+z+a    (19)
【００５７】
【数１５】
0005689844.tif 000016
【００５８】
【数１６】
0005689844.tif 000017
【００５９】
したがって、（１５）式を最大にするρｎ，ｋは、（１９）式においてｆ（ｚ）＝０となるｚを求めた後に、（２０）式から求めることができる。一方、（１９）式は、スカラー定数ａのみで形状が定まる１変数凸関数であり、効率的にｆ（ｚ）＝０の解を求める方法が存在する。例えば、ａの値毎にｆ（ｚ）＝０を与える解をあらかじめ求めておき、解の参照表を用意しておけば、参照表を見るだけで近似解を得ることができる。また、（１９）式を詳しく調べると、ａ＞−１／２でｆ（ｚ）≒ｅｘｐ（ｚ）＋ａ，ａ≦−１／２でｆ（ｚ）≒ｚ＋ａと荒く近似できることがわかる。これより、以下の近似解を得ることもできる。
【００６０】
【数１７】
0005689844.tif 000018
【００６１】
さらに、ｆ（ｚ）＝０を与える解の初期推定値としてこれらの近似解を用い、ニュートン法などの勾配法を用いて数値的な探索を行うことで、解の精度を上げることができる。しかも、このとき、ｆ（ｚ）は、１変数凸関数であるため、非常に効率的かつ効果的に勾配法による探索を実現できる。
【００６２】
（対数スペクトルの事前確率密度関数ｐ（ρｎ；Θθ，Θρ））
本実施形態では、対数スペクトルの事前確率密度関数ｐ（ρｎ；Θθ，Θρ）を、混合ガウス分布でモデル化する。状態パラメータθｎは、各短時間フレームｎにおいて、１からＮθで番号付されたＮθ個の有限状態の何れかの状態ｉをとるとする。事前確率密度関数ｐ（ρｎ；Θθ，Θρ）は、以下で定義される。
【００６３】
【数１８】
0005689844.tif 000019
【００６４】
ただし、
【００６５】
【数１９】
0005689844.tif 000020
【００６６】
p(θn=i;Θθ)=βi    (25)
スペクトル状態モデルΘθは、全ての状態ｉに関する混合比βｉからなり、状態依存スペクトルモデルΘρは、全ての状態ｉ、全ての周波数ｋに関する平均μｉｋと共分散行列ξｉｋとからなる。これらのモデルは、スペクトル推定の対象となる信号に関する学習データを用いて、事前に学習されているとする。混合ガウス分布のモデルパラメータの学習には、ＥＭアルゴリズムを用いる方法などが知られている。
【００６７】
以下、各部の詳細を説明する。
【００６８】
＜事後確率最大化スペクトル推定部１０６及びスペクトル状態推定部１０４の詳細＞
事後確率最大化スペクトル推定部１０６が、一つのスカラー変数ｚとそのスカラー変数に関する指数関数ｅｘｐ（ｚ）と一つのスカラー定数ａとの和によって規定される非線形方程式（例えば（１９）式）に関して、各短時間フレームｎにおける周波数ｋ毎の周波数信号ｘｎ，ｋと対数尤度重みｗθｎと状態依存スペクトルモデルΘρに依存してスカラー定数ａを定めるとともに（例えば（２１）式）、非線形方程式が０に一致するスカラー変数ｚの値を求め（例えば（１９）式、（２２）式）、その求めたスカラー変数ｚと周波数信号ｘｎ，ｋと対数尤度重みｗθｎと状態依存スペクトルモデルΘρとに基づき、対数スペクトルの推定値ρ＾ｎを更新する（例えば（２０）式）。
【００６９】
図５は事後確率最大化スペクトル推定部１０６の機能ブロック図を、図６はスペクトル状態推定部１０４及び事後確率最大化スペクトル推定部１０６の処理フローを表す。
【００７０】
事後確率最大化スペクトル推定部１０６は、初期値設定部１０６ａ、スカラー定数算出部１０６ｂ、スカラー変数算出部１０６ｃ、対数スペクトル算出部１０６ｄ、収束判定部１０６ｅ及びスペクトル算出部１０６ｆを備える。
【００７１】
初期値設定部１０６ａは、周波数信号ｘｎを受け取り、式（３’）のように対数スペクトルの推定値ρ＾ｎの初期値を、従来の最尤法により求める（ｓ２１）。
【００７２】
【数２０】
0005689844.tif 000021
【００７３】
スペクトル状態推定部１０４が、対数スペクトルの推定値ρ＾ｎに加えて、（２４）式と（２５）式のそれぞれで定義されるスペクトル状態モデルΘθである混合比βｉと状態依存スペクトルモデルΘρである平均μｉｋ及び共分散行列ξｉｋを受け取り、対数尤度重みｗｉを（１８）式に基づき以下のように求める（ｓ１）。
【００７４】
【数２１】
0005689844.tif 000022
【００７５】
さらに、周波数ｋ毎に、以下の手順により、対数スペクトルの推定値ρ＾ｎ，ｋを更新する。
【００７６】
スカラー定数算出部１０６ｂは、周波数信号ｘｎと、対数尤度重みｗｉと、状態依存スペクトルモデルΘρである全ての状態ｉ、全ての周波数ｋに関する平均μｉｋと共分散行列ξｉｋとを受け取り、（２１）式によりスカラー定数ａを求める（ｓ２２）。
【００７７】
スカラー変数算出部１０６ｃは、スカラー定数ａを受け取り、（１９）式に関して、ｆ（ｚ）＝０となるスカラー変数ｚを（近似的に）求める（ｓ２３）。
【００７８】
対数スペクトル算出部１０６ｄは、周波数信号ｘｎと、対数尤度重みｗｉと、状態依存スペクトルモデルΘρである全ての状態ｉ、全ての周波数ｋに関する共分散行列ξｉｋと、スカラー変数ｚとを受け取り、（２０）式を満たす対数スペクトルρｎ，ｋを求め、その推定値ρ＾ｎ，ｋとする（ｓ２４）。
【００７９】
スペクトル状態推定部１０４における処理をＥ−ｓｔｅｐとし、事後確率最大化スペクトル推定部１０６における処理をＭ−ｓｔｅｐとし、ＥＭアルゴリズムに基づき、ｓ１〜ｓ２４を収束条件を満たすまで繰り返す。そのため、収束判定部１０６ｅは、対数スペクトルの推定値ρ＾ｎ，ｋを受け取り、収束条件を満たすか否かを判定する（ｓ３）。収束条件を満たさない場合には、対数スペクトルの推定値ρ＾ｎをスペクトル状態推定部１０４に出力し、各部に対し、処理を繰り返すように制御信号を出力する。収束条件を満たす場合には、対数スペクトルの推定値ρ＾ｎをスペクトル算出部１０６ｆに出力する。
【００８０】
スペクトル算出部１０６ｆは、対数スペクトルの推定値ρ＾ｎを受け取り、各周波数ｋにおけるスペクトルの推定値σ＾ｎ，ｋを、σ＾ｎ，ｋ＝ｅｘｐ（ρ＾ｎ，ｋ）として求め（ｓ２６）、スペクトルの推定値σ＾ｎをスペクトル推定装置１０の出力値として出力する。
【００８１】
＜効果＞
このような構成により、スペクトルの分布を高精度に表現可能な対数スペクトルに関する潜在変数依存型ガウス分布をスペクトルの事前確率密度関数として用いて、効率的にスペクトルの値を推定できる。その結果、周波数信号が誤差を含むような場合でも、効率的かつ高精度に、そのスペクトルの推定が可能になる。
【００８２】
＜変形例＞
第一実施形態の変形例として、ＥＭアルゴリズムのＥ−ｓｔｅｐにおいて、各状態の事後確率ｐ（θｎ｜ρ＾ｎ；Θρ）を求める代わりに、最大の事後確率を与える状態を選択する場合の例を説明する。これは、混合ガウス分布や隠れマルコフモデルを用いた推定において、計算量削減のためにしばしば導入される近似計算である。この変形は、より厳密には、第一実施形態ではθｎを隠れ変数として扱っていたのに対し、ρｎと一緒にθｎも事後確率最大化推定で求めることに相当する。すなわち、以下の問題を解くことに相当する。
【００８３】
【数２２】
0005689844.tif 000023
【００８４】
具体的な処理手順としては、第一実施形態の処理手順の中のｓ１の処理が以下のように修正されるのみで、それ以外は、第一実施形態と同じである。
【００８５】
以下、変形例のスペクトル状態推定部１０４における処理（Ｅ−ｓｔｅｐ、ｓ１）を説明する。図７は変形例のスペクトル状態推定部１０４の機能ブロック図を、図８はその処理フローを示す。スペクトル状態推定部１０４は状態番号推定部１０４ａと対数尤度重み設定部１０４ｂとを含む。
【００８６】
状態番号推定部１０４ａは、対数スペクトルの推定値ρ＾ｎ、と、スペクトル状態モデルΘθである混合比βｉと、状態依存スペクトルモデルΘρである全ての状態ｉ、全ての周波数ｋに関する平均μｉｋ及び共分散行列ξｉｋとを受け取り、事後確率最大となる状態番号の推定値ｉ＾を
【００８７】
【数２３】
0005689844.tif 000024
【００８８】
として求める（ｓ１０４ａ）。
対数尤度重み設定部１０４ｂは推定値ｉ＾を受け取り、対数尤度重みｗｉを
【００８９】
【数２４】
0005689844.tif 000025
【００９０】
として定める（ｓ１０４ｂ）。
【００９１】
＜第二実施形態＞
第一実施形態と異なる部分についてのみ説明する。第二実施形態として、状態パラメータθｎが連続値をとる場合の実施形態について説明する。
【００９２】
スペクトル状態モデル記憶部１０１に記憶されているスペクトル状態モデル、状態依存スペクトルモデル記憶部１０２に記憶されている状態依存スペクトルモデル、及び各部の処理等が、第一実施形態とは異なる。
【００９３】
（状態パラメータの定義）
本実施形態では、状態パラメータθｎとして、周波数信号に対応するメル周波数ケプストラム係数（Mel-frequency cepstral coefficient、以下「ＭＦＣＣ」という）ｃｎを用いる。ＭＦＣＣｃｎは、各次数に対応するＮｃ個の要素ｃｎ，ｍを持つベクトルとして表現されているとする。よって、ｃｎ＝［ｃｎ，１，ｃｎ，２，…，ｃｎ，Ｎｃ］Ｔ、ただし下付添え字ＮｃはＮｃを表す。いま、ｃｎ＝Ｈ（ρｎ）を信号の対数スペクトルρｎをＭＦＣＣに変換する関数とする。すると、Ｈ（ρｎ）は、まず、対数スペクトルρｎの各要素に対数変換の逆変換（ｅｘｐ（・））を適用し、メルフィルタバンク処理（ｍｆｂ（・）と表記）を施し、個々のベクトル要素に対数変換（ｌｏｇ（・））を適用したのち、離散コサイン変換（Ｄ（・））を適応することに対応する。すなわち、Ｈ（ρｎ）は、以下の変換過程で表現される。
【００９４】
【数２５】
0005689844.tif 000026
【００９５】
（スペクトル状態モデルの定義）
本実施形態では、状態パラメータθｎの事前確率密度関数ｐ（θｎ；Θθ）としてＭＦＣＣの混合ガウス分布を用いるとする。これは、ｊをガウス分布の番号とすると、以下でモデル化される。
【００９６】
【数２６】
0005689844.tif 000027
【００９７】
ここで、γｊは分布番号ｊに対応する混合比、μｊとΣｊは、分布番号ｊに対応するガウス分布の平均と共分散行列である。したがって、スペクトル状態モデルΘθは、全てのｊに関するγｊとμｊとΣｊの集合とする。
【００９８】
（状態依存スペクトルモデルの定義）
本実施形態では、状態パラメータであるＭＦＣＣｃｎが既知の場合の対数スペクトルρｎの条件付き確率密度関数ｐ（ρｎ｜ｃｎ；Θρ）は、上記のｃｎ＝Ｈ（ρｎ）の逆変換過程としてモデル化する。一般に、ｃｎ＝Ｈ（ｓｎ）は多対一の変換となるため、その逆変換はユニークには定められない。したがって、その定め方には任意性がある。ここでは、一例を挙げる。まず、以下のように、線形回帰を用いて、ｃｎ＝Ｈ（ρｎ）の疑似逆変換であるρ＾ｎ＝Ｇ（ｃｎ）を定義する。
G(c)=Ac+b        (35)
ただし、Ａは行列（Ｎｋ×Ｎｃ）、ｂはベクトル（Ｎｋ×１）を表す。行列Ａとベクトルｂの値は、事前に音響信号のデータベースにより学習されるか、観測信号を用いて学習されるものとする。すなわち、いま学習用のデータベース（もしくは、観測信号）から、複数の周波数信号ｘｎにそれぞれ対応する複数の対数スペクトルρｎと、それに対応するＭＦＣＣｃｎ＝Ｈ（ρｎ）の組合せが与えられているときに、行列Ａとベクトルｂは、以下のように定められるものとする。
【００９９】
【数２７】
0005689844.tif 000028
【０１００】
また、逆変換誤差ｅ＝ρｎ−ρ＾ｎ＝ρｎ−Ｇ（Ｈ（ρｎ））は、平均０と共分散行列Ξのガウス分布に従うと仮定する。すなわち、
p(e)=N(e;0,Ξ)        (37)
これにより、条件付き確率密度関数ｐ（ρｎ｜ｃｎ；Θρ）は、以下のように定義される。
p(ρn|cn;Θρ)=N(ρn;G(cn),Ξ)        (38)
【０１０１】
本実施形態では、上記条件付き確率密度関数ｐ（ρｎ｜ｃｎ；Θρ）は、周波数毎の要素の積に分解可能と仮定されているので、共分散行列Ξは、対角要素にξｋをもつ対角行列になる。よって、Ξ＝ｄｉａｇ（ξｋ）と表すことができる。Ｇ（ｃ）のｋ番目の要素をＧｋ（ｃ）と書くとすると、ξｋは平均自乗回帰誤差Ｅ｛｜ρｎ，ｋ−Ｇｋ（ｃｎ）｜２｝として、事前に学習されるとする。すると、上記条件付き確率密度関数は、以下のように書くことができる。
【０１０２】
【数２８】
0005689844.tif 000029
【０１０３】
したがって、状態依存スペクトルモデルΘρとして、（３６）式の係数である行列Ａとベクトルｂ、及び全ての周波数ｋにおける逆変換誤差の分散ξｋを含んでいれば、上記条件付き確率密度関数ｐ（ρｎ｜ｃｎ；Θρ）は規定されることになる。
【０１０４】
（最適化関数）
本実施形態では、第一実施形態の変形例と同様に、対数スペクトルρｎと状態パラメータであるＭＦＣＣｃｎの両方をＭＡＰ推定により推定する場合を考える。
【０１０５】
【数２９】
0005689844.tif 000030
【０１０６】
したがって、第一実施形態の変形例と同様に、ρ＾ｎとｃ＾ｎを交互に更新することで、上式を最大化するρ＾ｎとｃ＾ｎを求める。
【０１０７】
スペクトル状態推定部１０４は、対数スペクトルの推定値ρ＾ｎが固定された下で、上式を最大化するＭＦＣＣｃｎの推定値ｃ＾ｎを求める。これは、例えば、ｐ（ｃｎ；Θθ）の混合ガウス分布の分布番号ｊを隠れ変数としたＥＭアルゴリズムで求めることができる（つまり、スペクトル状態推定部１０４と事後確率最大化スペクトル推定部１０６において行われるＥＭアルゴリズムのＥ−ｓｔｅｐ内で、後述する期待値算出部２０４ｂと状態パラメータ算出部２０４ｃとにおいてＥＭアルゴリズムを行う）。このための補助関数は以下のように定めることができる。
【０１０８】
【数３０】
0005689844.tif 000031
【０１０９】
ただし、
【０１１０】
【数３１】
0005689844.tif 000032
【０１１１】
したがって、ＥＭアルゴリズムでは、以下の処理を収束するまで繰り返すことで、（４３）式を最大化するｃｎを求める。これをＭＦＣＣである状態パラメータの推定値ｃ＾ｎとする。
１．Ｅ−ｓｔｅｐ：（４４）式により、Ｅ｛ｊ｜ｃ＾ｎ｝の値を更新する。
２．Ｍ−ｓｔｅｐ：（４３）式を最大化するｃｎの値として、ｃ＾ｎを更新する。具体的には、以下の式を計算する。
【０１１２】
【数３２】
0005689844.tif 000033
【０１１３】
そして、上記のように、本実施形態では、状態パラメータｃｎを潜在変数ではなく、ＭＡＰ推定により求めるべきパラメータとして扱う。このため、確定値として求めた上記の状態パラメータの推定値ｃ＾ｎに関する対数尤度重みは、以下のようにディラックデルタ関数δ（・）を用いて表現される。
【０１１４】
【数３３】
0005689844.tif 000034
【０１１５】
一方、事後確率最大化スペクトル推定部１０６は、ｗｃｎ（ただし、下付添え字ｃｎはｃｎを表す）を受け取り、（１５）式を最大化するρｎを求め、ρ＾ｎとする。（１５）式は、以下のように書き換えられる。
【０１１６】
【数３４】
0005689844.tif 000035
【０１１７】
上式は、（１６）式と同じ形をしているので、本実施形態により効率的に最大化することができる。例えば、本実施形態に基づくスペクトル状態推定部１０４の手順は以下のようになる。
【０１１８】
＜スペクトル状態推定部１０４の詳細＞
図９は第二実施形態のスペクトル状態推定部１０４の機能ブロック図を、図１０はその処理フローを示す。
【０１１９】
スペクトル状態推定部１０４は、初期値算出部２０４ａ、期待値算出部２０４ｂ、状態パラメータ算出部２０４ｃ、収束判定部２０４ｅ及び対数尤度重み算出部２０４ｆを含む。
【０１２０】
初期値算出部２０４ａは、対数スペクトルの推定値ρ＾ｎを受け取り、状態パラメータの初期値をｃ＾ｎ＝Ｈ（ρ＾ｎ）として定める（ｓ２０４ａ）（（３１）式参照）。
【０１２１】
期待値算出部２０４ｂは、状態パラメータの推定値ｃ＾ｎとスペクトル状態モデルΘθである混合比γｊ、平均μｊ及び共分散行列Σｊを受け取り、（４４）式により、期待値Ｅ｛ｊ｜ｃ＾ｎ｝を求める（ｓ２０４ｂ、Ｅ−ｓｔｅｐ）。
【０１２２】
状態パラメータ算出部２０４ｃは、対数スペクトルの推定値ρ＾ｎと、期待値Ｅ｛ｊ｜ｃ＾ｎ｝と、スペクトル状態モデルΘθである平均μｊ及び共分散行列Σｊと、状態依存スペクトルモデルΘρである行列Ａ、ベクトルｂ及び共分散行列Ξとを受け取り、（４５）式により、ＭＦＣＣである状態パラメータの推定値ｃ＾ｎを求める（ｓ２０４ｃ、Ｍ−ｓｔｅｐ）。
【０１２３】
期待値算出部２０４ｂにおける処理をＥ−ｓｔｅｐとし、状態パラメータ算出部２０４ｃにおける処理をＭ−ｓｔｅｐとし、ＥＭアルゴリズムに基づき、収束条件を満たすまでｓ２０４ｂ及びｓ２０４ｃを繰り返す。そのため、収束判定部２０４ｅは、状態パラメータの推定値ｃ＾ｎを受け取り、収束条件を満たすか否かを判定する（ｓ２０４ｅ）。収束条件を満たさない場合には、状態パラメータの推定値ｃ＾ｎを期待値算出部２０４ｂに出力し、各部に対し、処理を繰り返すように制御信号を出力する。収束条件を満たす場合には、状態パラメータの推定値ｃ＾ｎを対数尤度重み算出部２０４ｆに出力する。収束条件としては、例えば、（１）繰り返し回数が所定の回数を超えることや、（２）一つ前の繰り返し時に得られた状態パラメータの推定値と現在の繰り返し時に得られた状態パラメータの推定値との差分が閾値以下であること等が挙げられる。
【０１２４】
対数尤度重み算出部２０４ｆは、状態パラメータの推定値ｃ＾ｎを受け取り、式（４６）により、対数尤度重みｗｃｎを求め（ｓ２０４ｆ）、事後確率最大化スペクトル推定部１０６に出力する。
【０１２５】
なお、事後確率最大化スペクトル推定部１０６は、対数尤度重みｗｃｎと周波数信号ｘｎと状態依存スペクトルモデルΘρとを受け取り、（４７）式を最大化する各周波数ｋにおける対数スペクトルρｎ，ｋを求め、対数スペクトルの推定値ρ＾ｎ，ｋを更新する。
【０１２６】
なお、（４７）式の最大化は、前述までの例と同様、（４７）式を（１９）式の形に書き換えてｆ（ｚ）＝０となるスカラー変数ｚを求めたのち、求めたスカラー変数ｚに対応する対数スペクトルρｎを求めることで実現できる。
【０１２７】
最後に、事後確率最大化スペクトル推定部１０６のスペクトル算出部１０６ｆが、各周波数ｋにおけるスペクトルの推定値σ＾ｎ，ｋを、σ＾ｎ，ｋ＝ｅｘｐ（ρ＾ｎ，ｋ）として求め、スペクトルの推定値σ＾ｎをスペクトル推定装置１０の出力値として出力する。
【０１２８】
＜効果＞
このような構成により、第一実施形態と同様の効果を奏する。
【０１２９】
（シミュレーション結果）
この発明のスペクトル推定装置１０を評価する目的で確認実験を行った。このため、非特許文献１に記載されている残響除去法のなかで、残響除去された周波数信号の推定値からスペクトルを推定する処理の部分で第一実施形態及び第二実施形態を用いた実験を行った。
【０１３０】
非特許文献１による残響除去アルゴリズムは、以下になる。
１．残響除去された周波数信号の推定値ｘ＾ｎを観測信号とする。
２．周波数信号の推定値ｘ＾ｎからそのスペクトルの推定値σ＾ｎを最尤法により求める。
３．以下を収束するまで繰り返す。
（ａ）観測信号とスペクトルの推定値σ＾ｎから残響の予測係数を更新する。
（ｂ）観測信号と残響の予測係数から残響除去した信号の周波数信号の推定値ｘ＾ｎを求める。
（ｃ）周波数信号の推定値ｘ＾ｎからそのスペクトルの推定値σ＾ｎを最尤法により求める。
４．求められた周波数信号の推定値ｘ＾ｎを時間領域信号に変換し、残響除去された信号として出力する。
【０１３１】
本実験では、上記の３（ｃ）の処理において、最尤法の代わりに、第一実施形態及び第二実施形態を用いる場合と用いない場合の比較を行った。図１１は、その結果を示す。３つのグラフのそれぞれは、左から順に、長さの異なる３種類の観測信号（平均長さは、それぞれ１．１５秒、２．３秒、４．６秒）を用いた場合の結果を示している。各グラフの横軸は、上記の残響除去アルゴリズムの繰り返し回数を表している。繰り返し回数０は、観測信号を表す。縦軸は、残響除去された信号のケプストラム歪（ＣＤ）を示す。二点鎖線が非特許文献１の残響除去法で、一点鎖線が第一実施形態の方法でスペクトル推定を行った場合、実線が第二実施形態の方法でスペクトル推定を行った場合を示す。全ての場合において、２回目以降の繰り返しにおいて、非特許文献１の残響除去法よりも第一実施形態及び第二実施形態によるスペクトル推定を用いた場合の方が、ケプストラム歪を小さくできている。なお、上記の残響除去アルゴリズムにおいて、第一実施形態及び第二実施形態により推定されたスペクトルに基づき残響除去が行われるのは、２回目以降の繰り返しにおいてである。このため、一回目の繰り返しでは、第一実施形態及び第二実施形態を用いる場合と用いない場合で、ケプストラム歪の値に差は生じない。
【０１３２】
以上の結果より、第一実施形態及び第二実施形態により、潜在変数依存型ガウス分布を対数スペクトルの事前分布として導入し、対数スペクトルを事後確率最大化推定により求めることで、スペクトル推定精度を改善できることが確認された。
【０１３３】
＜その他の変形例＞
本発明は上記の実施形態及び変形例に限定されるものではない。例えば、上述の各種の処理は、記載に従って時系列に実行されるのみならず、処理を実行する装置の処理能力あるいは必要に応じて並列的にあるいは個別に実行されてもよい。その他、本発明の趣旨を逸脱しない範囲で適宜変更が可能である。
【０１３４】
＜プログラム及び記録媒体＞
上述したスペクトル推定装置は、コンピュータにより機能させることもできる。この場合はコンピュータに、目的とする装置（各種実施形態で図に示した機能ブロック図をもつ装置）として機能させるためのプログラム、またはその処理手順（各実施形態で示したもの）の各過程をコンピュータに実行させるためのプログラムを、ＣＤ−ＲＯＭ、磁気ディスク、半導体記憶装置などの記録媒体から、あるいは通信回線を介してそのコンピュータ内にダウンロードし、そのプログラムを実行させればよい。
【産業上の利用可能性】
【０１３５】
本発明は、各短時間フレームにおける周波数信号のスペクトル値を用いて行う様々な処理に利用することができる。
【符号の説明】
【０１３６】
１０   スペクトル推定装置
１０１ スペクトル状態モデル記憶部
１０２ 状態依存スペクトルモデル記憶部
１０４ スペクトル状態推定部
１０４ａ      状態番号推定部
１０４ｂ      設定部
１０６ 事後確率最大化スペクトル推定部
１０６ａ      初期値設定部
１０６ｂ      スカラー定数算出部
１０６ｃ      スカラー変数算出部
１０６ｄ      対数スペクトル算出部
１０６ｅ      収束判定部
１０６ｆ      スペクトル算出部
２０４ａ      初期値算出部
２０４ｂ      期待値算出部
２０４ｃ      状態パラメータ算出部
２０４ｅ      収束判定部
２０４ｆ      算出部
【図１】
0005689844.tif 000046
【図２】
0005689844.tif 000047
【図３】
0005689844.tif 000048
【図４】
0005689844.tif 000049
【図５】
0005689844.tif 000050
【図６】
0005689844.tif 000051
【図７】
0005689844.tif 000052
【図８】
0005689844.tif 000053
【図９】
0005689844.tif 000054
【図１０】
0005689844.tif 000055
【図１１】
0005689844.tif 000056