社団法人 電子情報通信学会
THE INSTITUTE OF ELECTRONICS,
INFORMATION AND COMMUNICATION ENGINEERS

信学技報
TECHNICAL REPORT OF IEICE.

全極型声道モデルと F0 パターン生成過程モデルを内部にもつ統計的音声スペクトルモデル
亀岡弘和†† 日本電信電話株式会社 NTT コミュニケーション科学基礎研究所〒 243-0198 神奈川県厚木市森の里若宮 3-1
E-mail: †kameoka@cs.brl.ntt.co.jp
あらまし

音声情報処理研究の歴史において，音声生成過程に関する重要な数理モデルとして，声帯振動の制御機構を模した藤崎モデルと，縦続連結した等長音響管で声道断面積関数をモデル化したことと等価な全極型声道モデルが有名である。本発表では，ある特定条件のウェーブレット変換のスペクトル領域において，これら両生成過程モデル
を同時に内包する確率モデルが立てられることを示し，そのパラメータの推論法について述べる。提案法は，全極型
声道モデルと藤崎モデルのパラメータをスペクトルから同時に直接推定することが可能な初めての音声分析手法であり，様々な興味深い応用展開が期待される。
キーワード 全極型声道モデル，F0 パターン生成過程モデル (藤崎モデル)，ウェーブレット変換，EM アルゴリズム

Statistical speech spectrum model incorporating all-pole vocal tract
model and F0 contour generating process model
Hirokazu KAMEOKA†† NTT Communication Science Laboratories
3-1 Morinosato Wakamiya, Atsugi, Kanagawa, 243-0198 Japan
E-mail: †kameoka@cs.brl.ntt.co.jp
Abstract In this paper, we propose to introduce a statistical speech spectrum model, simultaneously incorporating
the all-pole vocal tract model and the F0 contour generating process model.
Key words All-pole vocal tract model, F0 contour generating process model, wavelet transform, EM algorithm

1. ま え が き
音声は，肺から供給された空気が声帯振動によって周期的な疎密波となり，さらにその音波が声道で共振されることによって「声」となり，口から発せられる。人は音声対話において，声帯の固有振動数を巧みに制御しながら音声に適切な抑揚を与え，同時に，声道の形状を巧みに変化させながら異なる種類の母音を生成することで，聴き手に意図した情報を伝達している。
音声情報処理研究の歴史において，音声生成過程に関する重要な数理モデルとして，声帯振動の制御機構を模した藤崎モデルと，縦続連結した等長音響管で声道断面積関数をモデル化したことと等価な全極型声道モデルが有名である。本稿では，全極型声道モデルと藤崎モデルを同時に内包する統計的音声スペク
トルモデルを提案する。
声道スペクトルと基本周波数 (F0 ) は，分析，合成，認識，強調，符号化などの幅広い音声アプリケーションにおいて根幹となる音声の重要な特徴量であり，音声信号からこれらを自動的に推定する技術は長年にわたって研究され，発展を遂げてきた。
声道スペクトルを推定する問題と F0 を推定する問題は，一見すると別の問題のようで実は密接に関係している。例えば，以下，音声信号が声道のインパルス応答と声帯による周期的なパルス波との畳み込みで表現できる場合を考えると，声道スペク
トルは調波成分のピークをつないだ包絡線に対応するので，F0がもし既知であればどの成分が調波成分であるかが特定でき，
声道スペクトルを精度良く推定することができる。一方で，F0
推定では，実際の F0 の整数分の一倍と推定し誤ってしまうことがしばしば問題となるが，この半ピッチ誤りの問題は，声道ス
ペクトルが滑らかであるという先験的な知識を使えば効果的に回避することができる (半ピッチ誤りの状況では滑らかでない
スペクトル包絡を仮定してしまっていることに相当するため)。
また，音声の F0 の時間軌跡 (F0 パターンと呼ぶ。) も通常は滑らかであり，これを制約とすれば同様の推定誤りを回避できることが多い。以上のように，声道スペクトルと F0 パターン
は各々の推定において互いに有用な情報であるため，別個に求めるより，同時最適化問題として推定できるようにした方が多分に合理的である。以上の動機に基づき，スペクトル包絡構造と F0 パターンを同時モデル化したことが提案モデルにおける第一のポイントである。

—1—

