Vol.2012-MUS-96 No.8
2012/8/10

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

多重音基本周波数解析のための無限複合自己回帰モデル
吉井 和佳1,a)

後藤 真孝1,b)

概要：本稿では，音楽音響信号に含まれる複数の基本周波数 (F0) を推定するための確率的生成モデルについて述べる．従来，ソース・フィルタ理論に基づく非負値行列因子分解 (NMF) では，非負値の観測スペ
クトログラムを行列とみなし，二種類のスペクトル（ソースのスペクトルおよびフィルタの周波数応答）と，それらの各組み合わせに対する時間方向のアクティベーション（音量変化）との積に分解することが行われていた．このとき，楽器の発音機構は自己回帰 (AR) 系でモデル化できると仮定すれば，F0 をもつ
周期信号（くし型スペクトル）や白色雑音（平坦スペクトル）などの音源により，楽器の反響特性（音色）を表現する全極型フィルタが駆動されることで音響信号が生成される過程を考えることができる．しかし，この種の複合自己回帰モデルにおいては，本来未知であるはずのソースやフィルタの個数を事前に指定しておく必要があるという問題があった．この問題を解決するため，ガンマ過程を用いることで，上限のない個数のソースやフィルタを理論上許容可能なノンパラメトリックベイズモデルである無限複合自己回帰モデルを提案する．本モデルの学習のため，変分ベイズ法 (VB) と乗法更新則 (MU) を組み合わせた効率
的な反復最適化アルゴリズムを導出し，スパースな学習ができることを実験により確認した．

1. はじめにに増加させてしまうと，最適化アルゴリズムを用いて性質のよい局所解を求めるのが難しくなる．

多重音基本周波数解析は，自動採譜やコード認識などさ( 2 ) NMF を実行したあとで，各基底スペクトルから F0 を

まざまな音楽解析タスクの基礎をなす重要な技術である．

推定する処理が別途必要になる．しかし，基底スペク

近年，多重音の分離や採譜を行ううえで，非負値行列因子

トルの形状に制約を与えずに NMF を実行すると，調分解 (Nonnegative Matrix Factorization: NMF) に基づく波構造をもつ楽器音として不自然なものが基底として手法 [1–13] が多数提案されている．標準的な NMF [14,15]推定される可能性がある．このことは，F0 推定を難

は，時間・周波数平面上の振幅スペクトログラムあるいはしくするだけではなく，各基底に対して F0 が存在す

パワースペクトログラムを非負値行列とみなし，少数の基るかどうかの判定処理も追加で必要になる．

底スペクトルとそれらの各時刻における音量（アクティ( 3 ) 基底スペクトルの個数（モデルの複雑さ）は NMF にベーション）との積に分解する．このような行列の低ラン

よる分解結果に大きな影響を与えるにもかかわらず，

ク近似は，一見複雑に見える音楽音響信号であっても，構事前に手動で指定する必要がある．この問題は，第一成音の音高や音色の種類は限られているという事実から，の問題にも関連していることに注意されたい．単純にある程度妥当であると考えられている．また，NMF にはは，あらゆる複雑さのモデルを別々に試したあとで最乗法更新則 (Multiplicative Update: MU) と呼ばれる収束

善のものをひとつ選べばよいが，そのようなモデル選

保証付きの反復解法が存在し，基底と音量とを交互に更新

択は計算量の観点で現実的ではない．

することで局所最適解を効率的に求めることができる利点

上記で議論した通り，制約なしの NMF は柔軟すぎるうえ，がある．しかし，通常の NMF を多重音基本周波数解析に音楽的に意味のある分解結果を得られるよう手動で制御す利用するうえでの主要な問題として，以下の３点がある．

ることは困難であるといえる．しかし，これら３つの問題( 1 ) 基底スペクトル形状は時間変化しないため，楽器音のはこれまで部分的に扱われていただけであった [1–13]．

1
a)
b)

音色のバリエーションを十分に表現するためには，多統一的な枠組みのもとで上記３つの問題を解決するため，数の基底スペクトルが必要になる．例えば，同じ基本

本研究では無限複合自己回帰モデル (Inﬁnite Composite

周波数 (F0) であっても，スペクトル包絡の形状が異Autoregressive Model: iCARM) と呼ぶ音楽音響信号に対なる楽器音に対しては，それぞれ独立した基底スペク

する新しいベイズモデルを提案する（図 1）．本モデルのトルが必要になる．このようにモデルの複雑さを単純主要な特長は，以下の３点である．

産業技術総合研究所
Umezono 1-1-1, Tsukuba, Ibaraki 305–8568, Japank.yoshii(at)aist.go.jp
m.goto(at)aist.go.jp

c 2012 Information Processing Society of Japan( 1 ) ソース・フィルタ分解 [1]
基底スペクトルをソーススペクトルと全極型フィルタとの積にさらに分解する．これは，人間の発声機構を

1

Vol.2012-MUS-96 No.8
2012/8/10

情報処理学会研究報告
IPSJ SIG Technical Report
くし型スペクトル密度組み合わせ

全極型スペクトル密度⊗

IJ 個の基底…

…
白色雑音スペクトル⊗

IJ 個の音量⊕

J 個のフィルタ

I 個のソース
観測スペクトル

図 1

無限複合自己回帰モデル：I 個のソースと J 個のフィルタを組み合わせることで IJ 個の基底スペクトルが得られ，それぞれが各時刻において異なる音量でアクティベートされる．本研
究では，I および J を無限に発散させたときの極限を考える．

自己回帰系とみなすアプローチに着想を得ている．人

表 1 ソース・フィルタ理論に基づく NMF

基準
ソース #
フィルタ #
IS
I
AR
J
IS
H
I
MA
J
IS
(H)
I
J
KL
I
J
KL
H
I
J
KL
I∗
J
KL
H + N I∗
AR
J
Beta (0.5)
I
ARMA J ∗KL or IS H + N ∞AR∞(H：調波構造，N：ノイズ，-：その他，∗：時変化)

