JP5885210B2 - 基本周波数モデルパラメータ推定装置、方法、及びプログラム - Google Patents

Description

本発明は、基本周波数モデルパラメータ推定装置、方法、及びプログラムに係り、特に、音声信号から、観測基本周波数系列のパラメータを推定する基本周波数モデルパラメータ推定装置、方法、及びプログラムに関する。

＜藤崎モデル＞
音声のイントネーションを解析する手法に、藤崎の基本周波数（F₀）パターン生成過程モデル（藤崎モデル）が知られている（非特許文献１）。藤崎モデルとは、甲状軟骨の運動に注目してF₀パターンの生成過程を説明した、力学的モデルである。藤崎モデルでは、甲状軟骨の二つの独立な運動（平行移動運動と回転運動）にそれぞれ伴う声帯の伸びの合計がF₀の時間的変化をもたらすと解釈され、声帯の伸びとF₀パターンの対数値y(t)が比例関係にあるという仮定に基づいてF₀パターンがモデル化される。甲状軟骨の平行移動運動によって生じるF₀パターンy_p(t)をフレーズ成分、回転運動によって生じるF₀パターンy_a(t)をアクセント成分と呼ぶ。藤崎モデルでは、音声のF₀パターンy(t)は、これらの成分に声帯の物理的制約によって決まるベースライン成分y_bを足し合わせたものとして、

と表現される。これら二つの成分は二次の臨界制動系の出力であるとしてモデル化されており、

と表される（*は時刻tに関する畳み込み演算）。ここでu_p(t)はフレーズ指令関数と呼ばれ、デルタ関数（フレーズ指令）の列からなり、u_a(t)はアクセント指令関数と呼ばれ、矩形波（アクセント指令）の列からなる。これらの指令列には、発話の最初にはフレーズ指令が生起する、フレーズ指令は二連続で生起しない、異なる二つの指令（フレーズ指令とアクセント指令）は同時刻に生起しない、という制約条件がある。またαとβはそれぞれフレーズ制御機構、アクセント制御機構の固有角周波数であり、話者や発話内容によらず、おおよそα=3 rad/s、β=20 rad/s 程度であることが経験的に知られている。

＜藤崎モデルのパラメータ推定法１＞
従来、音声信号のF₀パターンから藤崎モデルのパラメータを推定する手法として、非特許文献２に記載の手法が知られている。この手法ではまず、観測F₀パターンに対して平滑化のための前処理を行う。具体的にはgross errorの除去、microprosodyの修正、短い無音区間と無声区間の補間を行った後、F₀パターンを至るところで連続かつ微分可能な区分的3次曲線で近似する。次に、そうして得られた平滑化F₀パターンの微分値の極大値・極小値を手掛かりに、アクセント指令列の位置と大きさを推定する。さらに観測F₀パターンから推定アクセント成分を差し引いたパターンをもとにleft-to-rightにフレーズ指令を挿入していく。最後に、推定指令列から生成したF₀パターンと観測F₀パターンの平均二乗誤差を最小とするように、指令列を微小変化させ、こうして得られた指令列を藤崎モデルの推定パラメータとする。

＜藤崎モデルのパラメータ推定法２＞
従来、音声信号のF₀パターンから藤崎モデルのパラメータを推定する手法として、他にも次のような手法がある（非特許文献３〜５）。この手法では、離散化した藤崎モデルをベースにして定式化したF₀パターン生成過程の確率モデルを使い、そのモデルに従ってP(y|θ)の最適化問題を解くことによって適切なパラメータを推定する（yは観測F₀パターン、θは藤崎モデルのパラメータ）。このモデル中では、制約条件があって扱いづらいフレーズ指令とアクセント指令のペアを、隠れマルコフモデル（HMM）から確率的に出力される値として扱う。また推定アルゴリズムにおいては、各成分について周辺化を行ったうえで、EMアルゴリズムによる反復解法で適切なパラメータを推定している。

Hiroya Fujisaki, Sumio Ohno and Wentao Gu, "Physiological and physical mechanisms for fundamental frequency control In some tone languages and a command-response model for generation of their F0 contours," Proceedings of International Symposium on Tonal Aspects of Languages: Emphasis on Tone Languages, Beijing, pp. 61-64 (2004-3). S. Narusawa, N. Minematsu, K. Hirose, and H. Fujisaki, "A method for automatic extraction of model parameters from fundamental frequency contours of speech," in Proc. ICASSP, 2002, pp. 509−512. H. Kameoka, J. L. Roux, and Y. Ohishi, "A statistical model of speech F0contours," in Proc.SAPA, 2010, pp. 43−48. 吉里幸太, 亀岡弘和, 齋藤大輔, 嵯峨山茂樹,"F0パターン生成過程の統計的モデルによる音声信号からのフレーズ・アクセント指令の推定," 日本音響学会春季研究発表会講演集, 2012, no. 1-11-9,

pp.311-314.
K. Yoshizato, H. Kameoka, D. Saito, and S. Sagayama, "Statistical approach to fujisaki-model parameter estimation from speech signals and its quantitative evaluation," in Proc. Speech

Prosody 2012, 2012, pp. 175−178.

本発明は、音声のF₀パターンから藤崎モデルのパラメータを推定する方法に関するものである。

この推定問題は不良設定の逆問題であるため、解析的に解くことは困難である。ここで、日本語のような非声調言語においては、フレーズ指令とアクセント指令の大きさは非負でなければならないという制約がある。この非負性は解を絞り込むための重要な制約になるにも関わらず、上記非特許文献３〜５に記載の藤崎モデルのパラメータ推定法では、この制約を最適化問題として直接導入することができなかった。上記の非特許文献５に記載の技術では、フレーズ・アクセント成分から指令列を逆算する際に非負制約付き逆畳み込み問題を解くというアドホックな方法で解決を図ったが、アルゴリズムの収束性が保証されなくなり、推定パラメータから生成されたF₀パターンと観測F₀パターンとの間の誤差が大きいという問題があった。この「誤差が大きい」という問題は、上記の非特許文献２に記載の、藤崎モデルのパラメータ推定法にも見られる問題である。

本発明は、上記の事情を鑑みてなされたもので、フレーズ指令とアクセント指令の非負性に関する制約を用いて、藤崎モデルのパラメータを推定することができる基本周波数モデルパラメータ装置、方法、及びプログラムを提供することを目的とする。

