3-P-31
複数振動基底に基づく歌声の F0 動特性の統計的モデリング∗◎大石康智, 亀岡弘和, 持橋大地, 永野秀尚, 柏野邦夫 (NTT CS 研)

1

はじめに歌声の声の高さを表す基本周波数（F0）系列には，階段状に変化する旋律成分と，ビブラートやオーバーシュートのような歌唱者の歌唱動作を表す動的な変動成分が複雑に重ね合わされている（図 1）．この F0 動的変動成分は，歌声知覚および個人性知覚に寄与する特徴であることが知られており [1, 2]，歌唱スタイルの記述やそれを利用した歌声検索，歌声評価 [4, 3, 8]のための有用な尺度となる．また，表情豊かかつ多様な歌声合成 [5, 6, 7] を実現するためにも必要不可欠な成分と言える．
これまで多くの研究では，2 次系モデルを利用して，歌唱の F0 動的変動成分を表現した [1, 9, 10]．これらの研究では，日本語の話声の F0 パターンを表現する
藤崎モデル [11] を参考にする．藤崎モデルは，臨界制動 2 次系のインパルス応答とステップ応答を利用して，日本語の句頭から句末に向けて緩やかに下降するフレーズ成分と，語句に対応して上昇下降するアクセント成分を表現し，これらを重畳することで，F0 系列を記述する．ただし，歌声の旋律に伴った急激な F0 の上昇・下降の制御及び，ビブラートのような周期的な振動は，臨界制動系では表現できないた
め，歌声の F0 制御モデルでは 2 次系の伝達関数

H(s) =Ω2
s2 + 2ζΩs + Ω2

観測F0系列2

旋律成分

4

ノート間の音高差
時間

68[sec]

歌唱動的変動成分
線形2次系
ガウス性
白色雑音

00時間
時間
入力ステップ信号： インパルス応答：観測信号：
系の出力信号：

00時間
残差信号：

Fig. 1 想定する歌声の F0 生成過程の概略図
ングの問題を解消するために，2 つの方策を導入した
歌声の F0 生成過程の新しい確率モデルを提案する．

1. モデルの自由度を効果的に下げる目的で，2 次系のインパルス応答 h(t) を，あらかじめ用意したいくつかの振動基底の疎（スパース）な線形和によって構成する．
2. ノートが切り替わる時点を始点終点と考え，F0系列をいくつかのセグメントに分割し，セグメントごとに，図 1 に示す信号の生成過程を仮定
する．h(t) と ²(t) はどちらも動的変動成分を表す信号であるが，h(t) はノートの切り替わりやオーバーシュートを表現する信号の大局的動特性を，²(t) はビブラートのような音高が安定するときの局所的動特性を表現する．

