JP6420781B2 - 声道スペクトル推定装置、声道スペクトル推定方法、及びプログラム - Google Patents

Description

本発明は、声道スペクトル推定装置、声道スペクトル推定方法、及びプログラムに係り、特に、音声信号から声道スペクトルを推定するための声道スペクトル推定装置、声道スペクトル推定方法、及びプログラムに関する。

音声合成、或いは音声変換をはじめ音声処理全般において、音声信号から声道スペクトルを推定する技術は多くの場面で用いられている。短区間ごとの音声信号を、周期デルタ関数（パルス列）を入力とした線形時不変系の出力としてモデル化できると仮定すれば、この線形系の入力とインパルス応答がそれぞれ声帯音源信号と声道特性に対応する。この仮定は、周波数領域では周期デルタ関数で表される声帯音源スペクトルと声道スペクトルの積で音声スペクトルが表されることに相当する。従って、音声スペクトルは、声道スペクトルを基本周波数間隔で周期的にサンプリングしたものであると見なすことができる。

この観点に基づいて、音声スペクトルから声道スペクトルを推定する方法がこれまで提案されている。代表的な方法の一つとして広く知られる“STRAIGHT ” は、音声信号を基本周期の幅で切り出し、切り出した各々の信号である切り出し信号のスペクトルを声道スペクトルの推定値とする方法である（非特許文献１）。これは周波数領域では、各調波成分のピークをsinc 補間したものを声道スペクトルと見なしていることに相当する。

しかし、STRAIGHTによって得られる声道スペクトル推定値は、定常な音声が対象であっても、音声信号を基本周期の幅で切り出す切り出しフレームのオフセットに依存して周期的に時間変化することが知られている。これは各調波成分が互いに干渉し合うからである。こうした周期的に時間変化するスペクトルの変動成分は、周期信号に対する有限窓を用いた周波数分析により不可避的に生じるものであり、声道スペクトル推定値に本来含めるべきものではない。したがって、STRAIGHTを用いた声道スペクトルの推定において、声道スペクトル推定値から変動成分を除くように改良された手法が提案されている（非特許文献２）。

前述したように、音声スペクトルは声道スペクトルを基本周波数(F₀)間隔でサンプリングしたものと見なせるため、音声のF₀が高いときほど声道スペクトル推定の手がかりは少なくなる。このことは、１フレームごとに独立な処理に本質的な限界があることを示唆している。

一方で、音声信号には同一の音素が繰り返し出現するため、類似した声道スペクトルが複数の異なる時刻で現れることも、声道スペクトル推定の手がかりとなる。複数のフレームが共通の声道スペクトルを持つと仮定でき、当該複数のフレームでF₀ が異なれば、実際に観測可能な声道スペクトルのサンプル点が単一のフレームの場合よりも増加するため、声道スペクトルの推定精度が向上すると考えられる。

こうした考えに基づき、同時に収録された調音運動データを用いて、複数フレームから声道スペクトルを推定する手法 が提案されている（非特許文献３）。また、同様の手法として、因子分析トラジェクトリ隠れマルコフモデルによる声道スペクトル推定法が提案されている（非特許文献４）。因子分析トラジェクトリ隠れマルコフモデルによる声道スペクトル推定法では、音声信号の各フレームに付与されているコンテキストラベルを用い、同一のコンテキストが付与された複数のフレームにおける調波成分の情報に加え動的特徴量を手がかりにすることで、声道スペクトルを推定する。

H. Kawahara, I. Masuda-Katsuse, A. de Cheveigne, "Restructuring speech representations using a pitch-adaptive time-frequency smoothing and an instantaneous-frequency-based F0 extraction:Possible role of a repetitive structure in sounds," Speech Commun., 27, pp. 187-207, 1999. H. Kawahara, M. Morise, T. Takahashi, R. Nisimura, T. Irino and H. Banno, "Tandem-STRAIGHT: A temporally stable power spectral representation for periodic signals and applications to interference-free spectrum, F0, and aperiodicity estimation," Proc. ICASSP, pp. 3933-3936, 2008. Y. Shiga and S. King, "Estimating the spectral envelope of voiced speech using multi-frame analysis," Proc. EUROSPEECH, pp. 1737-1740, 2003. T. Toda, "Statistical approach to vocal tract transfer function estimation based on factor analyzed trajectory HMM," Proc. ICASSP, pp. 3925-3928, 2008.

しかしながら、因子分析トラジェクトリ隠れマルコフモデルによる声道スペクトル推定法では、複数のフレームにおける調波成分の情報を手がかりにすることで、声道スペクトルを精度良く推定することは可能であるが，音声信号に対するコンテキストラベルの付与には膨大な労力を要するという問題がある。

本発明は、上記の事情を鑑みて成されたものであり、音声信号から声道スペクトルを精度良く推定することができる声道スペクトル推定装置、声道スペクトル推定方法、及びプログラムを提供することを目的とする。

上記の目的を達成するために本発明に係る声道スペクトル推定装置は、音声信号の時系列データを、基本周期の幅で切り出し、各々切り出した音声信号のスペクトルから、各時刻及び各正規化角周波数の観測時間周波数成分を表す観測スペクトログラムを出力する観測スペクトログラム推定部と、前記観測スペクトログラム推定部により出力された前記観測スペクトログラムと、各スペクトルパターンにおける各正規化角周波数のパワースペクトルを表す声道スペクトル、及び各スペクトルパターンの各時刻における重みから求められる各時刻及び各正規化角周波数の時間周波数成分を表す声道スペクトログラムとの距離を用いて表される目的関数を小さくするように、各スペクトルパターンにおける各正規化角周波数の前記声道スペクトル、及び各スペクトルパターンの各時刻における重みを推定する推定部と、を含んで構成されている。

本発明に係る声道スペクトル推定方法は、観測スペクトログラム推定部と推定部とを含む声道スペクトル推定装置における声道スペクトル推定方法であって、前記観測スペクトログラム推定部が、音声信号の時系列データを、基本周期の幅で切り出し、各々切り出した音声信号のスペクトルから、各時刻及び各正規化角周波数の観測時間周波数成分を表す観測スペクトログラムを出力し、前記推定部が、前記観測スペクトログラム推定部により出力された前記観測スペクトログラムと、各スペクトルパターンにおける各正規化角周波数のパワースペクトルを表す声道スペクトル、及び各スペクトルパターンの各時刻における重みから求められる各時刻及び各正規化角周波数の時間周波数成分を表す声道スペクトログラムとの距離を用いて表される目的関数を小さくするように、各スペクトルパターンにおける各正規化角周波数の前記声道スペクトル、及び各スペクトルパターンの各時刻における重みを推定する。