上記の目的を達成するために本発明に係る基本周波数モデルパラメータ推定装置は、音声信号を入力として、隠れマルコフモデルの各時刻ｋの状態ｓ_kからなる指令状態系列＾ｓと、各時刻ｋにおける甲状軟骨の平行移動運動によって生じる基本周波数パターンを表すフレーズ指令ｕ_p［ｋ］及び甲状軟骨の回転運動によって生じる基本周波数パターンを表すアクセント指令ｕ_a［ｋ］のペア＾ｏ[ｋ]からなる指令関数＾ｏと、各時刻ｋにおける状態ｓ_kに応じたフレーズ指令の振幅Ａ_p[ｋ]及び各アクセント指令ｎの振幅Ａ_a ⁽ⁿ⁾を表すパラメータ群θとを推定する基本周波数モデルパラメータ推定装置であって、前記音声信号の時系列データから、前記音声信号の各時刻ｋの基本周波数を表す観測基本周波数系列＾ｙを抽出する基本周波数抽出手段と、前記音声信号の時系列データについて、有声区間及び無声区間の何れであるかに応じて、各時刻ｋにおける前記基本周波数の不確かさの度合いを推定する有声無声区間推定手段と、前記指令関数＾ｏの初期値、及び前記パラメータ群θの初期値を設定する初期値設定手段と、前回更新された前記指令関数＾ｏ’または前記指令関数＾ｏの初期値＾ｏ’に基づいて、前記観測基本周波数系列＾ｙ、前記指令関数＾ｏ’、及び前記パラメータ群θ’が与えられたときの指令状態系列＾ｓの事後確率Ｐ（＾ｓ｜＾ｙ，＾ｏ’、θ’）を計算する指令状態系列事後確率更新手段と、前回更新された前記指令関数＾ｏ’または前記指令関数＾ｏの初期値＾ｏ’、前記観測基本周波数系列＾ｙ、各時刻ｋにおける前記不確かさの度合い、及び前記事後確率Ｐ（＾ｓ｜＾ｙ，＾ｏ’、θ’）に基づいて、前記観測基本周波数系列＾ｙが与えられたときの前記指令関数＾ｏ及び前記パラメータ群θの対数事後確率ｌｏｇＰ（＾ｏ，θ｜＾ｙ）を目的関数として、前記目的関数を増加させるように、各々非負値である前記指令関数＾ｏ、及び前記パラメータ群θを更新するモデルパラメータ更新手段と、予め定められた収束条件を満たすまで、前記指令状態系列事後確率更新手段による計算、及び前記モデルパラメータ更新手段による更新を繰り返し行う第１収束判定手段と、を含んで構成されている。

本発明に係る基本周波数モデルパラメータ推定方法は、音声信号を入力として、隠れマルコフモデルの各時刻ｋの状態ｓ_kからなる指令状態系列＾ｓと、各時刻ｋにおける甲状軟骨の平行移動運動によって生じる基本周波数パターンを表すフレーズ指令ｕ_p［ｋ］及び甲状軟骨の回転運動によって生じる基本周波数パターンを表すアクセント指令ｕ_a［ｋ］のペア＾ｏ[ｋ]からなる指令関数＾ｏと、各時刻ｋにおける状態ｓ_kに応じたフレーズ指令の振幅Ａ_p[ｋ]及び各アクセント指令ｎの振幅Ａ_a ⁽ⁿ⁾を表すパラメータ群θとを推定する基本周波数モデルパラメータ推定方法であって、基本周波数抽出手段によって、前記音声信号の時系列データから、前記音声信号の各時刻ｋの基本周波数を表す観測基本周波数系列＾ｙを抽出し、有声無声区間推定手段によって、前記音声信号の時系列データについて、有声区間及び無声区間の何れであるかに応じて、各時刻ｋにおける前記基本周波数の不確かさの度合いを推定し、初期値設定手段によって、前記指令関数＾ｏの初期値、及び前記パラメータ群θの初期値を設定し、指令状態系列事後確率更新手段によって、前回更新された前記指令関数＾ｏ’または前記指令関数＾ｏの初期値＾ｏ’に基づいて、前記観測基本周波数系列＾ｙ、前記指令関数＾ｏ’、及び前記パラメータ群θ’が与えられたときの指令状態系列＾ｓの事後確率Ｐ（＾ｓ｜＾ｙ，＾ｏ’、θ’）を計算し、モデルパラメータ更新手段によって、前回更新された前記指令関数＾ｏ’または前記指令関数＾ｏの初期値＾ｏ’、前記観測基本周波数系列＾ｙ、各時刻ｋにおける前記不確かさの度合い、及び前記事後確率Ｐ（＾ｓ｜＾ｙ，＾ｏ’、θ’）に基づいて、前記観測基本周波数系列＾ｙが与えられたときの前記指令関数＾ｏ及び前記パラメータ群θの対数事後確率ｌｏｇＰ（＾ｏ，θ｜＾ｙ）を目的関数として、前記目的関数を増加させるように、各々非負値である前記指令関数＾ｏ、及び前記パラメータ群θを更新し、第１収束判定手段によって、予め定められた収束条件を満たすまで、前記指令状態系列事後確率更新手段による計算、及び前記モデルパラメータ更新手段による更新を繰り返し行う。

本発明に係るプログラムは、上記の基本周波数モデルパラメータ推定装置の各手段としてコンピュータを機能させるためのプログラムである。

以上説明したように、本発明の基本周波数モデルパラメータ推定装置、方法、及びプログラムによれば、観測基本周波数系列＾ｙが与えられたときの指令関数＾ｏ及びパラメータ群θの対数事後確率ｌｏｇＰ（＾ｏ，θ｜＾ｙ）を目的関数として、目的関数を増加させるように、各々非負値である指令関数＾ｏ、及びパラメータ群θを更新することにより、フレーズ指令とアクセント指令の非負性に関する制約を用いて、藤崎モデルのパラメータを推定することができる、という効果が得られる。

ＨＭＭを説明するための図である。 状態の分割を説明するための図である。 本発明の実施の形態に係る基本周波数モデルパラメータ推定装置の構成を示す概略図である。 本発明の実施の形態に係る基本周波数モデルパラメータ推定装置における基本周波数モデルパラメータ推定処理ルーチンの内容を示すフローチャートである。 指令関数のマッチングを説明するための図である。 実験結果を示す図である。 実験結果を示す図である。