(1)の減衰率 ζ によって表現される，指数減衰 (ζ > 1)，減衰振動 (0 < ζ < 1)，臨界制動 (ζ = 1)，定常振動(ζ = 0) を利用する [1]．齋藤らは，減衰率 ζ と固有周波数 Ω を聴取実験に基づいて手動で調整し，それらによって得られる式 (1) のインパルス応答を，階段状に変化する信号に畳み込んで生成される F0 系列を利用して，表情豊かな歌声合成音を実現した [1]．
これに対し，我々は旋律成分および 2 次系の制御パラメータ ζ, Ω がいずれも未知の下で，観測されるF0 系列だけから，これらを同時に推定する統計的手法の枠組を提唱した [12]．これは，旋律成分を表す隠れマルコフモデル（HMM）と，差分近似に基づく式(1) のパラメトリックな表現によって，最尤のモデルパラメータを反復推定するアルゴリズムであった．しかし，このアルゴリズムでは，モデルパラメータの推定性能が悪かった．この理由は，図 1 の上部に示す
ように F0 系列には各ノート（音符）の切り替わりやオーバーシュート，ビブラートのような動的変動成分が所々で混在するのに対し，自由度が高すぎるモデルで，これらの動的変動成分を学習しようとしていたためである．本報告は，このようなオーバーフィッティ∗6000]t 5600
ne5200[c 48000F
4400
4000 0このような信号表現方法に基づき，パラメータ最適化アルゴリズムを導く．そして，パラメータの推定性能を評価し，推定されたパラメータの歌唱者ごとの違いから歌唱スタイルの違いについて分析する．

2複数の振動基底を利用した
線形 2 次系のインパルス応答表現

式 (1) のラプラス逆変換によって得られるインパルス応答は，ζ の値によって場合分けされる．
 −ζΩt ¡ √√ 2¢2Ωe√ 2e ζ −1Ωt − e− ζ −1Ωt , (ζ > 1)
2 ζ −1p¢−ζΩt ¡ Ωe√ 2 sin( 1 − ζ 2 Ωt ,(0 < ζ < 1)1−ζh(t) =Ω2 te−Ωt ,(ζ = 1)Ω sin(Ωt),(ζ = 0)
これらのインパルス応答をサンプリング周期 ∆ で離散化すると，図 1 の系の入出力関係は y = Φf と記述できる．ここで，y = [y1 , y2 , . . . , yN ]T , f =[f1 , f2 , . . . , fN ]T は，出力信号 y(t) と入力ステップ信

Statistical Modeling of F0 Dynamics in Singing Voices Based on Multiple Oscillation Bases. by OHISHI,Yasunori, KAMEOKA, Hirokazu, MOCHIHASHI, Daichi, NAGANO Hidehisa, KASHINO Kunio (NTTCommunication Science Laboratories)

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

号 f (t) をサンプリング周期 ∆ で離散化した時系列信号のベクトルを表す（N は信号の長さ）．この Φ が，系のインパルス応答を表し，例えば，ζ = 1 の場合，
Φ は下三角行列
 2 −Ω∆Ω ∆e2

 2Ω ∆e
Φ=.
.

−2Ω∆.
N Ω2 ∆e−N Ω∆Ω2 ∆e−Ω∆
..
.
...



..

.
2Ω2 ∆e−2Ω∆Ω2 ∆e−Ω∆Φ−1 ' w1 Υ(1) + w2 Υ(2) + . . . + wI Υ(I)(2)

ここでは，予め手動で ζ, Ω を決定し，I 個の振動現象を表すインパルス応答 {Φ(1) , Φ(2) , . . . , Φ(I) } を計算する．そして，これらの逆行列 Υ(i) := (Φ(i) )−1（逆フィルタのインパルス応答を表す．以後，振動基底と呼ぶ）の重み付き和で Φ−1 を近似する．ただし，
この重みパラメータ w := {w1 , w2 , . . . , wI } は疎（スパース）となるように正則化する．これは，Φ−1 が，ある限られた種類の振動基底のみによって表現されることを意味し，モデルの自由度を効果的に下げる目的として，このような操作を行う．後に説明するが，これは w の事前確率を課することで実現される．したがって，系の入出力関係は以下のように表現される．

(w1 Υ(1) + w2 Υ(2) + . . . + wI Υ(I) )y = f(3)便宜上 Ψ := w1 Υ(1) + w2 Υ(2) + . . . + wI Υ(I) とおく．

2 次系 F0 生成過程の統計的モデリング
2 次系の入出力関係を統計的にモデル化する．

3.1

入力信号と出力信号の確率モデル

入力信号 f はノート間の音高差を表すよう，ステップ信号を想定する．そのために，常に同じ値をもつベクトル u = [u1 , . . . , uN ]T = u[1, 1, . . . , 1]T = u1を用意する．ここで，スカラー値 u は音高差を表す
パラメータ，1 は N 個の 1 の値が並ぶベクトルとする．このベクトル u を平均とする多次元ガウス分布
N (u, αI N ) から生成される確率変数として，入力信号 f を表現する．α は分散を表す超パラメータであり，あらかじめ手動で値を設定する．I N は N × Nの単位行列を示す．
出力信号 y は，ガウス分布に従う変数集合 f の線形結合（y = Ψ−1 f ）であるから，y 自身もガウス分¡¢布に従い，y ∼ N Ψ−1 u, αΨ−1 (Ψ−1 )T と書ける．

3.2

尤度関数と事前確率

残差信号 ² = [²1 , ²2 , . . . , ²N ]T は，ガウス性白色雑音 N (0, βI N ) に従うとし，観測 F0 系列 o =[o1 , o2 , . . . , oN ]T は，系の出力信号 y に残差信号 ² が加わった信号 o = y + ² と仮定する．ここで，β は残差信号の分散を表す超パラメータである．y と ² は互いに独立であると仮定すると，観測信号 o が与えら

日本音響学会講演論文集

N
1log(2π) − log |Σ|22
1T −1− (o − µ) Σ (o − µ)
(4)
2
¡
¢
µ = Ψ−1 u, Σ = αΨ−1 (Ψ−1 )T + βI Nlog P (o|Θ) = −0となる．しかし，h(t) は，上式のように複数の場合からなるので，行列 Φ を以下のように構成する．

3れたときのモデルパラメータ Θ := {w, u, β} の対数尤度関数は，以下のように書ける．

Θ の事前確率 P (Θ) は，各要素の独立性 P (Θ) =P (w)P (u)P (β) を仮定し，u と β はそれぞれ一様分布に従うとする．パラメータ w の要素には疎な制約をもたせるため，その事前確率は一般化正規分布
P (w) =IYppλp
exp−λ |wi |2Γ(1/p)
i=1(5)に従うとする．p, λ は一般化正規分布の形状を規定する定数であり，0 < p < 2 のとき p(w) は優ガウス的となり，スパースネスを測るための尺度となる．

4

パラメータ最適化アルゴリズム

観 測 F0 系 列 o が 与 え ら れ た と き ，事 後 確 率
P (Θ|o) ∝ P (o|Θ)P (Θ) を最大化するパラメータ Θを決定したいが，°1 o が出力信号 y と残差信号 ² の和で構成される，°2 尤度関数が w に関して非線形となるため解析的に求めることは難しい．そこで，°1EM 法 [13] を適用して，その E-step で，o を y と ²に分離する，°2 M-step に補助関数法 [14] を適用して，Q 関数の補助関数を設計する，からなる 2 つの方策に基づいて，パラメータ最適化アルゴリズムを導出する．

4.1

完全データの定義

この事後推定問題に EM 法を適用する際の最初のステップは完全データを定義することである．ここでは，y と ² を完全データ x = [y T ²T ]T と見なして，EM 法を適用する．不完全データと完全データの関係は，o = Hx となり，ここで，H := [I N I N ]とする．x と現在のパラメータ Θ0 が与えられたときの，対数尤度関数の条件付き期待値を計算し，さらにlog P (Θ) を加算すると，以下のような Q 関数を得る．

Q(Θ, Θ0 ) =c¡¢1hlog |Λ−1 | − tr Λ−1 E[xxT |o; Θ0 ]2i+ 2mT Λ−1 E[x|o; Θ0 ] − mT Λ−1 m + log P (Θ)"#!
· −1 ¸T1ΨΨ0Ψ um :=, Λ−1 := α(6)100β INÃとなる．ここで，tr(·) は行列のトレースを表し，
E[x|o; Θ0 ] と E[xxT |o; Θ0 ] は，条件付きガウス分布の性質より，- 402 -

E[x|o; Θ0 ] = m + ΛH T (HΛH T )−1 (o − Hm)
E[xxT |o; Θ0 ] = Λ − ΛH T (HΛH T )−1 HΛ+ E[x|o; Θ0 ]E[x|o; Θ0 ]T (7)2010年9月

と書ける．後の計算のため，y, ² に対応するように
E[x|o; Θ0 ] と E[xxT |o; Θ0 ] を区分表現する．
· ¸·¸x̄yRy ∗0T0E[x|o; Θ ] =, E[xx |o; Θ ] =x̄²
∗ R²

一方，f + (Θ, w̄, γ) を u, β それぞれに関して微分
して 0 とおくと，ここで，x̄y と x̄² はどちらも長さ N のベクトルであり，Ry と R² はどちらも N × N の正方行列を表す．

が得られ，u と β の更新式となる．パラメータ最適化
アルゴリズムの流れを以下にまとめる．
初期化: パラメータ Θ に初期値を与える．
E-step: 条件付き期待値 E[x|o; Θ0 ], E[xxT |o; Θ0 ]，補助変数 w̄, γ の更新．
M-step: 式 (11)，(12) から，パラメータ Θ の更新．
収束判定: 式 (10) の値が収束していなければ，Θ0 = Θとして E-step に戻る．

4.2

M-step 更新式

式 (6) からパラメータ集合 Θ に関連する項を取り出し，最大化する目的関数を，

f (Θ) := −NI³X´XNlog αβ +log
wi Υ(i)n,n2n=1
i=151 T1u Ψx̄y −tr(ΨT ΨRy )α2αIX
11 T|wi |p−tr(R² ) −u u − λp2β2αi=1+(8)(i)と改めて定義する．ここで，Υn,n は Υ(i) の n 行 n 列目の対角要素を表す．この最大化問題を解く更新式を補助関数法 [14] により導く．式 (8) で与えられる目的関数の補助関数を以下の 2 つの不等式（証明略）を用いて導く．
NXlogI³X

N XI(i)´ X
wi Υn,nwi Υ(i)≥γlog
i,nn,nγi,nn=1 i=1n=1ip|wi | ≤ p|w̄i |p−1 wi + |w̄i |p − p|w̄i |p ,(0 < p ≤ 1)
(9)

ここで，補助変数 w̄ := {w̄1 , w̄2 , . . . , w̄I }, γ :={γ1,1 , . . . , γI,N } を定義する．式 (9) の二つの不等式を式 (8) に代入すると，

f + (Θ, w̄, γ) := −N X
I(i)XNwi Υn,nlog αβ +γi,n log
2γi,nn=1 i=11 T111 Tu Ψx̄y −tr(ΨT ΨRy ) −tr(R² ) −u uα2α2β2αI ³´X− λpp|w̄i |p−1 wi + |w̄i |p − p|w̄i |p(10)+i=1を得る．このとき，f (Θ) ≥ f + (Θ, w̄, γ) が成り立ち，等号成立は，w̄i = wi , γi,n = PIw̄i Υ(i)n,ni0 =1(i0 )w̄i0 Υn,n, (i =1, 2, . . . , I, n = 1, 2, . . . , N ) のときであるため，式(10) は補助関数の定義を満たす．
式 (10) を wi0 に関して微分して 0 とおくと，
I01 X ¡ T (i) T (i0 ) ¢1tr Ry ΥΥwi − uT Υ(i ) x̄yα i=1αp+ λ p|w̄ |i0p−1−NXγi0 ,nn=1

wi0= 0 (11)を得る．ただし，式 (11) は，i0 = 1, 2, . . . , I として，w1 , w2 , . . . , wI に関する非線形な連立方程式となるた
め，Coordinate descent 法 [15] を利用して解く．

日本音響学会講演論文集

u=1 T
1 Ψx̄y ,Nβ=1 ¡ ¢tr R²
N(12)

評価実験

提案するパラメータ最適化アルゴリズムの有効性を評価する．まず，オーバーフィッティングの問題を改善できるかどうか確認するために，提案法によってF0 系列から分解されたステップ信号とインパルス応答信号を評価する．そのために，ランダムに決定された ζ, Ω, u（値の範囲は，0 < ζ < 2，0.05 < Ω < 0.3，0 < u < 200 とする）を使って，人工的に 100 個の擬似 F0 系列を作成した．ここで，N, α, β はそれぞれ 300，2，100 とした．そして，F0 系列ごとにアルゴリズムを適用して，パラメータを推定する．前処理として，まず Υ を作成するために，ζ は 0 から 2までの間を 0.02 刻みで，Ω は 0.05 から 0.3 までの間を 0.005 刻みで変化させた．その結果，I = 3100 となる．したがって，wi = 1/I を初期値とした．また，
β = 100 を初期値とした．u の初期値は，F0 系列の中央値とした．α = 2，λ = 10000，p = 0.8 に固定する．これらの初期値は実験的に決定した．
評価方法として，疑似 F0 系列のステップ信号の作成に用いた u と推定された u の二乗平均平方根誤差(RMSE) を計算する．同様に，疑似 F0 系列の作成に用いた ζ, Ω から計算されるインパルス応答と，推定された w から計算されるインパルス応答の RMSE を計算する．ステップ信号とインパルス応答への分解性能が高ければ，オーバーフィッティングを解消できていると考える．表 1 は，それぞれの RMSE のすべての疑似 F0 系列に対する平均値を従来法 [12] と比較する．従来法に比べて値が小さくなり，モデルの自由度を下げた提案法の有効性を確認できた．
次に，実際の F0 系列を利用してパラメータの推定性能を定性的に評価する．声楽家，プロのポップス歌手，音楽的な訓練を受けていない素人 (それぞれ男女 1 名ずつの計 6 名) の歌声からなるデータベースを利用する．歌唱者は，きらきら星と喜びの歌を譜面を見ながら伴奏なしで歌唱した．2 パターンの日本語歌詞とハミングによる歌声を収録した．歌声信号は
計 36 サンプルとなる．F0 は，YIN[17] を利用して，5ms ごとに推定される (∆ = 5ms)．そして，半音が100cent となるように，Hz を cent の単位に変換する．
また，300ms 以内の F0 が推定されない無声音の区間は近傍の F0 値で線形補間を行う．本実験ではモデル- 403 -2010年9月

Table 1

ステップ信号に関する RMSE

インパルス応答に関する RMSE

21.7
0.67

5.29×10−2
2.98×10−4

従来法
提案法

ポップス歌手（女性）
声楽家（男性）0.16
素人（女性）
0.14 素人（男性）ポップス歌手
（男性）
0.12
0.46
0.50
0.54
0.58
0.62
過減衰
減衰振動
0.18

声楽家（女性）

6200

t]n5800
ce[0F5400

Fig. 3 歌唱者ごとに推定された ζ と Ω の平均値
60

5000

素人（男性）

5000t]nec4600[0F4200
3800
0

声楽家（女性）

0.20

二乗平均平方根誤差 (RMSE) による評価

24

時間 [sec]6

声楽家（女性） ポップス歌手（女性）素人（女性）

40度頻
200-400 -200 0 200 400 -400 -200 0 200 400 -400 -200 0 200 400[cent][cent][cent]8

Fig. 2 観測 F0 系列 o と生成 F0 系列 µ の比較

Fig. 4 歌唱者ごとに推定された u の頻度分布の性能を評価するために，図 1 に示すように，F0 系列に対してセグメントを手作業でラベル付けし，各セグメントごとにパラメータを推定する．歌声データ
ベース全体で，1323 個のセグメントをラベル付けした．前処理として，セグメントの最初の F0 値を，そのセグメントにおけるすべての F0 値から減算した．
図 2 は，各セグメントで推定されたパラメータw, β, u から計算される式 (4) の µ（生成 F0 系列）と観測 F0 系列 o を比較する．声楽家と素人のどちらも，生成 F0 系列は，オーバーシュートやポルタメントのようなノートが切り替わるときの大局的な動特性を表現できている．素人に比べて，声楽家の観測F0 系列には所々にビブラートが観測される．提案モ
デルでは，このような変動成分はすべて，ガウス性白色雑音に従うとし，分散パラメータ β によって表現される．この残差信号 ² をより精緻にモデル化して，ビブラートのような周期的な変動成分をも特徴抽出することは今後の課題である．例えば，近年機械学習の分野で注目を集める Multiple Kernel Learning を利用したガウス過程 [16] に基づく信号表現との類推から，² が従う確率分布の共分散行列を周期カーネルによって表現する方法が考えられる．
図 3 は，各セグメントで推定される w の要素の最も大きい値 wi に着目し，それに対応する ζ と Ω の，すべてのセグメントにわたる平均値を歌唱者ごとに計算した結果である．ζ が小さな値であることは，その振動現象がオーバーシュートのような減衰振動であることを意味する．一方で，Ω が小さな値であることは，ノートの立ち上がり時間が長いことを意味する．
したがって、素人歌唱者は歌唱技術が乏しいため，ζと Ω が他の歌唱者に比べて，小さい値になったと考えられる．図 4 は各セグメントで推定された u の歌唱者ごとの頻度分布を示す．声楽家やポップス歌手は，半音（100cent）の整数倍の位置にピークが観測されるが，素人は，そのピークが不鮮明となった．したがって，素人は正確な半音単位の音階で歌唱することが難しいと見ることができる．将来は，大規模歌声データベースを利用して，推定される ζ, Ω, u から歌唱スタイルについて詳細に分析する予定である．

日本音響学会講演論文集

6

おわりに歌声の F0 動特性のモデリングとそのモデルパラメータ推定に関する新しい解法を検討した．これを応用した歌声合成インタフェース [18]，さらに話声のF0 生成過程の確率モデル [19] について現在検討中である．今後は，提案法を多変量化して MFCC 信号に適用し，声質の動特性の制御についても検討する予定である．

参考文献[1] Saitou, T. et al., Proc. WASSPA 2007, pp. 215–218,
2007.
[2] Saitou, T. et al., Proc. EUROSPEECH 2009, pp. 832–835, 2009.
[3] Kako, T. et al, Proc. ISMIR 2009, pp. 393–397, 2009.
[4] Nakano, T. et al., Proc. ICSLP 2006, pp. 1706–1709,2006.
[5] Sundberg, J., Advances in Cognitive Psychology. Special
issue on Music Performance, Vol. 2, No. 2-3, pp. 131–143, 2006.
[6] Bonada, J. et al., Proc. SMAC 2003, 2003.
[7] Nakano, T. et al., Proc. SMC 2009, pp. 343–348, 2009.
[8] Proutskova, P. et al., Proc. ISMIR 2009, pp. 759–764,2009.
[9] 柏野邦夫ほか, 音講論集, 2-9-1，pp. 625–626 1998.
[10] Minematsu, N. et al., Proc. SpeechProsody 2004, pp.
487–490, 2004.
[11] Fujisaki, H., Vocal Physiology: Voice Production, Mechanisms and Functions, (O. Fujimura, ed.), Raven Press,
pp. 347–355, 1988.
[12] Ohishi, Y. et al., Proc. ICSLP 2008, pp. 139–142, 2008.
[13] Feder, M. et al., IEEE Trans. ASSP, Vol. 36, No. 4, pp.
477–489, 1988.
[14] Kameoka, H. et al., Proc. ICASSP 2009, pp. 3437–3440,2009.
[15] Meng, X. L. et al., Biometrika, Vol. 80, pp. 267–278,1993.
[16] Rasmussen, C. E. et al., Gaussian Processes for Machine
Learning, MIT Press, Cambridge, Mass, USA, 2006.
[17] de Cheveigné, A. et al., JASA, Vol. 111, No. 4, pp. 1917–1930, 2002.
[18] 大石康智ほか, 情処研報音楽情報科学, 2010–MUS–86–9,2010.
[19] 亀岡弘和ほか, 音講論集, 1-1-3, 2010.

- 404 -2010年9月