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

全極スペクトルモデルを用いた調波時間因子分解による多重音解析
中村 友彦1,a)

亀岡 弘和1,2,b)

概要：多重音解析では，観測振幅スペクトログラムを非負値行列とみなし非負値行列因子分解を適用するアプローチと，計算論的聴覚情景分析に基づくアプローチが主に用いられている．我々は，この２つのア
プローチの利点を兼ね備えた新たなスペクトログラムモデルを導出し，それに基づく多重音解析手法である調波時間因子分解を以前提案した．これらに加え，ソース・フィルタモデルに代表される楽音の物理的な生成過程も，多重音解析の性能向上には重要である．そこで本報告では，離散時間信号領域で定義されたソース・フィルタモデルが調波時間因子分解の音源スペクトログラムモデルに導入できることを示し，
そのパラメータ推論法について述べる．定量評価実験により，この導入によってモノラル音源分離性能が向上することを示す．

1. はじめに当該制約を満たすように観測信号の時間周波数成分を時間
周波数平面上でクラスタリングする．

多重音解析は，複数の音源信号が混合された観測信号か

一方，2 つ目のアプローチの代表的な手法は，観測振幅

ら個々の音源の情報（基本周波数，発音開始時刻，パワー

スペクトログラムを非負値行列とみなして非負値行列因子など）を得る処理であり，音楽情報検索や自動採譜，音響分解（Non-negative matrix factorization, NMF）を適用する

信号の自動編集など様々なアプリケーションの基礎技術でものである [5]．この手法は，限られた種類の音高の楽音ある．

がそれぞれ異なるタイミングで繰り返し生起するという音

多重音解析を行う際に，観測信号が多チャンネル信号で楽特有の性質に着眼し，各フレームのスペクトルが限られあれば音源の空間的な手がかりを利用することができる

た種類のスペクトルテンプレートの適当な重み付き和で表が，モノラル信号であるときはこの手がかりに代わる何ら現できるという仮定を用いている．スペクトルテンプレー

かの仮定が必要となる．モノラル信号を入力とする多重音

トと重みは非負値であるため観測信号は 2 つの非負値行列解析では，主に以下の２つのアプローチが取られてきた．1の積として記述される．したがって，観測振幅スペクトロ

つ目のアプローチは，計算論的聴覚情景分析のコンセプト

グラムをこの 2 つの行列の積に分解することによって，スに基づくものである（例，[1]）．計算論的聴覚情景分析とペクトルテンプレートと重みを同時推定し，観測振幅スペ

は，Bregman によって提唱された聴覚情景分析 [2] で示さクトログラムを音高毎のスペクトログラムに分離できる．

れた聴覚機能を計算機で実現しようという試みである．ここれら 2 つのアプローチは異なる手がかりを元に分離をのアプローチに則った手法である調波時間構造化クラスタ行っている．前者のアプローチは時間周波数表現の局所的リング（Harmonic-temporal clustering, HTC）[3, 4] は，同一な性質，後者は音楽信号特有の大域的な性質を利用しておの音源に由来する時間周波数成分が音脈と呼ばれる知覚的り，これらは相容れない関係ではなくいずれも高精度な多な構造に群化されるプロセスを模倣したものである．HTC

重音解析を実現するためには有用である．この考えに基づでは，群化の要件（調波性，連続性，同時性，同期性など）

き，我々はこれら 2 つの性質を同時に取り入れたスペクトが時間周波数領域の局所的な制約として記述されており，

ログラムモデルと，それに基づく多重音解析手法である調波時間因子分解（Harmonic-temporal factor decomposition,

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

c 1959 Information Processing Society of Japan⃝

HTFD）を提案した [6]．
HTFD の多重音解析性能をさらに向上させるためには，音源の物理的な生成過程に基づく手がかりが有用である．
管楽器や弦楽器から発せられた各音高の楽音は，管楽器に1

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

吹き込んだ呼気や弦の振動が楽器内部で共鳴して生成されak,n (u) の値のみに依存する．そのため，以下のように θk (t)る．前者は主に音高に対応し，後者は音色に対応することと ak,n (t) を 1 次の Taylor 展開で近似する．

が知られている．音高に対応する物理量である基本周波数（F0 ）はビブラートやポルタメントなどで時間的に変動する

ak,n (u) ≃ ak,n (t),θk (u) ≃ θk (t) + θ̇k (t)(u − t).

(4)ことも多いが，音色は比較的時不変とみなせる．したがっここで，θ̇k (u) は瞬時 F0 である．上述の近似を用いつつて，これらの性質を導入した音源スペクトログラムモデル

