3-1-15
拍長の連続性を考慮した潜在的調波配分法に基づくスコアアライメント手法∗☆前澤陽 (京大), 後藤真孝 (産総研), 尾形哲也, 奥乃博 (京大)
近年、計算機を用いて楽譜表現を援用した音楽音響信号の新たな楽しみ方が提唱されている。例えば、過去のヴァイオリン名演奏の指使いを推定する「演奏法の耳コピ」[1] や、自分の嗜好に合致した演奏者の検索 [2] や、特定の楽器を増幅し [3, 4]、市販の CDから、カラオケ音源を作成するといったことが可能となってきている。これらのアプリケーションでは、音楽音響信号の分析のために楽譜情報を用いる。楽譜というシンボリックな情報と音響信号という波形情報の橋渡しするためには、音響信号の位置と楽譜の位置の時間的対応付け (スコアアライメント、以下アライメント) を求めることが必須である。
アライメントに必要な要件は、音色や音量の変化に対するロバストネスと、音符の時系列の適切なモデル化である。特にクラシック音楽では、繰り返しの省略を検出する機構が必要である。というのは、クラシック音楽の楽譜に記載されている繰り返し指示は、しばし演奏者の解釈により、無視されることがあるためである。
従来、音量と音色のロバストネスを実現するため、アドホックな特徴設計 [5–7] や楽器音データベースを用いた音色の学習 [8] を行っていた。しかし、前者には、緻密なパラメータチューニングが必要であり、設計者のチューニングや音源の選定に性能が依存する問題がある。後者には、アライメントの品質が、楽器音データベースの良し悪しに関連する問題がある。
また、楽譜の時系列モデル化には、隠れマルコフモデ
ルや、線形動的システム (LDS) がある。前者は、繰り返し構造といった、楽譜上の状態遷移をうまく記述できる。しかし、モデルに暗黙に仮定される音長の独立性は、音楽的に妥当ではなく、これに起因する精度低下が問題となる。後者は、拍の連続性を考慮しているので、このような問題は起こりづらい。しかし、繰り返し構造のような離れた楽譜位置への遷移が扱えないという問題がある。
本稿では、音源の選定が不要であり、かつ音量と音色のロバストネスを実現し、繰り返し構造などを許容し、かつ拍長の連続性を保つアライメント手法を提案する。楽譜時系列のモデル化には、隠れセミマルコフモデル (HSMM) を、LDS による拍長モデルに条件づける。これにより、連続的なテンポと複雑な楽譜構造に対する許容を同時に実現する。音のモデルには、楽器音の混合音スペクトルをベイズ的に扱う音源モデル LHA[9] を用いる。音色と音量に無情報事前分布を置くことにより、これらに対するロバストネスを実現する。また、音色が無情報であるため、音源の選定が不要である。

1

成されると仮定する。 ただし、LHA の定式化と違い、調波構造は楽器音高ペア内で共有されているとし、また音量バランスは音符内で一貫していると仮定する。さらに、ある楽器の状態内に置ける周波数ビンは単一の楽器の、単一の倍音から生成されるとする。
(i)
Zi (f, d) を、状態 d において楽器音高ペア i が周波数 f が占拠している場合 1 でそれ以外は 0 の二値行(h)列とし、Zj (f, i) を、周波数 f が、楽器音高ペア iの第 j 倍音から生成される場合 1 の二値行列とする。
(s)Zl,d (t) は時刻 t が、状態 d で次の状態に遷移するまでのフレーム数が l のとき 1 の値をとる二値行列とする。i 番目の楽器音高ペアの基本周波数が μi であり、−1/2で隣接する周波窓関数の影響などにより分散 λi数でパワーが観測されるとする。以上より、観測信号の尤度は次のように表すことができる:
p(X|Z (i,h,s) , μ, λ) = 
Z (s) (t)X(f,t)Zi(i) (f,d)Zj(h) (f,i)
N log f /j|μi , λi −1 l,d調波構造と音量バランスは多項分布に従うと仮定する。
p(Z (i) |E, Z (s) ) =

(s)
(i)
ei (d)Zl,d (t)X(f,t)Zi (f,d)(2)

