1-1-4
全極型声道モデルと F0 パターン生成過程モデルを内部にもつ統一的音声生成モデル∗○亀岡弘和
日本電信電話 (株) NTT コミュニケーション科学基礎研究所

1

はじめにのように 0 次近似および 1 次近似する。µ は θ の導関数
(時間微分) であり，s の瞬時基本周波数を表す。式 (4)を式 (1) に代入すると，時刻 t 周辺の s を近似した信号

人は，音声対話において，声帯の振動周期と声道の共振特性を巧みに時間変化させながらそれぞれに抑揚や意N∑図に関係する情報と音韻に関係する情報を時間に載せてst (u) =wn (t)ej(nµ(t)(u−t)+cn (t))(5)
人に伝達している。音声情報処理研究の歴史において，
n=1
声道および声帯振動に関する重要な数理モデルとして，全極型モデルによる声道スペクトルモデル [1] と藤崎の を得る。ただし，cn (t) := nθ(t) + φn である。s と stF0 パターン生成過程モデル [2] が双璧をなしている。本 は大域的にはまったく異なる信号であるが，時刻 t 周辺
稿では，これらのモデルを同時に内包する初めての音 においてのみ似通っている点がポイントである。通常
声モデルを提案し，これを用いて観測信号から声道スペ ウェーブレット基底関数 ψα,t は時刻 t 周辺でのみ優勢
クトルモデルのパラメータと F0 パターンモデルのパラ であるため，時刻 t におけるウェーブレット変換スペク
トルは，式 (5) を式 (3) に代入して得られる st のウェーメータを一挙に推定できることを示す。
ブレット変換2

音声パワースペクトルの統計モデル

Sw (log α1 , t) ≃ ⟨st , ψα,t ⟩

本稿では，音声信号モデルとして，n 次調波成分の瞬時振幅が wn (u)，瞬時位相が nθ(u) + φn の擬似周期信号の解析信号表現s(u) =N∑j(nθ(u)+φn )wn (u)e, wn (u) ≥ 0(1)n=1を考え (u は時刻，N は調波成分の最大次数)，この信号のウェーブレット変換スペクトルを導く。まず，周波数が1 の，アドミッシブル条件を満たす適当なマザーウェーブレットを ψ(u) とし，ウェーブレット基底関数を()
1u−tψα,t (u) := √ψ
, α>0(2)α2παと定義する。α はスケールパラメータ，t はシフトパラメータであり，ψα,t (u) は時刻 t 周辺に局在する周期 αの成分を測るためのウェーブレットである。これを用いて，s の連続ウェーブレット変換を∫ ∞∗Sw (log α1 , t) = ⟨s, ψα,t ⟩ =s(u)ψα,t(u)du (3)−∞と定義する。ただし，⟨·, ·⟩ は内積を表す。
ここで，時刻 t におけるウェーブレット変換スペクトルについて着目する。音声信号 s において，θ と wn の時間変化はいずれも十分緩やかであると仮定し，まずこれらを時刻 t 周辺でそれぞれwn (u) ≃ wn (t), θ(u) ≃ θ(t) + µ(t)(u − t)
∗N
∑

=(6)

wn (t)ejcn (t) Ψ∗ (αnµ(t))(7)

n=1により近似できる。ただし，Ψ はマザーウェーブレット
ψ の Fourier 変換である。α は周期を表すので対数周波数 x = log α1 の関数として Sw を表すと，

Sw (x, t) =N∑

wn (t)ejcn (t) Ψ∗ (ne−x+Ω(t) )(8)n=1と書ける。ただし，Ω(t) := log µ(t) であり，Ω(t) は対数瞬時基本周波数を表す。
ここで，周波数特性 Ψ(ω) が，次のような ω = 1 で最大値をとる単峰的な実関数
e− 4σ12 (log ω)2 (ω > 0)
(9)Ψ(ω) :=
0(ω ≤ 0)となるようなマザーウェーブレットを選ぶと，Sw はSw (x, t) =N∑