提案モデルの第二の重要なポイントは，音声生成過程の確率モデル化を通して，与えられた任意の信号がどれだけ音声のものらしいかを定量的に評価することが可能になる点である。音
声生成過程から生成しえない信号には低いスコアを，音声生成過程から生成しうる信号には高いスコアを与えるような規準がもしあれば，音声分離，雑音除去，残響除去などの音声を対象としたあらゆる逆問題において音声を抽出する手掛かりとして大いに活用できることが期待されるからである。
以上より，スペクトル包絡構造と F0 パターンをいかにして同時モデル化するかという点と，いかに音声の生成過程を正確に統計モデルの形で表現できるかという点に主眼を置き，以下
で音声スペクトルのモデル化を行う。

wn (u)ej(nθ(u)+ρn ) , wn (u) >
=0

(1)

n=1

を考え (u は時刻，N は調波成分の最大次数)，この信号のウェー
ブレット変換スペクトルを導く。周波数が 1 の，アドミッシブル
条件を満たす適当なマザーウェーブレットを ψ(u) とし，ウェー
ブレット基底関数を„«
1u−t
ψx,t (u) := √
ψ
(2)
e−x
2πe−x
と定義する。e−x と t はスケールとシフトパラメータに該当し，
ψx,t (u) は時刻 t 周辺に局在する周期 e−x の成分を測るためのウェーブレットである。x は対数周波数に対応する。これを用いて，s の連続ウェーブレット変換を
Z ∞∗Sw (x, t) = ⟨s, ψx,t ⟩ =s(u)ψx,t
(u)du
(3)
−∞

と定義する。ただし，⟨·, ·⟩ は内積を表す。
ここで，時刻 t におけるウェーブレット変換スペクトルについて着目する。音声信号 s において，θ と wn の時間変化はい
ずれも十分緩やかであると仮定し，これらを時刻 t 周辺で

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

st (u) =NX

wn (t)ej(nθ̇(t)(u−t)+ξn (t))

(5)

n=1

を得る。ただし，ξn (t) := nθ(t) + ρn である。s と st は大域的にはまったく異なる信号であるが，時刻 t 周辺においてのみ
似通っている点がポイントである。通常ウェーブレット基底関
数 ψx,t は時刻 t 周辺でのみ優勢であるため，時刻 t におけるウェーブレット変換スペクトルは，式 (5) を式 (3) に代入して得られる st のウェーブレット変換

Sw (x, t) ≃ ⟨st , ψx,t ⟩
=

N
X
n=1

