マルチカーネル線形予測モデルによる音声分析∗○亀岡弘和, 大石康智, 持橋大地, ルルージョナトン
日本電信電話 (株) NTT コミュニケーション科学基礎研究所と仮定する。x = (x1 , · · · , xI )T ，
音声信号から声道伝達特性と声帯駆動信号を推定する
10.


ための古典的な手法として線形予測法が有名である [1]。
 −a1 . .

 . .

線形予測法は，声帯による駆動音源信号が白色ガウス過.. ...
 ..


Ψ := (3).. ..
程に従うという仮定に基づき，全極モデルによって表さ.
.
−aP



..
.. ..
れた声道伝達特性を推定する手法である。しかし，実際
.
.
.の音声における有声音の声帯振動波形は局所的に周期的
0−aP · · · −a1 1であり，その周期が短い (基本周波数が高い) ほど白色と置くと式 (1) は性の仮定からの乖離が顕著となる。このため，従来の線形予測法では，声道伝達特性と声帯駆動信号の情報を必Ψx = ϵ(4)ずしも適切に分解することができず，声道伝達特性の推と行列表現でき，式 (2) より定結果が基本周波数の影響を陽に受けてしまうという問題があった。線形予測法における以上の基本周波数の影
Ψx ∼ N (0, σ 2 I)
(5)
響を回避する目的で，これまで多数の改良アプローチがが言える。よって，提案されている。その代表的なアプローチとして，声帯()駆動信号が従う分布をラプラス分布と仮定したもの [2]，x ∼ N 0, σ 2 (ΨT Ψ)−1(6)基本周波数推定を予め行い，抽出した調波成分のスペクトルピークだけに対して全極モデルをフィッティングす となるため，|Ψ| = 1 であることに注意すると，x の確るもの [3]，声帯駆動信号を隠れマルコフモデルにより 率密度関数の対数は具体的にI1
確率モデル化したもの [4]，推定声門閉鎖区間を利用しlog p(x; a) = − log 2π − 2 xT ΨT Ψx(7)22σて全極モデルを推定するもの [5]，などが挙げられる。
1I
本稿では，上述の改良アプローチと同様に従来の線形= − log 2π − 2 (x − Xa)T (x − Xa)22σ
予測モデルによる声道スペクトル推定において問題とTなっていた基本周波数の影響を解消することを目的とす と書ける。ただし，a = (a1 , · · · , aP ) ，
るだけでなく，有声無声判定と基本周波数推定も同時に
0 ··· 0.. ..

行うことができる新しい音声分析法を提案する。
.
 x1.  . .
.. 0  ..
 .
X:=(8) .
2 マルチカーネル線形予測法
 .
x1  .
..  ..
. 
2.1 線形予測法 [1]xI−1 · · · xI−P
ここではまず従来の線形予測法を概説する。線形予測法では，短時間信号 x1 , · · · , xI を，P 次の自己回帰過程 である。よって，最尤な予測係数 a は，正規方程式

1序論

xi =P∑X T Xa = X T x
ap xi−p + ϵiを解くことによって得られる。

p=1からの標本値系列と仮定する。a := (a1 , · · · , aP )T は全極モデルのパラメータであり，線形予測法においてはこれを予測係数と呼ぶ。このモデルは，ϵi を入力，A(z) =1/(1 − a1 z −1 · · · − aP z −P ) を伝達関数，xi を出力とする線形系と見ることができ，声帯駆動信号に対応する入力 ϵ = (ϵ1 , · · · , ϵI )T を，定常な白色ガウス性雑音

ϵ ∼ N (0, σ 2 I)
∗(9)

(1)

(2)

2.2

ガウシアンプロセスによる声帯駆動のモデル化

以上のように線形予測法による音声生成モデルでは声帯駆動信号 {ϵi }Ii=1 を白色ガウス雑音と仮定したが，ここでは {ϵi }Ii=1 と標本時刻 {ti }Ii=1 との関係を表す回帰式をモデル化する見地から線形予測法を捉え直す。まず，
1 (t = c )hϕh (t) :=(10)0 (t ̸= ch )Speech analysis with multi-kernel linear prediction. by Hirokazu Kameoka, Yasunori Ohishi, Daichi Mochihashi, Jonathan Le Roux (NTT Communication Science Laboratories)と定義される基底関数を用いて表される線形回帰モデル

