3-6-18
ガウス過程回帰の混合エキスパートモデルを用いた
歌声 F0 軌跡の予測と生成∗◎大石康智 (NTT), 持橋大地 (統数研), 亀岡弘和, 柏野邦夫 (NTT)

1

はじめに

本研究の目的は，歌声の基本周波数の時間変化（F0軌跡）に表れる歌唱表現を特徴抽出し，任意の楽譜に対して，その歌唱表現を反映した F0 軌跡を予測可能な生成過程モデルを構築することである。歌唱表現の定義はまだ十分になされていないが，歌唱表現が歌唱フォルマントや F0 動的変動成分のような信号特徴と関連付けられることが明らかにされた [1]。特に，ビブラートやポルタメントのような F0 動的変動成分は知覚的な観点からも，歌唱表現と密接に関係する [2]。
Fig. 1 メロディの楽譜とそれを歌った歌声の F0 軌
ビブラートは，持続音において，F0 が準周期的に変跡: MIDI 信号からノートナンバと音長を算出してメ調される成分であり，速度，振幅，継続長，波形，規ロディの楽譜とした。
則性などによって特徴付けられる。一方，ポルタメン
トは，音符の切り替わりにおいて，F0 が緩やかに変つの拡張であり，混合係数それ自身を入力変数の関数化する成分である [3]。歌唱者は，楽譜のコンテキス（ゲート関数と呼ぶ）で表現する。すなわち，入力空ト（音符の高さや長さなど）が与えられた下で，これ間（楽譜のコンテキストからなる空間）の異なる領域らの動的変動成分を巧みに制御して，歌声に表情をごとに，個々の混合要素の密度関数（エキスパートと付ける。このような生成過程を計算機で実現できれ呼ばれ，ここではガウス過程回帰に相当する）が F0
ば，歌唱者の個性や表現力を学習することにつなが軌跡の予測を行う。ゲート関数はいずれの混合要素がり，歌声の認識や合成 [4–6] への応用を期待できる。
どの領域において優勢であるかを判定する。これによ
例えば，ある歌唱者の大量の歌声から歌唱表現を学習り，歌唱者がコンテキストに依存して，歌唱表現に起することによって，どんな新しい楽譜に対しても，そ因する F0 動的変動成分を多様に使い分ける動作を表の表現を付与して合成することが可能となるだろう。
現できると考える。さらに，動的変動成分を精緻に特これまで個々の F0 動的変動成分を特徴抽出して， 徴づけるために，個々のガウス過程回帰におけるカー歌声の認識や合成，歌唱力評価に利用する試みはな
ネル関数の設計を工夫する。評価実験では，新規の楽されてきたものの [2, 3, 7–9]，動的変動成分を楽譜の譜に対する F0 軌跡の予測の観点から提案法を評価すコンテキストと結びつけて，体系的にその動特性をる。また，モデルの学習データ量について検討する。
学習する方法は十分に検討されていなかった。隠れ
マルコフモデル (HMM) に基づく歌声合成のように， 2 歌声 F0 軌跡の生成過程モデル
コンテキストラベルが付与された区間毎に，例えば，メロディの楽譜とその歌声が同期して与えられた下5 状態の HMM を用いて F0 の動特性を学習することで，楽譜に含まれる様々なコンテキストからF0 軌跡もできるが，状態内で出力確率分布が一定であるとへの回帰問題を考える。入力変数は楽譜に基づいて，いう HMM の制約のために，短時間に細かく変動す例えば下記のように構成される。
る歌声の F0 を，固定された状態数と局所的な動的特徴で表現することは難しいとされる [10]。また，未知xt = [ 音符内位置（発音開始時刻からの時間）, 音長,のコンテキストに対して頑健なモデルを構築するた先行音符との音高差, 後続音符との音高差 ]Tめに，木構造に基づくクラスタリングを行うが，実際に生成される F0 軌跡は過剰に平滑化され，本来の動もちろん，歌詞の音素情報やクレッシェンドのような特性を再現できないという問題もある。
演奏情報の有無などを入力変数に加えることも可能
本稿では，楽譜のコンテキストから歌声の F0 軌跡である。一方で，出力変数 {yt }Tt=1 には，F0 軌跡かが生成される入出力関係を，ガウス過程回帰の混合ら音符情報を取り除いた後退差分信号を利用する。差エキスパートモデル [11] を利用して直接的にモデル分信号を扱うことで，Fig. 1 に示すように，動的変動化する。混合エキスパートモデルは混合モデルの一成分を顕著に観測することができる。F0 は YIN[12]∗Mixture of Gaussian process experts for predicting and generating sung melodic contour. by OHISHI,Yasunori (NTT), MOCHIHASHI, Daichi (ISM), KAMEOKA, Hirokazu, KASHINO, Kunio (NTT)

