Vol.2015-MUS-107 No.50
2015/5/24

情報処理学会研究報告
IPSJ SIG Technical Report

全極スペクトルモデルと擬似周期信号モデルのウェーブレッ
ト変換を用いた多重音スペクトログラムの調波時間因子分解
中村 友彦1,a)

亀岡 弘和1,2,b)

概要：多重音解析では，観測振幅スペクトログラムを非負値行列とみなし非負値行列因子分解を適用するアプローチと，計算論的聴覚情景分析に基づくアプローチが主に用いられてきた．本報告では，この２つのアプローチの利点を兼ね備えた新たなスペクトログラムモデルを導出し，それに基づく多重音解析手法
である調波時間因子分解について述べる．音源の連続時間信号モデルとして擬似周期信号から出発し，各音源のウェーブレット変換を導出した後，それらを重畳したもので音楽音響信号のスペクトログラムを表現する．さらに，離散時間信号領域で定義された自己回帰モデルをスペクトログラムモデルに導入にできることを示す．このモデルに対し，ビブラートやポルタメントで起こる基本周波数の時間変動や，調による音高の出現頻度の偏りなどの音楽特有の性質を補助情報として組み込んだ，効率的なパラメータ推論ア
ルゴリズムを導出する．

1. はじめに
多重音解析は，複数の音源信号が混合された観測信号から個々の音源の情報（基本周波数，発音開始時刻，音量な周波数成分をクラスタリングすることにより多重音解析を行う手法である．このアプローチでは，楽音のスペクトロ
グラムの局所的な構造に着目し楽音に詳細な仮定を課している．

ど）を得る処理である．これは音楽情報処理の重要課題の一方で，楽音に詳細な仮定を課さない手法として音楽の一つであり，音楽情報検索や自動採譜，音楽音響信号加工

大局的な構造に着目したアプローチも提案されている．こなど様々なアプリケーションの基礎技術である．

のアプローチでは，観測振幅スペクトログラムを非負値

多チャンネル音響信号を入力とする場合には音源の空間

行列とみなして非負値行列因子分解（Non-negative Matrix

的な手がかりを多重音解析に利用できるが，モノラル音

Factorization, NMF）[5] を適用する．NMF は，限られた種

響信号が入力である場合にはこれに代わる何らかの手が類の音高の楽音がそれぞれ異なるタイミングで繰り返し生

かりが必要である．モノラル音響信号を入力とする多重起するという音楽特有の性質に着目し，限られた種類のス

音解析のためのアプローチとして，計算論的聴覚情景分

ペクトルテンプレートの適当な重み付き和で各時刻の観測

析のコンセプトに基づく手法が提案されている [1–3]．計スペクトルを表現できるという仮定を用いている．すなわ

算論的聴覚情景分析とは，Bregman によって提唱された聴ち，観測振幅スペクトログラムは，2 つの非負値行列（各覚情景分析 [4] で示された人間の聴覚機能を計算機で実現列がスペクトルテンプレートを表す基底行列と，各行が対しようという試みである．調波時間構造化クラスタリン

応するスペクトルテンプレートの音量の時間発展を表すアグ（Harmonic-Temporal Clustering, HTC）[2, 3] はこのアプ

クティベーション行列）の積として表現される．したがっ

ローチに則り，人間が混在する複数の音の中で個々の音をて，観測振幅スペクトログラムをこの 2 つの行列に分解す聞き分けるために用いると考えられている要件（調波性，ることにより，スペクトルテンプレートと重みを同時推定

連続性，同時性，同期性）を時間周波数成分の局所的な制

し，各楽音の振幅スペクトログラムに分離できる．

約として記述し，当該要件を満たすように観測信号の時間

前者のアプローチは音源のスペクトログラムの局所的な構造，後者のアプローチは音楽音響信号のスペクトログラ

1

2

a)
b)

東京大学大学院情報理工学系研究科，東京都文京区本郷 7-3-1，
113-0033
日本電信電話株式会社 NTT コミュニケーション科学基礎研究所，
神奈川県厚木市森の里若宮 3-1，243-0198
Tomohiko Nakamura@ipc.i.u-tokyo.ac.jp
kameoka@hil.t.u-tokyo.ac.jp/kameoka.hirokazu@lab.ntt.co.jpⓒ 2015 Information Processing Society of Japanムの大局的な構造を手がかりとしている．これらは互いに矛盾するわけではなく，いずれも高精度な多重音解析を実現するには有用な手がかりである．そこで，我々はこれら

2 つの手がかりを同時に取り入れたスペクトログラムモデ

1Vol.2015-MUS-107 No.50
2015/5/24

情報処理学会研究報告
IPSJ SIG Technical Report