wn (t)ejcn (t) e− 4σ2 (x−Ω(t)−log n)12(10)

n=1と具体的に書ける。よってそのパワースペクトルは|Sw (x, t)|2 =λ(x, t) + o(x, t)λ(x, t) :=(4)N∑

wn (t)2 e− 2σ2 (x−Ω(t)−log n)
1(11)
2(12)

n=1

Uniﬁed speech production model incorporating all-pole vocal-tract model and F0 contour generating processmodel. by KAMEOKA, Hirokazu (NTT Communication Science Laboratories)

日本音響学会講演論文集- 211 -2010年9月

と書ける。o(x, t) は調波間干渉を表す交差項であり，調波成分同士の重なりが少ないほど小さくなる。
周期 u0 でサンプリングされた実際の観測音声信号に対して上記と同条件で算出したウェーブレット変換のパワースペクトルを Y [k, l] := Y (xk , tl )(ただし k は対数
周波数インデックス，l は離散時刻インデックス) とすると，wn (t) と Ω(t) が真の値と等しい場合には以上より

Y [k, l] ≃ λ[k, l]のような基本周波数 µ[l] := µ(tl ) の周期信号とする。また，上記声帯音源モデルにおける複素振幅 vl,n について

vl,n ∼ NC (0, 1), n = 1, · · · , Nを仮定する。ただし NC は複素正規分布を表す。式 (15) の両辺に離散時間 Fourier 変換 (DTFT) を行い，Al (z) :=al [0] − al [1]z −1 · · · − al [M ]z −M と置いて整理すると，√(13)
Sl (ω) =となることが期待される。ただし，λ[k, l] := λ(xk , tl )である。wn (t) と Ω(t) が真値と等しくともなお生じる
Y [k, l] と λ[k, l] との間の誤差にはさまざまな要因があり，例えば，式 (1) の形で定義した擬似周期性の仮定からの逸脱，式 (4) の近似，交差項 o の存在，などが挙げられる。提案法の枠組では，これらの複合的な誤差要因を詳細にモデル化することはせず，まとめて一挙に確率的な現象と見なし，

Y [k, l] ∼ Pois(λ[k, l])(14)と具体的に仮定する。ただし，Pois(c) は c をパラメータとする Poisson 分布を表す。この仮定の下での λ の最尤推定問題は，スペクトル間の近さを測る尺度の一つとして近年音響信号処理分野で多用される I ダイバージェンスと呼ぶ歪み尺度を規準とした Y と λ の最適フィッティング問題と等価となる。提案法の枠組におけるこの仮定の有難みはのちに明らかにする。

3

音声生成モデルに基づく事前分布設計

前述のパワースペクトルモデル λ[k, l] において推定すべき未知変数は，瞬時振幅 wn [l] := wn (tl ) および対数基本周波数 Ω[l] := Ω(tl ) である。提案法のポイントは，これら両未知変数の事前分布を，それぞれ音声生成プロセスのモデルとして古典的で有名な全極型声道モデル [1] と F0 パターン生成過程モデル [2] に基づいて設計する点にある。以下，その設計方法を順に述べる。

3.1

全極型声道モデル

時刻 tl における式 (5) の音声信号モデルの，観測信号と同じサンプリング周期 u0 の下での離散時間表現を
sl [i] := stl (iu0 ) とし，これを全極型システムからの出力

sl [i] =M∑

al [m]sl [i − m] + ϵl [i](15)m=0と考える。al [0], · · · , al [M ] は全極モデルのパラメータであり，線形予測モデルにおいて予測係数と呼ぶものに該当する。ここで，声帯音源に対応する入力 ϵl [i] を

ϵl [i] =N∑n=1

日本音響学会講演論文集vl,n ejnµ[l]iu0(16)(17)

N
2π ∑vl,n δ(ω − nµ[l]u0 )Al (ejω ) n=1(18)となる。ただし，ω は規格化角周波数である。上式に対し逆 DTFT を行うと sl [i] の陽な表現