Parseval の定理を適用すれば，対数周波数（ x := ln(1/α)）が適切に構築できれば，実際の楽器では起こりえないスペと対数瞬時 F0 （Ωk (t) = ln θ̇k (t)）を用いて，式 (3) をクトルの時間変動が抑制できるはずである．
本報告では，HTFD の多重音解析性能向上のため，上述の音源の物理的な生成過程に基づく手がかりを離散時間領Wk (x, t) =N∑ak,n (t)Ψ∗ (ne−x+Ωk (t) )ej(nθk (t)+φk,n ) ,

(5)

n=1

域で定義されたソース・フィルタモデルを用いて記述し，と変形できる．Ψ は ψ の Fourier 変換であり任意に選ぶこ

ある特定のウェーブレット変換領域で定義された HTFD のとができるので，以下の ω = 1 で最大値をとる対数正規分

スペクトログラムモデルに導入する．また，構築したスペ

布型の実関数を用いる．
 (ln ω)2−

e 4σ2
Ψ(ω) = 0

クトログラムモデルの効率的なパラメータ推論法を提案する．以下，実数集合と複素数集合，虚数単位をそれぞれ√
R, C, j := −1 と表記する．

2. 音楽音響信号のスペクトログラムの確率モ
デル
2.1 楽音信号のウェーブレット変換(ω ≤ 0).

(6)

σ は Ψ(ω) を ln ω 軸で見たときの標準偏差に対応する．これを用いて，式 (5) はWk (x, t) =本節では，[6] と同様に [4] に倣って楽音のウェーブレッ

(ω > 0)N∑

ak,n (t)e−(x−Ωk (t)−ln n)2
4σ2ej(nθk (t)+φk,n ) .

(7)

n=1

ト変換領域でのモデルを導出する．多くの楽音（特に調波と変形できる．ここで，n, n′ (n , n′ ) の指数項の重なりがほ楽器音）は，擬似周期信号（局所的に周期的とみなせ，周とんどない，すなわち調波性分が互いに重ならないと仮定

期や調波成分のパワーが時間的に滑らかに変化する信号）

できれば，|Wk (x, t)|2 は近似的にとみなせる．ここで，音高のインデックスを k，調波成分のインデックスを n = 1, 2, . . . , N − 1 とする．音高 k の楽音|Wk (x, t)|2 ≃の連続時間信号を，n 次調波成分の瞬時位相が nθk (u) ∈ R，
瞬時振幅が ak,n (u) ∈ C の擬似周期信号の解析信号表現

fk (u) =N∑N∑|ak,n (t)|2 e−(x−Ωk (t)−ln n)2
2σ2

(8)

n=1

と書ける．この仮定は調波成分のパワースペクトルが加法的であると近似していることに相当する．式 (8) で与えら

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

(1)

n=1

れるスペクトログラムモデルは HTC で用いられた調波時間構造化モデル [4] と同一であり，その時刻 t での断面は調として与える．ここで，u ∈ R は連続時間信号領域での時

波的に正規分布形の関数が並んだ混合正規分布モデルと同刻，φk,n ∈ R は初期位相である．この信号表現では，群化形の関数で表される．さらに，パワースペクトルの加法性の要件のうち調波性と周波数変調の同期性が暗に満たされを仮定すれば，K 個の楽音が重畳された音響信号のパワーている．ここで，アドミッシブル条件を満たす中心周波数

スペクトログラムは，式 (8) の k について和をとったもの

1 のアナライジングウェーブレットを ψ(u) ∈ C と書くと，として与えられる．

ウェーブレット基底関数 ψα,t (u) は1ψψα,t (u) = √2πα

ここまでスペクトログラムモデルを連続時間，連続対数(u − t)α(2)とかける．ここで，α > 0 はスケールパラメータ，t ∈ R

周波数領域で定義してきたが，我々が計算機によって実際に得ることのできる観測スペクトログラムは離散的な表現である．そのため本節以降は，等間隔に量子化された時間

は時間シフトパラメータである．このウェーブレット基底

tm (m = 0, 1, . . . , M − 1) と対数周波数 xl (l = 0, 1, . . . , L − 1)関数の複素共役 ψ∗α,t (u) を用いて fk (u) のウェーブレット変を用いて，観測スペクトログラム Yl,m := Y(xl , tm ) を表す．

換は，Wk (ln α1 , t) =∫ ∞∑N−∞ n=1ak,n (u)ej(nθk (u)+φk,n ) ψ∗α,t (u)du.

