1-8-19
歌声 F0 系列からの楽譜逸脱成分の抽出−動特性モデルに基づく楽譜との時間的対応付け−∗◎大石康智 (NTT), 亀岡弘和 (NTT/東大), 持橋大地 (統数研), 永野秀尚, 柏野邦夫 (NTT)1

歌唱者Aの歌声
F0系列
提案手法
楽譜

はじめに楽曲のメロディを歌った歌声の音高（F0）軌跡には，そのメロディを構成する楽譜の音符の並びだけでなく，楽譜に記載されない，“楽譜から逸脱した” 動的変動成分が含まれる [1]。これらは，発声器官の物理的制約に起因する成分（特に，オーバーシュートやプレパレーション，微細変動成分）と意図的表現による成分（特に，ビブラートやポルタメント）からなると考えられ，知覚的には，前者は人間らしさ・自然性に関係し，後者は巧拙感に関係することがわかってきている [2, 3]。さらに，後者は意図して意図通りにずれた場合と，意図通りに制御できなかったずれとに分かれ，習熟度に関連すると考えられる。ただし，これらの逸脱成分を F0 系列から特徴抽出することや歌唱者ごとにその多様性を学習することまではまだ十分に検討されていない。
本研究の目的は，このような物理的制約もしくは意図的表現による楽譜逸脱成分を F0 系列から楽譜と分離して抽出し，歌唱者ごとにどのようなパターンをもちうるのか，各パターンが文脈（楽譜の音符列）にどう依存するかを計算機に学習させることである(Fig. 1)。ここでは，うろ覚えの状態で歌った歌声ではなく，楽曲のメロディを既に知っているものとして，歌唱者なりに表情付けして歌った歌声を対象とする。
本研究は，歌唱者の歌い方や個性，癖を学習することにつながり，現在盛んに研究される歌声合成や歌声変換 [4, 5, 6] などへの応用が期待できる。例えば，ある歌声を別の歌唱者の歌い方に変換して合成することが可能となるだろう。文献 [5] と違って，事前に歌唱者の歌い方が学習されるため，どんなメロディにもその歌い方を転写できることを特長とする。
本稿では，F0 動特性モデル [7, 8] に基づいて，楽譜逸脱成分を特徴抽出するための，F0 系列と楽譜の音符列との時間的対応付け手法を提案する。まず歌声の F0 をy(t) = yt (t) + ys (t) = u(t) ∗ h(t) + ys (t)(1)と，連続時間領域 t で表現する。u(t) は楽譜の音符の並びが歌唱によって時間周波数軸上に伸縮した階段状の信号である。これに，減衰率 ζ と固有周波数 Ωをパラメータとする 2 次系伝達関数

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

∗

逸脱成分

オーバーシュート
ポルタメント

逸脱成分

ビブラート
プレパレーション
微細変動成分逸脱成分（物理的制約）
パターン
逸脱成分（意図表現）パターン
11パターン
パターン

Fig. 1 F0 系列と楽譜の音符列との時間的対応付けに基づく楽譜逸脱成分の抽出と歌い方の学習
別の音符に遷移するときの楽譜逸脱成分を含む信号
yt (t) を得る。ここで，u(t) を楽譜の音符の並びが既知の下で構成される隠れマルコフモデル（HMM）からの出力として表現する。一方，ys (t) はビブラートやプレパレーション，微細変動成分のような音符の音高が安定するときの楽譜逸脱成分を表す。このとき，y(t) である F0 と，そこに含まれる音符の並びが与えられた下で，これらを時間的に対応付けながら，音符の遷移における逸脱成分の因子 h(t) と音高が安定するときの逸脱成分 ys (t) を分離して特徴抽出することを考える。このモデルの表現方法とパラメータ推定アルゴリズムを提案し，評価実験結果を示す。

2

F0 系列と楽譜との時間的対応付け

動特性モデル [8] では，実測の F0 を手動で音符区間に区切り，その区間の動特性を 2 次系モデルで表現した。本手法は，音符の並びそのものを HMM で表現し，F0 の動特性を考慮しながら音符区間（状態系列）をも反復的に推定できることを特徴とする。