1sl [i] = √
2π=

∫ π√N
2π ∑vl,n δ(ω − nµ[l]u0 )ejiω dωjω
−π Al (e ) n=1

N∑vl,nejnµ[l]iu0jnµ[l]u0 )A(en=1 l(19)を得る。式 (5) を離散時間表現したもの

sl [i] =N∑

wn [l]ej(nµ[l](iu0 −tl )+cn (tl ))(20)n=1と式 (19) を比較すると，¯¯¯¯vl,n¯ , n = 1, · · · , N
wn [l] = ¯¯Al (ejnµ[l]u0 ) ¯(21)なる関係式を導ける。µ[l]u0 は基本周波数 µ[l] を規格化
角周波数の単位で表わしたものなので，1/|Al (ejµ[l]u0 )|,· · · , 1/|Al (ejN µ[l]u0 )| は全極スペクトルを基本周波数の間隔でサンプリングした値に他ならない。以上より wn [l]が従う分布を具体的に導くことができる。式 (17) より()vl,njnµ[l]u0 2∼
N0,1/|A(e)|ClAl (ejnµ[l]u0 )(22)が言え，実部と虚部が独立で，かつ等しい分散の正規分
布に従う複素数の絶対値は Rayleigh 分布に従うことが知られており，式 (21), (22) より最終的に()
wn [l] ∼ Rayleigh 1/|Al (ejnµ[l]u0 )|(23)が得られる。この事前分布において，Θw = {al [m] | 1 ≤l ≤ L, 0 ≤ m ≤ M } がハイパーパラメータとなる。

3.2

F0 パターン生成過程モデル

次に，対数基本周波数の時間軌跡 (以後，F0 パターン)Ω := (Ω[1], · · · , Ω[L])T の事前分布について述べる。音声 F0 パターンの代表的なモデルとして藤崎モデル [2] が有名である。藤崎モデルは喉頭の生理的，物理的特性に基づいて声帯振動制御機構をモデル化したものであり，その離散時間表現の確率過程のモデルが最近筆者らによって定式化された [4]。藤崎モデルの詳細およびその確率モデルの導出過程については [4] に委ねることとす- 212 -

2010年9月

るが，藤崎モデルパラメータ ΘΩ が与えられた下での Ωの条件つき分布が，多次元正規分布Ω ∼ N (µΩ , ΣΩ )(24)の形となることさえ踏まえておけば以後の議論においては十分である。なお，µΩ , ΣΩ は藤崎モデルパラメータΘΩ の値に応じて決まる変数であるため，ΘΩ がこの事前分布のハイパーパラメータとなる。

4パラメータ推定アルゴリズムと定数項を除いて等しくなる。L を最大化する Ω, w,
ΘΩ , Θw の解を解析的に解くことはできないが，補助関数法に基づく反復計算により局所最適解を探索すること∑ができる。L の補助関数は， n γn [k, l] = 1 を満たす
補助変数 γn [k, l] を用いて表される不等式∑21logwn [l]2 e− 2σ2 (xk −Ω[l]−log n) ≥(31)n{}wn [l]212γn [k, l] log−(xk − Ω[l] − log n)γn [k, l] 2σ 2n∑に基づいて設計される下限関数
以上のモデルにおいて，Ω, w := {wn [l] | 1 ≤ n ≤N, 1 ≤ l ≤ L}, ΘΩ , Θw が未知のパラメータであり，観測パワースペクトル Y := {Y [k, l] | 1 ≤ k ≤ K, 1 ≤ l ≤ L}を最も良く説明するような Ω, w, ΘΩ , Θw の値を推定したい。そこで，モデルパラメータの事後確率

L̃(γ) :=n,l1− (Ω − µΩ )T Σ−1Ω (Ω − µΩ )2)(∑wn [l] wn [l]2+ 2β−log
Pn [l]2Pn [l]を Ω, w, ΘΩ , Θw に関して最大化する最大事後確率(MAP) 推定問題を考える。2 章および 3 章より各確率は∏ λ[k, l]Y [k,l] e−λ[k,l](26)

Y [k, l]!

で与えられる。ただし，Ω̂ := (Ω̂[1], · · · , Ω̂[L])T は(27)∏ wn [l] − wn [l]2p(w|Θw ) =e 2Pn [l]Pn [l]Ω̂l =で与えられる。ただし，λ[k, l] は式 (12) のとおりであり，Pn [l] := 1/|Al (ejnµ[l]u0 )|2 である。ここで，以下を仮定する。

(A1) p(Θw ) を一様分布とする。
(A2) 対数周波数のサンプリング間隔 x0 := xk+1 − xkが十分小さく，以下の近似が成り立つ。
∫ ∞∑ 2 222ex /2σ dx ≃ x0exk /2σ(29)Y [k, l] logの対角行列である。
以上の準備により，補助関数法の各ステップにおける更新式が導ける。

1) 補助変数 γn [k, l]:wn [l]2 e−(xk −Ω[l]−log n) /2σ
γn [k, l] ← ∑22wn [l]2 e−(xk −Ω[l]−log n) /2σ−2πσ ∑x0n,l(35)2) F0 パターン Ω:Ω ← (Σ + ΣΩ )−1 (ΣΩ̂ + ΣΩ µΩ )2(36)3) n 次調波成分の瞬時パワー wn [l]:∑Y [k, l]γn [k, l] + β