ルを構築し，それに基づく多重音解析手法，調波時間因子と対数瞬時 F0 （Ωk (t) = ln θ̇k (t)）を用いて，Wk (x, t) は分解（Harmonic-Temporal Factor Decomposition, HTFD）をN∑提案する [6–8].

ak,n (t)Ψ∗ (ne−x+Ωk (t) )ej(nθk (t)+φk,n ) ,

(4)

n=1

ところで，近年の音楽情報検索に関する国際会議や国際コンテスト [9] では，調推定，和音推定などの研究が急速と表せる．ここで，θ̇k (u) は瞬時 F0 である．ψ の Fourier 変に進展している．調，和音などの情報が高精度に得られる換である Ψ は任意に選べるので，ω = 1 で最大値をとる対のであれば，多重音解析において有用な補助情報となりう数正規分布型の実関数 [3] (ln ω)2−

e 4σ2
Ψ(ω) = 0る．そこで，HTFD では調推定，和音推定により得られた情報（音楽事前情報）をスペクトログラムモデルに組み込
み，モデルパラメータの推定に活用するフレームワークも提案する．音源分離では，ユーザにより入力された情報を用いて分離精度を向上するユーザガイド付き音源分離手
法 [10, 11] や，楽譜を補助情報とする音源分離手法 [12, 13]などの研究が進められており，提案するフレームワークも.

(5)応する．式 (5) より，Wk (x, t) はWk (x, t) =N∑ak,n (t)e−(x−Ωk (t)−ln n)2
4σ2ej(nθk (t)+φk,n ) .

(6)

n=1

と変形でき，さらに調波成分のパワースペクトルが加法的であると仮定すれば，|Wk (x, t)|2 は

2. 音楽音響信号のスペクトログラムの確率モ
デル化
2.1 楽音信号のウェーブレット変換(ω ≤ 0)を用いる．σ は Ψ(ω) を ln ω 軸で見たときの標準偏差に対補助情報付き音源分離の一種として位置づけられる．
以下，実数集合と複素数集合，虚数単位をそれぞれ√
R, C, j := −1 と表記する．

(ω > 0)|Wk (x, t)|2 ≃N∑|ak,n (t)|2 e−(x−Ωk (t)−ln n)2
2σ2

(7)

n=1

と書ける．このモデルは，HTC で採用された調波時間構造

本節では，[3] に倣って楽音信号のウェーブレット変換化モデル [3] と同一であり，その時刻 t での断面は調波的領域でのモデルを導出する．多くの楽音は局所的に周期的に正規分布形の関数が並んだ混合正規分布モデルと同形のとみなせるので，楽音の解析的なモデルとして周期や調波関数で表される．

成分のパワーが時間的に滑らかに変化する擬似周期信号をここまでスペクトログラムモデルを連続時間，連続対数用いることができる．音高のインデックスを k，調波成分

周波数領域で定義してきたが，計算機で実際に得られる観測のインデックスを n = 1, 2, . . . , N − 1 とする．音高 k の楽音

スペクトログラムは時刻 t と対数周波数 x に関して離散的の連続時間信号を，n 次調波成分の瞬時位相が nθk (u) ∈ R，である．そのため本節以降は，等間隔に量子化された時刻

瞬時振幅が ak,n (u) ∈ C の擬似周期信号はtm (m = 0, 1, . . . , M − 1) と対数周波数 xl (l = 0, 1, . . . , L − 1)

fk (u) =N∑を用いて，観測スペクトログラムを Yl,m := Y(xl , tm ) と表す．

ak,n (u)ej(nθk (u)+φk,n )(1)

同様に Ωk,m := Ωk (tm )，ak,n,m := ak,n (tm ) とする．

n=1と記述できる．ここで，u ∈ R は連続時間信号領域での時
刻，φk,n ∈ R は初期位相である． fk (u) の連続ウェーブレッ
ト変換を得るための基底関数 ψα,t (u) を(u − t)1
ψα,t (u) = √ψ
α2πα2.2 ソース・フィルタモデルの導入
ここで，擬似周期信号について再考してみよう．擬似周期信号には局所的に周期的で滑らかに秋季や振幅が変化す(2)

ること以外は仮定が課されていない．そのため，実際の楽音とは異なるスペクトル形状も許容してしまっている．そと定義する．ここで，α > 0 はスケールパラメータ，t ∈ R

こで，楽音の物理的な生成過程に着目し適切に楽音のスペ

は時間シフトパラメータ，ψ(u) はアドミッシブル条件を満

クトル形状に制約を与えることを考えよう．

たす中心周波数 1 のアナライジングウェーブレットであく表現でき，離散時間信号に対する自己回帰過程として記る．これを用いて fk (u) のウェーブレット変換は，Wk (ln α1 , t) =∫ ∞∑N−∞ n=1