日本音響学会講演論文集- 365 -2014年3月

Fig. 2 MoGPEs における入力空間の割り当て: 各楕
円領域は全共分散型ガウス分布を表し，この領域では同じ GPR によって出力が予測される。
を利用して 10 ms ごとに推定され，下記のように Hzで表される周波数単位を cent に変換する。

ycent = 1200 log2yHz
440 × 2 12 −53(1)まず，ガウス過程回帰 (Gaussian process regression,GPR) を説明する。学習データ D = {xt , yt }Tt=1 が与えられ，出力変数を並べたベクトル y = [y1 , . . . , yT ]Tを GPR によって回帰する場合，その確率密度関数はy ∼ GP(y; 0, K + η 2 I) = N (y; 0, K + η 2 I)

(2)のような多次元ガウス分布で表現される。ここで，K は Ki,j = k(xi , xj ) を要素に持つグラム行列，k(xi , xj ) は 2 変数間の類似性を定義するカーネル関数である。また，η 2 と I はそれぞれ，観測雑音の分散値と単位行列を表す。GPR は新しい入力変数 x∗が与えられたときの出力変数 y∗ の予測分布を推論できる。予測分布は，

p(y∗ |D, x∗ ) = N (y∗ ; µ∗ , σ∗2 ),(3)2 −1µ∗ = k Ty,∗ (K + η I)
T2σ∗ = k(x∗ , x∗ ) − k∗ (K + η 2 I)−1 k∗

提案法のグラフィカル表現

はこれまで，2 つの構成方法 [11, 13] が提案されているが，文献 [11] は MoGPEs の完全な生成モデルを定義し，入力変数の欠損や出力変数から入力変数への逆予測が可能であるため，この構成方法を利用する。
生成過程の流れを以下に示す。

1. ディリクレ−多項分布モデルを用いて，T 個の入力変数を R 個の領域のいずれかに割り当てる。これは指示変数集合 {zt }Tt=1 によって表現される。
2. 指示変数の集合 {zt : zt = r} が与えられた下で，領域 r の入力変数の密度分布パラメータ θrx ≡{µr , Σr } が生成される。 ここで，密度分布は全共分散型ガウス分布を仮定する。
3. θrx が与えられた下で，領域 r に属する入力変数集合 Xr ≡ {xt : zt = r} が生成される。
4. カーネル関数のパラメータ θrGP ≡ {ψ r , wr2 , ηr2 }が領域毎に生成される。
5. 領域毎に Xr と θrGP を使ってグラム行列が計算され，出力変数を並べたベクトル y r ≡ {yt : zt =r} が生成される。
Fig. 3 はこのグラフィカル表現であり，同時分布は，x R
p({xt , yt }Tt=1 , {zt }Tt=1 , {θrGP }R
r=1 , {θr }r=1 |Ω)となり，k∗ は k(xt , x∗ ) (t = 1, . . . , T ) からなるベクトルを表す。学習データが “事例” となって，新たな入力変数に対する予測に貢献することがわかる。
カーネル関数の必要条件はグラム行列が半正定値対称行列となることである。出力信号が定常であることを仮定し，Squared exponential (SE) 共分散関数が一般的に利用されるが，Fig. 1 に示すように，差分信号は必ずしも定常ではない。歌唱者は音符の並びに依存して，様々な動的変動成分を使い分けながら歌唱するためである。このような非定常な動特性をモデル化するために，ガウス過程回帰の混合エキスパートモデル (Mixture of Gaussian process experts,MoGPEs) を利用する。
MoGPEs では，入力空間がゲート関数によって確率的に分割される (Fig. 2)。それぞれの領域では，異なる GPR によって，動的変動成分が特徴付けられて，出力が予測される。最終的に，出力は個々のゲート関数と GPR との混合モデルとして表現される。MoGPEs