2.1

歌声の F0 の離散時間表現

実測の F0 は離散時間信号であるので，式 (2) を後退差分変換（s ' (1 − z −1 )/t0 で近似する変換であり，t0 は離散時間表現におけるサンプリング周期）によって，離散時間領域で表現する。この変換により，H(s) の逆システムの伝達関数を z 領域でH−1 (z) = a2 z −2 + a1 z −1 + a0(3)と書ける。ただし，a2 = ϕ2 , a1 = −2ϕ(ψ + ϕ), a0 =1 + 2ϕψ + ϕ2 であり，ϕ と ψ は ϕ = 1/(Ωt0 ), ψ = ζと表現される。ここで，k を離散時刻インデックスとすると，階段状信号 u[k] はパラメータ ϕ と ψ によって特性が決まる拘束つき全極モデルからの出力

u[k] = a0 yt [k] + a1 yt [k − 1] + a2 yt [k − 2](2)のインパルス応答 h(t) を畳み込むことによって [1]，オーバーシュートやポルタメントのような，音符が歌唱者Aの歌い方(4)とみなすことができる。一方，ys (t) も離散時間表現
すると，歌声の F0 系列はこれらの成分の和 y[k] =Extraction of Note Deviation Manner in F0 contour of Singing Voice. by OHISHI, Yasunori (NTT),KAMEOKA, Hirokazu (NTT / University of Tokyo), MOCHIHASHI, Daichi (The Institute of StatisticalMathematics), NAGANO Hidehisa, KASHINO Kunio (NTT)

日本音響学会講演論文集- 279 -2011年9月

1st 音符

2nd 音符I th音符3Start・・・

21EndJ

Fig. 2 楽譜の音符の並びが既知の下で構成される隠れマルコフモデル
yt [k] + ys [k] と表現される。
2.2

歌声の F0 動特性モデルの統計モデル化

Left-to-Right 型 HMM を利用して，階段状に変化する信号 u[k] を表現する（Fig. 2）。u[k] を HMMの出力値系列として，u[k] ∼ N (µ[k] + d, σu2 ) のように正規分布する確率変数と見なし，平均 µ[k] + dが HMM の状態遷移に伴って階段状に変化するモデルを考える。Fig. 2 の点線枠で囲った部分は，メロディの楽譜の i 番目の音符を形成する状態集合{Si,1 , Si,2 , . . . , Si,J−1 , Si,J } であり，これらの状態の出力分布はすべて同じものとする。状態 Si,j (j =1, . . . , J − 1) は，Si,j+1 だけに遷移し (遷移確率 = 1)，
状態 Si,J は自己遷移，もしくは次の音符の初期状態Si+1,1 に遷移する。このような HMM を用いる理由は，短い音符区間が形成されることを抑制するためである。I は楽譜に含まれる音符の総数とする。HMMの構成は以下のとおりである。
出力値系列 : {u[k]}K
k=1
状態集合 : S := {S1,1 , S1,2 , . . . , S1,J , . . . , SI,J }状態系列 : {sk }K
k=1
状態出力分布 : P (u[k]|sk = Si,j ) = N (cSi,j [k] + d, σu2 )cSi,j [k] = A(i)
状態遷移確率 : φSi,j ,Si0 ,j 0:= log P (sk = Si0 ,j 0 |sk−1 = Si,j )ここで，状態出力分布の平均値の A(i) は楽譜から与えられる i 番目の音符の，対数周波数軸上の平均律の音高値とする。一方，d はその音高からのずれを表現し，キーシフトパラメータと呼ぶ。例えば，歌唱者が楽譜よりも少し高めに歌ったり，低めに歌ったりすることに対応でき，このパラメータも楽譜逸脱成分を表す。これは全体にわたる定数パラメータ dとしても，音符ごと（Fig. 2 の点線枠）に変化する
パラメータ {d1 , d2 , . . . , dI } としてもよい。以降は簡単のため，パラメータ d としてモデルを組み立てる。
φSi,j ,Si0 ,j0 を定数とすると，推定すべきパラメータは
2K
θu := {{sk }K
k=1 , d, σu } となる。また，{µ[k]}k=1 は，K{sk }k=1 が与えられたもとで µ[k] = csk [k] となる。
u[k] の統計モデルに基づき，F0 系列の確率密度関数を導く。まず，y := (y[1], . . . , y[K])T , u :=(u[1], . . . , u[K])T , µ := (µ[1], . . . , µ[K])T , d := d1Kとおくと，u|θu ∼ N (µ + d, Σu ),日本音響学会講演論文集