(3)

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

2.2 ソース・フィルタモデルの導入
1 節で述べたように，楽音の物理的な生成過程を適切にと得られる．典型的に ψ∗α,t (u) の優勢な部分が t の周りに局記述できれば，より詳細なスペクトログラムモデルが得ら在化しているので，式 (3) の積分結果は t 周りの θk (u) とれ多重音解析性能が向上するはずである．楽音の物理的なc 1959 Information Processing Society of Japan⃝2

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

ると，

fk,m [i] =N∑

vk,n,mn=1 Bk,mΩk,mΩ(ejne k,m u0 )

ejne

iu0

.

(13)

を得ることができ，これは fk,m [i] の別表現となっている（図 1 の灰色矢印）．式 (13) と式 (1) の離散時間表現

fk,m [i] =N∑

ak,n,m ej(neΩk,miu0 +φk,n )

(14)

n=1

を比較すると，全極システムのパラメータと 2.1 節で導入
図 1 自己回帰過程のパラメータと調波時間構造化モデルのパラメータの対応関係．

したパラメータの対応関係（図 1 の橙色矢印）を以下のよ
うに陽に記述することができる．

生成過程はソース・フィルタモデルでよく表現でき，このモデルは離散時間信号に対する自己回帰過程として記述|ak,n,m | =vk,n,mΩ

Bk,m (ejne k,m u0 )

.

(15)

できる．しかし，調波時間構造化モデルが定義されたのはウェーブレット変換領域であるため，式 (8) のパラメータ

2.3 モデルパラメータに対する拘束と自己回帰過程のパラメータとの対応関係を直接得ること

前節では全極システムをスペクトログラムモデルに導入は難しい．そこで本節では，[7] に従ってこれらのパラメー

したが，そのパラメータの値によっては楽音の性質とはかタの関係を式 (1) で定義された擬似周期信号の解析的な表

け離れたものも表しうる．そのため，モデルパラメータに現を介して得ることを目指す（図 1）．

適切な拘束を置く必要がある．NMF の重要な仮定は，各

時刻 tm における式 (8) の断面に対応する連続時間信号モ

音高のスペクトルを時変な成分と時不変な成分の積として

デルの離散時間表現を fk,m [i] とする．i は離散時間インデッ

表現したことにある．すなわち，どのスペクトルの構成要

クスである．この fk,m [i] が P 次の自己回帰過程によって素を時変とみなすか時不変とみなすががこの仮定の本質で

βk,m [0] fk,m [i] =P∑あり，ここではそれに従って考えてみよう．F0 はビブラー

βk,m [p] fk,m [i − p] + ϵk,m [i],(9)トやポルタメント中は時間的に変動すると仮定されるべきであるし，スペクトルのパワーについても同様である．

p=1と記述されると考えよう．βk,m [p] (p = 0, 1, . . . , P) は自己回

（NMF もパワーは時変としている．）しかし，楽器の音色

帰過程のパラメータ（線形予測係数とも呼ばれる）である．

については曲全体を通して比較的一定とみなせる．

この自己回帰過程は βk,m [p] をパラメータとして持つ全極

これらの仮定は以下のようにスペクトログラムモデルにシステムと等価であるため，ϵk,m [i] は全極システムの励起反映させられる．導出を簡単にするため，|ak,n,m | を以下の信号とみなせる．2.1 節で仮定したように， fk,m [i] の F0 は

ように 2 つの変数の積に分解する．

Ωk,m

e

であるため，励起信号 ϵk,m [i] の F0 も e

Ωk,m

である必要

がある．したがって，ϵk,m [i] はϵk,m [i] =N∑

vk,n,m ejne

Ωk,m
iu0

(10)

と記述できる．ただし，u0 > 0 は離散時間表現のサンプリング周期（観測信号のサンプリング周期と一致）であり，

vk,n,m ∈ C は n 番目の調波成分の複素振幅である．式 (9) に対して離散時間 Fourier 変換（discrete-time Fourier transform,

DTFT）を適用すると， fk,m の DTFT は√N
2π ∑vk,n,m δ(ω − neΩk,m u0 ),Fk,m (ω) =Bk,m (ejω ) n=1
Bk,m (z) :=βk,m [p]z−p(16)

wk,n,m は音高 k の楽音の時刻 tm における n 番目の調波成分
,

n=1

P
∑

√
|ak,n,m | = wk,n,m Uk,m .