日本音響学会講演論文集

Fig. 3=(4)

R
∏
[ x]
p(θr |Ω)p(Xr |θrx )p(y r |Xr , θrGP , Ω) × p({zt }Tt=1 |Ω)
r=1と書ける。R と Ω はそれぞれ，エキスパートの総数と超パラメータ集合を表す。本稿では，指示変数集合{zt }Tt=1 を積分消去せず直接表現して，これらの依存関係を考慮する。式 (4) における個々の分布を以下に示す。

p({zt }Tt=1 |Ω) =R
Γ(α) ∏ Γ(Tr + α/R)Γ(α + T ) r=1 Γ(α/R)

p(θrx |Ω) = N (µr ; m0 , Σr /β0 )W(Σ−1r ; W 0 , ν0 )p(Xr |θrx ) = N (Xr ; µr , Σr )()2p(y r |Xr , θ GPr , Ω) = GP y r ; 0, K r + ηr I r(5)

ここで，α, Tr , I r , ηr2 , W はそれぞれ，ディリクレ−多項分布の超パラメータ，Xr の要素数，Tr × Tr の単位行列，観測雑音の分散値，ウィシャート分布を示す。

- 366 -

2014年3月

となり，右辺の第 1 式が GPR とゲート関数の混合モデルとなることがわかる。ここで，µ∗ =R∑2−1cr kT
yrr,∗ (K r + ηr I r )r=1

Fig. 4

マルチカーネル学習に基づくグラム行列表現

本稿では，マルチカーネル学習に基づき，Fig. 4 のような複数のカーネルの線形結合

kr (xi , xj ) = wr2

M∑

ψr,m kr,m (xi , xj )

(6)

m=1

でグラム行列 K r を表現し，強度 wr2 と各カーネルの優勢度 ψr,m を推定すべき未知パラメータとみなす∑M[14]。ここで，xi , xj ∈ Xr ,m=1 ψr,m = 1 とする。
M は線形結合するカーネル関数の総数である。単位の異なる様々なコンテキスト（音符内位置や音高差，音長など）を扱うため，個々のカーネル関数 kr,m (xi , xj )を 2 つのカーネル関数の積で表現する [15]。
(p)

(p)

(c)

(c)

(p)
(c)
kr,m (xi , xj ) = kr,m(xi , xj )kr,m(xi , xj ) (7)
(p)

(p)

(p)

(c)

(c)

(c)

ここで， kr,m (xi , xj ) と kr,m (xi , xj ) はそれぞれ，音符内位置の類似を表す位置カーネルと音符の類似度を表すコンテキストカーネルと呼ぶ。入力変(p)(c)
数ベクトル xi を xi , xi のように位置に関する変数と音符のコンテキストに関する変数に分け，カーネル関数を計算する。位置カーネルは動的変動成分の時間的連続性や周期性を捉えるために，SE 共分散関数と周期共分散関数を利用する。
()(p)
(p)
(p)
(p)
(xi − xj )T (xi − xj )
(p)
(p)
(p)
kr,m (xi , xj ) = exp −(p) 22lm(())
(p)
lm(p)
(p)
(p)
(p)
2(p)
kr,m (xi , xj ) = exp −2 sin(x − xj )2π iσ∗2 =R∑

2−1c2r (kr (x∗ , x∗ ) − kTkr,∗ )r,∗ (K r + ηr I r )r=1

cr = p(z∗ = r|x∗ , θrx )であり，ゲート関数 cr はベイズの定理を利用して容易に計算できる。この予測分布に基づいて，F0 軌跡
Tm{ft }t=1を以下のように音符毎に再合成する。

ft =ot + c − ō, ot =t∑n=1µ∗,n , ō =Tm∑

ot /Tm (9)t=1