楽音の物理的な生成過程はソース・フィルタモデルでよ述できる．しかし，調波時間構造化モデルが定義されたの

ak,n (u)ej(nθk (u)+φk,n ) ψ∗α,t (u)du(3)

はウェーブレット変換領域であるため，式 (7) のパラメー
タと自己回帰過程のパラメータとの対応関係を直接得ること書ける．ψ∗α,t (u) の優勢な部分は時刻 t の周りにのみ存在すとは難しい．そこで本節では，[14] に従ってこれらのパラ

るので，Wk (ln α1 , t) は時刻 t 周りの θk (u) と ak,n (u) に強く依存メータの関係を式 (1) のパラメータを介して得ることを目する．そこで，θk (u) ≃ θk (t) + θ̇k (t)(u − t), ak,n (u) ≃ ak,n (t) と近指す．

似し Parseval の定理を適用すれば，対数周波数（ x := ln(1/α)）ⓒ 2015 Information Processing Society of Japan

離散時間インデックスを i とし，時刻 tm での式 (7) の断

2

Vol.2015-MUS-107 No.50
2015/5/24

情報処理学会研究報告
IPSJ SIG Technical Report面に対応する連続時間信号モデルの離散時間表現を fk,m [i]音高 k の励起信号の調波成分の相対的な複素振幅，Uk,mとする．この fk,m [i] が P 次の自己回帰過程によっては音高 k の楽音の時刻 tm における正規化振幅と解釈でき∑る．ただし，Uk,m の正規化条件は k,m Uk,m = 1 とする．全βk,m [0] fk,m [i] =P∑

βk,m [p] fk,m [i − p] + ϵk,m [i],(8)極スペクトルモデル 1/|Bk,m (ejω )|2 は音色に対応するので，

βk,m [p] と Bk,m (z) から時刻インデックス m を削除する．こp=1と記述できるとする．ここで，βk,m [p] (p = 0, 1, . . . , P) は自れにより式 (12) は

己回帰過程のパラメータ（線形予測係数とも呼ばれる）で極システムと等価であるため，ϵk,m [i] は全極システムの励
起信号とみなせる．2.1 節で仮定したように， fk,m [i] の F0
は eΩk,m であるから，励起信号 ϵk,m [i] の F0 も eΩk,m でなくてはならない．したがって，ϵk,m [i] はϵk,m [i] =N∑

以上をまとめると，音高 k のパワースペクトログラムモ

vk,n,m ejne

Ωk,m
iu0

,

(9)

デル Ck,l,m は

Ck,l,m =Hk,l,m Uk,m ,Hk,l,m =グ周期であり，vk,n,m ∈ C は n 番目の調波成分の複素振幅を表す．ここで詳細は省略するが， fk,m [i] の離散時間 Fourier
変換（discrete-time Fourier transform, DTFT）に逆 DTFT を適用すると，

vk,n,mn=1 Bk,m (e
P
∑

N
∑

jneΩk,m iu0jneΩk,m u0)eβk,m [p]z−pと記述され，HTC や NMF と同様にパワースペクトログラ
ムの加法性を仮定すれば，観測パワースペクトログラムモ

として， fk,m [i] の異なる表現が得られる．式 (10) と式 (1)の離散時間表現を比較すると，全極システムのパラメータと 2.1 節で導入したパラメータの対応関係が，vk,n,mΩ

Xl,m =

K∑

Ck,l,mと表せる（図 1）．Xl,m には，実際の音源信号の擬似周期性の仮定からの逸脱による誤差，調波間干渉を無視したことによる誤差，パワーの加法性に起因する誤差，背景雑音や残響に起因する誤差など，様々な要因による誤差が存在す
る．提案法では，一つ一つの誤差要因を詳細にモデル化す
る代わりに，まとめて一挙に確率的な現象と捉えることに(12)する．ここで，Yl,m の確率分布が平均 Xl,m の Poisson 分布から生成されたとすると，と陽に得られる．

2.3 モデルパラメータに対する拘束NMF における重要な仮定は，音源のスペクトルを時変な成分と時不変な成分の積として表現することであり，これによりモデルパラメータの空間を適切に制限できる．したがって，このモデルにおいてもスペクトルの構成要素がそれぞれ時変な成分または時不変な成分のどちらとみなされるべきかが重要となるはずである．F0 はビブラートやポルタメント中は大きく時間変動し，スペクトルのスケー
ルの時間変動も大きい．一方で，楽器の音色については曲全体を通して比較的一定とみなせることが多い．
簡単化のため，|ak,n,m | と vk,n,m を調波インデックスに依存する成分と依存しない成分に以下のように分解する．

(13)

wk,n,m は音高 k の楽音の調波成分の相対的な振幅，w̃k,n,m はⓒ 2015 Information Processing Society of Japan