の相対的なパワー，Uk,m は音高 k の楽音の時刻 tm における正規化振幅と解釈できる．ただし，Uk,m の正規化条件は∑
k,m U k,m = 1 である．同様に，vk,n,m についても√
vk,n,m = w̃k,n,m Uk,m .

(17)と分解できる．全極スペクトル 1/|Bk,m (ejω )|2 は音高 k の楽音の音色に対応するので，時不変なものとして扱いたい．

(11)

(12)

p=0

そのためには単純に βk,m [p] と Bk,m (z) から時刻インデック
ス m を削除すればよく，式 (15) は

wk,n,m =w̃k,n,mΩ

Bk (ejne k,m u0 )(18)

として得られる．ここで，ω は正規化角周波数，δ(ω) はと書き換えられる．HTFD では，物理的な生成過程を導入

Dirac のデルタ関数である．式 (11) に逆 DTFT を適用すしていなかったために全極スペクトルの時不変性を用いる

c 1959 Information Processing Society of Japan⃝3

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

確率モデルとして定式化する．これまでに用いた仮定や近似は現実には常に正確に成り立つとは限らないため，観測
パワースペクトログラム Yl,m は最適なパラメータが得られたとしても Xl,m と一致しないかもしれない．一つ一つの誤差要因を詳細にモデル化する代わりに，確率的生成モデル
によってこの逸脱を確率的現象と捉えることにする．ここで，Yl,m の確率分布が平均 Xl,m の Poisson 分布から生成されたとすると，Xl,ml,m e−Xl,mY

Yl,m ∼ Pois(Yl,m ; Xl,m ) =(21)

Γ(Yl,m )と記述でき，この尤度関数は
図 2 全極スペクトルモデルを導入したスペクトログラムモデル．緑∏

枠で囲まれた部分の点線が各音高に対するスペクトル包絡を表す．

Pois(Yl,m ; Xl,m ),

(22)

l,m

ことができず，解を制限するために wk,n,m を時不変なパラと書ける．ただし，Y := (Yl,m )l,m は観測パワースペクトロ

メータとして扱うことと等価な処理を行わざるをえなかっ

グラムの行列表現である．式 (22) の Xl,m に関する最大化た．一方で，提案モデルでは wk,n,m が時変なパラメータとは，I ダイバージェンス基準における Xl,m の Yl,m への最適して扱えるのが主要な違いであり，これにより励起信号のなフィッティングと等価である．

複素振幅の変化に伴うスペクトル形状の変化を記述できる．
まとめると，音高 k のパワースペクトログラムモデル

式 (18) から，条件付き確率分布 p(w|w̃, β, Ω) は Dirac のデルタ関数を用いて

Ck,l,m は音量軌跡 Uk,m ，調波成分の相対振幅 wk,m,n ，F0 軌跡

p(w|w̃, β, Ω) =Ωk,m によって N
(xl −Ωk,m −ln n)2 ∑ 2−2
2σCk,l,m =  wk,n,m e Uk,m ,∏ (δ wk,n,m −Ωu

Bk (ejne k,m 0 )

k,n,m

(19)

)

w̃k,n,m(23)と書ける．ただし，w := {wk,n,m }k,n,m , w̃ = {w̃k,n,m }k,n,m であり，β は (k, p) 番目の要素が βk [p] の行列，Ω は (k, m) 番目

n=1と記述され，観測信号のパワースペクトログラムモデルの要素が Ωk,m の行列を表す．これを元に条件付き確率分

Xl,m はそれらが重畳されたものとして以下のように表され

布 p(w|β, Ω) を導出するには， p(w̃) を定義し w̃ について周る（図 2）．

辺化すればよい．[7] のように，w̃k,n,m が等方的な複素正規

Xl,m =

K∑

Ck,l,m .

(20)

分布

k=1

提案モデル Xl,m は様々な仮定を置くことで従来のモデルと同一になる．wk,n,m を時不変にしたものを曲全体の音量を表す V > 0 と n 番目の調波成分の相対的なパワーを正∑規化した vk,n ( n vk,n = 1) の積に分解すれば，HTFD のス
ペクトログラムモデルと同一となる．また，式 (20) の括∑
弧内の項を Hk,l,m とおけば Xl,m = k Hk,l,m Uk,m と変形できるため，HTFD の場合と同様に NMF の時不変なスペクト
ルテンプレートが時変に拡張されたともみなせる．他にも

e−|w̃k,n,m | /νπν22w̃k,n,m ∼NC (w̃k,n,m ; 0, ν2 ) =