t,i,f,d,lp(Z (h) |A, Z (i,s) ) =

(s)

(i)

(h)

aj (i)Zl,d (t)X(f,t)Zi (f,d)Zj(f,i)

t,i,j,f,d,l(3)

e と a をそれぞれ音符生起確率と倍音生起確率と呼ぶ。これらは、音符の相対音量と倍音ピークの相対強度にそれぞれ対応すると考えることができる。これらを更に確率変数としてして扱い、事前分布を無情報にすることで、音色と音量の変化に対するロバストネ
スを実現できると考えられる。そこで、音符生起確率と倍音生起確率の事前分布としてディリクレ分布をおき、基本周波数の事前分布として Normal-Gamma
分布を置く:
p(μ, λ|ν, b, m, l) =p(E|E0 ) =p(A|A0 ) =(H) (H) (H) (H)N G(μi , λi |mi , bi , li , νii
D

Dir(e(d)|e0 (d))(5)

Dir(a(i)|a0 (i))(6)dIi) (4)

楽譜時系列 Z (s) の分布として HSMM を仮定する。
初期状態の確率分布を π とする。

モデルの定式化

本手法は、入力信号の定 Q 変換に対し、入力された楽譜表現とのアライメントを行う。以後、楽譜において、特定の楽器が奏でている特有の音高の対を楽器音高ペアと呼ぶ。すなわち、楽譜の特定の位置は複数の楽器音高ペアの集合であり、楽譜とはこれらを連結したものである。
音量と音色のロバストネスを実現するために、ス
ペクトルを潜在的調波配分法 (LHA) を用いてモデル化する。LHA の出力は、現在の楽譜位置に依存する。
各時間フレームにおけるスペクトルは LHA に従い生∗(1)t,i,j,f,d,l(s)p(Z (s) |T, π, τ ) = π Z (1)Z (s)  (t−1)Z (s) (t) l,d1,dlτd (d)N log|Td , σT2Ld
t=2,l,d,d =dp(π|π0 ) = Dir(π|π0 )p(τ |τ0 ) =Dir(τ (d)|τ0 (d))

(7)
(8)
(9)

d

Audio-to-Score alignment based on Latent Harmonic Allocation with smoothness of beat length. by Akira
MAEZAWA (Kyoto U.), Masataka GOTO (AIST), Tetsuya OGATA (Kyoto U.), Hiroshi G. OKUNO(Kyoto U.)

日本音響学会講演論文集- 1071 -2011年3月

式 (7) は、楽譜時系列を、拍長と楽譜上の状態遷移の組み合わせとして表すことを意味する。τ は、複雑な楽譜構造を HMM のように記述できる。Td は楽譜位置 d における対数拍長である。Td の連続性を保たせると、音楽的に妥当な拍長のモデル化が可能となる。そこで、Td を平滑化させるために、LDS をおく:p(T ) =

N

−1
Td |Td−1 , Ld−1 λ(T ) dp(λ)=

G(λ(T )(T )(T )d |ld , νd )(11)

本手法では、これらの事後分布を推定し、状態系列 Z (s) を音価 l に対して積分消去したものの事 (s)後確率を最大化させる状態系列 arg max l Zl,d (t)をスコアアライメントとする。しかし、事後分布の推定は困難であるため、変分近似に基づく EMアルゴリズム (VBEM) を用いて事後分布を推定する。VBEM では、事後分布 q(LDS, LHA, HSM M )が qLDS (LDS)qLHA (LHA)qHSMM (HSM M) と因子分解できると仮定する。このような分布を変分事 後 分 布 と 呼 ぶ 。こ こ で 、qHSMM (HSM M ) =qZ (s) (Z (s) )qπ (π)qτ (τ )と 因 子 分 解 で き 、qLHA (LHA) = qZ (h) (Z (h) )qZ (i) (Z (i) )qμ,λ (μ, λ)と因子分解でき、qLDS (LDS) = qT (T )qλ(T ) (λ(T ) )と因子分解できるとする。変分事後分布の推定は、同時分布との KL ダイバージェンスの最小化問題として定式化できる。すると、任意の因子 Z は、以下のように更新できる。
qZ (Z) ∝ exp log p(X, LDS, LHA, HSM M )¬Z(12)ただし、f (x)x は x の下での f (x) の期待値であり、¬y とは、y 以外のすべての確率変数のことを指す。推定は、KL ダイバージェンスが収束するまで、各確率変数の変分事後分布を交互に更新する。