以下、図面を参照して本発明の実施の形態を詳細に説明する。本発明で提案する手法では、観測F₀パターンの再現性が高い藤崎モデルのパラメータ推定を実現するために、藤崎モデルをベースにしたF₀パターン生成過程の確率モデルを定式化し、それに基づいて観測F₀パターンが生じたと仮定する。藤崎モデルのパラメータ推定アルゴリズムも、この確率モデルに基づく。

＜F₀パターン生成過程の確率モデル＞
本発明の原理について説明する。まず、F₀パターン生成過程の確率モデルについて説明する。

指令関数に付随する各種制約をモデルに組み入れるために、フレーズ指令、アクセント指令のペア＾o[k]=(u_p[k],u_a[k])^Tを出力するHMMを考える。このHMMでは、出力される指令関数{＾o[k]}_k=1 ^Kは、各時刻ごとにガウス分布に従い、

と確率的に表現される。ここで{s_k}_k=1 ^KはHMMの状態系列であり、平均ベクトル＾ν_sk[k]と分散共分散行列＾Υ_skはHMMの状態遷移の結果として定まる値である。具体的なHMMの構成は図１に示した。なお、行列又はベクトルを示す記号については、「＾」を付すこととする。

加えて、自己遷移の持続長をパラメータ化するために、それぞれの状態をいくつかの小状態に分割することを考える。なおこのとき、各々の小状態は全て同じ出力分布を持ち、小状態の数は十分大きな値となるようにしておく。図２に状態a_nを分割した例を示した。例えばこの図２のように全てのm≠0に対してa_n,mからa_n,m+1への状態遷移確率を1に設定することで、a_n,0からa_n,mへの遷移確率が状態a_nがmステップだけ持続する確率に対応し、アクセント指令の持続長を柔軟に制御できるようになる。同様にp₁とp₀とa₀も小状態に分割することで、フレーズ指令の持続長と指令間の間隔の長さの分布をパラメータ化することが可能になる。こうした分割をふまえて、以後は改めてp₀={p_0,0,p_0,1,...}、a₀={a_0,0,a_0,1,...}、a_n={a_n,0,a_n,1,...}と表記する。

提案するHMMの構成を定式化すると次のように書ける。

状態系列＾s={s_k}_k=1 ^Kが与えられたとき、このHMMはフレーズ指令関数u_p[k]とアクセント指令関数u_a[k]のペアを出力する。式(2)と式(4)で示した通り、u_p[k]とu_a[k]はそれぞれG_p[k]とG_a[k]というフィルタに畳み込まれてフレーズ成分x_p[k]とアクセント成分x_a[k]が出力される。これを式で表すと、

と書ける（*は離散時刻kに関する畳み込み演算）。このとき、F₀パターンx[k]は

と三種類の成分の重ね合わせで書ける。ただしu_bは時刻によらないベースライン成分である。

また、実音声においては、常に信頼のできる基本周波数F₀の値が観測できるとは限らない。例えば音声データからピッチ抽出によって得られた基本周波数F₀の推定値は、無声区間においては全く信頼できない値である。藤崎モデルのパラメータ推定を行うにあたっては、信頼のおける観測区間のF₀値のみを考慮に入れて、そうでない区間は無視することが望ましい。そこで、提案モデルに観測F₀値の時刻kにおける不確かさの程度v_n ²[k]を導入する。具体的には、観測F₀値y[k]を、真のF₀値x[k]とノイズ成分x_n[k]〜N(0,v_n ²[k])との重ね合わせで

と表現することで、信頼のおける区間かどうかに関わらず全ての観測区間を統一的に扱える。

φ_i′,i、u_b、v_p,sk ²、v_a,sk ²、v_n ²[k]、α、βを定数とみなし、{A_p[k]}_k=1 ^K,{A_a ⁽ⁿ⁾}_n=1 ^Nは一様に分布すると仮定する。するとx_n[k]を周辺化することで、出力値系列＾o={o[k]}_k=1 ^Kが与えられたときの＾y={y[k]}_k=1 ^Kの確率密度関数

が得られる。状態系列＾s = {s_k}_k=1 ^Kと指令の振幅を表すパラメータ群θ={{A_p[k]}_k=1 ^K,{A_a ⁽ⁿ⁾} _n=1 ^N}が与えられたとき、出力値系列＾oは

に従って生成される。また、P(＾s)は状態遷移確率の積として

と書ける。ただしφ_s1は初期状態がs₁である確率をあらわす。

＜藤崎モデルのパラメータ推定アルゴリズム＞
本発明では、藤崎モデルのパラメータ推定問題を、観測基本周波数系列＾yが与えられたときのパラメータ＾o,θの事後確率P(＾o,θ|＾y)を最大化する＾oとθを求める最大事後確率（Maximum A Posteriori;MAP）推定問題として定式化し、指令状態系列＾sを潜在変数と見なしてExpectation-Maximization(EM)アルゴリズムに基づく反復計算により＾oとθの局所最適解を探索する。EMアルゴリズムは、パラメータの対数事後確率ｌｏｇP(＾o,θ|＾y)の下限関数（Q関数と呼ぶ）を反復的に増加させることで、間接的にパラメータの対数事後確率ｌｏｇP(＾o,θ|＾y)を増加させていく方法であり、本問題におけるQ関数は

と書ける。ここで^c=は定数部分を除いて一致することを意味し、＾o′とθ′はそれぞれ＾oとθの直前の反復における値である。

＜システム構成＞
次に、観測された音声信号の時系列データを解析して、藤崎モデルのパラメータを推定する基本周波数モデルパラメータ推定装置に、本発明を適用した場合を例にして、本発明の実施の形態を説明する。

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

図３に示すように、基本周波数モデルパラメータ推定装置は、記憶部１と、基本周波数系列抽出部２と、有声無声区間推定部３と、初期値設定部４と、指令状態系列事後確率更新部５と、モデルパラメータ更新部６と、収束判定部７と、状態系列算出部８と、出力部９とを備えている。

記憶部１は、観測された音声信号の時系列データを記憶する。