1wn [l] − log |ΣΩ |221− (Ω − µΩ )T Σ−1Ω (Ω − µΩ )2()∑wn [l] wn [l]2+ 2βlog−Pn [l]2Pn [l]2nwn [l]2 e− 2σ2 (xk −Ω[l]−log n)1(34)k,n2nk,l√∑(33)
Y [k, l]γn [k, l]を要素とするベクトル，Σ は対角要素 Σl が/∑Y [k, l]γn [k, l]Σl = σ 2k以上より，対数事後確率 log p(Ω, w, ΘΩ , Θw |Y ) は∑∑k,n(A3) 事前分布の信頼度を調節できるようにする目的で，
式 (28) 右辺を 2β 乗したものを p(w|Θw ) とする。

L :=Y [k, l]γn [k, l](xk − log n)k,n(28)n,l(32)

n,l

∑

1/2
T −1|Σ−11Ω |p(Ω|ΘΩ ) =e− 2 (Ω−µΩ ) ΣΩ (Ω−µΩ )(2π)L/2−∞

wn [l]2γn [k, l]

1− (Ω − Ω̂)T Σ−1 (Ω − Ω̂)
2√
2πσ ∑
1−
wn [l]2 − log |ΣΩ |x02= p(Y |Ω, w)p(Ω|ΘΩ )p(w|Θw )p(ΘΩ )p(Θw ) (25)k,lY [k, l]γn [k, l] logn,k,lp(Ω, w, ΘΩ , Θw |Y ) ∝ p(Y |Ω, w)p(Ω, w, ΘΩ , Θw )p(Y |Ω, w) =∑wn [l]2 ←k

1 + β/Pn [l](37)

(30)

n,l

日本音響学会講演論文集- 213 -2010年9月

klim

β→∞

∑

1 + β/Pn [l]Y [k, l]γn [k, l] + βk

lim

β→0= Pn [l]1 + β/Pn [l]

(38)

8000

3200

Frequency (Hz)

この更新式は β → ∞ および β → 0 の極限をとると∑Y [k, l]γn [k, l] + β

1600
800
400200=∑100

Y [k, l]γn [k, l] (39)8000k
Frequency (Hz)3200となる。式 (38) は前段で推定された全極スペクトルを基本周波数の間隔でサンプリングした値，式 (39) は観測スペクトルに最も良くフィットする値をそれぞれ表しており，式 (37) は両者の中間的な値に相当する。