手法
Kameoka [1]Badeau [2]Durrieu [3]Virtanen [4]Carabias-Orti [5]Heittola [6]
安良岡 [7]Hennequin [8]本研究2. 関連研究
本章では，近年特に注目されている機械学習技術として，

NMF とノンパラメトリックベイズの関連研究を紹介する．
2.1 非負値行列因子分解

間は，同じ声の高さ（ソース）を維持しながら，声道のNMF は非負値行列に対してスパースな分解を行うため

形状（フィルタ）を変化させることで様々な母音を生の強力なアプローチである [14]．NMF は当初，与えられ成できる．楽器の発音機構も同様であるとみなせば，た顔画像の集合に対して，目や鼻といった少数の局所パーモデルの複雑さの増加を抑えつつ，F0 と音色というツへの分解を得る目的で利用されている．このような性質独立した観点から様々な楽器音を表現できる．

は，非負値の基底画像の加算しか許さないという制約から( 2 ) 調波構造モデリング [7, 9]

自然に導かれる．調波構造を持つ楽器音は局所的な周波数調波構造を持つ楽器音を表現するため，F0 をパラメー

帯域にエネルギーが集中しているため，NMF は音楽音響

タに持ち，周波数軸上で倍音を等間隔に並べたくし型

信号を局所パーツへ分解する目的に適している．

関数を用いてソーススペクトルを定式化する．このス

2.1.1 最適化のための基準

ペクトルが音響学的に妥当性のある全極型フィルタを

NMF を実行するうえでは，観測データ（スペクトログ

通過すると，各倍音の相対的な重みは，調波構造を持ラム）に対する再構成データ（基底スペクトルの線形結合）

つ楽器音として適切な値をとるようになる．

からの「乖離度」を測るための基準が必要である．良く利さらに，文献 [7, 16] における提案にしたがい，全周波用される基準として，カルバック・ライブラー (KL) ダイ

数帯域にわたって平坦なスペクトルをもつ白色雑音をバージェンス [14] や板倉・斎藤 (IS) ダイバージェンス [15]表す特別なソーススペクトルを導入する．これは，広が挙げられる（表 1）．両者を特別な場合として含むベータ

い周波数帯域にわたってなだらかなスペクトルをもつダイバージェンスに基づく NMF [17] も提案されている．

打楽器音（ギターの弦を指ではじいたり，ピアノのハ

音響信号のモデル化においては，IS-NMF の方が理論的

ンマーが弦を叩く際の打撃音なども含む）を生成する，KL-NMF の方がしばしには妥当であるが（3.2.1 節参照）

ための駆動音源である．調波構造を持つ楽器音と同様ば良い結果をもたらすことが知られている．この理由としに，音色は全極型フィルタによって特徴づけられる．

て，IS ダイバージェンスは凸性を持たないため，勾配降下( 3 ) ノンパラメトリックベイズ [12]型の最適化手法では，大域的な最適解を見つけることが困与えられた観測スペクトログラムを表現するのに必要難であることが考えられる．このような傾向は最尤推定になソースとフィルタの個数を自動的に調節するため，

おいて知られているが，ベイズ推定においても同様の傾向

ノンパラメトリックベイズモデルを構成する．本研究が見られるかは自明ではない．本研究では，各基準に基づでは，従来の複合自己回帰モデル [1] において，ソース

いて KL-iCARM と IS-iCARM を提案する．

とフィルタの個数を事前に指定することはせず，代わ

2.1.2 ソース・フィルタ分解

りにそれらを無限に発散させたときの極限を考える．

NMF の主要な拡張の一つに，ソース・フィルタ理論に基

このとき，無限次元の事前分布であるガンマ過程を導づくものがある（表 1）．例えば，亀岡ら [1] は，観測スペ

入することで，本来無限個あるうちで限られた個数のクトログラムを一定個数のソースと AR フィルタ（全極型

ソースとフィルタのみが実質的にアクティベートされ伝達関数）に分解するための複合自己回帰モデルを提案しるスパースな学習が実現できる．

ている．同様のモデルはいくつか提案されているが [3–5]，

本研究では，iCARM を学習するため，変分ベイズ法 (Vari-

フィルタは音響学的な視点に立って設計されたものではな

ational Bayes: VB) と乗法更新則 (MU) を組み合わせた新かった．一方，Badeau ら [2] は，移動平均 (MA) フィル

しいタイプの反復最適化アルゴリズムを提案する．

タ（全零型伝達関数）を用いることを提案している．

c 2012 Information Processing Society of Japan

2

Vol.2012-MUS-96 No.8
2012/8/10

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

いくつかの研究では，モデルの複雑さの増加と引き換えに，ソースやフィルタが時間変化することを許容するモデルを提案している．Heittola ら [6] や安良岡ら [7] は，フィルタは時不変のままであるが，ソースの時間変化を許容するモデルを提案している．一方，Hennequin ら [8] は，ソースを時不変として，自己回帰移動平均 (ARMA) フィルタの時間変化を許容するモデルを提案し，効率的な MU アル
ゴリズムに基づくパラメータ推定法を導出している．

2.1.3 調波構造モデリング
NMF の別の拡張として，基底スペクトルやソーススペ
クトルに対し，調波構造を持たせるような制約を与えるものがある．例えば，Vincent ら [10] は，ある周波数帯域内の調波構造スペクトルを表すテンプレートを複数用意し，各基底スペクトルをそれらの重み付き和で表現する手法を提案している．ソース・フィルタ型の NMF においては，