本発明に係るプログラムは、上記の声道スペクトル推定装置の各部としてコンピュータを機能させるためのプログラムである。

以上説明したように、本発明の声道スペクトル推定装置、声道スペクトル推定方法、及びプログラムによれば、観測スペクトログラムと、声道スペクトログラムとの距離を用いて表される目的関数を小さくするように、各スペクトルパターンにおける各正規化角周波数の声道スペクトル、及び各スペクトルパターンの各時刻における重みを推定することにより、音声信号から声道スペクトルを精度良く推定することができる、という効果が得られる。

本発明の実施の形態に係る声道スペクトル推定装置の構成を示す概略図である。 GMM-NMFに対する声道スペクトルの推定アルゴリズムを利用する声道スペクトル推定処理ルーチンの処理の一例を示すフローチャートである。 AR-NMFに対する声道スペクトルの推定アルゴリズムを利用する声道スペクトル推定処理ルーチンの処理の一例を示すフローチャートである。 声道スペクトル推定装置によって推定した声道スペクトルの評価結果を示す図である。

以下、図面を参照して本発明の実施の形態を詳細に説明する。
＜本発明の実施の形態の概要＞
声道スペクトログラムが低ランクな非負値行列で近似できるという仮定に基づいて、音声信号に付与されたコンテキストラベルを用いることなく、音声信号の複数のフレームにおける調波成分の情報を手がかりにして、声道スペクトルを精度良く推定する。

＜欠損データ補間による声道スペクトル推定＞
＜低ランク非負値行列による声道スペクトログラムのモデル化＞
まず、低ランク非負値行列による声道スペクトログラムのモデル化について説明する。時間インデックスをｔ(ｔ＝0,・・,T-1)とし、周波数インデックスｋ(ｋ＝0,・・,K-1)と対応する正規化角周波数をそれぞれω_kと表す。

音声における音素の種類は限られていることから、音声の一発話において、同じ声道スペクトルが複数の異なる時刻に現れる。この状況は、声道スペクトログラムを低ランクな非負値行列で表現可能であるということと同等であると解釈することができる。したがって、例えばＲ個の滑らかなスペクトルパターンを列方向に並べた非負値行列Ｈ＝(Ｈ_k,r)_k,rと、各スペクトルパターンの非負値の重みＵ＝(Ｕ_r,t)_r,tによって、声道スペクトログラムＸ_k,tは（１）式で表すことができる。

なお、ｒはスペクトルパターンのインデックスである。

周波数方向に滑らか、且つ、非負となるＨ_k,rは様々に設計可能であるが、本実施の形態では２種類のＨ_k,rを提案する。

１つ目の提案は、（２）式及び（３）式に示すように、Ｈ_k,rを混合正規分布型の関数で規定する。

ここで、ωは正規化角周波数、Ｎは混合数、ｎ(ｎ＝0,・・,N-1)は混合数Ｎのインデックスである。ｈ(ω)は周波数ワーピング関数であり、ｈ(ω)＝ωとすれば、Ｇ_n(ω)は線形周波数領域で平均ρ_n、分散ν² _nの正規分布と同形の関数となる。本実施の形態では、メル周波数領域で滑らかな声道スペクトルとなるように、周波数ワーピング関数ｈ(ω)を例えば（４）式のように設定する。

ここで、ｆ_sはサンプリング周波数であり、（４）式によって[0,π]の正規化周波数は[0,1]にマッピングされる。また、Ｗ＝(Ｗ_r,n)_r,n≧０は各正規分布の重みであり、非負値行列Ｈ及び重みＵに関するスケールの任意性を解消するため、各ｒに関してΣ_nＷ_r,n＝１とする。

２つ目の提案は、ソースフィルタモデルでよく用いられる全極フィルタを利用する方法である。具体的には、Ｐ次の全極フィルタの係数ａ_r:＝[ａ_r,0,ａ_r,1,・・,ａ_r,P]^Tを用いれば、全極フィルタの振幅スペクトルＨ_k,r ^(AR)、すなわち、声道スペクトルＨ_k,r ^(AR)は（５）式で表される。

ここで、Ｑ(ω)は、(ｐ,ｑ)成分がcos(ω(ｐ−ｑ))で表される(Ｐ＋１)×(Ｐ＋１)のToeplitz行列である。

＜欠損データに対する非負値行列因子分解アプローチ＞
STRAIGHTによって推定された声道スペクトログラムＹ＝(Ｙ_k,t)_k,tが与えられた場合、声道スペクトログラムの推定問題は、与えられた声道スペクトログラムＹ_k,tと推定した声道スペクトログラムＸ_k,tとの距離Ｄ_*(Ｙ_k,t;Ｘ_k,t)を用いて、（６）式のように定式化ができる。

ここで、Θはパラメータ集合であり、声道スペクトルＨ_k,rとして（２）式のＨ_k,r ^(GMM)を用いた場合は、Θ＝{Ｗ,Ｕ}であり、声道スペクトルＨ_k,rとして（５）式のＨ_k,r ^(AR)を用いた場合は、Θ＝{{ａ_r}_r,Ｕ}である。以降、説明の便宜上、Ｈ_k,r ^(GMM)を用いた場合の非負値行列因子分解(Nonnegative Matrix Factorization:NMF)を“GMM-NMF”、Ｈ_k,r ^(AR)を用いた場合のNMFを“AR-NMF”と称す。また、Ｚ_k,t∈[0,1]は、各時間周波数成分の信頼度を表すパラメータである。

なお、（６）式においてＺ_k,t＝0であれば、STRAIGHTによって推定された声道スペクトログラムＹ_k,tに対するコストは考慮されず、Ｚ_k,tが大きい時間周波数成分ほど重視される。

この信頼度Ｚ_k,tの単純な設計方法として、各時刻での基本周波数F₀とその高調波周波数に対応する時間周波数成分にはＺ_k,t=1、それ以外の時間周波数成分にはＺ_k,t=ξ（ただし、ξは0以上1以下の定数）とする方法が考えられる。ξの値は実験的、経験的に決定することもできるが，本実施の形態ではSTRAIGHTで得られる非周期性指標Ａ_k,t∈[0,1]を利用して設計する方法について述べる。非周期性指標Ａ_k,tは、各時間周波数成分に含まれる非周期成分の割合であるため、各ｋ,ｒに関してＺ_k,t＝1−Ａ_k,tとすれば、周期性成分を重視した声道スペクトログラムＸ_k,tの推定が可能になる。