ϵ = f (t) + ηf (t) =H∑

ϕh (t)wh = ϕ(t)T w(11)(12)は半正定値性を満たすカーネル関数の一つであるが，これは，基底関数として周期 T の周期関数を仮定していることに相当する [6]。本稿では，調波成分が等しいパワーをもつ駆動信号を仮定して基底関数をを考える。よって，ϵ = (ϵ1 , · · · , ϵI )T は，計画行列 Φ :=(ϕ(t1 ), · · · , ϕ(tI ))T を用いて

ϵ = Φw + η(13)と表される。ただし，η := (η1 , · · · , ηI )T である。ここで，重み w の事前確率を()2w ∼ N 0, σw
I(14)と仮定し，さらに η ∼ N (0, ση2 I) とすると，()2ϵ ∼ N 0, σw
ΦΦT + ση2 I(N∑h=1t − chϕh (t) :=sin 2πnTn=1)(21)と置き，そのもとで式 (18) に基づいて算出されるカーネル関数を用いる。このような周期カーネルを用いることで，ガウシアンプロセスにより周期性をもつランダム雑音を表現することが可能となる。

2.4

マルチカーネル線形予測モデル

式 (13) を式 (4) に代入することによって得られる

Ψx = Φw + η(22)(15)

が言える。H を十分大きくとり，{t1 , · · · , tI } ⊆ {c1 , · · · ,cH } となる状況を想定すると，

1 (i = j)ϕ(ti )T ϕ(tj ) =(16)0 (i ̸= j)

をカーネル線形予測モデルと呼ぶ。ここで式 (15) より，()2Ψx ∼ N 0, σw
ΦΦT + ση2 I(23)が言える。以後，K := ΦΦT と置く。Φ を，式 (21) で与えられる基底関数によって構成される計画行列とすると，このカーネル線形予測モデルは基本周期 T の音声2+ση2 と置くと， 信号を表現することができる。となるから，ΦΦT = I と書け，σ 2 := σw
さて，これまでは，声帯駆動信号の周期 T は既知とϵ ∼ N (0, σ 2 I)
(17)
見なしていたが，音声信号の基本周期は実際には未知でのように式 (2) の仮定に帰着する。以上のように，正規 ある。このため，予測係数だけでなく観測信号データか分布する結合係数による基底関数の線形和で表される回 らカーネル自体の推定も必要となる。そこで，ここでは帰モデルをガウシアンプロセス [6] といい，その確率分 マルチカーネル学習の考え方に基づき，複数の周期カー
布はカーネル関数
ネルの線形結合によって与えられるカーネル関数

k(t, t′ ) := ϕ(t)T ϕ(t′ )(18)Tk(t, t′ ) =M∑

θm km (t, t′ )(24)をもとに構成されるグラム行列 K = (Ki,j )I×I = ΦΦm=1

Ki,j = k(ti , tj )を導入し，周期 Tm の周期カーネル km (t, t′ ) の優勢度
θm を推定すべき未知パラメータと見なす。ただし，(19)によって規定される。本節では，従来の線形予測法における ϵ の仮定が，基底関数を式 (10) と定義した場合のガウシアンプロセスによる ϵ のモデル化に相当していることを見た。次節では，別のカーネル関数を導入して周期性をもつ雑音モデルを表現できることを示す。

2.3

周期カーネルの導入