2

(24)

から生成されると仮定すると，wk,n,m は以下の Rayleigh 分
布に従う．

(
wk,n,m ∼Rayleigh wk,n,m ;
=

)

ν
Ω

|Bk (ejne k,m u0 )|

wk,n,mΩ(ν/|Bk (ejne k,m u0 )|)2

e−wk,n,m /(2(ν/|Bk (e2Ωjne k,m u0 2
)|) )

.

(25)

HTFD と同じく様々な NMF の拡張や HTC との関連もある
が，既に [6] で記述しているため省略する．時間周波数表

これは全極スペクトルの効果が事前分布として確率モデル

現として短時間フーリエ変換を用いた場合に，Hk,l,m を時に導入できることを示している．

不変にした後励起信号のパワースペクトルを周波数点ごと

他のパラメータである Ω や U の事前分布についてはにパラメータとみなせば複合自己回帰系 [8] と一致し，全HTFD と同一のものを用いる．F0 に関しては 2 つの物理

極スペクトルモデル自体も周波数点ごとにパラメータとみ的な性質を用いて事前分布を設計している．1 つ目の性質

なせば [9] と一致する．

はビブラートやポルタメント中には F0 が時間に関して連
続的に変化しやすいことであり，2 つ目の性質は音高 k の

2.4 確率モデルとしての定式化
本節では，HTFD と同様に式 (20) で定義された Xl,m を

c 1959 Information Processing Society of Japan⃝

楽音はそれに対応する対数周波数 µk の周りに F0 が分布
しやすいことである．詳しい導出は省くが，この 2 つの4情報処理学会研究報告
IPSJ SIG Technical Report

性質は Ω の k 番目の行を転置した Ωk に対する確率分布

qlocal (Ωk ), qglobal (Ωk ) として以下のように記述できる．
qglobal (Ωk ) = N(Ωk ; µk 1 M , ξk2 I M ),(26)

qlocal (Ωk ) = N(Ω k ; 0 M , τ2k D−1 ), 1 −1 00−1 2 −1 0 0 −1 2 −1
D =  .
..
..
 ..
.
.
 0 . . . 0 −1
0 ... 0
0

(27)

...
...
...
..
.
2
−1

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

3. パラメータ推論アルゴリズム
観測パワースペクトログラム Y が与えられたときに，

2.4 節で構築した確率モデルに基づき事後確率 p(Θ|Y) ∝p(Y|Θ)p(Θ) を最大化するような w, Θ := {Ω, R, A} を求めた
い．したがって，目的関数

J(Θ) := ln p(Y|Θ) + ln p(Θ)
(28)

(33)

を最大化する問題として定式化できる．生成モデルを表す
項 ln p(Y|Θ)，事前分布を表す項 ln p(Θ) は w の定義域 W

N(Ωk ; µ, Σ) は平均 µ，分散 Σ をもつ K 次元正規分布であを用いて以下のように書ける．
∫ ∏ln p(Y|Θ) = ln
Pois(Yl,m ; Xl,m )p(w|β, Ω)dw

る．また，1 M は全要素が 1 の M 次元のベクトル，0 M は

M 次元の零ベクトル，I M は M × M の単位行列を表す．こln p(Θ) =ln p(Ωk ) + ln p(R) +k

れら 2 つの性質を同時に満たすような事前分布を設計する
ためは，Product-of-experts [10] と呼ばれるコンセプトを用

∑

W l,m

∑

ln p( Ak ).

(34)
(35)

k

式 (34) の右辺に，w に対する周辺化が入った対数関数があるため解析的に大域的最適解を導くことは困難であるが，いて，Ωk の事前分布をp(Ωk ) ∝ qlocal (Ωk )αlocal qglobal (Ωk )αglocal

この場合には補助関数法 [4, 7, 11] を用いることができる．

(29)

補助関数法は，目的関数値がちょうど上界になっているよとすればよい．αlocal , αglobal は，それぞれこの事前分布に対うな関数（補助関数と呼ぶ）を設計し，補助関数を補助変する qglobal (Ωk ), qglobal (Ωk ) の寄与を調節するハイパーパラ数と呼ぶ変数とパラメータ Θ 化することによって，J(Θ)

メータである．

を単調増加させる手法である．

U に関しても 2 つの性質を用いて事前分布を設計でき

対数関数は凹関数であるため，式 (34) に Jensen の不等る．ポピュラー音楽やクラシック音楽では調性があるため式が適用でき，以下のように下界が得られる．