基本周波数系列抽出部２は、音声信号の時系列データから、基本周波数の時系列データを抽出し、それらを離散時間ｋで表現するように変換して、音声信号の基本周波数の時系列データである観測基本周波数系列＾ｙ＝｛Ｆ₀［ｋ］｝（ｋ＝１,…,Ｋ）とする。この基本周波数の抽出処理は、周知技術により実現でき、例えば、非特許文献６（H. Kameoka, "Statistical speech spectrum model incorporating all-pole vocal tract model and F0 contour generating process model," in Tech. Rep. IEICE, 2010, in Japanese.）に記載の手法を利用して、８ｍｓごとに基本周波数を抽出する。

有声無声区間推定部３は、音声信号の時系列データから、有声区間と無声区間とを特定し、離散時間ｋ毎に、有声区間であるか無声区間であるかに応じて、観測Ｆ0［ｋ］値の不確かさの程度v_n ²[k]を推定する。無声区間では不確かさの程度を大きく推定し（例えば、v_n ²[k]=10¹⁵）、有声区間では不確かさの程度を小さく推定する（例えば、v_n ²[k]=0.22）。

初期値設定部４は、後述する処理で用いる各パラメータである、アクセント指令の数N、EMアルゴリズムの反復回数M、α、β、v_p ²[k]、v_a ²[k]、u_bを定数とみなし初期値を設定する。初期値として適当な値を設定する。また、初期値設定部４は、HMMの小状態の個数、遷移確率φ_i′,Iを、予め用意した正解データから学習して決定する。また、初期値設定部４は、上記非特許文献２に記載の藤崎モデルのパラメータ推定法を用いて、＾oの初期値（非負値）を設定する。また、初期値設定部４は、A_p[k]の初期値として、＾oのフレーズ指令関数の振幅を線形補間したものを設定し、A_a ⁽ⁿ⁾の初期値として適切な値を設定する。

本実施の形態では、上記の式(15)のQ関数にもとづき、藤崎モデルパラメータ＾oとθの局所最適解は、指令状態系列事後確率更新部５とモデルパラメータ更新部６の2つのステップを繰り返すことで得られる。

指令状態系列事後確率更新部５は、指令状態系列（潜在変数）の事後確率P(＾s|＾y,＾o′,θ′)を計算するステップであり、EMアルゴリズムではこれをEステップと呼ぶ。Forward-Backwardアルゴリズムを用いれば各k,tに対してP(s_k=t|＾y,＾o′,θ′)を効率的に求めることができる。具体的には、

と変形すると、各P(s_k=t|{y[l],＾o′[l],θ′[l]}_l=1 ^l=k)は、

という漸化式を順次（k=1,2,...,K）解くことによって計算でき、各P({y[l],o′[l],θ′[l]}_l=k+1 ^l=K|s_k=t)は、

という漸化式を順次（k=K,K−1,...,1）解くことによって計算できる。

このように、指令状態系列事後確率更新部５は、時刻ｋ、状態ｔの全ての組み合わせ（k,t）の各々に対して、前回更新された指令関数＾o′又は初期値＾o′に基づいて、事後確率P(s_k=t|＾y,＾o′,θ′)を算出することにより、観測基本周波数系列＾ｙ、指令関数＾ｏ’、及びパラメータ群θ’が与えられたときの指令状態系列＾ｓの事後確率P(＾s|＾y,＾o′,θ′)を計算する。

モデルパラメータ更新部６は、補助変数更新部６１、指令関数更新部６２、収束判定部６３、及び平均振幅更新部６４を備えている。

モデルパラメータ更新部６は、目的関数Q(＾o,θ,＾o′,θ′)を増加させるように、非負値である指令関数＾oとパラメータ群θを更新するステップであり、EMアルゴリズムではこれをMステップと呼ぶ。logP(＾y|＾o,θ)の項は、

と書ける。ただしG_b[k]=δ[k]（クロネッカーのデルタ）である。指令関数u_p[k],u_a[k]が非負であるという条件で式(21)を最大化する＾oを直接求めるのは難しいが、補助関数法に基づく反復計算により式(21)を局所的に最大化する＾oを求めることができる。補助関数法はEMアルゴリズムと同様最大化したい目的関数の下限関数を反復的に増加させていくことで目的関数を増加させる手法であるが、式(21)の下限関数は、ジェンセンの不等式

が成り立つことを利用して設計することができる。ただし、λ_i,k,l≧0を補助変数と呼び、Σ_iΣ_lλ_i,k,l=1を満たす。式(22)の等号成立条件は

である。また、Σ_^sP(^s|^y,^o′,θ′)logP(θ,^o|^s)の項は、

と書ける。よって、不等式

の右辺(補助関数と呼ぶ)は、Q(^o,θ,^o′,θ′)の下限関数となり、補助関数Q′(^o,θ,^o′,θ′)とする。補助関数Q′(^o,θ,^o′,θ′)をu_i[l]で微分すると、

となる。ゆえにMステップにおいてu_i[l]を求めるには、補助関数法の更新式

を用いてλ_i,k,lとu_i[l]を交互に更新することを十分な回数繰り返せばよい。

このように、補助変数更新部６１は、前回更新された各時刻ｋのフレーズ指令ｕ_p［ｋ］（又は初期値）に基づいて、時刻ｋ、ｌ（ｌ＜ｋ）の全ての組み合わせ（ｋ、ｌ）の各々について、上記の式（２８）に従って、補助変数λ_p,k,lを算出して更新する。また、補助変数更新部６１は、前回更新された各時刻ｋのアクセント指令ｕ_a［ｋ］（又は初期値）に基づいて、（ｋ、ｌ）の全ての組み合わせについて、上記の式（２８）に従って、補助変数λ_a,k,lを算出して更新する。

また、補助変数更新部６１は、ｕ_bに基づいて、（ｋ、ｌ）の全ての組み合わせについて、上記の式（２８）に従って、補助変数λ_b,k,lを算出して更新する。

指令関数更新部６２は、基本周波数系列＾ｙと、不確かさの程度v_n ² [k]と、指令状態系列事後確率更新部５によって更新された指令状態系列の事後確率P(＾s|＾y,＾o′,θ′)と、補助変数更新部６１によって更新された補助変数λ_p,k,lとに基づいて、上記式（２９）に従って、非負値である各時刻ｌのフレーズ指令ｕ_p［ｌ］を更新する。