Σu = σu2 I K(5)と書ける。ここで，1K は K 個の 1 からなるベクト
ル，I K は K × K の単位行列である。前節で述べたように，音符遷移における楽譜逸脱成分を含む信号
y t := (yt [1], . . . , yt [K])T と u の関係は，


a0O
a
 1 a0

aaa

10
Ψ= 2(6)


..
..
..
.
.
.

Oa2 a1 a0と置くと，u = Ψy t のように表現でき，式 (5) より，¡¢y t |θu , ϕ, ψ ∼ N Ψ−1 (µ + d), σu2 Ψ−1 (Ψ−1 )T (7)が導かれる。ここで，音符ごとに 2 次系の制御パラメータ ϕ, ψ が異なるとしてもよく，その場合は，Ψ =TT T(ΨT1 , Ψ2 , . . . , ΨI ) とし，各 Ψi は ϕi , ψi で表現される。以降は簡単のため，ϕ, ψ は全体で共通とする。
楽譜逸脱成分 y s := (ys [1], . . . , ys [K])T は Gauss 性¢¡
白色雑音を仮定し，y s ∼ N 0, σs2 I K と表現する。
これに含まれるビブラート，プレパレーション，微細変動成分の特徴抽出については今後の課題とする。
y t , y s は互いに独立と仮定し，Θ := {θu , ϕ, ψ, σs2 } の下での F0 系列 y = y t + y s の確率密度関数は，¡¢y|Θ ∼ N Ψ−1 (µ + d), σu2 Ψ−1 (Ψ−1 )T + σs2 I Kで与えられる。これが，パラメータ Θ に関する F0系列 y の尤度関数となる。Θ の事前確率は，各要素は独立で，状態系列 {s[k]}K
k=1 , ϕ, ψ, d 以外のパラメータは一様に分布すると仮定し，P (Θ) ∝QKP (s1 ) k=2 P (sk |sk−1 )P (ϕ)P (ψ)P (d) とする。
先の尤度関数は全区間で F0 が観測されることが暗に想定されているが，実際には F0 は有声音が発せられるときのみ観測可能であり，無声音の区間では F0は通常観測できない。従って，実測 F0 系列からのパラメータ推定を行うためには一般に欠損データの問題を扱う必要がある。統計的枠組では，この問題は不完全データ問題に他ならず，EM アルゴリズムにより扱うことができる。詳細は文献 [9] を参照されたい。

2.3

パラメータ推定アルゴリズム

実 測 F0 系 列 y が 与 え ら れ た と き ，事 後 確 率
P (Θ|y) ∝ P (y|Θ)P (Θ) を最大化するパラメータ Θ を推定する問題を考える。P (Θ|y) を最大化する問題は
T T
解析的に解くことはできないが，x := [y Tt , y s ] を完全データと見なすことで EM アルゴリズムによる不完全データ問題に帰着できる。不完全データと完全データの関係は，y = Hx となり，ここで，H := [I K I K ]とする。このとき，Q 関数は，
h¡¢c 1log |Λ−1 | − tr Λ−1 E[xxT |y; Θ0 ]Q(Θ, Θ0 ) =2i+ 2mT Λ−1 E[x|y; Θ0 ] − mT Λ−1 m + log P (Θ)"#"#Ψ−1 (µ + d)
ΨT Ψ/σu2O−1m :=, Λ :=0OI K /σs2- 280 -2011年9月