16008004002001005) 全極モデルパラメータ Θw : L̃ において Θw に関係する項は [3] で設定されている目的関数と同形となるので，Θw の更新方法には当該文献と同様のアイディアが適用できる。詳細は [3] に委ねるが，al :=(al [0], al [1], · · · , al [M ])T とすると，
al ← R−1 h(40)h ← R̂al(41)を繰り返すことで L̃ を増加させる al を得ることができる。ただし，Rl および R̂l はそれぞれ

1 ∑wn [l]2 cos nµ[l]u0 (m′ − m) (42)N n=1NRl (m′ − m) =1 ∑R̂l (m − m) =Pn [l] cos nµ[l]u0 (m′ − m)
N n=1
N′を m′ 行 m 列成分とする Toeplitz 行列である。

5動作実験

(43)

106KOGa U

2.7F0 (Hz)4) 藤崎モデルパラメータ ΘΩ : [4] の反復アルゴリズム
により L̃ を増加させる ΘΩ を得ることができる。

102.1



Time (s)図1

提案法の動作実験例。ウェーブレット変換による観測スペクトログラム Y [k, l] (上段) と λ[k, l] の推定値 (中段) とΩ の推定値 (下段・赤線) および —Ω の推定値 (下段・青線)であった。また，反復計算の初期段階において藤崎モデ
ルのパラメータ推定が不安定な動作をすることが多かった。これは，藤崎モデルのパラメータ更新ステップにおいては，有音区間の Ω[l] の推定値と本来は無視するべ
きはずの無音区間における Ω[l] の推定値も同等の価値を置いて推定を行っていることに原因があると考えられる。以上のことから，Ω と w を周辺化して Y から直接ΘΩ と Θw を推定できる方法をいずれ検討していく必要がありそうだという感触を得た。

6まとめ

本稿では，全極型声道スペクトルモデルと藤崎の F0パターン生成過程モデルを事前分布の形として内包する音声ウェーブレットスペクトルの統計モデルを提案し，パラメータの効率的な最適化方法を補助関数法に基づいて導出した。この枠組により，全極型声道モデルと藤崎モデルのパラメータをスペクトルからダイレクトに同時推定することが可能となるが，動作実験の結果，解決すべきいくつかの課題が明らかになった。

ATR 音声データベースの女性話者音声データ (サンプリング周波数 16kHz のディジタル信号) を用いて提案法の基本動作の確認を行った。図 1 にその動作実験例を示す。図 1 上段に上記の音声信号に対して x0 = 12cent,σ = 0.02, t0 = tk+1 − tk = 8ms の条件下で算出したウェーブレット変換スペクトログラム，中段に 4 章で述べた最適化アルゴリズムにより得た λ[k, l] の推定値を示 参考文献
す。また下段に F0 パラメータ Ω(赤線) および推定した藤崎モデルパラメータにより表現される F0 パターン (青 [1] F. Itakura, S. Saito, In Proc. 6th Int l Cong.
Acoust. (ICA 68), C-5-5, C17–20, 1968.
線) を示す。現段階では提案法の具体的な性能を定量的[2] H. Fujisaki, In Vocal Physiology: Voice Producに評価したわけではないが，動作実験の過程で既にいくtion, Mechanisms and Functions, (O. Fujimura,つかの課題が明らかになった。具体的には，Ω がまだ明ed.) Raven Press, pp. 347–355, 1988.
らかに正しい F0 の解に至っていない段階にもかかわら [3] A. El-Jaroudi, J. Makhoul, IEEE trans. Signal Proず収束してしまうケースがしばしば見られた。これは，cess., Vol. 39, No. 2, pp. 411–423, 1991.
Ω と µΩ は互いに近くなるようにいずれも更新されるた [4] 亀岡ら, 音講論 (秋), 1-1-3, 2010.
め，Ω の初期値に強い拘束を受けてしまうのが原因であると考えられる。この問題は w と Θw についても同様

日本音響学会講演論文集- 214 -

2010年9月