また、指令関数更新部６２は、基本周波数系列＾ｙと、不確かさの程度v_n ² [k]と、指令状態系列事後確率更新部５によって更新された指令状態系列の事後確率P(＾s|＾y,＾o′,θ′)と、補助変数更新部６１によって更新された補助変数λ_a,k,lとに基づいて、上記式（２９）に従って、非負値である各時刻ｌのアクセント指令ｕ_a［ｌ］を更新する。

また、指令関数更新部６２は、基本周波数系列＾ｙと、不確かさの程度v_n ² [k]と、指令状態系列事後確率更新部５によって更新された指令状態系列の事後確率P(＾s|＾y,＾o′,θ′)と、補助変数更新部６１によって更新された補助変数λ_b,k,lとに基づいて、上記式（２９）に従って、ベース成分ｕ_bを更新する。

収束判定部６３は、予め定められた収束条件を満足するか否かを判定し、収束条件を満足していない場合には、補助変数更新部６１及び指令関数更新部６２の各処理を繰り返す。収束判定部６３は、収束条件を満足したと判定した場合には、平均振幅更新部６４による処理に移行する。

収束条件としては、繰り返し回数ｓが予め定めた回数Ｓ（例えば、２０回）に達したことを用いればよい。なお、s-1回目のパラメータを用いたときの補助関数の値とs回目のパラメータを用いたときの補助関数の値との差が、予め定めた閾値よりも小さくなったことを、収束条件として用いてもよい。

平均振幅更新部６４は、Mステップとして、続けてθ={{A_p[k]}_k=1 ^K,{A_a ⁽ⁿ⁾}_n=1 ^N}を更新する。Q′(＾o,θ,＾o′,θ′)をA_p[k]とA_a ⁽ⁿ⁾で微分すると、

となる。ゆえにMステップにおいてA_p[k]とA_a ⁽ⁿ⁾を求める更新式は、

と書ける。

このように、平均振幅更新部６４は、指令関数更新部６２によって更新された各時刻ｋのフレーズ指令ｕ_p［ｋ］に基づいて、上記式（３２）に従って、各時刻ｋのフレーズ指令の振幅A_p[k]を更新すると共に、指令関数更新部６２によって更新された各時刻ｋのアクセント指令ｕ_a［ｋ］と、指令状態系列事後確率更新部５によって更新された指令状態系列の事後確率P(＾s|＾y,＾o′,θ′)とに基づいて、上記式（３２）に従って、各アクセント指令ｎの振幅A_a ⁽ⁿ⁾を更新することにより、パラメータ群θを更新する。

収束判定部７は、予め定められた収束条件を満足するか否かを判定し、収束条件を満足していない場合には、上記の更新値を改めてo′とθ′に代入して、反復アルゴリズム（指令状態系列事後確率更新部５及びモデルパラメータ更新部６の各処理）を繰り返す。収束判定部７は、収束条件を満足したと判定した場合には、状態系列算出部８による処理に移行する。

収束条件としては、繰り返し回数ｒが予め定めた回数Ｒ（例えば、２０回）に達したことを用いればよい。なお、ｒ-1回目のパラメータを用いたときの目的関数の値とｒ回目のパラメータを用いたときの目的関数の値との差が、予め定めた閾値よりも小さくなったことを、収束条件として用いてもよい。

状態系列算出部８は、最後に、Viterbiアルゴリズムを用いることで最適な状態系列^s^*を求める。具体的には、

という漸化式を順次（k=1,2,...,K）解くことによって求めたδ_t[k]とψ_t[k]を用いて、

という漸化式を順次（k=K,K−1,...,1）解くことによって計算できる。

このように、状態系列算出部８は、モデルパラメータ更新部６によって最終的に更新された指令関数＾ｏに基づいて、上記式（３３）〜式（３７）式に従って、状態系列＾ｓを算出する。そして、出力部９により、指令関数＾ｏ、パラメータ群θ、状態系列＾ｓを出力する。

＜基本周波数モデルパラメータ推定装置の作用＞
次に、本実施の形態に係る基本周波数モデルパラメータ推定装置１００の作用について説明する。まず、分析対象として、観測された音声信号の時系列データが、基本周波数モデルパラメータ推定装置１００に入力され、記憶部１に格納される。そして、基本周波数モデルパラメータ推定装置１００において、図４に示す基本周波数モデルパラメータ推定処理ルーチンが実行される。

まず、ステップＳ１０１において、記憶部１から、音声信号の時系列データを読み込み、各時刻ｋの基本周波数Ｆ₀からなる基本周波数系列ｙを抽出する。ステップＳ１０２において、音声信号の時系列データに基づいて、有声区間、無声区間を特定し、各時刻ｋの基本周波数の不確かさの程度v_n ² [ｋ]を推定する。

次のステップＳ１０３では、各パラメータN、M、α、β、v_p ²[k]、v_a ²[k]、u_bに対して適切な初期値を設定すると共に、HMMの小状態の個数、遷移確率φ_i′,Iを、予め用意した正解データから学習して決定する。また、従来手法により指令系列＾oを推定して、初期値として設定すると共に、A_p[k]の初期値及びA_a ⁽ⁿ⁾の初期値を設定する。

そして、ステップＳ１０４において、上記ステップＳ１０３で設定された指令系列＾oの初期値、または後述するステップＳ１０５で前回更新された指令系列＾oに基づいて、（k,t）の全ての組み合わせについて、事後確率P(s_k=t|＾y,＾o′,θ′)を更新することにより、指令状態系列の事後確率P(＾s|＾y,＾o′,θ′)を更新する。

ステップＳ１０５では、上記ステップＳ１０３で設定された指令系列＾oの初期値、または当該ステップＳ１０５で前回更新された指令系列＾oと、上記ステップＳ１０１で算出された基本周波数系列＾ｙと、上記ステップＳ１０２で算出された各時刻ｋの不確かさの程度v_n ²[ｋ]と、上記ステップＳ１０４で更新された指令状態系列の事後確率P(＾s|＾y,＾o′,θ′)とに基づいて、目的関数Q(＾o,θ,＾o′,θ′)を増加させるように、指令系列＾oと指令の振幅を表すパラメータ群θとを更新する

上記ステップＳ１０５は、以下のステップＳ１１１〜Ｓ１１４の各処理によって実現される。