モデルの推論lNf (i, j) = log N log f /j|μi , λi −1(13)
11 l¯i− log 2πν̄i + ψ(l̄i )(log f /j − m̄i )2+=−2 ν̄ib̄iel(d, l) = log p(log l|Td )T(14)(s)

Zd,l (t)llAj (i) = log aj (i)a(i)= ψ(ᾱj (i)) − ψ(15)Z (s)lE i (d) = log ei (d)e(d)= ψ(¯i (d)) − ψeτ d (d) = log τd (d)τ (d)= ψ(τ̄d (d)) − ψM
l=1

K
l=1

D
l=1


ᾱl (i)(i)qZ (h) (Z (h) ) =



(h)

j,i,f

ξj (f, i)Zj(i,f )(21)

φ (f,i)ただし ξj (f, i) =  jφk (f,i) であり、φ は次のようにk表される:
X(f, t)ηd (t) γi (f, d) ×t,dlNf (i, j) + lAj (i)log φj (f, i) =(22)

式 (22) も式 (20) と似たように、楽器音高ペア i 内の平均スペクトルを倍音毎に分配するものとみなせる。

2.1.1 LHA における変分 M ステップ

楽器音高ペアの独立性により qE = i qei と表せられる。 また、多項分布とディリクレ分布の共役性により、 ei の事後分布は次のように求められる:
qei ∼ Dir (ei | ¯i )(23)
ここで、¯i (d) = e0i (d) + t,f ηd (t) X(f, t)γi (f, d) と与えられる。同じく、倍音生起確率も qA = i qai と表す
qai ∼ Dir (ai |ᾱ)ᾱi = a0 (i) +¯l (d)γi (f, d)Zi (f,d)

X(f, t)γi (f, d)ηd (t) ξ(f, i)である。基本周波数とその分散の事後分布は、基本
周波数の独立性から qμ,λ = i qμi ,λi であり、また正規分布と N G 分布の共役性により、qμi ,λi は次のパ(H) (H) (H) (H)
ラメータで与えられる N G(m̄i , b̄i , l̄i , ν̄i ) である:
(H)

m̄i(H)b̄i(17):=

mi bi + Nγ,i log(f /j)ψi:= bi + Nγ,i1 bi Nγ,i(H)ν̄i := νi +2 bi + Nγ,i(18)(25)

bi + Nγ,i(H)l̄i := li + Nγ,iτ̄l (d)(24)

t,f,d

2.1 LHA の変分 E ステップ
HSMM の各状態 d における、楽器音高ペア i が周波数 f に占める割合 Z (i) を次のように更新する:
(20)

式 (20) 右辺第 1 項を状態 d で重み付けたスペクトルの周波数平均、また第 2 項を音符 i における音量の対数期待値と、音符 i の調波構造と倍音ピークの対数期待値による重み付けと考えると、ρi (f, d) は状態 d 内の平均スペクトルを、音符ごとに分配するとみなすことができる。
同じように、各楽器音高ペア i における、倍音 j が周波数 f に占める割合 Z (h) を次のように更新する:(16)

ここで、ψ(x) はディガンマ関数である。また、f (x)xは、確率変数 x の下での関数 f (x) の期待値である。

qZ (i) (Z (i) ) =j

ことができ、ai の事後分布は次のように求められる:

簡単のため、次の変数を定義する:

ηd (t) =X(f, t)ηd (t) ×
ξj (f, i) lNf (i, j) + lAj (i)

t(10)

d

2

log ρi (f, d) =lE i (d) +d(T )

ただし γi (f, d) = ρiρ(f,d)
であり、 ρ は次のようにi i (f,d)
表されるとする:

2+ Nγ,i

(26)
(27)
2



log(f /j)ψ(i) − mi
2log(f /j) − log(f /j)ψ(i)
(28)
ψ(i)