NMFで広く知られている距離Ｄ_*(Ｙ_k,t;Ｘ_k,t)として、例えば一般化Kullback-Leibler(KL)ダイバージェンスＤ_GKL、または２乗距離Ｄ_EUが挙げられる。一般化KLダイバージェンスＤ_GKLを用いた場合の目的関数Ｌ_GKL(Θ)を（７）式に示し、２乗距離Ｄ_EUを用いた場合の目的関数Ｌ_EU(Θ)を（８）式に示す。

＜補助関数法によるパラメータ推定アルゴリズム＞
次に、GMM-NMFに対するパラメータ推定アルゴリズム、及びAR-NMFに対するパラメータ推定アルゴリズムについて説明する。

＜GMM-NMFに対する反復アルゴリズムの導出＞
上記（７）式の一般化KLダイバージェンスＤ_GKLを用いた場合の目的関数Ｌ_GKL(Θ)における右辺括弧内の第2項は、対数関数の中に加算演算を含んでいるため、（７）式の目的関数を最小化する最適化問題を直接解くことは困難である。

しかし、多くのNMFを用いた研究で行われているように、最適化問題を直接解くことが困難である目的関数に、補助関数法と呼ばれる最適化原理を適用することによって、反復的に局所最適解を得ることができることが知られている。

補助関数法では、パラメータΘの目的関数Ｌ(Θ)に対して補助変数λを導入し、Ｌ(Θ)＝ｍｉｎ_λＬ⁺(Θ,λ)を満たす上界を規定する補助関数Ｌ⁺(Θ,λ)を導出する。補助関数Ｌ⁺(Θ,λ)をパラメータΘ、補助変数λに関して交互に最小化することによって、目的関数Ｌ(Θ)を広義単調減少させることができる。

対数関数は凹関数であるため、上記（７）式の右辺括弧内の第2項の上界は、Jensenの不等式を用いて（９）式で表される。

ここで、λ_k,t,r,n≧0は補助変数であり、各k,tに関してΣ_r,nλ_k,t,r,n＝1を満たす。なお、（９）式の等式成立条件は（１０）式を満たす場合となる。

したがって、目的関数Ｌ_GKL(Θ)に対する補助関数は（１１）式で表される。

ここで、λ:＝{λ_k,t,r,n}と定義した。（１１）式の補助関数のＷ_r,n、Ｕ_r,tに関する偏微分が0となる値を求めて（１０）式を代入することにより、（１２）式及び（１３）式に示す閉形式の更新式が得られる。

（１２）式及び（１３）式の更新式は全て非負値の項同士の積として計算されるため、初期値を非負値にすればＷ_r,n、Ｕ_r,tの非負値性は自然と保たれる。

次に、距離Ｄ_*(Ｙ_k,t;Ｘ_k,t)として２乗距離Ｄ_EUを用いた場合の目的関数Ｌ_EU(Θ)の更新式について検討する。（８）式も右辺括弧内の第2項に加算演算を含んでいるため、（８）式の目的関数を最小化する最適化問題を直接解くことは困難である。

しかし、２次関数は凸関数であるため、上記（８）式に対してJensenの不等式を用いることで、目的関数Ｌ_GKL(Θ)の場合と同様に補助関数を設計することができる。補助関数の設計は、目的関数Ｌ_GKL(Θ)における補助関数の設計方法と同様であるため省略するが、（８）式を最小化する更新式は、（１４）式及び（１５）式で表される。

＜AR-NMFに対する反復アルゴリズムの導出＞
AR-NMFに関しても、GMM-NMFに対する反復アルゴリズムの導出方法と同様にして、閉形式の更新式を導出することができる。

まず、一般化KLダイバージェンスの目的関数Ｌ_GKL(Θ)に対する補助関数Ｌ⁺ _GKL,AR(Θ,ξ)について検討する。補助関数Ｌ⁺ _GKL,AR(Θ,ξ)は、各ｋ,ｔに関してΣ_rξ_k,t,r＝１を満たす非負の補助変数ξ＝{ξ_k,t,r}_k,t,rを導入することで、（１６）式によって定義される。

また、（１６）式の等式成立条件はξ_k,t,r＝Ｈ^(AR) _k,r(Ｕ_r,t／Ｘ_k,t)となる。Ｕ_r,tの更新式は、（１３）式のＨ_k,r ^(GMM)をＨ_k,r ^(AR)に置換したものと同じになり、（１７）式で表される。

一方、AR-NMFに対する反復アルゴリズムを導出する場合、ａ_rの更新には、乗法更新型アルゴリズムを利用できる。

（７）式に示した一般化KLダイバージェンスの目的関数Ｌ_GKL(Θ)のａ_rに関する偏微分は、（１８）式〜（２０）式で表すことができる。

（１８）式の右辺括弧内の第１項及び第２項はどちらも正定値行列であり、（２１）式の乗法更新測を用いることで目的関数Ｌ_GKL(Θ,ξ)を広義単調減少させることができる。

詳細は省略するが、２乗距離Ｄ_EUを用いた場合の更新式についても、一般化KLダイバージェンスの場合と同様にして導出することができる。具体的には、Ｕ_r,tの更新式は、（１５）式のＨ_k,r ^(GMM)をＨ_k,r ^(AR)に置換した（２２）式で表され、ａ_rの更新式は（２１）式と同様に（２３）式となる。

＜第１の実施の形態＞
＜システム構成＞
次に、音声信号の複数のフレームにおける調波成分の情報を手がかりにして、声道スペクトルを推定する声道スペクトル推定装置に本発明を適用した場合を例にして、本発明の第１の実施の形態を説明する。なお、第１の実施の形態では距離Ｄ_*(Ｙ_k,t;Ｘ_k,t)として、一般化KLダイバージェンスＤ_GKLを用いた場合のGMM-NMFに対する声道スペクトルの推定アルゴリズムを利用する声道スペクトル推定装置の例について説明する。

図１に示すように、本発明の第１の実施の形態に係る声道スペクトル推定装置は、ＣＰＵと、ＲＡＭと、後述する声道スペクトル推定処理ルーチンを実行するためのプログラムを記憶したＲＯＭとを備えたコンピュータで構成され、機能的には次に示すように構成されている。

声道スペクトル推定装置１００は、入力部１０と、演算部２０と、記憶部３０と、出力部４０とを備えている。