wn (t)ejξn (t) Ψ∗ (ne−x+Ω(t) )( − 1 (log ω)2
e 4σ20(ω > 0)
(ω <
= 0)

(8)

となるようなマザーウェーブレットを選ぶと，Sw はN
X− 12 (x−Ω(t)−log n)24σ

wn (t)ejξn (t) e

(9)

n=1

本章では音声信号モデルとして，n 次調波成分の瞬時振幅がwn (u)，瞬時位相が nθ(u) + ρn の擬似周期信号の解析信号表現
NX

Ψ(ω) :=Sw (x, t) =2. 音声パワースペクトルモデル

s(u) =

により近似できる。ただし，Ψ はマザーウェーブレット ψ のフーリエ変換である。また，Ω(t) := log θ̇(t) であり，Ω(t) は対数瞬時基本周波数を表す。以上より，ウェーブレット変換のスペクトルはマザーウェーブレット ψ(t) を具体的にどう定めるかによって決まる。ここで，周波数特性 Ψ(ω) が，次のような ω = 1 で最大値をとる単峰的な実関数

(6)
(7)

と具体的に書ける。このマザーウェーブレットによって構成される定 Q フィルタバンクを定 Q 対数正規分布型フィルタバン
クと呼ぶことにする。ここで，|Sw (x, t)|2 は|Sw (x, t)|2 =λ(x, t) + o(x, t)

(10)

N
X

(11)

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

Y [k, l] ≃ λ[k, l](12)となることが期待される。ただし，λ[k, l] := λ(xk , tl ) である。
wn (t) と Ω(t) が真値と等しくともなお生じる Y [k, l] と λ[k, l]との間の誤差にはさまざまな要因があり，例えば，式 (1) の形で定義した擬似周期性の仮定からの逸脱，式 (4) の近似，交差
項 o の存在，などが挙げられる。提案法の枠組では，これらの複合的な誤差要因を詳細にモデル化することはせず，まとめて一挙に確率的な現象と見なし，

Y [k, l] ∼ Pois(λ[k, l])(13)と具体的に仮定する。ただし，Pois はポアソン分布を表す。この仮定の下での λ の最尤推定問題は，スペクトル間の近さを測る尺度の一つとして近年音響信号処理分野で多用される I ダイ
バージェンスと呼ぶ歪み尺度を規準とした Y と λ の最適フィッ
ティング問題と等価となる。
以上のパワースペクトルモデル λ[k, l] において推定すべき
未知変数は，瞬時振幅 wn [l] := wn (tl ) および対数基本周波数Ω[l] := Ω(tl ) である。提案法のポイントは，これら両未知変数の事前分布を，それぞれ音声生成過程のモデルとして古典的で有名な全極型声道モデル [1] と F0 パターン生成過程モデル [2]に基づいて設計し，スペクトルの生成過程を階層モデルにより表現する点にある。以下，各事前分布の設計方法を順に述べる。

3. 全極モデルに基づく瞬時振幅の確率分布
本章では，全極型声道モデルに基づき wn [l] の事前分布を導く。時刻 tl における式 (5) の音声信号モデルの，観測信号と同じ
サンプリング周期 u0 の下での離散時間表現を sl [i] := stl (iu0 )—2—とし，これを全極型システムからの出力

sl [i] =MX

4. 藤崎モデルに基づく F0 パターンの確率分布

al [m]sl [i − m] + ϵl [i]

(14)

m=0

と考えよう。al [0], · · · , al [M ] は全極モデルのパラメータであり，線形予測モデルにおいて予測係数と呼ぶものに該当する。
ここで，sl [i] の基本周波数は eΩ[l] であるため，声帯音源に対応する入力 ϵl [i] も基本周波数が eΩ[l] の周期信号

ϵl [i] =

N
X

Ω[l]

vl,n ejne

iu0

(15)

n=1

と仮定しておく必要がある (線形フィルタには入力の周波数成分の周波数値を変える効力はないため)。また，上記声帯音源
モデルにおける複素振幅 vl,n について

vl,n ∼ NC (0, 1), n = 1, · · · , N(16)を仮定する。ただし，NC は複素正規分布を表す。よって，この仮定は，駆動音源に関してすべての調波成分が等しいパワー
スペクトル密度をもち，位相を一様にとりうることを仮定していることに相当する。式 (14) の両辺に離散時間フーリエ変換(DTFT) を行い，Al (z) := al [0] − al [1]z −1 · · · − al [M ]z −M と置いて整理すると，√
N
2π X
vl,n δ(ω − neΩ[l] u0 )(17)
Sl (ω) =jω
Al (e ) n=1となる。ただし，Sl は sl [1], · · · , sl [T ] の DTFT，ω は規格化角周波数である。上式に対し逆 DTFT を行うと sl [i] の陽な表現

1
sl [i] = √
2π
=

N
X
n=1

Z π

√

N
2π X
vl,n δ(ω − neΩ[l] u0 )ejiω dωjω−π Al (e ) n=1vl,nAl (ejneΩ[l] u0 )

Ω[l]

ejne

iu0

(18)

を得る。式 (5) を離散時間表現したもの

sl [i] =

N
X

Ω[l]

wn [l]ej(ne(iu0 −tl )+cn (tl ))

(19)

n=1

と式 (18) を比較すると，
˛
˛
˛
˛
vl,n
˛ , n = 1, · · · , N
wn [l] = ˛˛Ω[l] u
˛
jne
0)
Al (e

(20)

なる関係式を導ける。eΩ[l] u0 は基本周波数 eΩ[l] を規格化角Ω[l]周波数の単位で表わしたものなので，1/|Al (eje u0 )|, · · · ,Ω[l]1/|Al (ejN e u0 )| は全極スペクトルを基本周波数の間隔でサンプリングした値に他ならない。以上より，全極モデルパラメー
タおよび Ω[l] が与えられた下での wn [l] が従う分布を具体的に導くことができる。式 (16) より

`
´
Ω[l]
vl,n
∼ NC 0, 1/|Al (ejne u0 )|2
Al (ejneΩ[l] u0 )(21)

が言え，実部と虚部が独立でかつ等しい分散の正規分布に従う
複素数の絶対値はレイリー分布に従うため，式 (20), (21) より

`
´
Ω[l]
wn [l] ∼ Rayleigh 1/|Al (ejne u0 )|(22)が得られる。この事前分布のパラメータ (ハイパーパラメータ)
は，πw = {al [m] | 1 <= M } である。この事前分布=m<= L, 0 <=l<は，Ω[l] に依存する点に注意が必要である。

本章では，藤崎モデルに基づいて Ω := (Ω[1], · · · , Ω[L])T の事前分布を導く。藤崎モデルは喉頭の声帯振動制御機構の物理モデルであるが，これを離散時間表現の確率過程として記述した確率モデルが最近筆者らによって提案された [4]。本章の内容
はこの新モデルがベースになっている。
藤崎モデルでは，甲状軟骨の二つの独立な運動 (平行移動と回転) に伴う声帯の長さの変化の合計が F0 の時間的変化をもたらすと解釈され，声帯の伸びと対数 F0 の変化が比例関係にあるという仮定に基づき F0 パターンがモデル化される。甲状軟骨の平行移動運動に関係する F0 パターンをフレーズ成分，回転運動に関係する F0 パターンをアクセント成分と呼び，それぞれ Ωp (t), Ωa (t) とする。Ωp (t) の生成過程 (フレーズ制御機構) はフレーズ指令と呼ぶパルス波を入力とした臨界制動の二次線形系，Ωa (t) の生成過程 (アクセント制御機構) はアクセン
ト指令と呼ぶ矩形波を入力とした臨界制動の二次線形系により表現される。以上の二つの成分と，声帯の物理的性質によって決まるベースライン成分 Ωb の和 Ωp (t) + Ωa (t) + Ωb で F0 パ
ターン Ω(t) を表したものが藤崎モデルである。フレーズ制御機構およびアクセント制御機構はGp (t) = α2 te−αt (t >
= 0)
2

−βt

Ga (t) = β te

(t >
= 0)

(23)
(24)

で与えられるインパルス応答 Gp (t) と Ga (t) によって特徴づけられる。α と β はそれぞれの制御機構の固有角周波数である。
フレーズ指令とアクセント指令をそれぞれ up (t) と ua (t) のよ
うに時間の関数として表せば，Ω(t) は Gp (t) と Ga (t) を用いてΩ(t) = Ωb + Ωp (t) + Ωa (t)(25)Ωp (t) = (Gp ∗ up )(t)(26)

Ωa (t) = (Ga ∗ ua )(t)(27)

と表される。ただし，∗ は畳み込みを表す。
藤崎モデルの離散時間表現を得るため，連続時間システムのフレーズ制御機構およびアクセント制御機構に対して後退差分変換により離散化を行う。まず，フレーズ制御機構の伝達ˆ˜
関数 (ラプラス変換) は Gp (s) = L Gp (t) = α2 /(s + α)2 である。後退差分変換は，時間微分演算子 s を z 領域における後退差分演算子 s ≃ (1 − z −1 )/t0 に置き換える変換であり (ただし，t0 = tl − tl−1 )，この変換により，Gp−1 (s) の逆システムの伝達関数を z 領域で Gp−1 (z) = gp2 z −2 + gp1 z −1 + gp0 と書ける。ただし，gp2 = (ϕ − 1)2 , gp1 = −2ϕ(ϕ − 1), gp0 = ϕ2 ,ϕ = 1 + 1/(αt0 ) である。ここで，時刻 tl におけるフレーズ指令関数およびフレーズ成分をそれぞれ up [l], Ωp [l] と表すと，Ωp [l] は単一のパラメータ ϕ(または α) によって特性が決まる
拘束つき全極モデルからの出力

up [l] = gp0 Ωp [l] + gp1 Ωp [l − 1] + gp2 Ωp [l − 2](28)と見なすことができる。同様に，アクセント指令関数 ua [l] とアクセント成分 Ωa [l] の関係も

ua [l] = ga0 Ωa [k] + ga1 Ωa [l − 1] + ga2 Ωa [l − 2](29)と書くことができる。ただし，ga2 = (φ−1)2 , ga1 = −2φ(φ−1),
ga0 = φ2 , φ = 1 + 1/(βt0 ) である。ベースライン成分 Ωb (t) の離散時間表現を Ωb [l] とすると，藤崎モデルの離散時間表現は
これらの三成分の和 Ω[l] = Ωp [l] + Ωa [l] + Ωb [l] で与えられる。

—3—

p1

a0

p0

a1

率密度関数を導く。式 (31) より，up := (up [1], · · · , up [L])T,ua := (ua [1], · · · , ua [L])T, —p := (µp [1], · · · , µp [L])T, —a :=(µa [1], · · · , µa [L])T とすると，(32)

Σa = diag(υa2 [1], · · · , υa2 [L])(33)ua ∼ N (—a , Σa ),a2Profile of Gaussian meanup ∼ N (—p , Σp ), Σp = diag(υp2 [1], · · · , υp2 [L])が言える。式 (28), (29) の関係式より，フレーズ成分 Ωp :=(Ωp [1], · · · , Ωp [L])T と up の関係，および，アクセント成分Ωa := (Ωa [1], · · · , Ωa [L])T と ua の関係は，

a32

図 1 Command function modeling with HMM.

次に，up [l] と ua [l] を確率モデル化する。そこで，藤崎モデ
ルにおいてフレーズ指令とアクセント指令に関して以下のような規則が定められていることに注意する。
(A1) フレーズ指令はパルス波，アクセント指令は矩形波である。(A2) 複数の指令が同時に発生することはない。
藤崎モデルのパラメータ推定における難しさの一つは，これらの制約の下で最適推定をいかにして行えるかという点にある。
提案法では，これを解決するため隠れマルコフモデル (HMM)
を用いて up [l] と ua [l] をモデル化する。up [l] と ua [l] を組にしたベクトル o[l] := (up [l], ua [l])T を考え，o[l] をo[l] = ‌[l] + ›[l](30)のようにフレーズ指令からなる時系列 µp [l] とアクセント指令
からなる時系列 µa [l] を組にしたベクトル ‌[l] = (µp [l], µa [l])Tに 雑 音 ›[l] が 加 わった も の と 見 な す。›[l] を 共 分 散 行 列 が
Υ[l] = diag(υp2 [l], υa2 [l]) のガウス性雑音とすると，以上より

o[l] ∼ N (‌[l], Υ[l])(31)が言える。ここで，‌[l] および Υ[l] が状態遷移に伴って変化する図 1 のような HMM を考える。HMM による o[l] のモデル化により，状態遷移の経路制限を通して両指令の組 ‌[l] に対して (A1), (A2) を満たすような制約を巧く与えることが可能となる。提案する HMM の構成は以下のとおりである。
出力系列: {o[l]}L
l=1
状態集合: S := {p0 , p1 , a0 , · · · , aN }状態系列: {sl }L
l=1
出力分布: P (o[l]|sl = j) = N (cj [l], Dj )cp0 [l] = ca0 [l] = (0, 0)T(n)cp1 [l] = (Ap [l], 0)T , can [l] = (0, Aa )T
Dp0 = Da0 = diag(σp20 , σa20 )Dp1 = diag(σp21 , σa20 ), Dan = diag(σp20 , σa2n )遷移確率: ϕj ′ ,j := log P (sl = j|sl−1 = j ′ )簡単のため状態遷移確率および分散 σp20 , σp21 σa20 , · · · , σa2Nを定数とすると，以上より指令入力モデルにおいて推定すべき未知パラメータは，状態遷移系列 sl ，フレーズ指令の(n)
大きさ Ap [l]，アクセント指令の大きさ {Aa }N
n=1 であり，
これらをまとめて πu と記す。また，平均値系列 {µp [l]}L
l=1 ,
L
L
{µa [l]}Lおよび分散値系列
{υ
[l]}
,
{υ
[l]}
は，状態
p
a
l=1
l=1
l=1
T
遷移系列 {sl }L
l=1 がひとたび決まれば (µp [l], µa [l]) ← csl [l],diag(υp [l], υa [l]) ← Dsl (l = 1, · · · , L) により一挙に決まる。
前述の指令入力モデルに基づき Ω = (Ω[1], · · · , Ω[L])T の確

2
3
3
gp0
O
ga0
O
6gp1 gp0
6ga1 ga0
7
7
6
6
7
7
gp2 gp1 gp0
7, Ga := 6ga2 ga1 ga0
7
Gp := 66
6
7
7.. .. ..
.. .. ..
4
4
5
5
. . .
. . .
O

gp2 gp1 gp0

O

ga2 ga1 ga0と置くと，それぞれ up = Gp Ωp , ua = Ga Ωa のように表現できることから，式 (32), (33) より`−1−1 T ´Ωp ∼ N G−1
p —p , Gp Σp (Gp )`−1−1 T ´Ωa ∼ N G−1
a —a , Ga Σa (Ga )

(34)
(35)

が 導 か れ る 。ベ ー ス 成 分 Ωb [l] に つ い て も ，同 様 に Ωb ∼N (µb 1, Σb ) とし，Ωp , Ωa とは独立と仮定する。ただし，Σbは定数とする。仮定より，Ωp , Ωa , Ωb は独立なので，藤崎モ
デルパラメータ πΩ := {πu , α, β, µb } が与えられたもとでのF0 パターン Ω = Ωp + Ωa + Ωb の確率密度関数は，`−1Ω ∼ N G−1
p —p + Ga —a + µb 1,
−1 T
−1
−1 T
G−1
p Σp (Gp ) + Ga Σa (Ga ) + Σb

´

(36)

で与えられる。以上より，以下を得る。

P (Ω|πΩ ) =ȷﬀ|Σ−1 |1/2
1
T −1
exp
−
(Ω
−
—
)
Σ
(Ω
−
—
)
Ω
Ω
Ω
2
(2π)T /2

−1—Ω = G−1(37)p —p + Ga —a + µb 1
`
´
`
´
−1
−1
T
T
+ Σb
+ G−1
ΣΩ = G−1p Σp Gpa Σa Ga

5. パラメータ推定アルゴリズム
以上のモデルをグラフィカルモデルで表わすと図 2 のようになる。Ω, w := {wn [l] | 1 <=n<= N, 1 <=l<= L}, πΩ , πw が未知のパラメータであり，観測スペクトル Y := {Y [k, l]}1<
=k <
=K,1<
=l<
=L
を最も良く説明するような Ω, w, πΩ , πw の値を推定することが目的である。そこで，モデルパラメータの事後確率p(Ω, w, πΩ , πw |Y ) ∝ p(Y |Ω, w)p(Ω, w, πΩ , πw )= p(Y |Ω, w)p(Ω|πΩ )p(w|Ω, πw )p(πΩ )(38)を Ω, w, πΩ , πw に関して最大化する最大事後確率 (MAP) 推定問題を考える。2.〜4. より各確率はp(Y |Ω, w) =Y λ[k, l]Y [k,l] e−λ[k,l]k,lY [k, l]!

1/2
T −1
1
|Σ−1
Ω |
e− 2 (Ω−—Ω ) ΣΩ (Ω−—Ω )L/2(2π)
Y wn [l] − wn [l]2p(w|Ω, πw ) =e Pn [l]Pn [l]p(Ω|πΩ ) =

(39)

(40)
(41)

n,l

—4—

で 与 え ら れ る 。た だ し ，λ[k, l] は 式 (11) の と お り で あ り，Ω[l]Pn [l] := 1/|Al (ejne u0 )|2 である。以上の MAP 推定問題の大域最適解は解析的に解くことはできないが，以下に示すように Expectation-Maximization(EM) アルゴリズムにより局所解を得ることができる。
EM アルゴリズムの導出では，潜在データと呼ぶ，観測不可な架空のデータをどのように定義するかがポイントとなる。ここでは，観測スペクトル Y [k, l] のうち n 番目の高調波に対応しているであろう成分 Cn [k, l] を潜在データと見なす。従って，XY [k, l] =Cn [k, l]
(42)

Ṡ

ṡ

ṟ
Ṟ
ṟ

Ṟ

Ṡ

n

を満たす。今，我々は Y [k, l] ∼ Pois(λ[k, l]) を仮定しているた
め，これと矛盾しないように Cn [k, l] が従う分布を仮定する必要がある。スペクトルモデル λ[k, l] が n 個の調波成分
1 2

λn [k, l] = wn [l]2 e− 2σ (xk −Ω[l]−log n)

2を最大化する Ω, w, πΩ , πw の更新則を導出できれば，以下のとおり EM アルゴリズムが実現できる。

(44)(E ステップ) 更新された Ω, w を Ω′ , w′ に代入する。
(M ステップ) I(Ω, w, πΩ , πw ) を最大化 (または増加) するよ
うに Ω, w, πΩ , πw を更新し，E ステップに戻る。

を仮定すれば，ポアソン分布の再生性より，矛盾することなく Y [k, l] ∼ Pois(λ[k, l]) が成り立つことが分かる。以上
より，観測データ Y := {Y [k, l]}1<と潜在データ
=k <
=K,1<
=l<
=L
C := {Cn [k, l]}1<を組にした完全データ
= n<
=N,1<
=k <
=K,1<
=l<
=L
{Y, C} の確率密度関数がp(Y, C|Ω, w) = p(Y |C)p(C|Ω, w)”
Y “X
p(Y |C) =δ Y [k, l] −Cn [k, l]p(C|Ω, w) =

(45)
(46)

n

k,l

Y Y λn [k, l]Cn [k,l] e−λn [k,l]n(47)Cn [k, l]!

k,lで与えられる。完全データの対数確率の，観測データとモデル
パラメータが与えられた下での条件つき期待値を Q 関数と呼び，
Q 関数をもとに EM アルゴリズムの各ステップの更新則が導かれる。以下で本問題における Q 関数を導く。log p(Y, C|Ω, w)の中で Ω, w に関係する項に対し，Cn [k, l] の期待値をとったものが Q 関数となるので，Q 関数はQ(Ω, w) =XX`⟨Cn [k, l]⟩p(Cn [k,l]|Y [k,l],Ω′ ,w′ )n´が言えるため，期待値は具体的に

(50)

X`
´
λ′ [k, l]log λn [k, l]−λn [k, l]Y [k, l] n′
λ [k, l]n,k,l+ log p(Ω|πΩ ) + log p(w|Ω, πw )k√P
よって， k λn [k, l] ≃ 2πσwn [l]2 /x0 である。
(B2) 事前分布の信頼度を調節できるようにする目的で，式(41) 右辺を 2γ 乗したものを p(w|Ω, Θw ) とする。
(B3) p(w|Ω, πw ) において Ω を 1 ステップ前の更新値 Ω′ に置き換える。すなわち，p(w|Ω, πw ) ≃ p(w|Ω′ , πw ) と近似する。
以上の準備により，M ステップの更新式が導ける。

1) F0 パターン Ω: I(Ω, w, πΩ , πw ) で Ω に関係する項は
1c
I(Ω, w, πΩ , πw ) = − (Ω − Ω̂)T Σ−1 (Ω − Ω̂)21− (Ω − —Ω )T Σ−1Ω (Ω − —Ω )2(53)

Ω = (Σ + ΣΩ )−1 (ΣΩ̂ + ΣΩ —Ω )(54)

で与えられる。ただし，Ω̂ := (Ω̂[1], · · · , Ω̂[L])T は
ﬃX
X Y [k, l] X ′λn [k, l](xk − log n)
Y [k, l] (55)Ω̂l =λ′ [k, l] nkkを要素とするベクトル，Σ は対角要素 Σl が.X
Σl = σ 2
Y [k, l](56)kと求まる。ただし，Binom は二項分布を表し，λ′n [k, l], λ′ [k, l]は Ω′ , w′ を代入した λn [k, l], λ[k, l] を表す。以上より Q 関数の具体形を得た。あとは，Q 関数に対数事前分布を加えたI(Ω, w, πΩ , πw ) :=

−∞

(48)

で与えられる。ただし ⟨f (x)⟩p(x|y) は条件つき期待値 Ex [f (x)|y]を表す。ここで，ポアソン分布の性質より，“λ′ [k, l] ”
Cn [k, l]|Y [k, l] ∼ Binom Y [k, l], n′(49)
λ [k, l]λ′ [k, l]⟨Cn [k, l]⟩p(Cn [k,l]|Y [k,l],Ω′ ,w′ ) = Y [k, l] n′
λ [k, l]M ステップの更新則を導出するにあたり，以下を仮定する。
(B1) 対数周波数のサンプリング間隔 x0 := xk+1 − xk が十分小さく，以下の近似が成り立つ。
Z ∞X −x2 /2σ222e−x /2σ dx ≃ x0
e k
(52)

であるので，I を最大化する Ω はk,llog λn [k, l] − λn [k, l]

図 2 Graphical model representation of the proposed model.

(43)の和で表されていることを思い出せば，Cn [k, l] ∼ Pois(λn [k, l]), (n = 1, · · · , N )

ṡ

(51)

の対角行列である。

2) n 次調波成分の瞬時パワー wn [l]: ∂I/∂wn [l] = 0 を解くと，I を最大化する wn [l] はXkwn [l]2 = √Y [k, l]λ′n [k, l]+γ
λ′ [k, l]2πσ/x0 + γ/Pn [l]

(57)

—5—

で与えられる。この更新式は γ → ∞, γ → 0 の極限をとると2πσ/x0 + γ/Pn [l]X
λ′ [k, l]Y [k, l] n′+γ
λ [k, l]λ′ [k, l]
x0 X
k
lim √
Y [k, l] n′= √γ→0λ [k, l]2πσ k2πσ/x0 + γ/Pn [l]となる。上段の式は直前に更新された全極スペクトルを基本周
波数の間隔でサンプリングした値，下段の式はスペクトルモデ
ルが観測スペクトルに最も良くフィットする wn [l]2 の値をそれぞれ表しており，式 (57) は両者の間を補間する値に相当する。
γ はその補間係数と見ることができる。

3) 藤崎モデルパラメータ πΩ : I において πΩ に関係する項は
1c
I(Ω, w, πΩ , πw ) = − log |ΣΩ |2
1− (Ω − —Ω )T Σ−1Ω (Ω − —Ω )
2

XX“
l

log

n

wn [l]2 ”
wn [l]2−Pn [l]2Pn [l](58)hl ← R̂l al(61)

を繰り返すことで R を増加させる al を得ることができる。ただし，Rl および R̂l はそれぞれN

2 X
Pn [l] cos neΩ[l] u0 (m′ − m)
N n=1

0

1.2

2.4

3.6

4.8

3.6

4.8

Time (s)
8000
3200
1600
800
400
200
100
0

1.2

2.4
Time (s)

400
300
200
1

2

3

4

図 3 観測スペクトログラム Y [k, l] (上段) と λ[k, l] の推定値 (中段)と Ω の推定値 (下段・赤線) および —Ω の推定値 (下段・青線)

0.4
0.2
0

1

2

3

4

Time (s)図 4 フレーズ指令 (青線) とアクセント指令 (緑線) の推定結果(60)

R̂l (m′ − m) =100

0al ← Rl−1 hl

1 X
wn [l]2 cos neΩ[l] u0 (m′ − m)
N n=1

200

(59)

であり，[3] で設定されている目的関数と類似した形となるので，
πw の更新方法には当該文献と同様のアイディアが適用できる。
詳細は省略するが，al := (al [0], al [1], · · · , al [M ])T とすると，

Rl (m′ − m) =400

Time (s)4) 全極モデルパラメータ πw : I において πw に関係する項は，
c

800

0

であり，[4] の目的関数と同形となるので，πΩ の更新方法には当該文献と同様の反復アルゴリズムが適用できる。詳細は [4]を参照されたい。

I(Ω, w, πΩ , πw ) = γ

Frequency (Hz)= Pn [l]Frequency (Hz)

γ→∞

√

1600

F0 (Hz)

lim

k

8000
3200

λ′ [k, l]Y [k, l] n′+γ
λ [k, l]Magnitude

X

(62)

N

(63)

を m′ 行 m 列成分とする Toeplitz 行列である。

6. 動 作 実 験
ATR 音声データベースの女性話者音声データ (サンプリング
周波数 16kHz のディジタル信号) を用いて提案法の基本動作の確認を行った。図 3 にその動作実験例を示す。図 3 上段に上記の音声信号に対して x0 = 12cent, σ = 0.02, t0 = tl − tl−1 = 8msの条件下で算出したウェーブレット変換スペクトログラム，中段に 5. で述べた最適化アルゴリズムにより得た λ[k, l] の推定値を示す。また，下段に F0 パラメータ Ω(赤線) および推定した藤崎モデルパラメータにより表現される F0 パターン (青線)を示す。F0 パラメータの推定結果を見ると，基本周波数推定においてしばしば問題となる倍ピッチ・半ピッチ誤りがほとん
どないことが分かる。これはスペクトル包絡の滑らかさと，F0
パターンの滑らかさの制約が効果的に働いていることを示している。また，図 4 に示すように藤崎モデルパラメータがスペク
トルから直接得られる点も提案法の大きな特徴の一つである。
現段階では提案法の具体的な性能を定量的に評価したわけではないが，動作実験の過程で既にいくつかの課題が明らかになった。具体的には，Ω がまだ明らかに正しい F0 の解に至っていない段階にもかかわらず収束してしまうケースがしばしば見られた。これは，Ω と —Ω は互いに近くなるようにいずれも更新されるため，Ω の初期値に強い拘束を受けてしまうのが原因であると考えられる。この問題は w と πw についても同様であった。また，反復計算の初期段階において藤崎モデルのパラメータ推定が不安定な動作をすることが多かった。これは，藤崎モデルのパラメータ更新ステップにおいては，有音区間の Ω[l] の推定値と本来は無視するべきはずの無音区間における Ω[l] の推定値も同等の価値を置いて推定を行っていることに原因があると考えられる。以上のことから，Ω と w を周辺化
して Y から直接 πΩ と πw を推定できる方法をいずれ検討していく必要がありそうだという感触を得た。筆者は既に，変分ベイズ法に基づき Ω と w の事後分布を近似計算するアルゴリズ
ムを導出しており，これをもとに，今後上記の課題解決に向けて検討を行っていく予定である。

7. ま と め
本稿では，全極型声道スペクトルモデルと藤崎の F0 パター
ン生成過程モデルを事前分布の形として内包する音声ウェーブ
レットスペクトルの統計モデルを提案し，パラメータの効率的な最適化方法を EM アルゴリズムに基づいて導出した。この枠
組により，全極型声道モデルと藤崎モデルのパラメータを直接スペクトルから同時推定することが可能となった。

—6—

文

献

[1] F. Itakura, S. Saito, “Analysis synthesis telephony based
upon the maximum likelihood method,” In Proc. 6th Int l
Cong. Acoust. (ICA 68), C-5-5, C17–20, 1968.
[2] H. Fujisaki, “A note on the physiological and physical basis for the phrase and accent components in the voice fundamental frequency contour,” In Vocal Physiology: Voice
Production, Mechanisms and Functions, (O. Fujimura, ed.)Raven Press, pp. 347–355, 1988.
[3] A. El-Jaroudi, J. Makhoul, “Discrete all-pole modeling,”
IEEE trans. Signal Process., 39(2), pp. 411–423, 1991.
[4] 亀岡, ルルー, 大石, “音声 F0 パターン生成過程の確率モデル,”
音講論 (秋), 1-1-3, pp.207-210, Sep. 2010.

—7—