ψi (f, j) は次のような多項分布である:(19)

ψf,j (i) =γi (f, d)ξj (f, i)ηd (t) X(f, t)/Nγ,i(29)γi (f, d)ξj (f, i)ηd (t) X(f, t)

(30)

d,t

Nγ,i =d,t,f,ji,f,d

日本音響学会講演論文集- 1072 -2011年3月

2.1.2 LDS の変分 E ステップ
LDS の更新には、カルマンスムーザーと同様、時系列の事後分布を前向き後ろ向きアルゴリズムを用い求める。
qT (T) は ψl (d) = t=2,d =d ζd (t, l, d) とすると、次のように、カルマンスムーザーと似た形で与えられる:

これらを用いて、拍長の事後分布を次のように得る：lψl (d) log N log|Td , σT2
Ld2
σTdl

−1+ log N Td |Td−1 , Ld−1 λ(T ) d(31)(T )2.2 HSMM の変分 EM ステップ
状態継続長の期待値を次のように得る:122el(d, l) = − 2 (log l/Ld − Td ) − log 2πσT2σTTd


21vd1
ud= − 2 log l/Ld −+−1−12σT
qsd
qd + sdd
1− log 2πσT2(39)− 22σT (qd −1 + sd −1 )log qT (T) =λd

通常の LDS では、ベクトルを出力するのに対し、本手法では l に対するヒストグラムを出力する点が異なる。前向きアルゴリズムは次のように表される:(L)

αd(Td ) = p(Td |ψ(1 : d))−1(L)∝ αd−1 (Td−1 ) p(Td |Td−1 , Ld−1 λ(L) d )p(logqτ =Dir (τd (d)|τ̄d (d))

ただし、τ̄d (d) = τ0d (d) + t,l ζd (t, l, d)qZ (s) (Z (s) ) は次のように求まる:とすると次のように求まる:1
ψl (d) 2 +l
σT

s−1d =⎛

Ld−1⎜− md ⎝λ(L) d
Ld−1

ud−1+sd−1

Ld−1l⎞2⎟⎠l ⎟ψl (d)log⎠
σT2Ld(34)

同様に後ろ向き変数を次のように求める:(L)

βd(Td ) = p(ψ(d + 1 : T )|Td )∝ p(ψ(d + 2 : T )|Td+1 )p(Td+1 |Td )l, Td+1 )ψl (i+1) dTd+1Ldlψl (d+1)
 l= βd+1 (Td+1 )N log|Td+1 , σT2×Ld+1l
−1N Td+1 |Td , Ld λ(L) d+1 dTd+1 = N (Td |vd , qd ) (35)×t=1log κd (f ) = γi (f, d)ξj (f, i) lNf (i, j) + lAj (i) + lEi (d)×κd (f ) は状態 d が出力する、正規化されていないスペクトルの期待値と解釈できる。 LHA の不確定さが高い状態 d の κd (f ) の周波数軸の累計は小さな値をとるため、 状態系列の期待値に、LHA がどれだけ信号を説明できるかの良し悪しが影響する。
また、これは通常の HSMM と同じ形をしているた
め、期待値の計算において前向き後ろ向きアルゴリズ
ムを使用できる。α(H) を HSMM の前向き変数、β (H)を後ろ向き変数とすると、次の漸化式が求まる:(H)

qd−1 =Ld⎜− nd ⎝

λ(T ) d+1
vd = nd qd

Ld

日本音響学会講演論文集

l

Ld

⎞2
⎟
⎠

(36)

l
vd+1
ψl (d + 1)
log+σT2
Ld+1
qd+1(37)(s)

αl,d (t) = p(Z(l,d) (t) = 1|X(1) · · · X(t))(H)∝X(f,t)αl ,d (t) exp log τl ,d (l, d)τκd (f )f (H)X(f,t) =κd (f )
× αl−1,d (t − 1)f(H)exp (eτ d (d ) + el(d, l)) α1,d (t − 1)(43)l ,d+(H)d(s)βl,d (t) = p(Xt+1 (f ) · · · XT (f )|Z(l,d) (t) = 1)

T