∫
∑
p(w|β, Ω)  dw
ln p(Y|Θ) ≥q(w)  ln Pois(Yl,m ; Xl,m ) + ln
q(w) 
Wl,m

曲中の調や和音によって音高の出現に偏りがあり，この偏りは各音高の相対的なエネルギーに対する確率分布として記述できる．また，特定の音高のみが曲中演奏されるわけではなく，一つの音高に着目したときにはその音高をもつ音符の分布が時間的にスパースになりやすい．これも各音高のパワーの時間発展に対する確率分布として記述できる．これらの確率分布を簡単に導入するため，Uk,m を音高∑∑方向に正規化された音量 Rk = m Uk,m ( k Rk = 1) と時刻∑
方向に正規化された音量 Ak,m = Uk,m /Rk ( m Ak,m = 1) に分解する．すなわち，Uk,m = Rk Ak,m .

(30)とする．これによって上記の 2 つの確率分布が，それぞれ

R := [R0 , R1 , . . . , RK−1 ]⊤ , Ak := [Ak,0 , Ak,1 , . . . , Ak,M−1 ]⊤ の事
前分布として以下のように記述できる．
∏ γ(R) −1
R ∼ Dir(R; γ(R) ) :=Rk k ,

Ak ∼ Dir( Ak ; γ(A)k ) :=

(A)
γk,m
−1
Ak,m

(31)

γ

:=

(32)

⊤
[γ1(R) , . . . , γ(R)K ] はは R の事前分布のハイパーパ

ただし，=c は定数を除いて一致することを示している．こ∫こで導入された補助変数 q(w) は W q(w)dw = 1, q(w) > 0 を満たす. 等号成立条件は q(w) が w の事後分布一致するときである．すなわち，この補助関数は EM アルゴリズムの Q 関∫数であるとも言える．以下では，Eq(w) [w2 ] := W q(w)w2 dwと表記する．

(37) の括弧内の第 1 項は対数関数の中に和を含む関数型にる．これは以下のように式 (37) の下界を設計することによって解決できる．式 (37) の括弧内の第 1 項に Jensen の不等式を適用すると，ラメータであり，これに音高の出現頻度を反映できる．
(A)
(A) ⊤
γ(A)
:= [γk,1, . . . , γk,M
] は Ak の事前分布のハイパーパラ
k

メータであり，この値を小さくすればよりスパースな Akの推定値に誘導できる．

c 1959 Information Processing Society of Japan⃝(37)

なっており，解析的に解くのが困難となる原因となってい

m
(R)

(36)

式 (20) を見ると Xl,m は k, n に関する和を含むため，式

k

∏


∑
∑
=
q(w)  Yl,m ln Xl,m −
Xl,m
c W
l,m
l,m
)
p(w|β, Ω)
+ ln
dw
q(w)
∫

Yl,m ln Xl,m ≥ Yl,m∑k,nλk,n,l,m ln

w2k,n,m e−(xl −Ωk,m −ln n)22σ2λk,n,l,m

Uk,m

,

(38)

5

情報処理学会研究報告
IPSJ SIG Technical Reportとして下界が得られる．ここで，補助変数 λ := {λk,n,l,m }k,n,l,m∑は正の実数であり k,n λk,n,l,m = 1 を満たす．不等式 (38) の等号成立条件は

λk,n,l,m =w2k,n,m e−(xl −Ωk,m −ln n)2
2σ2

Uk,m

Xl,m

(39)

である．したがって，J(Θ) の補助関数 J + (λ, q(w), Θ) はJ + (λ, q(w), Θ)(xl −Ωk,m −ln n)2∑∑
w2k,n,m e− 2σ2Uk,m
λk,n,l,m ln=Eq(w)  Yl,mc
λk,n,l,m
l,m

k,n

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

(40)

∑

∑

(R)
n λk,n,l,m + γk − 1,
2m,n Ak,m Eq(w) [wk,m,n ]∑∑(A)l Yl,m n λk,n,l,m + γk,m − 1,
Ak,m ∝∑
Rk n Eq(w) [w2k,m,n ]Rk ∝l,m Yl,m

∑

(44)

(45)

と導出できる．ln p(w|β, Ω) に含まれる Ω を無視すれば，Ω
の更新式は

−1
∑
 αlocal
αglobalΩk ←  2 D +IM +diag( pk,n,l )
2
τk
ξk
n,l


∑
 αglobal× µk 2 1 M +(xl − ln n)pk,n,l  ,