入力部１０により、声道スペクトルの推定対象である音声信号の時系列データが入力される。記憶部３０は、入力部１０により入力された音声信号の時系列データを記憶する。また、記憶部３０は、後述する各処理での結果を記憶すると共に、本処理ルーチンで用いる各パラメータの初期値を記憶している。

演算部２０は、観測スペクトログラム推定部２１と、初期設定部２２と、推定部２３と、終了判定部２４と、出力部２５とを備えている。

観測スペクトログラム推定部２１は、入力部１０で収集した音声信号の時系列データを入力として、公知の声道スペクトルの推定手法であるSTRAIGHTを用いて推定した声道スペクトログラムＹ_k,tを計算する。また、計算した声道スペクトログラムＹ_k,tを、記憶部３０に記憶しておく。より詳細には、観測スペクトログラム推定部２１は、入力部１０で収集した音声信号の時系列データを基本周期の幅で切り出し、各々切り出した音声信号から、各時刻ｔ及び各正規化角周波数ω_kの観測時間周波数成分を表す声道スペクトログラムＹ_k,tを計算する。なお、声道スペクトログラムＹ_k,tは、本実施の形態に係る観測スペクトログラムの一例である。

また、観測スペクトログラム推定部２１は、声道スペクトログラムＹ_k,tを計算する際、観測時間周波数成分毎に非周期性指標Ａ_k,tを計算し、Ｚ_k,t＝1−Ａ_k,tにより各時間周波数成分の信頼度Ｚ_k,tを算出する。算出した各時間周波数成分の信頼度Ｚ_k,tは、記憶部３０に記憶される。

初期設定部２２は、後述する処理で用いるパラメータＷ_r,n、Ｕ_r,t、及びＧ_n(ω_k)の各初期値を設定する。なお、Ｗ_r,nは非負値であり、Σ_nＷ_r,n＝１を満たすように初期値を設定する。Ｕ_r,tも非負値であり、例えば乱数を用いて適当な値に初期値を設定する。なお、ｒ(r=0,・・,R-1)は、Ｒ個のスペクトルパターンを指し示すスペクトルインデックスである。Ｇ_n(ω_k)に対しては、上記（３）式を満たすように初期値を設定する。この際、平均ρ_n、分散ν² _nには適当な値を用いればよい。設定したパラメータＷ_r,n、Ｕ_r,t、及びＧ_n(ω_k)の各初期値は、記憶部３０に記憶される。

推定部２３は、（ｋ、ｒ）の全ての組み合わせの各々について、記憶部３０に記憶されているＷ_r,n及びＧ_n(ω_k)に基づいて、上記（２）式に従って、声道スペクトルＨ_k,r ^(GMM)（以降、単に「Ｈ_k,r」と記載する）を計算し、記憶部３０に格納する。

推定部２３は、（ｋ、ｔ）の全ての組み合わせの各々について、記憶部３０に記憶されているＨ_k,r及びＵ_r,tに基づいて、上記（１）式に従って、声道スペクトログラムＸ_k,tを計算し、記憶部３０に格納する。

推定部２３は、（ｒ、ｎ）の全ての組み合わせの各々について、記憶部３０に記憶されているＷ_r,n、Ｚ_k,t、Ｘ_k,t、Ｙ_k,t、Ｇ_n(ω_k)、Ｕ_r,tに基づいて、上記（１２）式に従って、正規分布の重みＷ_r,nを更新し、記憶部３０に格納する。

推定部２３は、正規分布の重みＷ_r,nの更新に伴い、（ｋ、ｒ）の全ての組み合わせの各々について、Ｗ_r,n及びＧ_n(ω_k)に基づいて、上記（２）式に従って、声道スペクトルＨ_k,rを計算し、記憶部３０に格納する。

推定部２３は、声道スペクトルＨ_k,rの更新に伴い、（ｋ、ｔ）の全ての組み合わせの各々について、記憶部３０に記憶されているＨ_k,r及びＵ_r,tに基づいて、上記（１）式に従って、声道スペクトログラムＸ_k,tを計算し、記憶部３０に格納する。

推定部２３は、（ｒ、ｔ）の全ての組み合わせの各々について、記憶部３０に記憶されているＵ_r,t、Ｚ_k,t、Ｘ_k,t、Ｙ_k,t、Ｈ^(GMM) _k,r、すなわちＨ_k,rに基づいて、上記（１３）式に従って、声道スペクトログラムＸ_k,tのスペクトルパターンの重みＵ_r,tを更新し、記憶部３０に格納する。

更に、推定部２３は、声道スペクトログラムＸ_k,tのスペクトルパターンの重みＵ_r,tの更新に伴い、（ｋ、ｔ）の全ての組み合わせの各々について、記憶部３０に記憶されているＨ_k,r及びＵ_r,tに基づいて、上記（１）式に従って、声道スペクトログラムＸ_k,tを計算し、記憶部３０に格納する。

終了判定部２４は、予め定められた終了条件を満足するか否かを判定し、終了条件を満足していない場合には、推定部２３の各処理を繰り返す。終了判定部２４は、終了条件を満足したと判定した場合には、出力部２５による処理に移行する。

出力部２５は、記憶部３０に記憶されている声道スペクトログラムＸ_k,tを、入力部１０に入力された音声信号から推定される声道スペクトログラムとして出力する。

なお、終了条件としては、繰り返し回数がＬ-1回目の目的関数（７）式の値と、繰り返し回数がＬ回目の目的関数（７）式の値との差が、予め定めた閾値よりも小さくなったことを用いればよい。あるいは、終了条件として、繰り返し回数が、予め定められた上限回数に到達したことを用いてもよい。

＜声道スペクトル推定装置の作用＞
次に、GMM-NMFに対する声道スペクトルの推定アルゴリズムを利用する第１の実施の形態に係る声道スペクトル推定装置１００の作用について説明する。

マイクロホンで取得された音声信号の時系列データが声道スペクトル推定装置１００に入力され、記憶部３０に格納される。そして、声道スペクトル推定装置１００において、図２に示す声道スペクトル推定処理ルーチンが実行される。

まず、ステップＳ１００において、記憶部３０から、音声信号の時系列データを読み込み、当該音声信号の時系列データに対して、STRAIGHTによる声道スペクトルの推定を行い、各時間ｔ及び各正規化角周波数ω_k(k=0,・・,K-1)のの観測時間周波数成分を表す声道スペクトログラムＹ_k,tを計算し、得られた声道スペクトログラムＹ_k,tを、記憶部３０に格納する。