Badeau ら [2] は，各ソーススペクトルを，倍音に対応する
周波数で１を取りやすくなる確率的なバイナリベクトルで表現する手法を提案している．安良岡ら [7] や Hennequinら [9] は，各ソーススペクトルを，パラメトリックな関数（例：倍音周波数を中心として配置した急峻なガウス関数の和）で表現している．Carabias-Orti ら [5] は，倍音の重みをいくつかのテンプレートの重み付き和としてさらに分解している．文献 [5, 9] では，効率的な MU アルゴリズムに基づくパラメータ推定法を導出している．
安良岡ら [7] の手法の特徴は，全てのソースに調波構造を持たせるのではなく，平坦なスペクトルをもつ特別なソースをひとつ用意しておくことである．これは，人間の発声機構において，音源は周期信号（くし型スペクトル）か白
色雑音（平坦スペクトル）に大別でき，それら音源信号は，
声道（AR フィルタ）によって調音されるという仮定に基づいている．この仮定は楽器音のモデル化においても，ある程度妥当であると考えられる．本研究では，ベイズ的な枠組みのもとで楽器音の生成過程を記述する．

2.2 ノンパラメトリックベイズ
ノンパラメトリックベイズ理論 [18] は，古典的なベイズ理論を拡張したものである．従来のベイズ推定の枠組みでは，未知のパラメータに対して事前分布を導入することで，観測データが与えられたもとでの事後分布を求めることができる．しかし，本来は未知であるはずのモデルの複雑さ（本研究においてはソースやフィルタの個数）は，単なる超パラメータとして事前に指定してしておく必要があった．
このとき，網羅的に探索を行おうとしても，指定すべき複
雑さが複数種類ある場合には計算量が組み合わせ爆発を起
こしてしまう．ノンパラメトリックベイズモデルでは，事

表 2

M, N
I, J
Xmn
Ymn
θi
φj
Wim
Ajm
Hijn

無限複合自己回帰モデルで利用されるシンボル

周波数ビン数，フレーム数
ソース数，フィルタ数（いずれも無限大に発散）
ビン m・フレーム n における振幅またはパワー
ビン m・フレーム n における再構成された値
ソース i のグローバルな重みフィルタ j のグローバルな重み
ソース i におけるビン m の振幅またはパワーフィルタ j におけるビン m のゲイン
フレーム n におけるソース i とフィルタ j の音量きモデルに，Hoﬀman らによるノンパラメトリックベイズ
拡張された NMF であるガンマ過程 NMF (GaP-NMF) [12]が挙げられる．このモデルは，理論上無限個の基底スペク
トルが存在することを仮定するが，観測データに合わせて
少数の基底が実質的にアクティベートされる機構を有している．この機構は，各基底の時間変化を許容するような複雑なモデル [13] を構成する際にも利用されている．

3. 無限複合自己回帰モデル
本章では，無限複合自己回帰モデル (iCARM) と呼ぶ，
多重音基本周波数解析のためのノンパラメトリックベイズ
モデルについて述べる（図 1 および図 2）．本モデルは，IS
ダイバージェンス基準のもとで，パワースペクトログラムを指定された個数のソースと AR フィルタの積に分解する
複合自己回帰モデル [1] を拡張したものである．本研究では，KL ダイバージェンス基準のもとで，振幅スペクトログ
ラムを分解する別の定式化も行う（3.2.1 節参照）
．ここで，
文献 [7] に提案されている通り，各ソースに調波構造を持たせるため，くし型関数を用いることにする．ただし，平坦なスペクトル（ホワイトノイズ）を持ったソースを別に一つだけ用意しておく．さらに，文献 [12] に提案されている通り，ガンマ過程を用いることで，理論上無限個のソースとフィルタが存在することを許容するように拡張する．

3.1 スペクトログラムの分解
最初に，提案モデルの定式化に必要な数学的記号を定義しておく（表 2）．M を周波数ビン数，N をフレーム数として，M × N の複素スペクトログラムを行列 X とみなす．
また，I をソース数，J をフィルタ数とし，これらの値が無限に発散するときを考える．小文字の m, n, i, j はそれぞれインデックスを表すものとする．
本研究では，X の非負値表現（振幅あるいはパワースペ
クトログラム）を 3 種類の独立した因子 W ，A および Hに分解することを試みる．

|Xmn | or |Xmn |2 ≈I,J →∞



θi φj Wim Ajm Hijn

(1)

i,j

後分布推定時にそれらを一度に最適化することができる．

ここで，Wim ，Ajm および Hijn はそれぞれ，ソース i に近年，ベイズモデルは音楽情報処理においても盛んに利おけるビン m の振幅またはパワー，フィルタ j におけるビ

用され，ベイズ拡張された NMF [19] は音源分離タスクにおン m のゲイン，フレーム n におけるソース i とフィルタ j

いて良い結果を与えることが知られている．特に注目すべの音量を表す．また，2 種類の変数 θi および φj が導入さ

c 2012 Information Processing Society of Japan

3

Vol.2012-MUS-96 No.8
2012/8/10

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

れており，全フレームにわたるソース i あるいはフィルタに発散しても，θ や φ のうちで限られた要素のみがある程
度以上の大きさを持ち，その他の要素は無視できるくらい

ソース i

Wim

フィルタ jA jmに小さい値をとる（厳密にはゼロではない）ようにできる．
基底 ijすなわち，ソースやフィルタの「実質的な」個数 I + および J + がデータから自動的に決まることを示している．

σ2µi 2µi 3µi

Wim A jm⊗m＝

j のグローバルな重みを表す．このとき，I や J が無限大m

図 2 ソース W i とフィルタ Aj との組み合わせ

我々の目的は，観測データ X が与えられたもとで，未知のm確率変数の事後分布 p(θ, φ, H|X; W , A) を計算し，W お

して，θi >  となる要素の個数 I + はほとんど確実に有限よび A を構成するパラメータを推定することである．以降，であることが証明されている．実用上は，打ち切りレベル

事前分布 p(θ)，p(φ)，p(H)，尤度関数 p(X|θ, φ, H; W , A)，I を α より十分に大きくとれば，θ の I 個の要素のうちいパラメトリックな関数 W ，A の設計について議論する．

くつかがある程度大きくなるようなスパースな学習が実現

3.2 確率モデルの定式化できる．同様の議論は，φ に対しても適用できる．