Xl,ml,m e−Xl,mY

Yl,m ∼ Pois(Yl,m ; Xl,m ) =√√|ak,n,m | = wk,n,m Uk,m , vk,n,m = w̃k,n,m Uk,m .

(17)

k=1

(11)

Bk,m (ejne k,m u0 )

(16)

n=1

(10)

p=0

|ak,n,m | =

(15)
(x −Ω
−ln n)2− l k,m22σw2k,n,m e

デル Xl,m は

N
∑

(14)

2.4 確率モデルとしての定式化と記述できる．ただし，u0 > 0 は離散時間表現のサンプリン

Bk,m (z) :=Ω

Bk (ejne k,m u0 )と書き直せる．

n=1

fk,m [i] =w̃k,n,mwk,n,m =ある．この自己回帰過程は βk,m [p] をパラメータとする全Γ(Yl,m )(18)と記述できる．この生成モデルの尤度関数を Xl,m に関し最大化する問題は，I ダイバージェンス基準における Xl,m のYl,m への最適フィッティング問題と等価である．
詳細は省略するが，[14] のように w̃k,n,m が平均 0，分散

ν の等方的な複素正規分布に従うとすれば，式 (14) より
2

wk,n,m は以下の Rayleigh 分布に従う．
()
ν
wk,n,m ∼Rayleigh wk,n,m ;Ω|Bk (ejne k,m u0 )|Ωwk,n,m2
jne k,m u0 2)|) )=e−wk,n,m /(2(ν/|Bk (e.
Ωk,m(ν/|Bk (ejne u0 )|)2(19)
2.5 従来のスペクトログラムモデルとの関連提案モデル Xl,m は，様々な仮定を置くことで従来のスペ
クトログラムモデルと同一になる．式 (16) のように Hk,l,m

3

Vol.2015-MUS-107 No.50
2015/5/24

情報処理学会研究報告
IPSJ SIG Technical Report

メントなどによって F0 が局所的には時間に関して連続的に変化する傾向がある．この 2 つの性質は，Ω の k 番目の行を転置した Ωk に対する確率分布 ql (Ωk ), qg (Ωk ) として

qg (Ωk ) = N(Ωk ; µk 1 M , ξk2 I M ),(20)

ql (Ωk ) = N(Ω k ; 0 M , τ2k D−1 ),図 1 HTFD のスペクトログラムモデル．緑枠で囲まれた部分の点線が各音高に対するスペクトル包絡を表す．

をパラメトリックな関数で表さず，それ自体をパラメータとして扱えば，Xl,m は可変基底 NMF [15] のスペクトログ


 1

−1

 0

D =  .
 ..
 0


0

−1

0

0

...

2

−1

0

...

−1

2
..
.

−1
..
.

...
..
.

...

0

−1

2

...

0

0

−1


0 

0 

0 ..  .
. 
−1

1

(21)

(22)

と記述できる．ここで，N(Ωk ; µ, Σ) は平均 µ，分散 Σ をもつ M 次元正規分布，1 M は全要素が 1 の M 次元のベクト
ル，0 M は M 次元の零ベクトル，I M は M × M の単位行列

ラムと一致する．さらに，Hk,l,m を時不変として扱えば通を表す．この 2 つの確率分布を用いて Ωk の事前分布は常の NMF のスペクトログラムモデルと一致する．ここで，p(Ωk ) ∝ ql (Ωk )αl qg (Ωk )αgl

時不変にした Hk,l,m が調波構造を持つことを仮定すれば調(23)

波 NMF [16, 17] のスペクトログラムモデルと同一である．

と設計できる．ハイパーパラメータ αl , αg により，事前分

一方で，式 (16) において Ωk,m を時不変にすれば，[18] のス

布に対する qg (Ωk ), qg (Ωk ) の寄与を調節可能である．

ペクトログラムモデルと一致する．さらに，Uk,m に拘束付

ポピュラー音楽やクラシック音楽では調性があるため，曲きの混合正規分布モデルと同形の関数を仮定すれば，HTC

中の調や和音によって音高の出現頻度に偏りがある．このにおける [2, 3] のスペクトログラムモデルと同一となる．

偏りは，Uk,m に対する確率分布として記述できる．また，時3. 音楽事前情報の組み込み
3.1 音楽事前情報の確率モデルへの導入方針

刻に関する Uk,m のスパース制約が NMF において効果的であることが知られており，このスパース制約も確率分布として記述でき，音高の偏りを表す確率分布と統合的に扱える．

現できるはずである．音楽事前情報には推定誤りを含むこ

これらの確率分布の導入を簡単にするため，Uk,m = Rk Ak,m∑∑と分解する．Rk := m Uk,m ( k Rk = 1) は音高方向に正規化∑された振幅であり，Ak,m := Uk,m /Rk ( m Ak,m = 1) は時刻方ともあるが，補助情報としての信頼度を確率と捉えれば，向に正規化された振幅である．これらにより音高の出現頻

各パラメータの事前分布として推論に組み込める．

度の偏りと時刻に関するスパース性を表す確率分布は，それ1 節で述べたように，前段処理で得られた音楽事前情報を，補助情報として活用できれば高精度な多重音解析を実また，従来有効であることが知られている事前分布や複数の推定結果を同時に活用することもできる．例えば，調と和音の情報が得られたときには，その調で出現しやすく
かつその和音で出現しやすい音高に，高い事前確率が割りぞれ R := [R0 , R1 , . . . , RK−1 ]⊤ , Ak := [Ak,0 , Ak,1 , . . . , Ak,M−1 ]⊤の事前分布として以下のように記述できる．

R ∼ Dir(R; γ(R) ),Ak ∼ Dir( Ak ; γ(A)
k )

(24)

当てられるべきである．この「かつ」に相当する演算は，⊤ここで，γ(R) := [γ1(R) , . . . , γ(R)K ] は R の事前分布のハイパー両方の条件を表す確率分布の積で表現できる．そのため，

パラメータである．出現頻度の高い音高に対応する γk(R) に複数の条件を同時に成立させるように事前分布を設計する

大きな値を与えれば，その音高が出現しやすいように R のには，各条件を表す確率分布の積をとり，それを正規化す

(A)
(A) ⊤
推定値を誘導できる．γk(A) := [γk,1, . . . , γk,M
] は Ak の事前ればよい．この考え方は Products of Experts (PoE) [19] と分布のハイパーパラメータであり，この値を小さくすれば呼ばれる．

よりスパースな Ak の推定値に誘導できる．

3.2 音楽事前情報を反映した事前分布の設計例

4. パラメータ推論アルゴリズム

具体例として，Ωk,m と Uk,m の事前分布を設計する．演

ここまでで構築した確率モデルのもとで，与えられ奏中の各音符の F0 を考えると，大域的にはその音符に対た 観 測 パ ワ ー ス ペ ク ト ロ グ ラ ム Y に 対 し ，事 後 確 率応する F0 の近辺に存在する．一方で，特にバイオリンな

p(Θ|Y) ∝ p(Y|Θ)p(Θ) を最大化するような w, Θ := {Ω, R, A}

どの弦楽器や管楽器による演奏では，ビブラートやポルタを求めたい．このとき最大化したい目的関数は，ⓒ 2015 Information Processing Society of Japan4

Vol.2015-MUS-107 No.50
2015/5/24

情報処理学会研究報告
IPSJ SIG Technical ReportJ(Θ) := ln p(Y|Θ) + ln p(Θ)(25)

である．ここで，w := {wk,n,m }k,n,m の定義域 W を用いて，∫ln p(Y|Θ) = ln
×

∏
W l,m

∏

l,m

(

Rayleigh wk,n,m ;∑p(w|β, Ω)  + ln p(Θ).
−
Xl,m + ln
q(w) 
l,m

ν

)
dw

Ω
|Bk (ejne k,m u0 )|ln p(Ωk ) + ln p(R) +

∑

k

k,n

∑

Pois(Yl,m ; Xl,m )k,n,mln p(Θ) =J + (λ, q(w), Θ)(xl −Ωk,m −ln n)2
∑
∑
Uk,m
w2k,n,m e− 2σ2=Eq(w)  Yl,m
λk,n,l,m lnc
λk,n,l,mと書ける．ただし，=c は定数を除いて一致することを示(26)
ln p( Ak ).

(31)

(27)

k

と書ける．
式 (26) の右辺に w に対する周辺化が入った対数関数があるため解析的に大域最適解を導くことは難しいが，補助関数法を用いることによって閉形式の更新則を導出できる．
補助関数法は，目的関数値 J(Θ) の上界となる関数（補助関数）を設計し，補助関数を補助変数と呼ぶ変数とパラメー

す．さらに，式 (31) の括弧内の第 2 項に対して，x に関す
る区分求積法による近似を用い，x < x0 と xL−1 < x の範囲の値が無視できるものとすると，∫ ∞∑1
Xl,m ≃X(x, tm )dx∆ x −∞
l
√
∑
2πσ ∑
=
Rk Ak,mw2k,n,m .
∆x
n
k

(32)

が得られる．この近似は，式 (31) の括弧内の第 2 項が Ωk,mにほとんど依存しないことを仮定することに相当する．
この近似を適用した J + (λ, q(w), Θ) を J ++ (λ, q(w), Θ) とタ Θ に対して交互に最大化することにより，J(Θ) を単調置き，これを用いて補助変数とパラメータに関する更新式

増加させる手法である．

を導出する．補助変数の更新式は等号成立条件であり，モ対数関数は凹関数であるため，式 (26) に Jensen の不等式を適用すると，

デルパラメータ Θ の更新式はそれぞれ J ++ (λ, q(w), Θ) に関する各パラメータでの偏微分が 0 となる値を求めること


∫
∑
∑
Xl,m ∑
ln p(Y|Θ) ≥q(w)  Yl,m ln−
Xl,m +
Yl,m
Yl,m
W
l,m
l,m
l,m
∑
Ωk,m
+
ln Rayleigh(wk,n,m ; ν/Bk (ejne u0 )|)k,n,mにより導出できる．詳しい導出や更新式は本稿では省略する（詳細は [7] 参照）．

5. 多重音解析の動作確認実験と定量評価実験
5.1 時間変化する F0 の推定実験− ln q(w)) dw(28)

HTFD が時間的に変化する F0 を推定できるかを確認す∫
補助変数 q(w) は W q(w)dw = 1, q(w) ≥ 0 を満たす. 等号成るため，RWC 楽器音データベース [20] から D♭4，F4，A♭4

立条件は q(w) が w の事後分布と一致するときである．以
∫下では，Eq(w) [w2 ] := W q(w)w2 dw と表記する．

したサンプリング周波数 16 kHz の音響信号を分離する実験のバイオリン音源のビブラート音を用いて，人工的に合成を行った．簡単のため，本節と次節の実験では全極スペク

式 (17) の Xl,m は k, n に関する和を含むため，式 (28) のトルモデルを省いて wk,m,n の時不変性を仮定した HTFD [6]

括弧内の第 1 項は対数関数の中に和を含むため，解析的にを用いた．スペクトログラムを求める際には，時間シフ解くのが困難である．そこでさらに，式 (28) の括弧内の第ト間隔を 14.6 ms，解析周波数を 55 Hz から 7040 Hz まで

1 項に Jensen の不等式を適用すると，

10 cent 間隔に設定し，高速近似ウェーブレット変換 [21]Yl,m ln Xl,m ≥ Yl,m∑

λk,n,l,m ln

w2k,n,m e−k,n(xl −Ωk,m −ln n)22σ2λk,n,l,mを行った．アナライジングウェーブレットとして，式 (5)

Uk,m

,

(29)

で定義された対数正規分布型のウェーブレットを用いた（σ = 0.02）．調波成分の個数を N = 8，音高数を K = 73 とし，µk として A1（55 Hz）から A♯7 までの基本周波数にそれぞれ対応させた．他のパラメータは，全ての k に対しとして 下界 が 得 られ る ．補 助 変数 λ := {λk,n,l,m }k,n,l,m は∑
λk,n,l,m ≥ 0, k,n λk,n,l,m = 1 を満たし，等号成立条件は−6γ(A)k = (1 − 3.96 × 10 )1I ，τk = 0.83，vk = 1.25，αg = α s = 1，

γ(R) = (1 − 2.4 × 10−3 )1K とした．
図 2 に推定されたスペクトログラムモデルを示す．F0 のλk,n,l,m =w2k,n,m e−(xl −Ωk,m −ln n)2
2σ2

Uk,m

Xl,m

時間変化を記述できない NMF（基底数 3，I ダイバージェ(30)ンス基準）では，ビブラート時の F0 の時間変化が平均化されてしまっていることが確認できる（図 2 (a)）
．一方で，

+

である．したがって，J(Θ) の補助関数 J (λ, q(w), Θ) は
ⓒ 2015 Information Processing Society of JapanHTFD により得られたスペクトログラムモデルでは，図5

Vol.2015-MUS-107 No.50
2015/5/24

情報処理学会研究報告
IPSJ SIG Technical Report

554

311
2

6

0

Frequency [Hz]

554

311
0

311
2

4

Time [s]

554

311

6

0

2

4

6

4

6

Time [s](a) HTFD により推定されたスペクトログラムモデルと F0 軌跡（橙色の曲線）．
988
988

988

Frequency [Hz]4Time [s]

554

2

4

Time [s]6Frequency [Hz]0988Frequency [Hz]988

Frequency [Hz]Frequency [Hz]

988

554

311
0

2

4

Time [s]

554

311

6

0

2

Time [s](b) NMF により推定されたスペクトログラムモデル
図2

HTFD [6] と NMF により推定されたスペクトログラムモデル（ [8] より抜粋）．左から順に，それぞれ D♭4，F4，A♭4 に対応．

A♭4

Pitch

F4
D4
D♭4

Time
(a) 調の情報なし
図4A♭4

F0 が時間変化しない拘束をおいた HTFD（Harmonic NMF）と，
全極スペクトルモデルを導入していない HTFD（HTFD）[6]，
全極スペクトルモデルを導入した HTFD（HTFD+）[7] で得ら

Pitch

F4

れた各楽曲に対する平均 SNR 改善量と標準誤差．

D4D♭4のアクティベーションを抑制できるかを確認する．この音Time(b) 調の情報あり
図 3 調の情報を用いない場合 (a) と用いた場合 (b) の HTFD により

響信号の調を変ニ長調と仮定し，音階内の音高に対応する γk(R) は 1 − 2.4 × 10−3 ，音階外の音高に対応する γk(R) は

1 − 3.0 × 10−3 とした．他の条件は前節と同一である．

推定された A3 から A♭4 に対応するアクティベーション（ [8] よ調の情報を用いない場合（図 3 (a)）には，実際に音響信り抜粋）．赤色の曲線は D4 のアクティベーションに対応する．

号に含まれていないにも関わらず，D4 に対応するアクティ
ベーションが高くなってしまっている．この場合に比べ，

2 (b) のように適切に F0 軌跡が推定された．

調の情報を用いた上述の事前分布を用いると，D4 は音階外であるため，この音高に対応するアクティベーションが5.2 音楽事前情報の効果確認実験
次に，調の情報を用いて事前分布を設計し音階外の音高ⓒ 2015 Information Processing Society of Japan抑制された．したがって，γk(R) に調の情報を反映し多重音解析性能を改善できることが確認できた．

6

Vol.2015-MUS-107 No.50
2015/5/24

情報処理学会研究報告
IPSJ SIG Technical Report

5.3 モノラル音源分離性能の定量評価実験

らが行った研究の成果をまとめたものであり，実装や実最後に，全極スペクトルモデルの効果を確認するためモ

験，議論を通して本研究に貢献していただいた両氏に感謝

ノラル音楽音響信号の各音高への分離実験を行った．比較する．本研究は JSPS 科研費 26730100 の助成を受けたもの

手法として，全極スペクトルモデルを用いない HTFD [6]である．

と，さらに F0 を時不変とした HTFD （スペクトルテンプ
レートを，調波的に正規分布形の関数が並んだ混合正規分

参考文献

布モデルと同形の関数に拘束した場合の調波 NMF [16, 17][1]に対応．以後，Harmonic NMF）を用いた．分離では時間∑
周波数マスクを Ck,l,m / k Ck,l,m と設計して用いた．ただし，[2]

提案法では w の代わりに事後分布 q(w) が推定されるため，

w2 の推定値を Eq(w) [w2 ] として計算した．各音高ごとの演奏の録音を用意するのは困難であったため，RWC クラシッ

[3]

ク音楽データベース [20] の RM-C001 から RM-C005 の最初の 30 秒を MIDI シンセサイザー FluidSynth [22] で合成[4]した音響信号（サンプリング周波数 16 kHz）を入力とし[5]て用いた．スペクトログラムの計算には，時間シフトを14.6 ms とした高速近似ウェーブレット変換を用いた．調波成分の個数は N = 20 とし，全ての k に対して τk = 1.0，

[6]

−4
(R)
γ(A)
= 0.8 × 1K ，P = 20, ν = 1 とk = (1 − 1.0 × 10 )1I ，γした．他のパラメータやアナライジングウェーブレットは[7]

5.1 節と同一とした．また，Harmonic NMF は 100 イテレー
ション，HTFD と HTFD+は 20 イテレーション時点での推[8]定値を用いた．
楽曲毎の分離音の signal-to-noise ratio（SNR）改善量の平均と標準誤差を図 4 に示す．Harmonic NMF と HTFD を比べると SNR の差の平均値は約 0.02 dB であり，分離性能に大きな差は見られなかったが，MIDI で音源を作成し

[9]
[10]

たため F0 の変動が少なかったからと考えられる．一方，

HTFD+は全曲で Harmonic NMF と HTFD よりも SNR が約
0.80 dB，約 0.78 dB 平均的に改善した．この結果より，全[11]

極スペクトルモデルの導入によってモノラル音源分離性能が向上することを確認できた．

6. 結論

[12]

本報告では，NMF と HTC の利点を兼ね備えた新たなス
ペクトログラムモデルを提案し，そのモデルに基づく多重[13]

音解析手法 HTFD について述べた．各音高の音響信号の連
続時間におけるモデルとして擬似周期信号から出発し，そのウェーブレット変換を近似的に導出した．また，離散時[14]間信号領域で定義された全極システムをウェーブレット変換で定義されたスペクトルモデルに組み込んだ．また，音楽事前情報を事前分布として導入できることを示し，複数[15]の音楽事前情報を用いる場合や従来用いられてきた事前分
布と併用するときにも，PoE を用いて事前分布を適切に設[16]計する方法を示した．提案した確率モデルに対して，閉形式の更新式からなるパラメータ推論アルゴリズムを導いた．
モノラル音源分離実験により提案法の有効性を確認した．
謝辞 本稿は，四方紘太郎氏，高宗典玄氏らと共に著者ⓒ 2015 Information Processing Society of Japan[17]

Hu, G. and Wang, D. L.: An auditory scene analysis approach
to monaural speech segregation, Topics in Acoust. Echo andNoise Contr., pp. 485–515 (2006).
Kameoka, H., Nishimoto, T. and Sagayama, S.: A Multipitch
Analyzer Based on Harmonic Temporal Structured Clustering, IEEE Trans. Acoust., Speech, and Language Process.,
Vol. 15, No. 3, pp. 982–994 (2007).
Kameoka, H.: Statistical Approach to Multipitch Analysis,
PhD Thesis, The University of Tokyo (2007).
Bregman, A. S.: Auditory scene analysis: The perceptual organization of sound, MIT press (1994).
Smaragdis, P. and Brown, J. C.: Non-negative matrix factorization for polyphonic music transcription, Proc. IEEE Workshop Applications Signal Process. Audio Acoust., IEEE, pp.
177–180 (2003).
四方紘太郎，高宗典弘，中村友彦，亀岡弘和：調波時間
因子分解法に基づく事前情報付き多重音解析，情処研報，
No. 39 (2014).
中村友彦，亀岡弘和：全極スペクトルモデルを用いた調波
時間因子分解による多重音解析，情処研報，No. 26 (2015).
Nakamura, T., Shikata, K., Takamune, N. and Kameoka,
H.: Harmonic-Temporal Factor Decomposition Incorporating Music Prior Information for Informed Monaural Source
Separation, Proc. Int. Symposium Music Info. Retrieval, pp.
623–628 (2014).
[Online: 18, Apr. 2015], http://www.music-ir.org/
mirex/wiki/MIREX_HOME.
Smaragdis, P. and Mysore, G. J.: Separation by ”humming”:
User-guided sound extraction from monophonic mixtures,
Proc. IEEE Workshop Applications Signal Process. Audio
Acoust., IEEE, pp. 69–72 (2009).
Ozerov, A., Févotte, C., Blouet, R. and Durrieu, J. L.: Multichannel nonnegative tensor factorization with structured constraints for user-guided audio source separation, Proc. Int.
Conf. Acoust. Speech Signal Process., IEEE, pp. 257–260(2011).
Hennequin, R., David, B. and Badeau, R.: Score informed
audio source separation using a parametric model of nonnegative spectrogram, Proc. Int. Conf. Acoust. Speech Signal
Process., pp. 45–48 (2011).
Simsekli, U. and Cemgil, A. T.: Score guided musicalsource separation using generalized coupled tensor factorization, Proc. Eur. Signal Process. Conf., IEEE, pp. 2639–2643(2012).
亀岡弘和：全極型声道モデルと F0 パターン生成過程モデ
ルを内部にもつ統計的音声スペクトルモデル，信学技報，Vol. SP2010-74, pp. 29–34 (2010).
Nakano, M., Le Roux, J., Kameoka, H., Ono, N. andSagayama, S.: Infinite-state spectrum model for music signal analysis, Proc. Int. Conf. Acoust. Speech Signal Process.,
pp. 1972–1975 (2011).
Raczyński, S. A., Ono, N. and Sagayama, S.: Multipitch analysis with harmonic nonnegative matrix approximation, Proc.
Int. Conf. Music Info. Retrieval, pp. 381–386 (2007).
Vincent, E., Bertin, N. and Badeau, R.: Harmonic and inharmonic Nonnegative Matrix Factorization for Polyphonic
Pitch transcription, Proc. Int. Conf. Acoust. Speech Signal

7

情報処理学会研究報告
IPSJ SIG Technical Report

[18]

[19]

[20]
[21]
[22]

Vol.2015-MUS-107 No.50
2015/5/24

Process., pp. 109–112 (2008).
Yoshii, K. and Goto, M.: Infinite Latent Harmonic Allocation: A Nonparametric Bayesian Approach to Multipitch
Analysis, Proc. Int. Soc. Music Info. Retrieval, pp. 309–314(2010).
Hinton, G. E.: Training products of experts by minimizing
contrastive divergence, Neural Comput., Vol. 14, No. 8, pp.
1771–1800 (2002).
Goto, M.: Development of the RWC Music Database, Proc.
Int. Congress Acoust., pp. l–553–556 (2004).
亀岡弘和，田原鉄也，西本卓也，嵯峨山茂樹：信号処理
方法及び装置. 特開 2008-281898, (20. Nov. 2008).
[Online: 21, Apr. 2015], http://www.fluidsynth.org/.

ⓒ 2015 Information Processing Society of Japan

8