ステップＳ１０２において、ステップＳ１００で得られた声道スペクトログラムＹ_k,tに対応する非周期性指標Ａ_k,tを計算し、Ｚ_k,t＝1−Ａ_k,tにより各時間周波数成分の信頼度Ｚ_k,tを計算し、得られた信頼度Ｚ_k,tを記憶部３０に格納する。

ステップＳ１０４において、乱数を用いてＷ_r,n及びＵ_r,tの初期値を設定する。なお、Ｗ_r,n及びＵ_r,tの初期値は非負値とし、Ｗ_r,nに対しては、Σ_nＷ_r,n＝１を満たすように初期値を設定する。また、Ｇ_n(ω_k)については、上記（３）式を満たすように初期値を設定する。この際、平均ρ_n、分散ν² _nには適当な値を用いればよい。

こうして設定されたＷ_r,n、Ｕ_r,t、及びＧ_n(ω_k)の各初期値は、記憶部３０に記憶される。

次に、ステップＳ１０６では、ステップＳ１０４で設定されたＷ_r,n及びＧ_n(ω_k)に基づいて、上記（２）式に従って、声道スペクトルＨ_k,rを各（ｋ、ｒ）の組み合わせについて算出して、記憶部３０に格納する。また、ステップＳ１０６では、ステップＳ１０４で設定されたＵ_r,t、及び本ステップで算出されたＨ_k,rに基づいて、上記（１）式に従って、声道スペクトログラムＸ_k,tを各（ｋ、ｔ）の組み合わせについて算出して、記憶部３０に格納する。

ステップＳ１０８では、ステップＳ１０４で設定されたＷ_r,n、Ｇ_n(ω_k) 及びＵ_r,t、ステップＳ１０２で算出されたＺ_k,t、ステップＳ１０６で算出されたＸ_k,t、並びにステップＳ１００で算出されたＹ_k,t、すなわち、記憶部３０に記憶されている最新の各パラメータＷ_r,n、Ｚ_k,t、Ｘ_k,t、Ｙ_k,t、Ｇ_n(ω_k)、Ｕ_r,tに基づいて、上記（１２）式に従って、正規分布の重みＷ_r,nを各（ｒ、ｎ）の組み合わせについて算出して、記憶部３０に格納する。

ステップＳ１１０では、ステップＳ１０６と同様に、記憶部３０に記憶されている最新の各パラメータＷ_r,n及びＧ_n(ω_k)、すなわち、ステップＳ１０８で算出されたＷ_r,n及びステップＳ１０４で設定されたＧ_n(ω_k)に基づいて、上記（２）式に従って、声道スペクトルＨ_k,rを各（ｋ、ｒ）の組み合わせについて算出して、記憶部３０に格納する。また、ステップＳ１１０では、記憶部３０に記憶されている最新の各パラメータＵ_r,t及びＨ_k,r、すなわち、ステップＳ１０４で設定されたＵ_r,t、及び本ステップで算出されたＨ_k,rに基づいて、上記（１）式に従って、声道スペクトログラムＸ_k,tを各（ｋ、ｔ）の組み合わせについて算出して、記憶部３０に格納する。

ステップＳ１１２では、記憶部３０に記憶されている最新の各パラメータＺ_k,t、Ｘ_k,t、Ｙ_k,t、Ｈ_k,r、Ｕ_r,t、すなわち、ステップＳ１０４で設定されたＵ_r,t、ステップＳ１０２で算出されたＺ_k,t、ステップＳ１１０で算出したＸ_k,t及びＨ_k,r、ステップＳ１００で算出したＹ_k,tに基づいて、上記（１３）式に従って、声道スペクトログラムＸ_k,tのスペクトルパターンの重みＵ_r,tを各（ｒ、ｔ）の組み合わせについて算出して、記憶部３０に格納する。

ステップＳ１１４では、記憶部３０に記憶されている最新の各パラメータＵ_r,t及びＨ_k,r、すなわち、ステップＳ１１２で算出されたＵ_r,t、及びステップＳ１１０で算出されたＨ_k,rに基づいて、上記（１）式に従って、声道スペクトログラムＸ_k,tを各（ｋ、ｔ）の組み合わせについて算出して、記憶部３０に格納する。

次のステップＳ１１６では、ステップＳ１００で算出したＹ_k,tと、ステップＳ１０２で算出したＺ_k,tと、ステップＳ１１４で算出したＸ_k,tに基づいて、（７）式に従って目的関数Ｌ_GKL(Θ)の値を算出して、記憶部３０に記憶する。そして、前回のステップＳ１１６で算出した目的関数Ｌ_GKL(Θ)の値を記憶部３０から読み込み、今回のステップＳ１１６で算出した目的関数Ｌ_GKL(Θ)の値と、前回のステップＳ１１６で算出した目的関数Ｌ_GKL(Θ)の値との差分が、予め記憶部３０に記憶されている予め定められた閾値よりも小さいか否かを判定し、差分が予め定められた閾値以上の場合には、終了条件を満足していないと判断して、上記ステップＳ１０８へ戻り、上記ステップＳ１０８〜ステップＳ１１６の処理を繰り返す。

一方、差分が予め定められた閾値未満の場合には、終了条件を満足したと判断して、ステップＳ１１８で、ステップＳ１１４で算出した最新の声道スペクトログラムＸ_k,tを出力して、声道スペクトル推定処理ルーチンを終了する。

なお、上記では、距離Ｄ_*(Ｙ_k,t;Ｘ_k,t)として、一般化KLダイバージェンスＤ_GKLを用いた場合のGMM-NMFに対する声道スペクトルの推定アルゴリズムを利用する声道スペクトル推定装置の例について説明したが、距離Ｄ_*(Ｙ_k,t;Ｘ_k,t)として、２乗距離Ｄ_EUを用いてもよいことは言うまでもない。

この場合、例えばステップＳ１０８が（１４）式に示したＷ_r,nの更新式の計算に置き換わり、ステップＳ１１２が（１５）式に示したＵ_r,tの更新式の計算に置き換わる。

＜第２の実施の形態＞
＜システム構成＞
次に、第２の実施の形態に係る声道スペクトル推定装置について説明する。第２の実施の形態では、距離Ｄ_*(Ｙ_k,t;Ｘ_k,t)として、一般化KLダイバージェンスＤ_GKLを用いた場合のAR-NMFに対する声道スペクトルの推定アルゴリズムを利用する声道スペクトル推定装置の例について説明する。本発明の第２の実施の形態に係る声道スペクトル推定装置は、図１に示した第１の実施の形態に係る声道スペクトル推定装置のシステム構成と同様に、入力部１０と、演算部２０と、記憶部３０と、出力部４０とを備えている。また、演算部２０は、観測スペクトログラム推定部２１と、初期設定部２２と、推定部２３と、終了判定部２４と、出力部２５とを備えている。