本節では，KL ダイバージェンスあるいは IS ダイバー3.2.3 H に対するガンマチェイン事前分布
各ソース・フィルタの組み合わせの音量が時間方向にな

ジェンスに基づく iCARM の設計方法について述べる．

めらかに変化するようにするため，H に対してガンマチェ

3.2.1 X に対する尤度関数
本研究では，KL ダイバージェンスあるいは IS ダイバーイン事前分布 [20] をおく．具体的には，以下を仮定する．

Hij1 ∼ Gamma (β, β/d)

ジェンスを乖離度基準として用いる．いま，再構成データ
ij
ij
Y を Ymn = ij Ymn
および Ymn= θi φj Wim Ajm Hijn とGijn ∼ Gamma (β, βHijn−1 )定義すると，式 (1) から Xmn と Ymn との乖離度を最小化するような Ymn を求めたい．これは，Xmn が与えられた
ときに，Ymn をパラメータにもつポアソン分布あるいは指数分布の最尤推定を行うことと等価である．ここで，Xmnijに含まれるソース i とフィルタ j の組み合わせ成分を Xmnijとしておく．ただし，Xmn = ij Xmnが成り立つ．

まず，KL-iCARM では，振幅領域での加算性 |Xmn | =ij
ij |Xmn | が成り立つことを仮定する．これは数学的には正しくないが，実用上はうまく働くことが知られている．

Hijn ∼ Gamma (β, βGijn )(5)

ここで，β は，事前分布の強度（音量の時間変化のなめらかさ）を調節する超パラメータである．Gijn は，隣り合う音量 Hijn−1 および Hijn が正の相関を持つようにするための，両者と負の相関をもつ補助変数である−1（Eprior [Gijn ] = Hijn−1，Eprior [Hijn ] = G−1．このとき，ijn ）Gijn を積分消去すると，正の相関をもつマルコフ連鎖確率
H

(H

)β

−1
ijn−1 ijnp(Hijn |Hijn−1 ) = Γ(2β)2Γ(β) (Hijn−1 +Hijn )2β Hijn が得られる.

ij
ij
いま，|Xmn| ∼ Poisson(Ymn) を仮定すると，ポアソン分


ij
ij
布の再生性から ij |Xmn | ∼ Poisson( ij Ymn) が成立す3.2.4 W に対するくし型関数

るので，振幅 |Xmn | の尤度は次式で与えられる．

成されるくし型関数で表現する（図 2）．ここで，H は考|Xmn | ∼ Poisson (Ymn )(2)

一方，IS-iCARM では，複素領域での加算性 Xmn =ij
ij Xmn に基づいて，数学的には正確な定式化を行う．い

ij
ij
ま，Xmn
∼ Nc (0, Ymn) を仮定すると，複素ガウス分布の


ij
ij
再生性から ij Xmn ∼ Nc (0, ij Ymn) が成立する．したがって，パワー |Xmn |2 の尤度は次式で与えられる．

|Xmn |2 ∼ Exponential (Ymn )(3)

3.2.2 θ および φ に対するガンマ過程事前分布
無限次元の非負値ベクトルである θ および φ に対して

各ソーススペクトル W i を H 個のガウス関数の和で構慮すべき倍音の個数である．具体的には，以下を考える．
H
(m − hμi )2
(6)
Wim =
exp −
2σ 2
h=1

ここで，μi *1 は基本周波数を，σ は倍音周辺のエネルギーの拡散度合いを表す．また，白色雑音を表す特別なソースに対しては，全周波数帯域にわたって WIm = 1 とする．

3.2.5 A に対する全極型伝達関数
いま，あるフレームにおいて，ソース i とフィルタ j のij*2 が P 次の組み合わせで生成された信号 xij ≡ {xt }2M
t=1

自己回帰過程に従うと仮定すると，以下が成り立つ．

は，ガンマ過程事前分布をおくのが都合がよい．具体的には，θ および φ の各要素に対し，以下を仮定する．
α 
γ 
θi ∼ Gamma, α , φj ∼ Gamma
,γ
(4)
I
J
こうすると，I が無限大に近づくに連れて，ベクトル θ は，
形状パラメータ α を持つガンマ過程から生成される無限個の系列を良く近似するようになる．任意の実数  > 0 に対

c 2012 Information Processing Society of Japan


xij
t = −

P


i
ajp xij
t−p + st(7)

p=1

ここで，si ≡ {sit }2M
t=1 はソース i に対応する信号であり，

aj ≡ {aj0 , · · · , ajp }T はフィルタ j の係数である（aj0 = 1）．
*1*2基本周波数を μ̃i [Hz] とすると，μi = μ̃i /(r/2M ) [bins] で変換できる．ここで，r は音響信号のサンプリング周波数である．
2M はフーリエ変換に用いるサンプル数（窓幅）である．

4

Vol.2012-MUS-96 No.8
2012/8/10

情報処理学会研究報告
IPSJ SIG Technical Reporti2M
iwi ≡ {wti }2M
t=1 を s の自己相関関数，{Wim }m=1 を w の離また，W および A（μ，σ および a）を最適化するため，

散フーリエ変換 (DFT) で得られる複素スペクトル（ウィーVB に乗法更新則 (MU) [8, 9] を組み合わせる手法を提案

ナー・ヒンチンの定理から si のパワースペクトルと等価）する．一般的に，MU における更新式は，最小化すべきコ

ij 2Mとする．また，{Xm}m=1 を xij の DFT で得られる複素スiスト関数（本研究では −L）を偏微分することで得られる．

ペクトルとする．このとき，音源信号 s が定常なガウス性

例えば，基本周波数 μi の更新式は，コスト関数の偏微分がij
雑音であると仮定すると，Xmは複素ガウス分布 Nc (0, Σijm)

２つの正の項の差 −∂L∂μi = Gμi − Fμi として記述できるとに従う．ここで，Σij
m = Wim Ajm であり，Ajm は