λ(T ) d+1(42)i,j同じく平方完成を行うと次を得る。ただし
nd =
−1
λ(L) d+1   ψl (d+1)
1++ l σ2とする:
qd+1
Ld⎛

(s)

Zl,d (t)X(f, t) log κd (f ) (41)ただしp(logλ(T ) d+1(s)Zl,d (1) log ππ+(33)⎞λ(L) d⎜ud = sd ⎝md⎛(40)(s)(s)+Z1,d (t − 1)Zl,d (t)(eτ d (d) + el(d, l))t=2,d =d−1

λ(L) d(38)l,d非積分項の指数の中において Td−1 を積分消去
Td に 対 し 平方 完 成 す る と 、 ud , sd は, md =し(L)dlog q(Z (s) ) =l

λ d 1+ Ld−1
sd−1(L)

状態遷移確率 τ の期待値は次のように求まる:

l|Td , σT2 )ψl (d) dTd−1lLd−1= N (Td−1 |ud−1 , sd−1 ) N Td |Td−1 , Ld−1 λ(L) d dTd−1ψ (d)
N log l/Ld |Td , σT2 l = N (Td |ud , sd ) (32)×(L)

q(Td |l1:T ) = αd (Td ) βd (Td ) =1
udvd1+N Td | −1,qdsd
qd + s−1
qd−1 + s−1dd (H)=βl ,d (t + 1) elog τl,d (l ,d )τκd (f )X(f,t+1)

l ,df⎧ (H)X(f,t+1)⎪βl−1,d (t + 1)
l>1⎪f κd (f )⎨ X(f,t+1)=exp (eτ d (d))
df κd (f )
⎪ ⎪
⎩(H)
× l βl ,d (t + 1) exp (el(d , l ))l=1(44)- 1073 -2011年3月

これらを用い、次の期待値を求める:

ηd (t) ∝(H)

lαt(H)(l, d) βt(l, d)ξd,l (d , t) ∝ αt−1 (1, d )eeτ d (d)+el(d) βt (l, d)

3(45)
(46)

評価実験

実験では、(1) 現状で多用されているシステムとの性能差 (2)LDS を用いた拍長モデルの有用性、(3) 音
色と音量に不確定性を持たせる LHA を用いることの有用性、の三点を評価する。(1) は、クロマベクトルの総コサイン距離最小化基準に基づく DTW を使用する。近年高性能である手法は、クロマベクトル同士の距離を DTW を用いて最小化するものが多い [6]。
(2) を評価するために、タイミングモデルに LDS を用いない手法を用意する。音価に比例するような音長の期待値を持った HSMM を用意した。固定されたテンポに依存するという意味では、このタイミングモデルは [8] と同等である。(3) を評価するために、調波構造と音量バランスに事前分布を持たせないものを用意する。スペクトルモデルは [5] と同等になる。
調波構造のモデルは [5] で用いられた値を使った。サンプリング周波数 8kHz、 分析フレームレート 20 E0と A0 は無情報に設定し、調波構造の事前分布は楽譜に記載された音高を平均とし標準偏差を 20 cent とした。CQT は 0.25 半音毎に評価した。
まず、RWC クラシック音楽データベース [10] 60 曲の楽譜表現 (SMF) に対し、シンセサイザーを用いて合成した音響信号を用意する。この音響信号を用いてスコアアライメントを行った結果の拍位置と、SMFから算出される拍位置の絶対誤差のパーセンタイルを評価基準として用いる。このような評価方法は、タイミング情報が正確に取れるというメリットがある。
また、実際に人間が演奏した録音でも同じような性能を発揮することが示唆されている [6]。
結果を表 1 に示す1 。人間の拍位置指定精度がおお
よそ 100 ミリ秒であることを踏まえると、オーケスト
ラのような複雑な楽器構成をもち音符が密である楽曲
でも、人間の拍位置精度と同程度の性能を 7 割方発揮する。また、現状多く使用されている手法 (Chroma)より、はるかに性能が高いことが分かる。LH と LHLを比較すると、タイミングモデルの有効性が示唆される。MLHL の結果から、音色と音量を固定した場合は、スペクトルをモデル化するアライメント手法は破綻することが分かる。これは、音色と音量に多様性を持たせることの重要性を表している。