ステップＳ１１１では、上記ステップＳ１０３で設定された指令系列＾oの初期値、または後述するステップＳ１１２で前回更新された指令系列＾oに基づいて、（ｋ、ｌ）の全ての組み合わせについて、上記の式（２８）に従って、補助変数λ_p,k,l、λ_a,k,l、λ_b,k,lを算出して更新する。

次のステップＳ１１２では、上記ステップＳ１０１で算出された基本周波数系列＾ｙと、上記ステップＳ１０２で算出された各時刻ｋの不確かさの程度v_n ² [ｋ]と、上記ステップＳ１０４で更新された指令状態系列の事後確率P(＾s|＾y,＾o′,θ′)と、上記ステップＳ１１１で更新された補助変数λ_p,k,l、λ_a,k,l、λ_b,k,lとに基づいて、上記式（２９）に従って、非負値である各時刻ｌのフレーズ指令ｕ_p［ｌ］及びアクセント指令ｕ_a［ｌ］からなる指令系列＾oとベース成分ｕ_bとを更新する。

次のステップＳ１１３では、収束条件として、繰り返し回数ｓが、Ｓに到達したか否かを判定し、繰り返し回数ｓがＳに到達していない場合には、収束条件を満足していないと判断して、上記ステップＳ１１１へ戻り、上記ステップＳ１１１〜ステップＳ１１２の処理を繰り返す。一方、繰り返し回数ｓがＳに到達した場合には、収束条件を満足したと判断し、ステップＳ１１４で、上記ステップＳ１１２で更新された各時刻ｋのフレーズ指令ｕ_p［ｋ］及びアクセント指令ｕ_a［ｋ］と、上記ステップＳ１０４で更新された指令状態系列の事後確率P(＾s|＾y,＾o′,θ′)とに基づいて、上記式（３２）に従って、各時刻ｋのフレーズ指令の振幅A_p[k]、及び各位置ｎのアクセント指令の振幅A_a[k]を更新することにより、パラメータ群θを更新する。

そして、ステップＳ１０６において、収束条件として、繰り返し回数ｒが、Ｒに到達したか否かを判定し、繰り返し回数ｒがＲに到達していない場合には、収束条件を満足していないと判断して、ステップＳ１０７で、上記ステップＳ１０５で更新された指令関数＾o,パラメータ群θを、＾o′,θ′に代入して、上記ステップＳ１０４へ戻り、上記ステップＳ１０４〜ステップＳ１０５の処理を繰り返す。一方、繰り返し回数ｒがＲに到達した場合には、収束条件を満足したと判断し、ステップＳ１０８で、上記ステップＳ１０５で最終的に更新された指令関数＾ｏに基づいて、上記式（３３）〜式（３７）式に従って、状態系列＾ｓを算出し、出力部９により、指令関数＾ｏ、指令の振幅を表すパラメータ群θ、状態系列＾ｓを出力して、基本周波数モデルパラメータ推定処理ルーチンを終了する。

＜実験＞
本実施の形態における重要な成果は、藤崎モデルを確率モデルとして表現することに成功したことである。本発明者らは、数多くの統計的手法に基づく音声アプリケーションに、本実施の形態で提案したモデルを組み込むことによって、将来的には韻律を扱う強力な手法が得られると考えている。そのためには、スペクトル特徴量と同じようにして、藤崎モデルのパラメータであるフレーズ、アクセント指令関数が音声コーパスから自動的に学習できると非常に便利である。この点において、確率モデルとして定式化した、本実施の形態の提案モデルと提案アルゴリズムは、たとえば、上記非特許文献２のような統計的でない手法よりも優れていると言える。しかし、提案アルゴリズムを用いた実音声からの藤崎モデルパラメータの推定性能が、既存手法の性能を上回っているかどうかはまだ明らかでない。そこで、本実施の形態で提案した手法のパラメータ推定性能を定量的に評価するための実験を行った。

詳しい実験条件を以下に記す。本実験で実音声データとして用いたのは、ATR 日本語音声データベースのB セット（非特許文献７（A. Kurematsu, K. Takeda, Y. Sagisaka, S. Katagiri, H. Kuwabara, and K. Shikano, "ATR japanese speech database as a tool of speech recognition and synthesis," Speech Communication, vol. 27, pp. 187-207, 1999.）を参照）である。これは503 文の音素バランス文からなる音声データベースであり、その中から一人の男性話者(MHT) を選択した。また、その音声データベースに対して、ある韻律研究の専門家が手動で求めたフレーズ、アクセント指令関数を正解データとして用いた。提案手法の入力として与える観測F₀ パターンを音声データから抽出する手法には、本発明者らが以前提案した上記非特許文献６に記載のアルゴリズムを用いた。

定数のパラメータについては、N=10、離散時刻のサンプリング間隔t₀=8ms、EMアルゴリズムの反復回数M=20回、α=3.0 rad/s、β=20.0 rad/s、v_p ²[k]=0.2²、v_a ²[k]=0.1²、無声区間ではv_n ²[k]=10¹⁵、有声区間ではv_n ²[k]=0.22、そしてu_bは有声区間におけるlogF₀の最小値に、それぞれ設定した。HMMの小状態の個数や遷移確率φ_i′,iについては、ATR日本語音声データベースのBセットのNo.1からNo.200までの200文の正解データから学習して決定した。＾oの初期値を、上記非特許文献２に記載した方法を用いて定めた。A_p[k]の初期値は＾oのフレーズ指令関数の振幅を線形補間したものとし、A_a ⁽ⁿ⁾の初期値を0.1nとした。

パラメータ推定実験は、No.201 からNo.503 までの303 文を対象にして行った。推定パラメータを評価する方法として、観測F₀ パターンの再現性と、言語学的な妥当性の二つを考慮した。これらは一般にトレードオフの関係にある。例えば、短い区間に細かく大量の指令を立てれば観測F₀ パターンを非常によく再現することができるが、そうして作った指令関数は言語学的に妥当なものであるとは言えない。そこで本実験では、本実施の形態で提案した手法によって得られた推定パラメータが言語学的に十分妥当なものでありつつ、観測F₀ パターンの再現性が非常に高いことを確認することを目的とする。