ここで，c と Tm はそれぞれ，楽譜に記載される音符の音高と音長を表す。
GPR はこれまで音声合成 [15] や F0 生成 [16]，声質変換 [17] などに適用されており，合成音の観点から性能の向上が確認されている。提案法はこれらの研究と関連するが，動的変動成分の動特性を精緻に特徴づけるために，MoGPEs と多様なカーネル関数を利用することが特徴的である。また，MoGPEs を利用すれば，通常の GPR よりも計算コストを低減できることも特筆すべき点である [11]。

3パラメータの推論

MCMC-EM アルゴリズムを用いて，すべてのパラメータを推論する。具体的には，{zt }Tt=1 と {θrx }Rr=1の推論のためにギブスサンプリングを利用する [11]。
また，式 (5) の多次元ガウス分布を独立なガウス分布の和とみなして，{θrGP }Rr=1 の推論のために EM アルゴリズムを利用する [14]。

コンテキストカーネルには SE 共分散関数を利用する。
4 評価実験
()
(c)
(c) T
(c)
(c)
(xi − xj ) Λ(xi − xj )
(c)
(c)
(c)
新たな入力（楽譜）に対する出力（F0 軌跡の差分信kr,m(xi , xj ) = exp −2号）の予測性能の観点から提案法を評価する。
「RWC
2 (c) 22(c)(c)
研究用音楽データベース: ポピュラー音楽」(RWCΛ−1 = diag(lm,1 , lm,2 , . . . , lm,Dc )MDB-P-2001) における 10 楽曲 (No. 38, 39, 42, 44,
まとめると，パラメータと超パラメータはそれぞれ， 45, 46, 64, 72, 74, 76) の歌声の F 軌跡を手作業でラ0xGPΘ = {z1 , . . . , zT , θ1x , . . . , θR, θ1GP , . . . , θR}， Ω =ベル付けした結果[18]を利用する。また，そのMIDI(p)(p) (c)(c){α, m0 , W 0 , β0 , ν0 , l1 , . . . , lM , l1,1 , . . . , lM,Dc } と信号を楽譜として利用する。これらの楽曲はすべて同(c)なる。Dc は xi の次元数を表す。
一歌手によって歌われたものである。本来ならばこれ新たな入力変数 x∗ に対する出力変数 y∗ の予測分
らの楽曲の音響信号から F0 を推定すべきであるが，布を考える。式 (3) を参考にすると，MoGPEs では， 今回は提案法の性能の上限を調べるためにこのようなデータを用いた。xt と F0 の後退差分値 yt は 10 msp(y∗ |D, x∗ , Θ, Ω)(8)ごとに計算され，無声音のため F0 が観測されない区
R∑GPx
間は除去した。楽曲 No. 38, 39, 42, 44, 45 はモデル
=p(y∗ |y r , Xr , x∗ , z∗ = r, θr )p(z∗ = r|x∗ , θr )r=1
パラメータの学習に，楽曲 No. 46, 64, 72, 74, 76 は= N (y∗ ; µ∗ , σ∗2 )予測性能の評価に利用した。学習データと評価データ

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

2014年3月

RRMSE に基づく予測性能の比較
1020304050

Table 2 学習データの量に対する RMSE
学習データの曲数
135MoPRs
MoGPEs

71.8
25.0

MoGPEsTable 1

59.7
24.0

73.1
23.9

79.6
23.1

56.5
22.3

Fig. 6

Fig. 5