入力部１０及び観測スペクトログラム推定部２１については、第１の実施の形態と同様であるため、説明を省略する。

初期設定部２２は、後述する処理で用いるパラメータａ_r、及びＵ_r,tの各初期値を設定する。ここで、Ｐ次の全極フィルタ係数ａ_rは、例えば乱数を用いて適当な値に初期値を設定する。Ｕ_r,tの各初期値の設定は、（システム構成その１）での説明と同様である。設定したパラメータａ_r及びＵ_r,tの各初期値は、記憶部３０に記憶される。

推定部２３は、（ｋ、ｒ）の全ての組み合わせの各々について、記憶部３０に記憶されているａ_rに基づいて、上記（５）式に従って、声道スペクトルＨ_k,r ^(AR)（以降、単に「Ｈ_k,r」と記載する）を計算し、記憶部３０に格納する。

推定部２３は、（ｋ、ｔ）の全ての組み合わせの各々について、記憶部３０に記憶されているＨ_k,r及びＵ_r,tに基づいて、上記（１）式に従って、声道スペクトログラムＸ_k,tを計算し、記憶部３０に格納する。

推定部２３は、記憶部３０に記憶されているａ_r、Ｚ_k,t、Ｈ_k,r、Ｘ_k,t、Ｙ_k,t、Ｕ_r,tに基づいて、上記（２０）式に従って、全極フィルタ係数ａ_rを更新し、記憶部３０に格納する。

推定部２３は、全極フィルタ係数ａ_rの更新に伴い、（ｋ、ｒ）の全ての組み合わせの各々について、ａ_rに基づいて、上記（５）式に従って、声道スペクトルＨ_k,rを計算し、記憶部３０に格納する。

推定部２３は、声道スペクトルＨ_k,rの更新に伴い、（ｋ、ｔ）の全ての組み合わせの各々について、記憶部３０に記憶されているＨ_k,r及びＵ_r,tに基づいて、上記（１）式に従って、声道スペクトログラムＸ_k,tを計算し、記憶部３０に格納する。

推定部２３は、（ｒ、ｔ）の全ての組み合わせの各々について、記憶部３０に記憶されているＵ_r,t、Ｚ_k,t、Ｘ_k,t、Ｙ_k,t、Ｈ^(AR) _k,r、すなわちＨ_k,rに基づいて、上記（１７）式に従って、声道スペクトログラムＸ_k,tのスペクトルパターンの重みＵ_r,tを更新し、記憶部３０に格納する。

更に、推定部２３は、声道スペクトログラムＸ_k,tのスペクトルパターンの重みＵ_r,tの更新に伴い、（ｋ、ｔ）の全ての組み合わせの各々について、記憶部３０に記憶されているＨ_k,r及びＵ_r,tに基づいて、上記（１）式に従って、声道スペクトログラムＸ_k,tを計算し、記憶部３０に格納する。

終了判定部２４は、予め定められた終了条件を満足するか否かを判定し、終了条件を満足していない場合には、推定部２３の各処理を繰り返す。終了判定部２４は、終了条件を満足したと判定した場合には、出力部２５による処理に移行する。

出力部２５は、記憶部３０に記憶されている声道スペクトログラムＸ_k,tを、入力部１０に入力された音声信号から推定される声道スペクトログラムとして出力する。

＜声道スペクトル推定装置の作用＞
次に、AR-NMFに対する声道スペクトルの推定アルゴリズムを利用する第２の実施の形態に係る声道スペクトル推定装置１００の作用について説明する。

マイクロホンで取得された音声信号の時系列データが声道スペクトル推定装置１００に入力され、記憶部３０に格納される。そして、声道スペクトル推定装置１００において、図３に示す声道スペクトル推定処理ルーチンが実行される。

図３に示す声道スペクトル推定処理ルーチンのステップＳ１００、Ｓ１０２、及びＳ１１８は、既に説明した図２の声道スペクトル推定処理ルーチンの対応するステップと同様であるため、説明を省略する。

ステップＳ１０５において、乱数を用いてａ_r及びＵ_r,tの初期値を設定する。この際、ａ_r及びＵ_r,tの初期値は非負値を設定する。

こうして設定されたａ_r及びＵ_r,tの各初期値は、記憶部３０に記憶される。

次に、ステップＳ１０７では、ステップＳ１０５で設定されたａ_rに基づいて、上記（５）式に従って、声道スペクトルＨ_k,rを各（ｋ、ｒ）の組み合わせについて算出して、記憶部３０に格納する。また、ステップＳ１０７では、ステップＳ１０５で設定されたＵ_r,t、及び本ステップで算出されたＨ_k,rに基づいて、上記（１）式に従って、声道スペクトログラムＸ_k,tを各（ｋ、ｔ）の組み合わせについて算出して、記憶部３０に格納する。

ステップＳ１０９では、ステップＳ１００で算出されたＹ_k,t、ステップＳ１０２で算出されたＺ_k,t、ステップＳ１０５で設定されたａ_r及びＵ_r,t、ステップＳ１０７で算出されたＨ_k,r及びＸ_k,t、すなわち、記憶部３０に記憶されている最新の各パラメータａ_r、Ｚ_k,t、Ｘ_k,t、Ｙ_k,t、Ｈ_k,r、Ｕ_r,tに基づいて、上記（２１）式に従って、全極フィルタ係数ａ_rを算出して、記憶部３０に格納する。

ステップＳ１１１では、記憶部３０に記憶されている最新のａ_r、すなわち、ステップＳ１０９で算出されたａ_rに基づいて、上記（５）式に従って、声道スペクトルＨ_k,rを各（ｋ、ｒ）の組み合わせについて算出して、記憶部３０に格納する。また、ステップＳ１１１では、記憶部３０に記憶されている最新の各パラメータＵ_r,t及びＨ_k,r、すなわち、ステップＳ１０５で設定されたＵ_r,t、及び本ステップで算出されたＨ_k,rに基づいて、上記（１）式に従って、声道スペクトログラムＸ_k,tを各（ｋ、ｔ）の組み合わせについて算出して、記憶部３０に格納する。