観測F₀パターンの再現性の評価基準には、観測F₀ パターンと推定指令関数から再構成されたF₀ パターンとの平均二乗誤差（log F₀[Hz] RMSE）を用い、この値が小さいほど再現性が高いとした。言語学的な妥当性の評価基準には、検出率という値を用い、これが大きいほど言語学的に妥当なパラメータであるとした。検出率は以下のように定義される。図５に例を示したように、推定パラメータ列と正解パラメータ列を比較して、指令単位でのマッチングをとる。指令と指令のマッチングがとれる条件は、二つの指令が同種の指令であること（フレーズ指令同士またはアクセント指令同士）と、二つの指令の時間のずれがS = 0.3 秒以下であることとした。ただし、アクセント指令に関しては生起時刻と終了時刻の平均を基準にした。また、二つのマッチングは時刻に関して交差していてはならない。マッチングがとれた指令同士の距離を1、そうでないときの距離を0 として、これらの条件を満たしなおかつ距離最大になるようなマッチングは、動的計画法によって求めることができる。推定実験に用いた303 文全てに対してこのマッチングをとったとき、マッチングの総数をN_Mとする。また、推定パラメータ列における指令の総数をN_E、正解パラメータ列における指令の総数をN_A とおく。ここで、挿入エラーE_Iを(N_E ‐N_M)/N_Aと定義し、脱落エラーE_D を(N_A‐N_M)/N_A と定義し、最終的な検出率Dは1‐E_I‐E_Dであると定義した。なお、この検出率の定義では指令の振幅を考慮に入れていない。これは、フレーズ、アクセント指令の振幅はベースライン成分の値に強く依存するが、このベースライン成分の値が提案手法と正解データで大きく異なるためである。具体的には、提案手法ではベースライン成分の値u_b を有声区間におけるlog F₀の最小値に設定しているが、正解データでは常にlog 60 Hzに固定しており、提案手法でu_b の値を固定すると推定性能が落ちることが確認されたためである。

提案手法を用いたパラメータ推定結果と、比較手法として選んだ、上記の非特許文献２に記載のパラメータ推定アルゴリズム（非統計的手法）を用いた推定結果とを図６にまとめた。この結果を見れば分かる通り、提案手法の検出率は比較手法と同程度である一方で、提案手法のlog F₀ RMSE の値は比較手法を大きく下回っている。つまり、提案手法を用いた実音声からのフレーズ、アクセント指令関数の推定は、既存手法に匹敵する言語学的な妥当性を満たしつつ、観測F₀ パターンの再現性では既存手法を上回る性能を持っていることが確認できた。

図７に、提案手法のパラメータ推定結果を示した。この図の上のグラフは、有声区間の観測F₀パターン（実線）と推定パラメータから再構築したF₀パターン（点線）であり、下のグラフは推定フレーズ指令関数と推定アクセント指令関数を示したものである。入力F₀パターンは、例として、ATR日本語音声データベースのBセットのNo.353から得られたものを用いた。この例で示したように、本発明は、観測F₀パターンと推定パラメータから再構築したF₀パターンが非常によく一致するようなパラメータ推定が可能である。

以上説明したように、本発明の実施の形態に係る基本周波数モデルパラメータ推定装置によれば、観測基本周波数系列＾ｙが与えられたときの藤崎モデルのパラメータ＾ｏ及びθの対数事後確率ｌｏｇＰ（＾ｏ，θ｜＾ｙ）の下限関数Ｑ（＾ｏ，θ、＾ｏ’，θ’）を目的関数として、目的関数を増加させるように、各々非負値である指令関数＾ｏ、及びパラメータ群θを更新することにより、フレーズ指令とアクセント指令の非負性に関する制約を用いて、藤崎モデルのパラメータを推定することができる。

本実施の形態では、EMアルゴリズムにおけるMステップが補助関数法による反復計算（λの更新ステップと＾oとθの更新ステップの反復計算）により構成され、この反復計算によりQ(＾o,θ,＾o',θ')が減少しないことが保証されているため、目的関数値（＝p(＾o,θ|y)）の収束性が保証される。

また、音声のF₀パターンを入力として藤崎モデルのパラメータを推定する、本実施の形態に係る推定アルゴリズムでは、フレーズ成分とアクセント成分の非負性を直接導入することが可能で、なおかつ収束性が保証されている。具体的には、藤崎モデルの確率モデル表現を畳み込み混合隠れマルコフモデルに基づいて定式化することによって、フレーズ成分とアクセント成分の非負性を直接導入することと共に、収束性を保証することができる。その結果、観測F₀パターンの再現性が非常に高い藤崎モデルのパラメータを推定することができる。

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

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

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

１ 記憶部
２ 基本周波数系列抽出部
３ 有声無声区間推定部
４ 初期値設定部
５ 指令状態系列事後確率更新部
６ モデルパラメータ更新部
７ 収束判定部
８ 状態系列算出部
６１ 補助変数更新部
６２ 指令関数更新部
６３ 収束判定部
６４ 平均振幅更新部
１００ 基本周波数モデルパラメータ推定装置

Claims (5)

1. 音声信号を入力として、隠れマルコフモデルの各時刻ｋの状態ｓ_kからなる指令状態系列＾ｓと、各時刻ｋにおける甲状軟骨の平行移動運動によって生じる基本周波数パターンを表すフレーズ指令ｕ_p［ｋ］及び甲状軟骨の回転運動によって生じる基本周波数パターンを表すアクセント指令ｕ_a［ｋ］のペア＾ｏ[ｋ]からなる指令関数＾ｏと、各時刻ｋにおける状態ｓ_kに応じたフレーズ指令の振幅Ａ_p[ｋ]及び各アクセント指令ｎの振幅Ａ_a ⁽ⁿ⁾を表すパラメータ群θとを推定する基本周波数モデルパラメータ推定装置であって、
   前記音声信号の時系列データから、前記音声信号の各時刻ｋの基本周波数を表す観測基本周波数系列＾ｙを抽出する基本周波数抽出手段と、
   前記音声信号の時系列データについて、有声区間及び無声区間の何れであるかに応じて、各時刻ｋにおける前記基本周波数の不確かさの度合いを推定する有声無声区間推定手段と、
   前記指令関数＾ｏの初期値、及び前記パラメータ群θの初期値を設定する初期値設定手段と、
   前回更新された前記指令関数＾ｏ’または前記指令関数＾ｏの初期値＾ｏ’に基づいて、前記観測基本周波数系列＾ｙ、前記指令関数＾ｏ’、及び前記パラメータ群θ’が与えられたときの指令状態系列＾ｓの事後確率Ｐ（＾ｓ｜＾ｙ，＾ｏ’、θ’）を計算する指令状態系列事後確率更新手段と、
   前回更新された前記指令関数＾ｏ’または前記指令関数＾ｏの初期値＾ｏ’、前記観測基本周波数系列＾ｙ、各時刻ｋにおける前記不確かさの度合い、及び前記事後確率Ｐ（＾ｓ｜＾ｙ，＾ｏ’、θ’）に基づいて、前記観測基本周波数系列＾ｙが与えられたときの前記指令関数＾ｏ及び前記パラメータ群θの対数事後確率ｌｏｇＰ（＾ｏ，θ｜＾ｙ）を目的関数として、前記目的関数を増加させるように、各々非負値である前記指令関数＾ｏ、及び前記パラメータ群θを更新するモデルパラメータ更新手段と、
   予め定められた収束条件を満たすまで、前記指令状態系列事後確率更新手段による計算、及び前記モデルパラメータ更新手段による更新を繰り返し行う第１収束判定手段と、
   を含む基本周波数モデルパラメータ推定装置。