1

1
Ajm = 
2 = T P
aU
j −2π m pi 
m aj
j
2M
 p=0 ap e


i

(8)

m
行列で，[U m ]pq = cos(2π 2M(p − q)) で与えられる．した
ij 2がって，|Xm| は，Wim Ajm をパラメータにもつ指数分布に従うことが分かる．すなわち，xij が与えられたもとで
ij 2 Mの aj の最尤推定は，{|Xm| }m=1 と {Wim Ajm }M
m=1 の間の IS ダイバージェンスを最小化することと等価である*3 ．
これまでの議論から，KL ダイバージェンスあるいは IS
ダイバージェンスに基づく iCARM では次式を仮定する．

Ajm =

1
T
aj U m ajor1
Ajm = T
aj U m aj
IS(9)

ここで，KL-iCARM の場合に平方根を計算する理由は，パワースペクトログラムではなく振幅スペクトログラムを分解するからである（3.2.1 節参照）．

3.3 確率モデルの学習
我々の目的は，対数周辺尤度 log p(X; W , A) が最大化されるように，事後分布 p(θ, φ, H|X; W , A) と関数 W お
よび A （パラメータ μ，σ ，a）を最適化することである．
しかし，対数周辺尤度や真の事後分布は解析的に計算することができないため，本研究では変分ベイズ法 (VB) を用
いて近似計算を行う．まず，因子分解が可能なように限定された関数形をもつ事後分布を考える．

q(θ, φ, H) =q(θi )iq(φj )jq(Hijn )分がゼロであれば，μi は変化せず，そうでなければ勾配とは逆の方向に更新されるため，コスト関数は減少する．

である（図 2）
．U m は (P + 1) × (P + 1) の Toeplitz 型のKLF

すると，μi ← μi × Gμμi で与えられる．したがって，偏微3.3.1 KL-iCARM における変分推論
本節では，KL-iCARM における変分事後分布の更新式を導く．まず，式 (11) に含まれる対数尤度の期待値

E[log p(X|θ, φ, H; W , A)] は次式で与えられる．
E[log p(X|θ, φ, H; W , A)]
log Poisson(|Xmn ||Ymn )
=
mn

=



|Xmn |E log

mn

−





θi φj Wim Ajm Hijnij E θi φj Wim Ajm Hijn −log |Xmn |! (13)

mn ij

mn

しかし，第一項の「和の対数関数の期待値」は解析的に計算できない．そこで，対数関数 log(x) は凹関数であることから，λk ≥ 0 かつ k λk = 1 を満たす任意の補助変数

λ = {λ1 , · · · , λK } に対してイェンセンの不等式
 x  
 
x 
k
k
≥
(14)
log
xk = log
λk
λk logλkλkkkが成り立つことを用いると，式 (13) の変分下限はE[log p(X|θ, φ, H; W , A)]θi φj Wim Ajm Hijn≥
|Xmn |
λmnij E log
λmnijmn  ij E θi φj Wim Ajm Hijn −log |Xmn |! (15)
−
mn ij

(10)

k

mn


で与えられる．ここで， ij λmnij = 1 となる補助変数 λを導入した．第一項は「対数関数の和の期待値」に変換さijn

次に，この変分事後分布と真の事後分布との KL ダイバー

れ，解析的に計算可能になっている．

ジェンスを最小化するように反復最適化を行う．これは，式 (15) の十分統計量は log(x) および x であり，式 (4),(5)対数周辺尤度の変分下限 L を最大化することと等価であで与えられる θ, φ, H に対するガンマ事前分布と共役性がる*4 ．ここで，L および更新式は次式で与えられる．

成立していることから，変分事後分布をq(θi ) = Gamma(aθi , bθi ),log p(X; W , A) ≥ E[log p(X|θ, φ, H; W , A)]H
q(Hijn ) = Gamma(aHijn , bijn )+ E[log p(θ)] + E[log p(φ)] + E[log p(H)]− E[log q(θ)] − E[log q(φ)] − E[log q(H)] ≡ L(11)q(θ) ∝ exp(Eq(φ,H ) [log p(X, θ, φ, H; W , A)])q(φ) ∝ exp(Eq(θ ,H ) [log p(X, θ, φ, H; W , A)])q(H) ∝ exp(Eq(θ ,φ) [log p(X, θ, φ, H; W , A)])
*3

(12)

線形予測分析 (Linear Predictive Coding: LPC) においては，
一般的に，音源信号 si は白色雑音であることが仮定されている．本研究では，周期信号の場合を主に考えている．
（Wim = 1）*4 より正確にいえば，L は依然として解析的に計算できないので，
さらなる変分下限を計算して，それを最大化することになる．

c 2012 Information Processing Society of Japanq(φj ) = Gamma(aφj , bφj )(16)と仮定すると，変分パラメータは次式で求められる．
α aθi = + mnj |Xmn |λmnijIbθi = α + mnj E[φj Wim Ajm Hijn ]γ 
aφj = + mni |Xmn |λmnijJbφj = γ + mni E[θi Wim Ajm Hijn ]

(17)
aH
ijn = 2β +m |Xmn |λmnijHbijn = βE[Gijn +Gijn+1 ] + m E[θi φj Wim Ajm ]5

Vol.2012-MUS-96 No.8
2012/8/10

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

また，式 (15) を微分してゼロとおくことで，変分下限を

2.0

6.0

x が大きければ0.5程度最大化する λmnij は次式で与えられる．

λmnij ∝ exp(E[log(θi φj Wim Ajm Hijn )])

1.5

(18)

すなわち，λmnij は，ビン m・フレーム n において，ソース i・フィルタ j の組み合わせに対応する成分 θi φj Wim Ajm Hijn

y = exp(log( x ))E[θi φj Wim Ajm Hijn ] を計算する場合とは異なり，小さな成分は消滅に向かう効果が得られる．例えば，θ に着目すると，