ステップＳ１１３では、記憶部３０に記憶されている最新の各パラメータＺ_k,t、Ｘ_k,t、Ｙ_k,t、Ｈ_k,r、Ｕ_r,t、すなわち、ステップＳ１０５で設定されたＵ_r,t、ステップＳ１０２で算出されたＺ_k,t、ステップＳ１１１で算出したＸ_k,t及びＨ_k,r、ステップＳ１００で算出したＹ_k,tに基づいて、上記（１７）式に従って、声道スペクトログラムＸ_k,tのスペクトルパターンの重みＵ_r,tを各（ｒ、ｔ）の組み合わせについて算出して、記憶部３０に格納する。

ステップＳ１１５では、記憶部３０に記憶されている最新の各パラメータＵ_r,t及びＨ_k,r、すなわち、ステップＳ１１３で算出されたＵ_r,t、及びステップＳ１１１で算出されたＨ_k,rに基づいて、上記（１）式に従って、声道スペクトログラムＸ_k,tを各（ｋ、ｔ）の組み合わせについて算出して、記憶部３０に格納する。

次のステップＳ１１７では、ステップＳ１００で算出したＹ_k,tと、ステップＳ１０２で算出したＺ_k,tと、ステップＳ１１５で算出したＸ_k,tに基づいて、（７）式に従って目的関数Ｌ_GKL(Θ)の値を算出して、記憶部３０に記憶する。そして、前回のステップＳ１１７で算出した目的関数Ｌ_GKL(Θ)の値を記憶部３０から読み込み、今回のステップＳ１１７で算出した目的関数Ｌ_GKL(Θ)の値と、前回のステップＳ１１７で算出した目的関数Ｌ_GKL(Θ)の値との差分が、予め記憶部３０に記憶されている予め定められた閾値よりも小さいか否かを判定し、差分が予め定められた閾値以上の場合には、終了条件を満足していないと判断して、上記ステップＳ１０９へ戻り、上記ステップＳ１０９〜ステップＳ１１７の処理を繰り返す。

一方、差分が予め定められた閾値未満の場合には、終了条件を満足したと判断して、ステップＳ１１８で、ステップＳ１１５で算出した最新の声道スペクトログラムＸ_k,tを出力して、声道スペクトル推定処理ルーチンを終了する。

なお、上記では、距離Ｄ_*(Ｙ_k,t;Ｘ_k,t)として、一般化KLダイバージェンスＤ_GKLを用いた場合のAR-NMFに対する声道スペクトルの推定アルゴリズムを利用する声道スペクトル推定装置の例について説明したが、距離Ｄ_*(Ｙ_k,t;Ｘ_k,t)として、２乗距離Ｄ_EUを用いてもよいことは言うまでもない。

この場合、例えばステップＳ１０９が、（２３）式に示したａ_rの更新式の計算に置き換わり、ステップＳ１１３が、（２２）式に示したＵ_r,tの更新式の計算に置き換わる。

＜声道スペクトル推定精度評価実験＞
＜評価実験の条件＞
次に、第１の実施の形態に係る声道スペクトルの推定方法（以降、「提案法」という）の有効性を示す目的で、提案法で推定した声道スペクトルと、STRAIGHTで推定した声道スペクトルの推定精度を比較する評価実験を行った。

ATRデジタル音声データベースのAセットから、日本人女性話者1名による20文の音声信号（サンプリング周波数は16kHz）をSTRAIGHTの手法で分析し、基本周波数F₀、声道スペクトル、非周期性指標Ａ_k,tを抽出した。なお、ここで得られたスペクトルを正解の声道スペクトルとみなす。

そして、正解の声道スペクトルと2^-1.0、2^-0.5、2^0.0、2^0.5、2^1.0、2^1.5倍したF₀を用いて、音声信号をそれぞれSTRAIGHTで再合成した。

そして、再合成音声信号からSTRAIGHT及び提案法で声道スペクトルを推定し、声道スペクトル推定値と正解の声道スペクトルとのメルケプストラム歪みを用いて、各々の手法の声道スペクトル推定性能を比較した。なお、メルケプストラム歪みは、1次から24次のメルケプストラム係数を用いて計算し、STRAIGHTによる声道スペクトルの推定では、フレームシフトを5ms(Ｔ＝81761)、声道スペクトルの次元をＫ＝513とした。

また、提案法として、GMM-NMFに対するパラメータ推定アルゴリズムにおいて、一般化KLダイバージェンスＤ_GKLと2乗距離Ｄ_EUを用いた場合の評価を行った。具体的には、STRAIGHTで推定された再合成音声信号の声道スペクトルを用い、F₀の各定数倍毎に20発話の声道スペクトログラムＹ_k,tを同時に用いて、Ｗ_r,n及びＵ_r,tを推定した。また、スペクトルパターンＲ＝90、混合数Ｎ＝100とし、ｎ＝0,・・,Ｎ−１に対して平均ρ_n＝ｎ／(Ｎ−１)、標準偏差ν_n＝１／(Ｎ−１)とした。Ｗ_r,n及びＵ_r,tは非負の乱数で初期化し、提案法における声道スペクトルの推定アルゴリズムの反復回数は100回とした。

＜評価実験結果＞
図４に、提案法とSTRAIGHTによる声道スペクトルの推定結果のメルケプストラム歪みを示す。なお、“GKL”の列は、Ｆ₀の倍率ｘに対する一般化KLダイバージェンスを用いた場合の評価結果、“EU”の列は、Ｆ₀の倍率ｘに対する2乗距離を用いた場合の評価結果、及び“STRAIGHT”の列は、STRAIGHTによる評価結果を表している。各々の評価結果は、［平均値±標準偏差］[dB]の形式で記載されており、括弧内の値は、非周期性指標Ａ_k,tを用いなかった場合の評価結果を示している。

Ｆ₀の倍率ｘが高くなるにしたがって、観測できる調波成分も少なくなるため、当該フレーム以外の調波成分を利用することによる効果が現れると考えられるが、図４に示すように、Ｆ₀の倍率ｘが高くなるほど、STRAIGHTを用いた評価結果に比べて、GKLの評価結果の方がメルケプストラム歪みが小さくなり、当該フレーム以外の調波成分が声道スペクトルの推定に有効であることが確認できる。

EUの評価結果も、Ｆ₀の倍率ｘが高くなるほど、STRAIGHTを用いた評価結果に比べてメルケプストラム歪みが小さくなる。しかし、GKLの評価結果に比べると、平均的にメルケプストラム歪みが大きくなる傾向が見られ、一般化KLダイバージェンスを用いたGMM-NMFに対するパラメータ推定アルゴリズムの方が、声道スペクトルの推定に適しているということができる。