2. 前記モデルパラメータ更新手段は、
   前回更新された各時刻ｌのフレーズ指令ｕ_p［ｌ］又は各時刻ｌのフレーズ指令ｕ_p［ｌ］の初期値に基づいて、時刻ｋ、ｌの組み合わせ（ｋ、ｌ）の各々について、補助変数λ_p,k,lを算出して更新し、かつ、前回更新された各時刻ｋのアクセント指令ｕ_a［ｋ］又は各時刻ｋのアクセント指令ｕ_a［ｋ］の初期値に基づいて、時刻ｋ、ｌの組み合わせ（ｋ、ｌ）の各々について、補助変数λ_a,k,lを算出して更新する補助変数更新手段と、
   前記観測基本周波数系列＾ｙと、各時刻ｋの前記不確かさの程度と、前記計算された指令状態系列の事後確率P(＾s|＾y,＾o′,θ′)と、前記補助変数更新手段によって更新された前記補助変数λ_p,k,l、λ_a,k,lとに基づいて、前記目的関数の下限関数Ｑ（＾ｏ，θ、＾ｏ’，θ’）のさらに下限となる関数を補助関数として、前記補助関数を増加させるように、各時刻ｌのフレーズ指令ｕ_p［ｌ］及びアクセント指令ｕ_a［ｌ］を更新する指令関数更新手段と、
   予め定められた収束条件を満たすまで、前記補助変数更新手段による更新、及び前記指令関数更新手段による更新を繰り返し行う第２収束判定手段と、
   前記指令関数更新手段によって更新された各時刻ｌのフレーズ指令ｕ_p［ｌ］に基づいて、各時刻ｋのフレーズ指令の振幅A_p[ｋ]を更新し、かつ、前記指令関数更新手段よって更新された各時刻ｌのアクセント指令ｕ_a［ｌ］と、前記計算された指令状態系列の事後確率P(＾s|＾y,＾o′,θ′)とに基づいて、各アクセント指令ｎの振幅Ａ_a ⁽ⁿ⁾を更新することにより、前記パラメータ群θを更新する平均振幅更新手段と、
   を含む請求項１記載の基本周波数モデルパラメータ推定装置。
3. モデルパラメータ更新手段によって最終的に更新された指令関数＾ｏに基づいて、前記状態系列＾ｓを算出する状態系列算出手段を更に含む請求項１又は２記載の基本周波数モデルパラメータ推定装置。
4. 音声信号を入力として、隠れマルコフモデルの各時刻ｋの状態ｓ_kからなる指令状態系列＾ｓと、各時刻ｋにおける甲状軟骨の平行移動運動によって生じる基本周波数パターンを表すフレーズ指令ｕ_p［ｋ］及び甲状軟骨の回転運動によって生じる基本周波数パターンを表すアクセント指令ｕ_a［ｋ］のペア＾ｏ[ｋ]からなる指令関数＾ｏと、各時刻ｋにおける状態ｓ_kに応じたフレーズ指令の振幅Ａ_p[ｋ]及び各アクセント指令ｎの振幅Ａ_a ⁽ⁿ⁾を表すパラメータ群θとを推定する基本周波数モデルパラメータ推定方法であって、
   基本周波数抽出手段によって、前記音声信号の時系列データから、前記音声信号の各時刻ｋの基本周波数を表す観測基本周波数系列＾ｙを抽出し、
   有声無声区間推定手段によって、前記音声信号の時系列データについて、有声区間及び無声区間の何れであるかに応じて、各時刻ｋにおける前記基本周波数の不確かさの度合いを推定し、
   初期値設定手段によって、前記指令関数＾ｏの初期値、及び前記パラメータ群θの初期値を設定し、
   指令状態系列事後確率更新手段によって、前回更新された前記指令関数＾ｏ’または前記指令関数＾ｏの初期値＾ｏ’に基づいて、前記観測基本周波数系列＾ｙ、前記指令関数＾ｏ’、及び前記パラメータ群θ’が与えられたときの指令状態系列＾ｓの事後確率Ｐ（＾ｓ｜＾ｙ，＾ｏ’、θ’）を計算し、
   モデルパラメータ更新手段によって、前回更新された前記指令関数＾ｏ’または前記指令関数＾ｏの初期値＾ｏ’、前記観測基本周波数系列＾ｙ、各時刻ｋにおける前記不確かさの度合い、及び前記事後確率Ｐ（＾ｓ｜＾ｙ，＾ｏ’、θ’）に基づいて、前記観測基本周波数系列＾ｙが与えられたときの前記指令関数＾ｏ及び前記パラメータ群θの対数事後確率ｌｏｇＰ（＾ｏ，θ｜＾ｙ）を目的関数として、前記目的関数を増加させるように、各々非負値である前記指令関数＾ｏ、及び前記パラメータ群θを更新し、
   第１収束判定手段によって、予め定められた収束条件を満たすまで、前記指令状態系列事後確率更新手段による計算、及び前記モデルパラメータ更新手段による更新を繰り返し行う
   基本周波数モデルパラメータ推定方法。
5. 請求項１〜請求項３の何れか１項に記載の基本周波数モデルパラメータ推定装置の各手段としてコンピュータを機能させるためのプログラム。