ξk
n,l
[
]⊤
1
pk,n,l := 2 Yl,1 λk,n,l,1 , Yl,2 λk,n,l,2 , . . . , Yl,M λk,n,l,M
σ

(46)
(47)

と得られる．diag(p) は p の要素を対角成分に並べた対角行と書ける．

列である．この更新則では一部の項を無視しているため，

式 (40) の括弧内の第 2 項は Ωk,m に対して非線形ではあ∑るが， l Xl,m は X(x0 , t0 ), . . . , X(xL−1 , t M−1 ) を等間隔 ∆ x でサ∑ンプルした点であるため， x に関する積分を用いて l Xl,mJ ++ (λ, q(w), Θ) に対する単調増加性が保証されないことには以下のように精度良く近似することができる．
∫ ∞∑1
Xl,m ≃X(x, tm )dx∆ x −∞l∫ ∞ (x−Ω −ln n)2
k,m
1 ∑ 2
wk,n,m Rk Ak,me− 2σ2 dx=∆ x k,n
−∞
√
∑
∑
2πσ
=
Rk Ak,mw2k,n,m .
(41)
∆x
n
k

れた目的関数と β に依存する補助関数の項が同型であるた

注意されたい．
全極システムのパラメータ β に関しては，[13] で定義さめ，[13] で提案された手法を修正して用いることができる．
紙面の都合上詳細は省くが，以下の更新則を反復的に繰り返すことによって J ++ (λ, q(w), Θ) は単調増加する．

hk ← Ĉk (βk )βk ,(48)

ただし，Ck と Ĉk (βk ) は (P + 1) × (P + 1) の Toeplitz 行列であり，それらの (p, q) 番目の要素はそれぞれ

この近似は，式 (40) の括弧内の第 2 項が Ωk,m にほとんど依存しないことを意味している．

βk ← Ck−1 hk ,Ck,p,q =この近似を適用した J + (λ, q(w), Θ) を J ++ (λ, q(w), Θ) と2
1 ∑ Eq(w) [wk,m,n ]cos[(p − q)neΩk,m u0 ],MN m,n2ν2るため，λ は式 (39) にしたがって更新し，q(w) は w の事後

(49)
∑
1
1
Ĉk,p,q (βk ) =cos[(p − q)neΩk,m u0 ]MN m,n |Bk (ejneΩk,m u0 )|2

分布を求めれば良い．λ 内の w に関する積分のため λ を直(50)置き，これを用いて補助変数とパラメータに関する更新式を導出する．まず，補助変数の更新式は等号成立条件であ接計算することは難しいが，w2k,n,m を Eq(w) [w2k,n,m ] として近似的に計算する．詳しい計算は紙面の都合上省くが，q(w)と書ける．

4. モノラル音源分離性能の定量評価

の更新則は


∑

q(wk,m,n ) =Nakagami wk,m,n ;Yl,m λk,n,l,m + 1,∑ll Yl,m λk,n,l,m + 1√Ω
2πRk Ak,m σ/∆ x + |Bk (e jne k,m u0 )|2 /(2ν2 )全極スペクトル導入の効果を確認するため，モノラル音 (42)2aa 2a−1 −aζ 2 /bζe
Γ(a)baて HTFD と F0 を時不変とした HTFD （スペクトルテンプ
レートを正規分布状に拘束した場合の調波 NMF [14, 15] に対応．以後，Harmonic NMF）を用いた．分離には時間周∑波数マスクを Ck,l,m / k Ck,l,m と設計して用いた．ただし，となる [12]．ただし Nakagami(ζ; a, b) は仲上分布

Nakagami(ζ; a, b) =楽音響信号の各音高への分離実験を行った．比較手法とし(43)

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

w2 の推定値を Eq(w) [w2 ] として計算した．各音高ごとに録を表す．Θ に関する更新式は J ++ (λ, q(w), Θ) に対する各パ

音された演奏を用意するのは困難であったため，RWC クラ

ラメータでの偏微分が 0 となる値を求めることによって導

シック音楽データベース [16] の RM-C001 から RM-C005

出できる． R, A に関する更新式は，の最初の 30 秒を MIDI シンセサイザー FluidSynth [17] で

c 1959 Information Processing Society of Japan⃝6

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

築したのち，そのパラメータ推論法を示した．定量評価実験により，提案法が HTFD よりも高いモノラル音源分離性能を持つことを確認した．今後は，調やオンセットに関する事前情報を事前分布に導入し，多重音解析に対する効果を定量的に評価する予定である．
謝辞 本研究の一部は，JSPS 科研費 26730100 の助成を受けたものである．
参考文献
[1]
図3