数波
周数対となる。ここで，tr(·) は行列のトレースを表し，
E[x|y; Θ0 ] と E[xxT |y; Θ0 ] は，条件付きガウス分布の性質より，

E[x|y; Θ0 ] = m + ΛH T (HΛH T )−1 (y − Hm)

音符列
時間

数波周数対

合成F0系列時間

Fig. 3 文献 [1] に基づく合成 F0 系列の作成

E[xxT |y; Θ0 ] = Λ − ΛH T (HΛH T )−1 HΛ+ E[x|y; Θ0 ]E[x|y; Θ0 ]T (8)と書ける。E ステップでは，直前のステップで更新されたパラメータを Θ0 に代入し，E[x|y; Θ0 ] とE[xxT |y; Θ0 ] を算出する。M ステップでは，各パラメータの更新を行う。便宜上，y t , y s に対応するように E[x|y; Θ0 ] および E[xxT |y; Θ0 ] を
" #"
#x̄t
Rt ∗0
T0
E[x|y; Θ ] =, E[xx |y; Θ ] =x̄s
∗ Rsのように区分表現する。各パラメータの更新式を以下に示す。
1) 状態系列: Q 関数で s := {sk }K
k=1 に関する項はTable 1

合成 F0 系列を利用した評価実験結果
音符区間の正解率 [%]

提案法
比較法

99.094.8

二乗平均平方根 (RMS)ζΩd0.0398
0.09720.00209
0.01612.2820.66

こで，M は各行に 1 が一個あり残りはすべて 0 であるような K × I のバイナリ行列であり，階段状に変化する µ の各要素に基づいて構成される。各 i に関して事前分布 di ∼ N (0, σd2i ) を導入し，Q 関数を diに関して最大化すると更新式が求められる。
4) 分散パラメータ: σu2 , σs2 の更新式は，

K1 X(o[k] − csk [k])2 + log P (s1 )I1 (s) := − 22σu¡¢σu2 = tr(ΨT ΨRt ) − 2r T Ψx̄t + r T r) /K(12)σs2 = tr(Rs )/K(13)k=1+KXlog P (sk |sk−1 )(9)k=23である。ただし，o[k] := [Ψx̄t ]k であり，[·]k はベクトルの k 番目の要素を表す。これを最大化する {sk }K
k=1は Viterbi 探索により効率的に解くことができる。
2) 2 次系制御パラメータ: ϕ の事前分布を ϕ ∼N (µϕ , σϕ2 )，ψ の事前分布を ψ ∼ N (µψ , σψ2 ) とする。
Q 関数より，ϕ, ψ に関係する項は11tr(ΨT ΨRt ) + 2 (µ + d)T Ψx̄t
2σu2σu11+ log |Ψ| − 2 (ϕ − µϕ )2 − 2 (ψ − µψ )2 (10)2σϕ2σψ

I2 (ϕ, ψ) := −と な る 。こ こ で ，式 (3) よ り，Ψ は 定 数 行 列U 0 , U 1 , U 2 を用いて Ψ = U 2 ϕ2 + U 1 ϕψ + U 0 と表せるので，式 (10) の ϕ と ψ に関する偏導関数を求めて，ニュートン・ラフソン法を利用すれば，数値的に ϕ と ψ が求められる。
TT T
Ψ = (ΨT1 , Ψ2 , . . . , ΨI ) と表現した場合も同様
2に，各 i に関して事前分布 ϕi ∼ N (µϕi , σϕ), ψi ∼i2N (µψi , σψi ) を導入する。そして，Q 関数の ϕi と ψiに関する偏導関数から，ニュートン・ラフソン法により数値的に ϕi と ψi が求められる。
3) キーシフトパラメータ: d の事前分布を d ∼N (0, σd2 ) として，Q 関数を d に関して最大化すると

d = (1T Ψx̄t − 1T µ)/(K + σu2 /σd2 )(11)の更新式を得る。
音符ごとに変化するキーシフトパラメータを導入した場合，d = M (d1 , d2 , . . . , dI )T と表現できる。こ