また、非周期性指標Ａ_k,tを全ての時間周波数成分で一様、すなわち、全ての（ｋ、ｔ）の組み合わせに対してＺ_k,t＝1とした場合、それぞれの音声スペクトルの推定手法において、括弧内の値が取得される。

GKLの評価結果において、何れのＦ₀の倍率ｘについても、非周期性指標Ａ_k,tを用いた方が、非周期性指標Ａ_k,tを一様にした場合のメルケプストラム歪みより小さくなっていることから、非周期性指標Ａ_k,tが声道スペクトルの推定に関する性能向上に寄与することが確認できる。

また、GKLにおける括弧内の評価結果と、STRAIGHTを用いた評価結果を比べると、何れのＦ₀の倍率ｘについても、GKLにおけるメルケプストラム歪みが小さいことがわかる。したがって、声道スペクトログラムＸ_k,tが低ランクな非負値行列で近似できるという仮定が、声道スペクトル推定に有用であることが示唆される。

このように、本発明に係る提案手法では、音声信号に付与されたコンテキストラベルを用いることなく、音声信号の複数のフレームにおける調波成分の情報を手がかりにして、音声信号から声道スペクトルを精度良く推定することができる。

なお、本発明は、上述した実施形態に限定されるものではなく、この発明の要旨を逸脱しない範囲内で様々な変形や応用が可能である。

例えば、上述の声道スペクトル推定装置は、内部にコンピュータシステムを有しているが、「コンピュータシステム」は、ＷＷＷシステムを利用している場合であれば、ホームページ提供環境（あるいは表示環境）も含むものとする。

また、本願明細書中において、プログラムが予めインストールされている実施形態として説明したが、当該プログラムを、コンピュータ読み取り可能な記録媒体に格納して提供することも可能である。

１０ 入力部
２０ 演算部
２１ 観測スペクトログラム推定部
２２ 初期設定部
２３ 推定部
２４ 終了判定部
２５ 出力部
３０ 記憶部
１００ 声道スペクトル推定装置

Claims (7)

1. 音声信号の時系列データを、基本周期の幅で切り出し、各々切り出した音声信号のスペクトルから、各時刻及び各正規化角周波数の観測時間周波数成分を表す観測スペクトログラムを出力する観測スペクトログラム推定部と、
   前記観測スペクトログラム推定部により出力された前記観測スペクトログラムと、各スペクトルパターンにおける各正規化角周波数のパワースペクトルを表す声道スペクトル、及び各スペクトルパターンの各時刻における重みから求められる各時刻及び各正規化角周波数の時間周波数成分を表す声道スペクトログラムとの距離を用いて表される目的関数を小さくするように、各スペクトルパターンにおける各正規化角周波数の前記声道スペクトル、及び各スペクトルパターンの各時刻における重みを推定する推定部と、
   を含む声道スペクトル推定装置。
2. ｒ番目のスペクトルパターンにおけるｋ番目の正規化角周波数ω_kの前記声道スペクトルを、以下の式で表わされるＨ_k,r ^(GMM)とし、
   前記推定部は、ｒ番目のスペクトルパターンにおけるｋ番目の正規化角周波数ω_kの前記声道スペクトルとして、ｎ番目の正規分布の各々に対する重みＷ_r,nを推定する請求項１記載の声道スペクトル推定装置。
   
   ただし、Ｇ_n（ω）は、平均ρ_n、分散ν_n ²の正規分布を表し、ｈ（ω）は、周波数ワーピング関数を表す。
3. ｒ番目のスペクトルパターンにおけるｋ番目の正規化角周波数ω_kの前記声道スペクトルを、以下の式で表わされるＨ_k,r ^(AR)とし、
   前記推定部は、ｒ番目のスペクトルパターンにおけるｋ番目の正規化角周波数ω_kの前記声道スペクトルとして、Ｐ次の全極フィルタの係数ａ_rを推定する請求項１記載の声道スペクトル推定装置。
   
   ただし、Ｑ（ω）は、（ｐ，ｑ）成分が、ｃｏｓ（ω（ｐ−ｑ））で表わされる（Ｐ＋１）×（Ｐ＋１）のToeplitz行列である。
4. 観測スペクトログラム推定部と推定部とを含む声道スペクトル推定装置における声道スペクトル推定方法であって、
   前記観測スペクトログラム推定部が、音声信号の時系列データを、基本周期の幅で切り出し、各々切り出した音声信号のスペクトルから、各時刻及び各正規化角周波数の観測時間周波数成分を表す観測スペクトログラムを出力し、
   前記推定部が、前記観測スペクトログラム推定部により出力された前記観測スペクトログラムと、各スペクトルパターンにおける各正規化角周波数のパワースペクトルを表す声道スペクトル、及び各スペクトルパターンの各時刻における重みから求められる各時刻及び各正規化角周波数の時間周波数成分を表す声道スペクトログラムとの距離を用いて表される目的関数を小さくするように、各スペクトルパターンにおける各正規化角周波数の前記声道スペクトル、及び各スペクトルパターンの各時刻における重みを推定する
   声道スペクトル推定方法。
5. ｒ番目のスペクトルパターンにおけるｋ番目の正規化角周波数ω_kの前記声道スペクトルを、以下の式で表わされるＨ_k,r ^(GMM)とし、
   前記推定部が推定することでは、ｒ番目のスペクトルパターンにおけるｋ番目の正規化角周波数ω_kの前記声道スペクトルとして、ｎ番目の正規分布の各々に対する重みＷ_r,nを推定する請求項４記載の声道スペクトル推定方法。
   
   ただし、Ｇ_n（ω）は、平均ρ_n、分散ν_n ²の正規分布を表し、ｈ（ω）は、周波数ワーピング関数を表す。
   声道スペクトル推定方法。
6. ｒ番目のスペクトルパターンにおけるｋ番目の正規化角周波数ω_kの前記声道スペクトルを、以下の式で表わされるＨ_k,r ^(AR)とし、
   前記推定部が推定することでは、ｒ番目のスペクトルパターンにおけるｋ番目の正規化角周波数ω_kの前記声道スペクトルとして、Ｐ次の全極フィルタの係数ａ_rを推定する請求項４記載の声道スペクトル推定方法。
   
   ただし、Ｑ（ω）は、（ｐ，ｑ）成分が、ｃｏｓ（ω（ｐ−ｑ））で表わされる（Ｐ＋１）×（Ｐ＋１）のToeplitz行列である。
   声道スペクトル推定方法。
7. 請求項１〜請求項３の何れか１項に記載の声道スペクトル推定装置の各部としてコンピュータを機能させるためのプログラム。