M∑θm = 1(25)m=1とする。従って，{θm }M
m=1 の推定を通して音声信号に含まれる優勢な基本周期を知ることができる。km (ti , tj )を i 行 j 列の要素にもつ行列を Km とすると，先に述べたとおり，有声音の場合には声帯駆動信号はM∑
K=
θm Km(26)周期性をもつことが多い。本節では，その周期が既知でm=1あると仮定して議論する。
2.2 節で見たとおり，ガウシアンプロセスの確率分布 となり，これと式 (23) より，はカーネル関数によって決まる。カーネル関数は，必ず()M∑
2−1−T2T−1
しも基底関数の定義を介して導き出す必要はなく，半正x ∼ N 0, σw
θm Ψ Km Ψ + ση (Ψ Ψ)
定値性を満たすものであれば自由に定義できる。例えば，
m=1)/ )((t − t′ℓ2(20) が言える。以後，モデルパラメータ a, θ, σ をまとめてk(t, t′ ) := exp −2 sin2 πTΘ で表す。以上より，音声信号データ x が与えられたとc

きの，モデルパラメータ Θ の対数尤度は具体的にlog p(x; Θ))1(= − I log 2π + log |Σ| + xT ΨT Σ−1 Ψx2M∑2
Σ = σwθm Km + ση2 Iで与えられる。ただし，= は定数項以外の等号を表す。
上式に対し，x, Θ = Θ′ が与えられたときの条件つき期待値をとり，事前分布項 log p(θ) を加えると，Q 関数(27)1(log |Λ|2( −1 [ T ¯]) )+ tr Λ E yy ¯x; Θ′(34)Q(Θ, Θ′ ) = log p(θ) −(28)m=1と書ける。ただし，Σ は声帯駆動信号 ϵ の分散共分散行列のモデルである。以上の音声信号モデルをマルチカー
ネル線形予測モデルと呼ぶ。

を得る。ここで，H := [I, · · · , I] と置くと，不完全デー
タ x と完全データ y との間には

x = Hy2.5

スパース正則化パラメータ学習

ここで，モデルの過適応を防ぐ目的で，θ をスパースに誘導する事前分布として一般化正規分布

p(θ) =M∏{}λβexp − λβ |θm |β2Γ(1/β)
m=1X T Σ−1 Xa = X T Σ−1 x¯][なる関係式が成り立つことから，E yy T ¯x; Θ は¯[]E yy T ¯x; Θ = Λ − ΛH T (HΛH T )−1 HΛ+ ΛH T (HΛH T )−1 xxT (HΛH T )−T HΛT(29)2を仮定する。ただし，0 < β < 2 とする。以上より，σwと ση2 が周期成分と白色成分のエネルギーの割合を，θmが周期成分に含まれる特定周期の優勢度をそれぞれ表すパラメータとなり，予測係数 a とともにこれらの推定も同時に行うことで，基本周波数の影響を軽減した予測係数推定が行えることが期待される。
以下では，式 (27), (29) によって構成されるパラメータ事後確率を最大化する問題を考える。Θ についての同時最適解は解析的に解くことはできないが，a と θ とσ それぞれに関して順番に目的関数を最大化する操作を繰り返すことで局所最適解を探索することができる。まず，θ と σ を固定した下での a の最適解は正規方程式(35)(36)と具体的に与えられる。式 (33) の各対角ブロックに対¯[]応する E yy T ¯x; Θ′ の対角ブロックを R1 , · · · , RM +1と置くと，Q 関数は′cQ(Θ, Θ ) = −M
∑λβ |θm |β −m=1−( 2M 2)Ilog σwση θ 1 · · · θ M2M)
1 ∑ 1 ( T −1tr Ψ Km ΨRm22σw m=1 θm()
1− 2 tr ΨT ΨRM +12ση(37)と書ける。ここで，|θm |β (β < 2) に関する以下の不等式′ β−1′′ β|θm |β ≤ β|θm|(θm − θm) + |θm| (θm > 0) (38)(30)

を用いて Q(Θ, Θ′ ) の下限関数を設計することができ，補を解くことにより解析的に求まる。θ と σ については， 助関数法 [7] に基づいて σ と θ の更新式が導ける。これEM(Expectation-Maximization) 法と補助関数法 [7] を により，θm の更新式は二次方程式組み合わせることで目的関数の増加を保証する更新式を導くことができる。まず，観測信号 x を，各周期の√√成分 θ1 Ψ−1 Φ1 w, · · · , θM Ψ−1 ΦM w および白色成分Ψ−1 η に対応する M + 1 個の独立な確率変数()
2ym ∼ N 0, σw
θm Ψ−1 Km Ψ−T , m = 1, · · ·, M()
yM +1 ∼ N 0, ση2 (ΨT Ψ)−1(31)()2 β′ β−1 22−12σw
λ β|θm|θm + Iσw
θm − tr ΨT Km
ΨRm = 0の正の根で与えられる。なお，ここでは式 (25) の制約を考慮しなかったため，上式に基づいて θ を更新後，強制的に式 (25) を満たすように規格化するステップを設2けることにする。σwおよび ση2 の更新式についても解析的に得ることができ，それぞれの和に分解し，これらを完全データと扱う。よって，完TT
全データ y = (y1T , · · · , yM+1 ) に対する対数尤度は，)
1(clog p(y; Θ) = − log |Λ| + y T Λ−1 y2 2σw θ1 Ψ−1K1 Ψ−T..
.
Λ=σw2 θM Ψ−1KM Ψ−TOM)
1 ∑ 1 ( T −1tr Ψ Km ΨRm
IM m=1 θm)
1 (ση2 = tr ΨT ΨRM +1I2σw=(32)
Oση2 (ΨT Ψ)−1

(39)
(40)

で与えられる。

(33)

2.6

従来法 [2] との関係

[2] では，予測係数の更新と各標本時刻における予測残差の分散の更新を反復的に繰り返すことで，予測残差8Frequency (kHz)64

2F0 (Hz)

0500
400
300
200

Periodic component ratio100

図110.50012367

消することを目的とするだけでなく，有声無声判定と基本周波数推定も同時に行うことができるマルチカーネル
線形予測分析法と呼ぶ新しい音声分析法を提案した。動
作実験により，提案法の基本周波数推定機能，および，周期・白色成分比の推定機能を確認した。
謝辞 本研究に関して，有益な議論と情報提供を頂いた廣谷定男氏 (NTT CS 研) に感謝する。

動作実験

ATR 音声データベースの音声データ (女性話者) を用いて提案法の基本動作の確認を行った。標本化周波数 16kHz のディジタル信号に対し，フレーム長 64ms のHanning 窓を 32ms シフトさせながら分析を行った。各周期カーネルは，100Hz から 600Hz までの間を 6cent 間隔で離散化した 518 個の基本周波数値 1/Tm をもとに構成し，予測次数 P は 13 とした。また，アルゴリズムの反復計算の回数は 600 とした。
図 1 に，分析した音声信号のスペクトログラムと，提案法により推定された各時刻における θ1 , · · · , θM の分布，および，周期成分と白色成分のエネルギー比を表す
22σw/(σw+ ση2 ) の値をグラフ化したものを示す。θ1 , · · · ,θM の分布は各時刻の短時間信号においてどの基本周期が優勢であるかを表し，音声の基本周波数軌跡にピークが鮮明に現われていることが確認できる。

45

観測信号のスペクトログラム (上段) と推定した θ の分布 (中段) と周期・白色成分比 (下段)の ℓ1 ノルムを最小化する線形予測法を実現している。この手法は式 (30) において Σ を反復的に更新する点では提案法と共通しているが，その更新方法が異なる (従来法では Σ は対角行列であり，その対角要素は評価規準として各予測誤差値をどれだけ重く見積もるかを表したパラメータとなっている)。提案法とこの従来法の性質の違いについては稿を改めて議論したい。

34Time (s)まとめ

本稿では，従来の線形予測モデルによる声道スペクト
ル推定において問題となっていた基本周波数の影響を解

参考文献[1] F. Itakura and 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] C.-H. Lee, “On robust linear prediction ofspeech,” IEEE trans. ASSP, 36(5), 642–650, 1988.
[3] A. El-Jaroudi and J. Makhoul, “Discrete all-polemodeling,” IEEE trans. Signal Process., 39(2),411–423, 1991.
[4] 佐宗, 田中, “HMM による音源のモデリングと高基本周波数に頑健な声道特性抽出,” 信学論 D-II,
J84-D-II(9), 1960–1969, 2001.
[5] Y. Miyoshi et al., “Analysis of speech signals of
short pitch period by a sample-selective linear prediction,” IEEE trans. ASSP, 35(9), 1233–1240,1987.
[6] C.E. Rasmussen and C.K.I. Williams, GaussianProcesses for Machine Learning, MIT Press, Cambridge, Mass, USA, 2006.
[7] 亀岡, 小野, 柏野, 嵯峨山, “複素 NMF: 新しいス
パース信号分解表現と基底系学習アルゴリズム,”
音講論 (秋), 2-8-13, 657–660, 2008.