4

Table 1 絶対推定誤差のパーセンタイル [ミリ秒]。
小さいほど高精度な推定。Chr は従来法、LH は時間長を独立に扱った本手法 ( p(Td ) = δ(Td − 10) に設定) 、MLHL は音量と音色を固定した本手法、LHL は提案手法。

まとめ

本稿では、音色や音量の不確定性を扱い、演奏のタイミングモデルを取り入れつつも、繰り返し構造といった、楽譜上の遷移を取り扱えるスコアアライメン
ト手法を提案した。また、音色音量モデルとタイミングモデルの有効性と、現状で多用されている手法の性能差を評価し、その有効性を確認した。今後の課題としては、単一パートのアライメントがある。今までの多くのアライメント手法は、楽譜位置と音響信号の対応付けを求めるが、実際には特定のパートが他より速く弾くといったことがある。単一パートのアライメントを、通常のアライメントから算出出来れば、音源分離や演奏分析といった、楽譜を援用した音楽音響信号分析の性能の向上が期待される。また、信号モデルにアタックや打楽器音も取り入れることにより、更なる精度の向上が期待できる。

歌声+ピアノ
伴奏楽器+ピアノ
伴奏
ピアノ
ソロ
小規模
アンサンブル
オーケ
ストラ

Chr
LH
MLHL
LHL
Chr
LH
MLHL
LHL
Chr
LH
MLHL
LHL
Chr
LH
MLHL
LHL
Chr
LH
MLHL
LHL

25%
88
13
749
7
68
14
863
8
90
17
1485
990
16
1927
10
123
38
3111
23

50%
289
37
2175
19
182
32
2549
21
304
48
4520
21
259
46
4296
22
394
104
10463
51

75%
831
184
4811
51
619
86
6437
45
1363
224
10468
50
891
131
8827
45
1384
574
21788
119

90%
2566
658
9973
119
2714
255
9373
93
6422
891
19415
126
2804
393
16260
886688479334275
805

95%
7319
1023
13737
220
9848
473
11219
163
11736
2040
26728
269
4710
816
25178
133
36550
16768
44847
2996

参考文献
[1] A. Maezawa et al. Violin Fingering Estimation
Based on Violin Pedagogical Fingering Model
Constrained by Bowed Sequence Estimationfrom Audio. In IEA/AIE, 2010.
[2] A. Maezawa et al. Query-By-Conducting: An
Interface to retrieve classical-music interpretations by real-time tempo input. In ISMIR,
pages 477–482, 2010.
[3] K. Itoyama et al. Integration and Adaptationof Harmonic and Inharmonic Models for Separating Polyphonic Musical Signals. In ICASSP,pages I–57–I–60, April 2007.
[4] Y. Han and C. Raphael. Desoloing Monaural
Audio Using Mixture Models. In ISMIR, pages
145–148, 2007.
[5] C. Raphael.
A Hybrid Graphical Modelfor Aligning Polyphonic Audio with MusicalScores. In ISMIR, pages 387–394, 2004.
[6] M. Muller and S. Ewert. Towards TimbreInvariant Audio Features for Harmony-BasedMusic. IEEE TASLP, 18(3):649–662, March2010.
[7] N. Hu et al. Polyphonic audio matching andalignment for music retrieval. In WASPAA,
pages 185–188, 2003.
[8] A.T. Peeling, P. Cemgil and S. Godsill. A Probabilistic Framework for Matching Music Representations. In ISMIR, pages 267–272, 2007.
[9] K. Yoshii and M. Goto. Infinite Latent Harmonic Allocation: A nonparametric Bayesian
approach to multipitch analysis. In ISMIR,
pages 309–314, 2010.
[10] M. Goto. Development of the RWC Music
Database. In Int’l Congress on Acoustics, volume I, pages 553–556, 2004.

1 本手法によるアライメントの推定結果に基づき拍位置を合成したオーディオデモを以下の URL にて公開している：
http://winnie.kuis.kyoto-u.ac.jp/~amaezaw1/alignment j

日本音響学会講演論文集- 1074 -2011年3月