F0 を固定した HTFD（Harmonic NMF）と HTFD，提案法（Proposed）で得られた各楽曲に対する平均 SNR 改善量と標準誤差．

[2]

合成したの音響信号（サンプリング周波数 16 kHz）を入[3]

力として用いた．スペクトログラムを得る際には，時間シフト間隔を 14.6 ms とし解析周波数を 55 Hz から 7040

Hz まで 10 cent 間隔でとった高速近似ウェーブレット変[4]

換 [18] を用いた．アナライジングウェーブレットは式 (6)[5]で定義された対数正規分布型のウェーブレットを用いた（σ = 0.02）．調波成分の個数は N = 20，音高数は K = 73とし，µk は A1（55 Hz）から A♯7（3322 Hz）までの音高に[6]半音単位で対応させた．Ωk に関するハイパーパラメータは(αlocal , αglobal , τk , vk ) = (1, 1, 1.0, 1.25) とし，R と Ak に関する

[7]

−4
ハイパーパラメータは γ(R) = 0.8×1K ，γ(A)k = (1−1.0×10 )1Iとした．全極システムの次数は P = 20 としハイパーパラ
メータは ν = 1 とした．また，Harmonic NMF は 100 イテ[8]

レーション，HTFD と提案法は 20 イテレーションとした．
図 3 に，分離音の signal-to-noise ratio（SNR）改善量の平[9]均と標準誤差を楽曲毎に示した．Harmonic NMF と HTFDで得られた SNR の差の平均値は約 0.02 dB であり，分離性能に大きな差は見られなかった．これは，MIDI で音源[10]を作成したため F0 の変動が少ないからであると考えられる．提案法は，全曲で Harmonic NMF と HTFD よりも平均

SNR 改善量が高く，Harmonic NMF と HTFD との SNR の[11]

平均的な差はそれぞれ約 0.80 dB，約 0.78 dB であった．この結果から，全極スペクトルの導入によるモノラル音源分[12]

離性能の向上が確認できる．本実験では全極スペクトルの[13]

パラメータも教師なしで推定したが，事前にこのパラメータを学習して分離に用いることもでき，その場合にはより[14]

高精度な分離が期待できる．この効果を確かめるため，今後は教師あり学習を用いて実音源に対する分離実験を行う
予定である．

[15]

5. 結論
本報告では，HTFD の多重音解析性能を向上させるために楽音の生成過程をソース・フィルタモデルを導入した多重音解析手法を提案した．ソース・フィルタモデルを全極

[16]
[17]
[18]

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).
Bregman, A. S.: Auditory scene analysis: The perceptual organization of sound, MIT press (1994).
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).
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).
亀岡弘和：全極型声道モデルと F0 パターン生成過程モデ
ルを内部にもつ統計的音声スペクトルモデル，信学技報，Vol. SP2010-74, pp. 29–34 (2010).
Kameoka, H. and Kashino, K.: Composite autoregressive
system for sparse source-filter representation of speech, Proc.
2009 IEEE International Symposium on Circuits and Systems, IEEE, pp. 2477–2480 (2009).
Virtanen, T. and Klapuri, A.: Analysis of polyphonic audio
using source-filter model and non-negative matrix factorization, Advances in Models for Acoust. Process., Neural Inf.
Process. Syst. Workshop (2006).
Hinton, G. E.: Training products of experts by minimizing
contrastive divergence, Neural Comput., Vol. 14, No. 8, pp.
1771–1800 (2002).
亀岡弘和：音声音響信号処理のための確率モデルと学習
アルゴリズム，第 17 回情報論的学習理論ワークショップ(2014).
亀岡弘和：音声生成過程の確率モデル，第 13 回情報論的
学習理論ワークショップ (2010).
El-Jaroudi, A. and Makhoul, J.: Discrete all-pole modeling,IEEE Trans. Signal Process., Vol. 39, No. 2, pp. 411–423(1991).
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
Process., pp. 109–112 (2008).
Goto, M.: Development of the RWC Music Database, Proc.
Int. Congress Acoust., pp. l–553–556 (2004).
: FluidSynth. http://www.fluidsynth.org/.
亀岡弘和，田原鉄也，西本卓也，嵯峨山茂樹：信号処理
方法及び装置. 特開 2008-281898, (20. Nov. 2008).

システムとして導入しパワースペクトログラムモデルを構

c 1959 Information Processing Society of Japan⃝

7