音声 F0 パターン生成過程の確率モデル∗○亀岡弘和,ルルージョナトン,大石康智

日本電信電話 (株) NTT コミュニケーション科学基礎研究所

モデル化される。甲状軟骨の平行移動運動に関係する F0
1 序論
パターンをフレーズ成分，回転運動に関係する F0 パター
音声の韻律的特徴を表現する物理量の一つである基本 ンをアクセント成分と呼び，それぞれ yp (t), ya (t) とす
周波数 (F0 ) の時間変化パターン (以後，F0 パターン) に る。ただし，t は時刻である。yp (t) の生成過程 (フレーは，構文や意図の伝達に関連する情報が豊富に含まれて ズ制御機構) はフレーズ指令と呼ぶパルス波を入力としおり，言語情報に関連する量であるフォルマントの時間 た臨界制動の二次線形系，ya (t) の生成過程 (アクセント変化パターンと並んで音声の情報伝達における重要な役 制御機構) はアクセント指令と呼ぶ矩形波を入力とした割を担っている。フォルマントは声道の形状の変化に付 臨界制動の二次線形系により表現される。以上の二つの随して変化する声道共振特性を特徴づける量であるのに 成分と，声帯の物理的性質によって決まるベースライン
対し，F0 パターンは肺からの空気供給量と甲状軟骨の運 成分 yb の和 yp (t) + ya (t) + yb で F0 パターン y(t) を表動に付随して変化する声帯の伸びに応じて決まる量であ したものが藤崎モデルである。フレーズ成分は短区間のる。それゆえ，音声からは音声生成メカニズムの物理的 上昇のあとに緩やかな下降をなす成分で，アクセント成制約の範囲内の F0 パターンしか生成され得ない。この 分は急激な上昇下降をなす成分であるため，多くの言語ように考えると，ある F0 パターンがどれだけ音声のもの に共通して前者が句単位の比較的緩やかな音調を，後者であるらしいかを客観的に評価するためには，それが音 が語または音節単位の比較的急激で局所的な音調を表現声生成メカニズムの物理的制約をどれだけ満たしている する役割を担っていると考えられている。フレーズ制御かを規準に評価すれば良い，ということになる。従って， 機構およびアクセント制御機構は{ 2 −αtF0 パターンの生成過程のモデル化は，音声の自然性を考
α te(t ≥ 0)
(1)
Gp (t) =慮に入れることで性能向上が望めるいかなる音声アプリ
0(t < 0){ケーションにおいて大変有益となる可能性がある。
1 − (1 + βt)e−βt (t ≥ 0)
(2)
Sa (t) =F0 パターン生成過程のモデルの代表格といえば藤崎の0(t < 0)
モデル [1] (以後，藤崎モデル) であろう。藤崎モデルは，発話の言語学的情報と密接な関係にあるパラメータを与 で与えられるインパルス応答 Gp (t) とステップ応答 Sa (t)えることによって，実測の音声 F0 パターンを極めて良く によって特徴づけられる。α と β はそれぞれの制御機構近似するパターンを生成することが可能であり，その有 の固有角周波数であるが，これらは話者ごとにほぼ一効性は音声合成に広く利用されていることからも実証さ 定値をとること，話者の個人差も言語による差も比較的れている。また，実測の F0 パターンから藤崎モデルのパ 小さいことが確かめられており，おおよそ α = 3rad/s,ラメータを推定する逆問題を扱った研究例も多く [2]，藤 β = 20rad/s 程度であることが経験的に知られている。
崎モデルのパラメータを音声認識や感情認識などの精度 これらを用いて F0 パターン y(t) は具体的に向上に役立てる試みも取り組まれている。しかし，逆問∑
題の解析的な複雑さゆえにいまだ有効なパラメータ推定
y(t) = yb +Ap,i Gp (t − T0,i )
i法は確立されていないのが現状である。
∑{}
F0 パターンのモデル化の重要性と藤崎モデルの (潜在Aa,j Sa (t − T1,j ) − Sa (t − T2,j ) (3)+j的な) 有効性を踏まえ，本稿では，藤崎モデルを確率過程に基づいて統計モデル化する。その目的は，(1) 統計的と表される。ただし，Ap,i と T0,i はそれぞれ i 番目のフ手法を駆使した強力なパラメータ推定法の枠組を確立す
レーズ指令の大きさと生起時刻であり，Aa,j と T1,j , とること，(2) 統計モデルに基づく多くの音声処理問題 (分T2,j はそれぞれ j 番目のアクセント指令の振幅と始端時
析，合成，分離，強調) に，F0 パターンの統計モデルを刻と終端時刻を表す。ところで，アクセント制御機構の新たな音声らしさの規準として組み込めるような下地をインパルス応答 Ga (t) は Sa (t) の時間微分ゆえ作ること，にある。

2

藤崎の F0 パターン生成過程モデル [1]藤崎モデルでは，甲状軟骨の二つの独立な運動 (平行移動と回転) に伴う声帯の長さの変化の合計が F0 の時間的変化をもたらすと解釈され，声帯の伸びと対数 F0 の変化が比例関係にあるという仮定に基づき F0 パターンが∗{ 2 −βtβ teGa (t) =0(t ≥ 0)
(t < 0)(4)となり，Gp (t) と同形であることが分かる。従って，アクセント成分 ya (t) は複数の矩形波からなるアクセント指令関数 ua (t) と Ga (t) の畳み込みによって表される。

A probabilistic model of speech F0 contours. by KAMEOKA Hirokazu, LE ROUX Jonathan and OHISHI
Yasunori (NTT Communication Science Laboratories)

3

藤崎モデルの離散時間表現

本章では，藤崎モデルの離散時間表現を得るため，連続時間システムのフレーズ制御機構およびアクセント制御機構に対して後退差分変換により離散化を行う。
まず，フレーズ制御機構の伝達関数 (Laplace 変換) は[]
Gp (s) = L Gp (t) = α2 /(s + α)2 で与えられる。後退差分変換は，時間微分演算子 s を z 領域における後退差分演算子 s ≃ (1 − z −1 )/t0 に置き換える変換であり (t0 は変換先の離散時間表現におけるサンプリング周期)，この変換により，Gp−1 (s) の逆システムの伝達関数を z 領域でHp−1 (z) = a2 z −2 + a1 z −1 + a0(5)と書くことができる。ただし，

a2 = (ψ − 1) , a1 = −2ψ(ψ − 1), a0 = ψ22(6)および，ψ = 1 + 1/(αt0 ) である。ここで，k を離散時刻インデックスとし，フレーズ指令関数およびフレーズ成分の離散時間表現をそれぞれ up [k] および yp [k] とすると，yp [k] は，単一のパラメータ ψ(あるいは α) によって特性が決まる拘束つき全極モデルからの出力

up [k] = a0 yp [k] + a1 yp [k − 1] + a2 yp [k − 2](7)と見なすことができる。同様に，アクセント指令関数ua [k] とアクセント成分 ya [k] の関係もua [k] = b0 ya [k] + b1 ya [k − 1] + b2 ya [k − 2](8)と 書 く こ と が で き る 。た だ し b2 = (φ − 1)2 , b1 =−2φ(φ − 1), b0 = φ2 , φ = 1 + 1/(βt0 ) で あ る 。
ベースライン成分 yb (t) の離散時間表現を yb [k] とすると，藤崎モデルの離散時間表現はこれらの三成分の和
y[k] = yp [k] + ya [k] + yb [k] で与えられる。

4

藤崎モデルの統計モデル化

まず，up [k] と ua [k] のモデル化について述べる。藤崎モデルにおいて，フレーズ指令とアクセント指令に関して以下のような規則が定められている。
(A1) フレーズ指令はパルス波，アクセント指令は矩形波である。(A2) 発話の開始または先行フレーズ内のアクセント指令終了後にフレーズ指令が発生する。また，フレーズ指令の後にアクセント開始指令が発生する。つまり，アクセント指令発生中はフレーズ指令は発生しない。
(A3) アクセント指令の開始した後には必ずアクセント終了指令が発生する。つまり，隣接するアクセント指令同士は重なり合うことはない。
藤崎モデルのパラメータ推定における難しさの一つは，これらの制約の下で最適推定をいかにして行えるかという点にあろう。特に，上記の (A2), (A3) より，up [k] とua [k] は相互に依存し合う関係がある点に注意が必要である。提案法では，これを解決するため隠れマルコフモデ
ル (HMM) を用いて指令入力信号を確率モデル化する。

p1a0

p0a1a2

Profile ofa3

p0 p1 a0

図1a1a0 a2

p0 p1 a0state

Command function modeling with HMM.

まず，o[k] := (up [k], ua [k])T を，o[k] ∼ N (ν[k], Υ)[][ 2σµ [k]ν[k] := µp[k] , Υ := p0a0σa2](9)(10)のように正規分布する確率変数と見なし，平均 ν[k] が
図 1 のような状態遷移に伴って変化するモデルを考える。
これは HMM に他ならず，このように o[k] を HMM でモデル化したことにより，状態遷移の経路制限 (状態遷移確率の設定) を通して ν[k] に対して上記の (A1)〜(A3)を満たすような制約を与えることが可能となる。提案する HMM の構成は以下のとおりである。
出力値系列: {o[k]}K
k=1
状態集合: S := {p0 , p1 , a0 , · · · , aN }状態系列: {sk }K
k=1
状態出力分布: P (o[k]|sk = i) = N (ci [k], Υ)( )T(i = p0 , a0 )(0, 0)Tci [k] = Ap [k], 0 (i = p1 )( 0, A(n) )T (i = a )an
状態遷移確率: ϕi′ ,i := log P (sk = i|sk−1 = i′ )簡単のため状態遷移確率を定数とすると，以上より指令入力モデルにおいて推定すべきパラメータは，フレーズ指令の大きさ Ap [k]，状態遷移系列 sk ，アクセント指(n)22
令の大きさ {Aa }N
n=1 ，指令入力信号の分散 σp , σa であり，これらをまとめて θu と記す。また，平均値系列KK{µp [k]}K
k=1 および {µa [k]}k=1 は，状態遷移系列 {sk }k=1が与えられたもとで (µp [k], µa [k])T ← csk [k] で与えられる。
前述の指令入力モデルに基づき y = (y[1], · · · , y[K])Tの確率密度関数を導く。式 (9), (10) より，up := (up [1],· · · , up [K])T, ua := (ua [1], · · · , ua [K])T, µp := (µp [1],· · · , µp [K])T, µa := (µa [1], · · · , µa [K])T とすると，
up |θu ∼ N (µp , Σp ), Σp = σp2 I(11)ua |θu ∼ N (µa , Σa ),(12)Σa = σa2 Iが言える。3 章で得た関係式より，フレーズ成分 yp :=(yp [1], · · · , yp [K])T と up の関係，および，アクセント

成分 ya := (ya [1], · · · , ya [K])T と ua の関係は，





a0Ob0
O
aabb
a1 a0 a
 1 0

2 1 0
, B := b2 b1 a0
 (13)A := 
 .. .. ..
 .. .. ..

. . .
. . .
Oa2 a1 a0
Ob2 b1 b0としよう。y と yobs との関係は，各行に 1 が一個あり残りはすべて 0 であるような K ′ × K のバイナリ行列 Mを用いて yobs = M y と表される。これを利用して EMの各ステップを導くことができる。紙面の都合上導出は省略せざるを得なかったが，E ステップでは y をy ← µ + ΣM T (M ΣM T )−1 (yobs − M µ)と置くと，それぞれ up = Ayp , ua = Bya のように表現できることから，式 (11), (12) より(−1−1−1 T)

yp |θu , α ∼ N A µp , A Σp (A )()ya |θu , β ∼ N B −1 µa , B −1 Σa (B −1 )T(16)が言える。ただし，Σb = σb2 I であり，θb := {µb , σb2 }と置く。仮定より，yp , ya , yb は独立なので，Θ := {θu ,α, β, θb } が与えられたもとでの F0 パターン y = yp +ya + yb の確率密度関数は，式 (14), (15) と式 (16) より，A−1 Σp (A−1 )T + B −1 Σa (B −1 )T + Σb)(17)で与えられる。以上より，{}
1|Σ |T −1exp − (y − µ) Σ (y − µ)P (y|Θ) =2(2π)T /2−1 1/2µ = A−1 µp + B −1 µa + µb 1(18)()()−1−1Σ = A−1 Σp AT+ B −1 Σa B T+ Σbが，F0 パターン y が与えられたときの藤崎モデルパラ
メータ Θ の尤度関数である。
Θ の 事 前 確 率 に つ い て は ，各 要 素 は 独 立 で ，状 態
遷移系列 {s[k]}Kk=1 と制御パラメータの ψ と φ 以外の パ ラ メ ー タ は 一 様 に 分 布 す る と 仮 定 し ，P (Θ) ∝∏KP (ψ)P (φ)P (s1 ) k=2 P (sk |sk−1 ) とする。

5Θ ← argmax log P (y|Θ)P (Θ)

提案モデルの応用場面

以上の F0 パターンの統計モデルの具体的な応用例について述べる。まず，提案モデルを用いて実測 F0 パターンからパラメータを推定する問題を考える。式 (18) では
全区間で F0 データが観測されていることが暗に想定されているが，実際には F0 は有声音が発せられるときのみ観測可能であり，無声音の区間では F0 データは通常観測できない。従って，実測 F0 パターンからのパラメータ推定を行うためには一般に欠損データの問題を扱う必要がある。統計的枠組では，この類の問題は不完全データ問題に他ならず，Expectation-Maximization(EM) アルゴリズムにより扱うことができる。y ∈ RK を観測 F0 データと欠損 F0 データからなる完全データとし，実際に観測さ′れた F0 データを並べたベクトルを yobs ∈ RK (K ′ ≤ K)

(20)

Θ(15)

yb |µb ∼ N (µb 1, Σb )(y|Θ ∼ N A−1 µp + B −1 µa + µb 1,により更新し，M ステップではその y を用いて(14)が導かれる。ベース成分 yb [k] についても，同様に白色
Gauss 性雑音 ϵb [k] に起因する確率変数 yb [k] = µb +ϵb [k] と仮定し，ϵb [k] ∼ N (0, σb2 ) とし，同様に，ϵξ [j] とϵξ′ [j ′ ] は (ξ, j) ̸= (ξ ′ , j ′ ) のとき独立とすると，(19)を実行すれば良い。なお自明ではあるが，欠損区間がない場合を表す M = I においては式 (19) は y ← yobs となり，yobs をそのまま完全データと見なして良いという意味にたしかになっていることが分かる。
第二の応用方法としては，F0 をパラメータにもつ何らかの音声モデルの事前分布として利用するやり方である。
具体的な応用方法については [3] で詳しく論じているのでそちらを参照されたい。

6

パラメータ推定アルゴリズム

P (Θ|y) を最大化する問題 (式 (20) に相当) は解析的に解くことはできないが，x := (ypT , yaT , ybT )T を完全データと見なすことで EM アルゴリズムによる不完全データ問題に帰着できる。この場合，完全データの対数尤度は，先に見たとおり，11log |Λ−1 | − (x − m)T Λ−1 (x − m)22
[ ][ −1 ]ypA µp
x := ya , m := B −1 µa(21)ybµb 1 T −1OOA Σp A
Λ−1 := OB T Σ−1O  (22)a BOO
Σ−1bclog P (x|Θ) =で与えられる。このとき，Q 関数 Q(Θ, Θ′ ) は，

Q(Θ, Θ′ ) =c1[log |Λ−1 | − tr(Λ−1 E[xxT |y; Θ′ ])2]+ 2mT Λ−1 E[x|y; Θ′ ] − mT Λ−1 m + log P (Θ) (23)となる。ここで，y = Hx(ただし，H = [I I I]) であるから，E[x|y; Θ] と E[xxT |y; Θ] は，具体的にE[x|y; Θ] = m + ΛH T (HΛH T )−1 (y−Hm) (24)
E[xxT |y; Θ] = Λ − ΛH T (HΛH T )−1 HΛ+ E[x|y; Θ]E[x|y; Θ]T(25)と書ける。E ステップでは，直前のステップで更新されたモデルパラメータを Θ′ に代入し，上記に基づいて
E[x|y; Θ′ ] と E[xxT |y; Θ′ ] が算出される。yp , ya , yb に対応するように E[x|y; Θ′ ] および E[xxT |y; Θ′ ] を]][Rp ∗ ∗x̄pT′
E[x|y; Θ ] = x̄a , E[xx |y; Θ ] = ∗ Ra ∗ (26)∗ ∗ Rbx̄b[′+ log P (Θ)(27)と書き直せて，これを用いて各パラメータについて M ス
テップの更新式を求めることができる。

1) 状態系列: Q 関数の中で s := {sk }K
k=1 に関係する項は2(o[k] − csk [k])T Υ−1 (o[k] − csk [k])k=1+ log P (s1 ) +K
∑log P (sk |sk−1 )(28)

k=2F0 (Hz)
Magnitude
Magnitude

6KOGa U
0.5

06KOGa U

0.50012345Time (s)6789Test data: F0 contour in solid line and phrase component in
dotted line (top), phrase (middle) and accent commands (bottom).

Magnitude

I1 (s) := −K1∑210

図2F0 (Hz)T −1T −1− tr(B T Σ−1a BRa ) + 2—a Σa B x̄a − —a Σa —aiT −1
T −1− tr(Σ−1b Rb ) + 2—b Σb x̄b − —b Σb —bMagnitudeのように区分表現すると，Q 関数は
hc 1
T −1−1log |AT Σ−1Q(Θ, Θ′ ) =p A| + log |B Σa B| + log |Σb |2T −1
T −1− tr(AT Σ−1p ARp ) + 2—p Σp Ax̄p − —p Σp —p

106KOGa U
0.506KOGa U

0.5001

図32345Time (s)6789Parameter estimates: µ in solid line and A−1 µp + µb indotted line (top), µp (middle), and µa (bottom).となる。ただし，o[k] := ([Ax̄p ]k , [B x̄a ]k )T であり，[·]k はベクトルの k 番目の要素を表す。これを最大化 7 実験する状態遷移系列 {sk }K
提案モデルを 6 章で提案した最適化法のアルゴリズム
k=1 は動的計画法により効率的に解くことができる。まず，すべての状態 i について としての純粋なふるまいを調べる目的で，元の藤崎モデル
δ1 (i) を δ1 (i) = − 12 (o[1] − ci [1])T Υ−1 (o[1] − ci [1]) と置 を用いて合成した人工の F0 パターン (図 2) を対象とし
くと，k = 2, · · · , K について逐次的に δk (i) を δk (i) = た動作実験を行った。F0 パターンの長さは 10s，サンプ[]maxi′ δk−1 (i′ )− 12 (o[k]−ci [k])T Υ−1 (o[k]−ci [k])+ϕi′,i
リング周期は 10ms，各パラメータは図 2 の中段と下段のにより計算していくことができる。各ステップで選択さ とおりとした。アルゴリズムの反復回数は 10 とし，N はれる状態番号 Ψk (i) = argmaxi′ [δk−1 (i′ ) + ϕi′,i ] を記憶し 10，状態遷移確率はそれぞれ ϕp ,p =log(0.999), ϕp ,p =0 00 1ておくことで，k = K まで到達後に sk−1 = Ψk (sk ) (k = log(0.001), ϕp ,a =log(1.0), ϕa ,a =log(0.999), ϕa ,a =1 00 0n 0K, · · · , 2) により選択された状態番号を辿っていくこと log(0.001), ϕa0 ,an = log(0.0001), ϕan ,an = log(0.899),ができ，最適経路 s1 , · · · , sK を得ることができる。
ϕa ,p = log(0.1), 1 ≤ n ≤ 10 とした。以上の条件で得n02) フ レ ー ズ 制 御 パ ラ メ ー タ: ψ の 事 前 分 布 を ψ ∼ た µ, µp ，µa の推定結果を図 3 に示す。図 3 と図 2 をN (µψ , 1/νψ2 ) とする。式 (6) より，A は定数行列 U2 , 比較すると，真の値に極めて近い値が得られていること
U1 , U0 を用いて A = U2 ψ 2 + U1 ψ + U0 と表せるので， が分かる。このことから式 (20) に相当する最適化問題は提案アルゴリズムにより効果的に解けることが分かった。
Q 関数の ψ に関する偏導関数は 4 次式となり，また，離散時間表現への変換による近似の影響はさほど
2tr(U2T U2 Rp )ψ 4 + 3tr(U2T U1 Rp )ψ 3
大きくないことも確認できた。
2 22+ {tr((2U2T U0 + U1T U1 )Rp ) − 2µT
p U2 x̄p + σp νψ }ψ
2 22+ {tr(U1T U0 Rp ) − µT
p U1 x̄p − 2σp νψ µψ }ψ − 2Kσp8まとめ

本稿では，藤崎モデルの離散時間表現への変換とそのの根を解くことで極値を求めることができる。求まった統計モデル化を通して F0 パターン生成過程の統計モデ
4 つの極値の中で I2 (ψ) を最大にする ψ が更新値となる。
ルを構築した。提案モデルの応用例について簡単に触れ，
3) アクセント制御パラメータ: 同様に φ の事前分布を パラメータ推定を EM アルゴリズムによって実現できるφ ∼ N (µφ , 1/νφ2 ) とする。更新値の導出過程は 4) と同 ことを示した。人工の F0 パターンに対しパラメータ推様なので省略する。
定アルゴリズムを実行し，その基本動作を確認した。歌
声の F0 パターンのモデル化についても [4] で検討してい
4) その他:るので興味のある読者は是非参照されたい。

Ap [k] = ûp [k], (k ∈ Tp1 )
(29)
∑1
A(n)
[B x̄a ]k , Tan = {k|sk = an } (30)a =|Tan |k∈TanTµb = 1 x̄b /T( ())
T
σp2 = tr AT ARp − 2µTp Ax̄p + µp µp /K( ())Tσa2 = tr B T BRa − 2µTa B x̄a + µa µa /K()
σb2 = tr(Rb ) − 2µb 1T x̄b /K + µ2b(31)
(32)
(33)
(34)

参考文献
[1] H. Fujisaki, In Vocal Physiology: Voice Production, Mechanisms and Functions, (O. Fujimura,ed.) Raven Press, pp. 347–355, 1988.
[2] S. Narusawa et al., In Proc. ICASSP’02, Vol. 1,pp. 509–512, 2002.
[3] 亀岡, 音講論 (秋)’10, 1-1-4, 2010.
[4] 大石ら, 音講論 (秋)’10, 3-P-31, 2010.