日本音響学会講演論文集となる。ここで，r := µ + d とした。

評価実験

提案法による F0 系列と楽譜の音符列との時間的対応付け，および楽譜逸脱成分を表現するモデルパラメータの推定性能を評価する。実測の F0 系列を用いた場合，これらの正解となる値はわからないため，まずは人工的に作成した合成 F0 系列で評価する。事前に童謡や歌謡曲から抜粋したメロディの楽譜（合計 16 種類，抜粋部分の長さは平均 11.4 秒）に基づいて，対数周波数 (単位は cent，半音は 100 cent に相当する) 軸上で階段状に変化する音符列を作成する（サンプリング周期 5 ms，休符区間はその前の音符の音高を休符区間も伸ばす）。これに，文献 [1] で調査されたプロ歌唱者とアマチュア歌唱者の F0 系列を合成するためのパラメータを用いて，楽譜逸脱成分を付加する (Fig. 3)。具体的には，オーバーシュート (プロ：ζ = 0.570, Ω = 0.0363 [rad/ms]，アマ：ζ = 0.518, Ω = 0.0312 [rad/ms])，プレパレーショ
ン (プロ：ζ = 0.675, Ω = 0.0363 [rad/ms]，アマ：ζ = 0.660, Ω = 0.0287 [rad/ms])，ビブラート (プロ：Vibrato extent = 83.0 cent, Vibrato rate = 6.25 Hz、
アマ：Vibrato extent = 33.5 cent, Vibrato rate =4.70 Hz)，微細変動成分 (白色雑音をカットオフ周波数 10 Hz のローパスフィルタに通し，振幅の平均値を 20 cent に調整)，キーシフト (+25 cent，0 cent，−25 cent の 3 パターン) とした。
作成された全 96 パターンの合成 F0 系列にパラメータ推定アルゴリズムを適用する。固定パラメー
2
タは，文献 [1] を参考に，J = 20, µϕ = 5.74, σϕ=2
20.1, µψ = 0.57, σψ = 0.02, σd = 400, φSi,J ,Si,J =- 281 -

2011年9月

C#B
A
G
F
D#C#
C#
B
A
G
F
D#
C#C#
B
A
G
F
D#C#F0系列

歌唱者A

0.036

歌唱者A

歌唱者B

歌唱者C

0.035

）名音（数波
周数対

(1)

(2)

(3)

0.034
0.55

歌唱者B
(1)

(2)

(3)

1

2345(2)78(3)0.6

0.65 0.550.60.65(1)0(2)

0-200200(3)0.99999, φSi+1,1 ,Si,J = 0.00001 とした。評価尺度は，(1) 正解の音符区間と，推定された状態系列に基づく音符区間の一致する割合（正解率），(2) オーバーシュー
ト作成の ζ, Ω と，推定されたパラメータ ψ, ϕ から計算される ζ, Ω との二乗平均平方根 (RMS)，(3) キーシフトの値と，推定された d との RMS，とする。
Tab. 1 に評価結果を示す。σs2 が 10−5 のオーダで変化しなくなるまで反復推定を繰り返した。比較法とは，2.2 節の HMM を用いた Viterbi 探索により，F0系列を前処理で音符区間に分割し，文献 [8] の動特性モデルによってその区間ごとにパラメータを推定して評価した結果である。提案法は F0 の音符遷移における動特性を考慮しながら Viterbi 探索を反復して行うため，前処理のみの比較法に比べて，音符区間（状態系列）の推定性能が高くなる。それによって，ζ, Ω, dの RMS が小さくなり，比較法に比べて，楽譜逸脱成分を抽出できていると言える。提案法の ζ, Ω の RMSが許容できる誤差であるか確認するため，音符列に，正解の ζ, Ω を代入した 2 次系インパルス応答を畳み込んだ信号と，推定された ζ, Ω を代入した 2 次系インパルス応答を畳み込んだ信号との RMS を計算したところ，平均 6.57 cent であった。これは半音 100 centの 6.57%の誤差であるので，提案法は楽譜逸脱成分を抽出するための妥当な方法であると考えている。
次に実測の F0 系列から楽譜逸脱成分の特徴抽出を試みる。音大出身の歌唱者 3 名（男性 1 名，女性 2名）がレミオロメンの「粉雪」の出だしの部分をヘッドフォンで伴奏を聞きながら歌詞付きで歌った歌声を収録した。F0 は YIN[10] を利用して，5 ms ごとに推定した (t0 = 5 ms)。ここでは，音符ごとに 2 次系パラメータおよびキーシフトパラメータが変化するものとして，アルゴリズムを実行した。従って，事前2分布のパラメータを µϕi = 5.74, σϕ= 0.1, µψi =i
220.57, σψi = 0.02, σdi = 400 (i = 1, . . . , I) とした。
Fig. 4 の点線は楽譜の音符列を時間軸上で伸縮させた結果であり，さらに各音符にキーシフトパラメータ di を周波数軸上に足し合わせた結果が赤実線（式(5) の平均値）となる。定性的にみると階段状信号がうまく推定され，どの歌唱者も楽譜の音符の音高よ0-200[s]Fig. 4 実測の F0 系列と楽譜の音符列との時間的対応付け結果