E[θi ] = exp(log(aθi ))/bθi ，exp(E[log θi ]) = exp(ψ(aθi ))/bθi

4.0

a=2

a が１減ると

a =1

ディスカウント

1.0
2.0
0.5

が占める割合の期待値を表している．ただし，式 (18) の右辺
は θi φj Wim Ajm Hijn の「幾何平均」であるので，通常の平均

K a +1 (2 x)x
K a (2 x)y=のディスカウント

0y = exp(ψ ( x ))x が小さいとほぼゼロに
0.5

1.0

1.5

2.0

a=0
a = −1

0

0.5

1.0

1.5

2.0

図 3 ディガンマ関数 ψ(x) や第二種変形ベッセル関数 Ka (x) の性質に起因するディスカウント効果

であり，後者は前者からほぼ一定のディスカウントがなさ

いま，式 (22) の十分統計量は log(x) および x1 であり，

れ，aθi が小さいほど影響が大きいことが分かる（図 3 左）．

θ, φ, H に対するガンマ事前分布の十分統計量は log(x) お

これは NMF のベイズ推定に特有の現象であり，通常の平よび x であるので，共役性が満たされていない．しかし，均を計算する最尤推定や事後確率最大化 (MAP) 推定では，log(x)，x および x1 を十分統計量にもつ一般化逆ガウス

このようなスパース化効果を得ることはできない．

(GIG) 分布（付録参照）を変分事後分布に用いることで，3.3.2 IS-iCARM における変分推論

閉形式の更新式導出が可能になる [12]．具体的には，本 節 で は ，IS-iCARM に お け る 変 分 事 後 分 布 の 更 新式を導く．まず，式 (11) に含まれる対数尤度の期待値

E[log p(X|θ, φ, H; W , A)] は次式で与えられる．
E[log p(X|θ, φ, H; W , A)]
log Exponential(|Xmn |2 |Ymn )
=
mn

=−


mn

−


mn

1
θ
φ
W
i
j
im Ajm Hijn
ij

|Xmn |2 E 
E log



θi φj Wim Ajm Hijn

(19)

ij

しかし，第一項の「和の逆数の期待値」も第二項の「和の対数関数の期待値」も解析的に計算できないため，文献 [12, 21] で提案されている手法を用いる．まず，第一項の下限関数について考える．関数 − x1 は凹関数であることから，ηk ≥ 0 かつ k ηk = 1 を満たす任意の補助変数

η = {η1 , · · · , ηK } に対してイェンセンの不等式


−1
1
1
−1

= 
ηk xk = −ηk2
(20)
xk ≥ −xk
k xk
k ηk ηkηkkkが成り立つことを利用する．次に，第二項の下限関数について考える．式 (13) と異なり，負の対数関数 − log(x) は凸関数であることから，ω 周りのテイラー展開を用いて
x− log(x) ≥ − log(ω) + 1 −(21)
ωとして ω における接線を得る．式 (19) の変分下限はE[log p(X|θ, φ, H; W , A)]


1
2
(22)
|Xmn |2
ηmnij
E
≥−
θ
φ
W
Ajm Hijn
i
j
im
mn ⎛
ij
⎞



⎝log(ξmn ) − 1 + 1−
E θi φj Wim Ajm Hijn ⎠ξmn ijmn
で与えられる．ここで， ij ηmnij = 1 となる補助変数 η
および任意の正の実数 ξmnij を導入した．この結果，全ての項が解析的に計算可能となっている．

c 2012 Information Processing Society of Japanq(θi ) = GIG(aθi , bθi , cθi ),q(φj ) = GIG(aφj , bφj , cφj )HHq(Hijn ) = GIG(aHijn , bijn , cijn )(23)と仮定すると，変分パラメータは次式で求められる．
αE[φj Wim Ajm Hijn ]aθi = , bθi = α + mnj
ξmn
I


θ
2 2
ci = mnj |Xmn | ηmnij E φj Wim A1jm Hijn
γ
E[θ W A H]
aφj = , bφj = γ + mni i imξmnjm ijn
J2
cφj = mni |Xmn |2 ηmnij
E θi Wim A1jm Hijn


H
1
2 2
aH
ijn = 2β, cijn =m |Xmn | ηmnij E θi φj Wim Ajm
 E[θi φj Wim Ajm ]bH(24)
ijn = βE[Gijn +Gijn+1 ] +mξmn
また，式 (22) を微分してゼロとおくことで，変分下限を最大化する ηmnij および ξmn は次式で与えられる．

ηmnij ∝ E[ θi φj Wim1Ajm Hijn ]−1ξmn = ij E[θi φj Wim Ajm Hijn ]

(25)
(26)

2
式 (22) において，ηmnij
は，KL-iCARM における λmnij と同様に，ビン m・フレーム n におけるパワー |Xmn |2 を各成分に分配する役割を果たしている．ただし，式 (25) の右辺
は θi φj Wim Ajm Hijn の「逆数の平均の逆数」であり，やはり単純な平均 E[θi φj Wim Ajm Hijn ] に比べてディスカウント

θに効果が発生し，スパースな解へと誘導される．例えば，√√√√
Ka+1 (2 bc) c
Ka (2 bc) c√ √ ，E[ 1 ]−1 =√ √θiKa (2 bc) b
Ka−1 (2 bc) b（添え字 θi は省略）であるので，常に E[θi ] より E[ θ1i ]−1 の方が小さくなる（図 3 右）．特にガンマ過程の場合は aθi → 0であるので，cθi がゼロに近いほど影響が大きい．

着目すると，E[θi ] =最後に，式 (5) の事前分布で導入された補助変数 G に関して考える．変分事後分布としてガンマ分布 q(Gijn ) =GGamma(aGijn , bijn ) を仮定すると，変分パラメータはaG
ijn = 2β,bG
ijn = βE[Hijn−1 +Hijn ](27)で計算できる．これは KL-iCARM でも同様である．

6

Vol.2012-MUS-96 No.8
2012/8/10

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