予測結果の比較の量はそれぞれ，649.3 秒，625.6 秒である。 指示変数の初期値は，学習データの入力ベクトル集合 {xt }を k-means クラスタリングすることによって得られた割当て結果を利用した。各領域に割当てられた入力ベクトルを利用して，θrx を初期化した。M = 216とし，θrGP の初期値はそれぞれ，wr2 = 100, ψr,1 =1/M, . . . , ψr,M = 1/M, ηr2 = 10 (r = 1, . . . , R) とした。超パラメータの設定を脚注に示す1 。パラメータの推論回数は 50 回とした。
評価尺度として，実際の出力値系列と式 (8) の予測平均値系列間の二乗平均平方根誤差 (Root meansquare error, RMSE) を利用した。
√∑Tt2RMSE =(10)t=1 (yt − µ∗,t ) /Tt ,ここで，Tt は評価データの長さを表す。比較として，多項式回帰の混合エキスパートモデル (Mixture ofpolynomial regression experts, MoPRs) を定義した。
式 (5) の GPR が n 次の多項式回帰に置き換わる。
Tab. 1 はエキスパートの数 R を変化させたときのRMSE の平均値を示す。この結果は楽曲 No. 35 だけを利用してパラメータを学習し，評価した結果である。MoPRs の多項式の次数は n = 8 とした。これより高次の多項式を利用した場合，多項式の重み係数の推定が安定しなかったためである。RMSE の観点か
ら，提案法は多項式回帰からなる単純な手法と比較して，予測性能が上回ることを確認した。また R を増やすにつれて予測性能は向上した。Fig. 5 は 2 つの手法の予測結果を比較する。動的変動成分の動特性を表現するために，周期カーネルを利用する提案法の有効性を確認できた。Tab. 2 は学習データの量に対1 入力ベクトルの次元数 D は 4 であり，α = 1, β = 0.1, ν =00
D + 1 に固定する。m0 と W 0 はそれぞれ，学習データのすべての入力ベクトルの平均ベクトルおよび共分散行列の逆行列を ν0で割算した行列に設定する。位置カーネルの超パラメータとして，(p)SE カーネルは {lm , m = 1, . . . , 108|0.05, 0.11, 0.23, 0.5}，周期(p)カーネルは {lm , m = 109, . . . , 216|0.13, 0.15, 0.17, 0.2} から選(c)択した。コンテキストカーネルの超パラメータは {lm,1 |1, 2.2, 5}，
(c)

22.3

20.7

20.5

メロディの楽譜とそれを歌った歌声の F0 軌跡する予測性能を検証する。R = 50 に固定した。実験
結果より，学習データを増やすに連れて RMSE が小さくなることを確認した。Fig. 6 は式 (9) によって再合成される F0 軌跡を示す。正解となる実際の F0 軌跡とほぼ同じ軌跡が得られるが，細かい動特性に関してはまだ予測が十分でない。予測性能を高めるた
めに，エキスパートの数や学習データ量，入力コンテキストの種類を今後さらに検討する予定である。

5

おわりに歌声の F0 軌跡に表れる様々な動的変動成分を特徴抽出して歌唱者の歌唱表現を学習し，任意の楽譜に対して，その歌唱表現を反映した F0 軌跡を予測する生成過程モデルを提案した。実験結果より，単純な多項式回帰を用いるよりも GPR の混合エキスパートモデルを利用することの有効性を確認した。予測性能を高めるために，学習データを増やすこと，エキスパートの数や超パラメータを調整すること，エキスパートの割当てにディリクレ過程を導入することが挙げられる。この枠組を歌唱者認識や歌唱表現認識，歌声合成に応用することも検討中である。

参考文献[1] Sundberg, “The Science of the Singing Voice”
Northern Illinois University Press, 1987.
[2] Saitou et al., in Proc. WASPAA 2007.
[3] Regnier, “Localization, Characterization andRecognition of Singing Voices” IRCAM / UPMC
in Paris, France, 2013.
[4] Oura et al., in Proc. ICASSP 2012.
[5] Nose et al., in Proc. INTERSPEECH 2013.
[6] Doi et al., in Proc. INTERSPEECH 2013.
[7] Nakano et al., in Proc. INTERSPEECH 2006.
[8] Migita et al., in Proc. ICA 2010.
[9] Ohishi et al., in Proc. INTERSPEECH 2012.
[10] Zen et al., Speech Communication, Vol. 51, No. 11,pp. 1039–1064, 2009.
[11] Meeds et al., in Proc. NIPS 2006.
[12] de Cheveigné et al., JASA, Vol. 111, No. 4, pp.
1917–1930, 2002.
[13] Rasmussen et al., in Proc. NIPS 2002.
[14] Yoshii et al., in Proc. ICASSP 2013.
[15] Koriyama et al., in Proc. INTERSPEECH 2013.
[16] Fernandez et al., in Proc. ICASSP 2013.
[17] Pilkington et al., in Proc. INTERSPEECH 2011.
[18] Goto, in Proc. ISMIR 2006.

(c){lm,2 |1, 2.2, 5}，{lm,3 |0.1, 0.55, 3} の組合せで設定した。

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

2014年3月