日本音響学会講演論文集

0.65 0.55-200200

歌唱者C
(1)
60.6[cent]2000.5

1.0 [s]0.5

1.0 [s]0.5

1.0 [s]Fig. 5 各歌唱者の F0 系列における楽譜逸脱成分の抽出結果：図中の (1)〜(3) は Fig. 4 の音符区間 (1)〜(3) の逸脱成分に対応する。
りも半音程度高めに歌う楽譜逸脱成分を観測できる。
Fig. 5 の上部は，各音符の遷移における楽譜逸脱成分を構成する 2 次系パラメータ ζ, Ω の分布を示す。
これは歌唱者ごとに大きな違いはなく，先行研究 [1]の結果と整合する。一方，Fig. 5 の下部には Fig. 4の (1)〜(3) の音符区間における音高が安定するときの楽譜逸脱成分を示す。ここでは，各音符の x̄s を図示した。同じ音符なのに，歌唱者ごとにその挙動が異なり，例えば，歌唱者 A はオンセットから周期的に振動するが，歌唱者 C はオンセットからしばらくして振動する傾向にある。この楽譜逸脱成分を周期的なビブラート，プレパレーション，微細変動成分に分離すること，そしてそれらが歌唱者ごとにどのようなパターンをもつのか，音符の並びにどう影響を受けるか学習することを今後の課題とする。

4

おわりに歌唱者の歌い方や個性，癖の学習を目指して，F0系列に含まれる楽譜逸脱成分を抽出するための，F0系列と楽譜の音符列との時間的対応付け手法を提案し，その評価実験を行った。今後は，抽出された楽譜逸脱成分のパターンを学習すること，音高とともに音色，音量に関する楽譜逸脱成分の特徴抽出に取り組むことも検討中である。

参考文献[1] Saitou, T. et al., Proc. EUROSPEECH 2009, pp.
832–835, 2009.
[2] 齋藤ほか, 音講論集, 2-Q-26, 2007.
[3] 右田ほか, 音講論集, 3-P-28, 2010.
[4] Bonada, J. et al., Proc. SMAC 2003, 2003.
[5] Nakano, T. et al., Proc. SMC 2009, pp. 343–348,2009.
[6] 大浦ほか, 情処研報音楽情報科学, 2010–MUS–86–1,2010.
[7] 柏野ほか, 音講論集, 2-9-1，pp. 625–626 1998.
[8] Ohishi, Y. et al., Proc. INTERSPEECH 2010, pp.
2598–2601, 2010.
[9] 亀岡ほか, 音講論集, 1-1-3, 2010.
[10] de Cheveigné, A. et al., JASA, Vol. 111, No. 4, pp.
1917–1930, 2002.

- 282 -2011年9月