3.3.3 KL-および IS-iCARM における乗法更新則
ソース W およびフィルタ A を構成するパラメータ μ，σおよび a は，乗法更新則を用いて最適化できる．具体的にるには小さすぎるが，ピアノの音色のみを扱う本実験においては十分である．基本周波数 {μi }88
i=1 は，標準的なピアノの 88 鍵に対応するように初期化した．

は，式 (15) あるいは式 (22) で与えられる対数尤度の期待値

複数の基本周波数は，閾値処理により各フレームごとに（または対数周辺尤度 log p(X; W , A)）の変分下限の符号検出した．もし，フレーム n においてソース i のゲインj θi φj Hijn が閾値よりも大きければ，そのフレームにはを反転したものをコスト関数とみなせば，μi ← G−1
μi Fμi μi ，−122
σ 2 ← G−1
σ 2 Fσ σ ，aj ← Gaj F aj aj とすると，コスト関数を減少させることができる．このとき，各項は


2

i)
F
G
Fμi = mnjh h(mVmnij+ hμi Vmnij) exp − (m−hμ
2
2σ



(m−hμi )2GF
Gμi = mnjh h(mVmnij + hμi Vmnij ) exp − 2σ2


2

i)
F
(m − hμi )2 exp − (m−hμFσ2 = mnijh Vmnij
2σ 2


2

i)
G
Gσ2 = mnijh Vmnij(28)(m − hμi )2 exp − (m−hμ2σ 2

F 値の平均が最大化されるように，全体でひとつ定めた．
4.2 実験結果
最初に，提案手法の動作確認のため，楽曲「alb se2」の冒頭 4.9 秒に対する分解を試みた．この音響信号では，D4,

C#4, C4, A3, F#3 の順にピアノ音が一つずつ重畳し，最終的に五重音を構成する．図 4 に示すように，KL-iCARM は
適切な個数のソース（5 個の調波構造と 1 個の白色雑音）をで与えられる．さらに，KL-iCARM の場合は，
−1
F
G
Vmnij
= E[θi φj Ajm Hijn ], Vmnij= |Xmn |λmnij Wim
F aj = mni θi φj Wim Hijn A3jm U m
Gaj = mni |Xmn |λmnij A2jm U m(29)となる．一方，IS-iCARM の場合は，
E[θi φj Ajm Hijn ]
F
Vmnij
=
ξmn


G
2 2
Vmnij = |Xmn | ηmnij E θi φj W 2 1Ajm HijnimE[θ φ Wim Hijn ] 2
Ajm U mF aj = mni i j ξmn2
Gaj = mni |Xmn |2 ηmnij
E θi φj W1im Hijn U m

μi が含まれていると判定した．閾値は，30 曲全体に対する見つけることができている．また，ビン m・フレーム n におi
ける各ソース i の成分 E[Ymn] = j E[θi φj Wim Ajm Hijn ]を計算することで，観測データ X を調波成分と雑音成
分に分離することができる．一方，図 5 に示すように，

IS-iCARM ではソースやフィルタの個数は過剰に見積もられ，オクターブエラー（正解の F0 の半分を検出）が多く発生する傾向があった．この理由として，IS ダイバージェ
ンスでは，再構成データ Ymn が観測データ Xmn を超えてもコスト関数があまり増加しないことが考えられる．

(30)

jとなる．ただし，aj ← G−1
aj F aj aj で更新すると，a0 = 1

は満たされなくなる．そのため，aj0 = 1 となるように各

次に，実験データ全体を評価に用いた．KL-iCARM お
よび IS-iCARM で得られた F 値は，全 30 曲平均でそれぞれ 48.4%および 35.1%であった．この結果は文献 [5, 11] で報告されている最新の成果には及ばないが，無限の複雑さ

フィルタ Aj をスケーリングする処理を導入した．これにを持つ空間においてスパースな学習を行うことができる性

より，コスト関数の単調減少性は保証されなくなるが，収質から，iCARM は依然有望であると考えている．現状の束性に実用上問題はないことを実験により確認している．

モデルの制約として，適切な個数の基底が推定できたとし

4. 評価実験ても，各時刻におけるオン・オフを判定するために閾値処理が必要なことが挙げられる．今後は，オン・オフに対応

本章では，KL ダイバージェンスあるいは IS ダイバージェ
ンスに基づく iCARM（KL-iCARM および IS-iCARM）の性能を評価するために行った比較実験について報告する．

4.1 実験条件

するバイナリ変数を導入することを検討したい．

5. おわりに本研究では，無限複合自己回帰モデル (iCARM) とよぶ

実験データは文献 [5,11] と同様であり，MAPS ピアノデー

ノンパラメトリックベイズモデルを提案し，観測スペクト

タベース [22] に含まれる「ENSTDkCl」サブセット 30 曲か

ログラムを３つの因子，すなわち，ソース，フィルタ，それら冒頭 30 秒を抜き出したものを用いた．本研究では，44.1

らの組み合わせの音量の積に分解することを試みた．実験[kHz]・ステレオのオリジナル音源を 16 [kHz]・モノラルに変の結果，ガンマ過程を用いることで，与えられたデータに換し，窓幅 2048 点，シフト長 160 点の短時間フーリエ変換

対して適切な個数のソースやフィルタを推定することがで(STFT) を行うことで各曲のスペクトログラムを得た．このきることが分かった．今後は，多重音基本周波数推定の精結果，M = 1024 かつ N = 3000 となった．振幅あるいはパ
ワースペクトログラムは，各曲ごとに M1N mn |Xmn | = 1

度を上げるため，ウェーブレット変換などで得られる対数あるいは maxmn |Xmn | = 1 となるように正規化した．超

また，フィルタ（スペクトル包絡）は楽器音の音色を表し

パラメータは，I = 88 + 1, J = 10, α = 1, β = γ = 0.1,ていると考えられるので，さらに，楽器パートの音源分離

H = 20, P = 4, d = Eemp [|Xmn |] or Eemp [|Xmn |2 ] としにも取り組んでいきたい．ソースやフィルタを構成するパ

た．最大フィルタ数 J = 10 はガンマ過程を正確に近似す

ラメータのベイズ的な取り扱いも今後の課題である．

2

c 2012 Information Processing Society of Japan

周波数軸をもつスペクトログラムの利用を検討している．

7

Vol.2012-MUS-96 No.8
2012/8/10

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

観測データ X

再構成データ Y

観測データ X

再構成データ Y

Y 中の調波成分

Y 中の雑音成分

Y 中の調波成分

Y 中の雑音成分

ソースの重み E[θ]優位なフィルタ A ([dB])

ソースの重み E[θ]優位なフィルタ A ([dB])

C#4
雑音
D4
F#3
A3 C4

必要なソースの個数が過大評価されている

小さな重みを持つ
不要なソースフィルタの重み E[φ]フィルタの重み E[φ]必要なフィルタの個数が過大評価されている

非常に小さな重みを持つ不要なフィルタ

図 4 KL-iCARM による音源分離結果

図 5

謝辞: 本研究の一部は，JSPS 科研費 23700184 および

JST OngaCREST プロジェクトの支援を受けた．
参考文献
[1]
[2]

[3]

[4]
[5]

[6]

[7]
[8]

[9]
[10]
[11]

H. Kameoka and K. Kashino. Composite autoregressive
system for sparse source-ﬁlter representation of speech.
ICASSP, pp. 2477–2480, 2009.
R. Badeau, V. Emiya, and B. David. Expectationmaximization algorithm for multi-pitch estimation andseparation of overlapping harmonic spectra. ICASSP,
pp. 3073–3076, 2009.
J.-L. Durrieu, G. Richard, B. David, and C. Févotte.
Source/Filter model for unsupervised main melody extraction from polyphonic audio signals. IEEE Trans.
on ASLP, 18(3):564–575, 2010.
T. Virtanen and A. Klapuri. Analysis of polyphonicaudio using source-ﬁlter model and non-negative matrix factorization. NIPS Workshop, 2009.
J. J. Carabias-Orti, T. Virtanen, P. Vera-Candeas,
N. Ruiz-Reyes, and F. J. Cañadas-Quesada. Musical instrument sound multi-excitation model for non-negative
spectrogram factorization. IEEE J. of Sel. Top. in Sig.
Proc., 5(6):1144–1158, 2011.
T. Heittola, A. Klapuri, and T. Virtanen. Musical instrument recognition in polyphonic audio using sourceﬁlter model for sound separation. ISMIR, pp. 327–332,
2009.
安良岡直希, 奥乃博. 調波・非調波・音色構造因子分解による音響信号分析と音源分離インターフェースへの応用. 情処研報, 2012-MUS-94, pp. 1–8, 2012.
R. Hennequin, R. Badeau, and B. David. NMF with
time-frequency activations to model nonstationary audio events. IEEE Trans. on ASLP, 19(4):744–753,
2011.
R. Hennequin, R. Badeau, and B. David. Timedependent parametric and harmonic templates in nonnegative matrix factorization. DAFx, pp. 1–8, 2010.
E. Vincent, N. Bertin, and R. Badeau. Adaptive harmonic spectral decomposition for multiple pitch estimation. IEEE Trans. on ASLP, 18(3):528–537, 2010.
N. Bertin, R. Badeau, and E. Vincent. Enforcing har-c 2012 Information Processing Society of Japan


[12]
[13]
[14]
[15]

[16]
[17]
[18]
[19]
[20]
[21]
[22]

IS-iCARM による音源分離結果

monicity and smoothness in Bayesian non-negative matrix factorization applied to polyphonic music transcription. IEEE Trans. on ASLP, 18(3):538–549, 2010.
M. Hoﬀman, D. Blei, and P. Cook. Bayesian nonparametric matrix factorization for recorded music. ICML,2010.
M. Nakano et al. Bayesian nonparametric spectrogram modeling based on inﬁnite factorial inﬁnite hidden
Markov model. WASPAA, pp. 325–328, 2011.
D. Lee and H. Seung. Algorithms for non-negative matrix factorization. NIPS, pp. 556–562, 2000.
C. Févotte, N. Bertin, and J.-L. Durrieu. Nonnegative
matrix factorization with the Itakura-Saito divergence:With application to music analysis. NECO, 21(3):793–
830, 2009.
B. Niedermayer. Improving accuracy of polyphonicmusic-to-score alignment. ISMIR, pp. 585–590, 2009.
C. Févotte and J. Idier. Algorithms for nonnegative
matrix factorization with the beta-divergence. NECO,
23(9):2421–2456, 2011.
P. Orbanz and Y. W. Teh. Bayesian nonparametricmodels. Encyclopedia of Machine Learning. Springer,
2010.
A. T. Cemgil. Bayesian inference for nonnegative matrix factorisation models. Computational Intelligence
and Neuroscience, 2009:Article ID 785152, 2009.
A. T. Cemgil and O. Dikmen. Conjugate gamma
Markov random ﬁelds for modelling nonstationary
sources. ICA, 2007.
亀岡弘和, 後藤真孝, 嵯峨山茂樹. スペクトル制御エンベ
ロープによる混合音中の周期および非周期成分の選択的イコライザ. 情処研報, 2006-MUS-66, pp. 77–84, 2006.
V. Emiya, R. Badeau, and B. David. Multipitch estimation of piano sounds using a new probabilistic spectral smoothness principle. IEEE Trans. on ASLP,
18(6):1643–1654, 2010.

付

Gamma(x|a, b) =ba a−1 −bx
e
,
x
Γ(a)

録
Exponential(x|λ) =

1 −x
e λ
λ

a

Poisson(x|λ) =λx −λ(b/c) 2 xa−1 −(bx+ c )x√ee , GIG(x|a, b, c) =x!
2Ka (2 bc)

8