Vol.2014-MUS-104 No.9
2014/8/26

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

多重音基本周波数解析のための無限重畳離散全極型モデル
吉井 和佳1,a)

糸山 克寿1,b)

後藤 真孝2,c)

概要：本稿では，多数の楽器音が重畳している音楽音響信号を，音の三要素である音高（基本周波数）
・音色（スペクトル包絡）・音量に分解するための確率的ソース・フィルタモデルについて述べる．ソース・
フィルタ理論は楽器音分析に広く利用されており，楽器音のフーリエ変換スペクトルは，音源信号の基本周波数に起因するスペクトル微細構造と楽器音の音色を表すスペクトル包絡との積に分解される．このとき，スペクトル包絡が全極型モデルで表現できると仮定すると，理論的には線形予測分析 (LPC) を用いて，線形周波数領域でスペクトル包絡を推定することができる．しかし，実際には，調波構造のピークのみがスペクトル包絡からの信頼できるサンプルであるとみなせるため，スペクトル包絡推定に全周波数帯域を利用することは適切ではない．この問題の解決法のひとつに離散全極型モデルが知られているが，多重音に対して適用することはできなかった．本研究では，離散全極型モデルを LPC の多重音拡張である
複合自己回帰モデルの枠組みに組み入れることで，調波構造が複数重畳した音響信号を扱うことができる
無限重畳離散全極型モデルを提案する．本モデルは，人間の聴覚特性に則した対数周波数領域で定式化されるノンパラメトリックベイズモデルであり，適切な個数のスペクトル包絡とそこからサンプルされた適切な個数の調波構造を推定することができる．実験の結果，提案手法の有効性を確認した．

1. はじめにす微細構造（調波構造）と楽器個体の共鳴特性を表すスペ
クトル包絡によってよく特徴づけられる．人間の聴覚系は音楽情報処理分野において，機械学習技術に基づく音楽

スペクトルのピーク（フォルマント）に対して敏感である

音響信号の統計的モデリングは，ホットなトピックのひとつ

ため，楽器音の各時間フレームにおけるスペクトル包絡を，である．特に，多重音に対する基本周波数推定や音源分離

全極型周波数伝達関数（自己回帰フィルタの周波数応答）における有用性から，非負値行列分解 (Nonnegative Matrixを用いて表現することが一般的である [2]．全極型スペク

Factorization: NMF) は大きな注目を集めている [1–13]．

トル包絡推定の古典的な方法である線形予測分析 (Linear

標準的な NMF では，多重音の振幅あるいはパワースペク

Predictive Coding: LPC) [16] は，音源信号がガウス性白トログラム（非負値行列）を二つの非負値行列，すなわち，

色雑音であるという強い仮定のもとで，ある特定の確率モ

周波数方向の基底スペクトルの集合と各基底スペクトルにデルの最尤推定に対応していることが知られている．

対応する時間方向の音量変化の集合とに分解することがで

LPC の問題点を解決する有望な統計的アプローチとしきる．このとき，観測スペクトログラムと再構成スペクトて，複合自己回帰モデル (Composite Autoregressive Model:

ログラムとの近似誤差を表すコスト関数を最小化するた

CAR) [5] と呼ばれるソース・フィルタ NMF が提案され

め，効率的な乗法更新アルゴリズムが利用できる [14]．こている．本モデルでは，観測音響信号のスペクトログラムのアルゴリズムは，ある特定の確率モデルの最尤推定と等は，複数個の微細構造（ソース）と複数個のスペクトル包価であることがあとになって示されている [15]．

絡（フィルタ）との組み合わせから構成されているとみな

近年，音声信号の生成過程を説明する目的で考案された

す．このアプローチの重要な特徴は，音源信号がガウス性

ソース・フィルタ理論を楽器音の統計的モデリングに援用

白色雑音であるとは仮定されておらず，全極型スペクトルすることがしばしば行われている [2–4]．周波数領域におい包絡推定と同時に音源信号のスペクトル自体が推定されるて，楽器音の音高と音色はそれぞれ，音源信号の性質を表点である．本モデルは確率的な解釈が可能であるため，ノンパラメトリックベイズ理論を用いて，適切な個数のソー

1
2
a)
b)
c)

京都大学 大学院情報学研究科 知能情報学専攻
Yoshida-honmachi, Sakyo-ku, Kyoto 606-8501, Japan
産業技術総合研究所 情報技術研究部門
1-1-1 Umezono, Tsukuba, Ibaraki 305-8568, Japan
yoshii(at)i.kyoto-u.ac.jp
itoyama(at)kuis.kyoto-u.ac.jp
m.goto(at)aist.go.jp

c 2014 Information Processing Society of Japanスとフィルタを推定できる無限複合自己回帰モデルに拡張できる [6]．また，音源信号のスペクトルが調波構造をもつように制約を加えることも可能であり，音源信号の基本
周波数 (F0) を最尤推定することができる [6]．

1Vol.2014-MUS-104 No.9
2014/8/26

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

複合自己回帰モデルを含む従来のソース・フィルタ NMFの主要な問題点として，以下の二点が挙げられる．

離散全極型モデル
倍音成分の重みをフィルタに従うように調節
全極型スペクトル包絡

調波構造

⊗

• 調波構造のピークのみがスペクトル包絡からの信頼できるサンプルとみなせるにもかかわらず，調波構造に基底⊗対してスペクトル包絡を推定する際に，すべての周波数帯域が等しく考慮に入れられていた．

音量

対数周波数 [cents]

ソース (音高)で得られる線形周波数領域では，ピアノ音やギター音に見られるような，高次倍音の非調波性（F0 の整数倍からのずれ）が，等間隔に配置された倍音成分からなる理想的な調波構造モデルと実際の観測スペクトルとの致命的な不整合を起こしていた．

…

性に則しているにもかかわらず，短時間フーリエ変換(STFT) を用いるのが一般的であった．フーリエ変換

対数周波数 [cents]⊕…

• ウェーブレット変換や定 Q 変換の方が人間の聴覚特フィルタ (音色)各フレームで観測される多重音スペクトル

図 1 無限重畳離散全極型モデル (iSDAP)：ソース数 I とフィルタ
数 J がいずれも無限に発散した場合の極限を考える．

ここで，a = [a0 , · · · , aP ]T は自己回帰フィルタの係数（たこれらの問題を解決するため，本研究では，無限重畳離だし a0 = 1）であり，{sm }M
m=1 は線形予測誤差である．

散全極型モデル (Inﬁnite Superimposed Discrete All-pole

式 (1) はソース・フィルタモデルの観点から説明すること

Model: iSDAP) とよぶ多重音解析のための新しいソース・ができる．すなわち，x が音声信号であるとすると，s はフィルタ NMF を提案する（図 1）．本モデルは，スペク

声帯（ソース）から生成される音源信号に対応し，a は声トル包絡が調波構造のピークで不要に鋭いピークを持つ道（フィルタ）の反響特性を表している．

ことを防ぐため，調波構造の離散的なピークのみを用いてスペクトル包絡を推定できる離散全極型モデル (Discrete

式 (1) は，s を入力にとり，x を出力する線形系とみなすことができ，その振る舞いはパラメータ a で決定される．

All-pole Model: DAP) [17] に着想を得ている．しかし，多式 (1) は a と x との畳み込みであることから

重音に対して DAP を適用するには，個々の調波構造を分A(z)X(z) = S(z) i.e., X(z) = S(z)F (z)離し，F0（スペクトル包絡推定に利用する離散的な倍音周波数）を推定する必要があった．本研究の重要な貢献の一つは，無限複合自己回帰モデルと同様に，重畳された複数の調波構造を扱えるように DAP を拡張したことである．
提案モデルでは，離散的な周波数を取り扱えるという DAPの副次的な効果として，「対数周波数」スペクトログラムを・音色（スペクトル包絡）・音量に適切な個数の音高（F0）
分解することができる．

が成立する．ここで，X(z) および S(z) はそれぞれ，x お
よび s の z 変換であり，X(z) =

本章では，音響信号に対するソース・フィルタモデルとNMF への応用について紹介する．従来のモデルは一般に，STFT で得られる線形周波数領域で定式化されている．次章で提案するモデルを定式化する際の基礎となるため，これらのモデルの確率的解釈について説明する．

M

xm z −mS(z) =

m=1

M

sm z −m

(3)

m=1

1
で与えられる．ここで，F (z) = A(z)は全極型モデルであdef

り，次式で与えられる．

F (z) =

2. 関連研究

(2)

1
1
= P
−p
A(z)
p=0 ap z(4)m
いま，2π M= ωm とし，z = eiωm を式 (2) に代入すると,

この線形系の伝達特性のフーリエ領域表現を得る．

X(eiωm ) = S(eiωm )F (eiωm )(5)

ここで，{X(eiωm )}M
m=1 は観測信号 x の複素スペクトルであり，{S(eiωm )}M
m=1 は音源信号 s の複素スペクトル，{F (eiωm )}M
m=1 は全極型伝達関数の周波数応答である．

2.1 線形予測分析（全極型モデル）
線形予測分析 (LPC) は観測スペクトルに対してスペク

LPC の目的は，s がガウス性白色雑音であるという強いトル包絡を推定するための音響信号モデリング手法であ仮定の下で，フィルタ係数 a を推定することである．まる．LPC では，与えられたスペクトルに対応する観測音響

ず，S(eiωm ) が，すべての周波数帯域で同じ分散を持つ複信号 x = {xm }M
m=1 （M は窓幅）が，P 次の自己回帰過程

素ガウス分布に従うことを仮定する．

に従うことを仮定している．

xm = −

P
p=1

ap xm−p + sm i.e.,S(eiωm ) ∼ Nc (0, σ 2 )P

ap xm−p = sm (1)p=0c 2014 Information Processing Society of Japan(6)

ここで，σ 2 は音源信号 s の白色スペクトルのパワーを表す．式 (5) および式 (6) を用いると，

2

Vol.2014-MUS-104 No.9
2014/8/26

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

[dB]

最尤推定における式 (9) の最大化は，次式で与えられる

DAPでは調波構造のピークの影響を受けにくい

0

IS ダイバージェンスの最小化と等価である．
  XmXm
DIS (Xm |σ 2 Fm ) =
−
log
−
1
(10)
σ 2 Fmσ 2 Fm離散全極型モデル (DAP)

-20
-40
-60

m∈Ω

-80
-100
-120
-140

図 2

線形予測分析 (LPC)

観測スペクトル（調波構造）
1

2

3

4

5

6

7

8 [kHz]

ここで，Fm は次式で与えられる．
1
1
Fm = 
2 = TaU P
ma p=0 ap e−iωm p (11)

調波構造をもつ観測スペクトルに対して線形予測分析 (LPC)

U m は (P + 1) × (P + 1) のテプリッツ行列であり，各要素

および離散全極型モデル (DAP) で推定されたスペクトル包絡は [U m ]pq = cos(ωm (p − q)) で与えられる．式 (10) を σ 2X(eiωm ) ∼ Nc (0, σ 2 |F (eiωm )|2 )(7)を得る．ここで，Xm = |X(eiωm )|2 および Fm = |F (eiωm )|2と定義すると，式 (7) は簡潔にに関して偏微分してゼロとおくことで，更新則1  Xm
σ2 =
|Ω|
Fm

(12)

m∈Ω

を得る．ここで，|Ω| は Ω に含まれる周波数の個数を表す.
一方，式 (10) を a に関して偏微分すると，

2

Xm ∼ Exponential(σ Fm )スペクトル包絡である．式 (8) は LPC の確率モデルであ∂DIS (Xm |σ 2 Fm )(13)
= 2(R − R )a∂aを得る．ここで，R および R は正定値行列であり，

1 
R= 2
Xm U m R =Fm U m(14)σ2るので，{Fm }M
m=1 （すなわち a）および σ は最尤推定をで与えられる．ベクトル a に関する乗法更新則は(8)と書ける．ここで，図 2 で示される通り，{Xm }M
m=1 は観
M
測信号 x のパワースペクトル，{Fm }M
m=1 は {Xm }m=1 の用いて求めることができる [16].

LPC の主な問題のひとつは，音源信号の周期性（例：バイオリンの弦振動）に起因する音高をもつ音響信号を解析すると，観測スペクトル {Xm }M
m=1 に含まれる調波構造に影響を受け，推定されるスペクトル包絡 {Fm }M
m=1 は倍音
周波数において不要に急峻なピークをもつことである．こ
の理由は，実際には調波構造を構成する倍音のみがスペク

m∈Ω

m∈Ω

a ← R−1 R a

(15)

で与えられる．R は正定値行列であることから常に逆行列をもち，式 (15) は安定的に動作する．式 (10) が収束するまで式 (12) および式 (15) を交互に反復することでパラメータの最適化を行う．ただし，σ 2 を調節することで，a0 = 1を満たすようスペクトル包絡 {Fm }m∈Ω を正規化する．

トル包絡からの信頼できるサンプルとみなせるにもかかわ調波構造を対象とした DAP の主な問題のひとつは，Ω

らず，M 個すべての周波数帯域が全極型モデルの推定に用を定義するため，観測スペクトルの F0 をあらかじめ与えいられるからである．

る必要があることである．したがって，複数の調波構造が重畳した多重音を解析する場合には，個々の調波構造を分離し，F0 を推定する必要がある．

2.2 離散全極型モデル
離散全極型モデル (DAP) は，LPC のもつ問題を解決するために提案されたスペクトル包絡推定手法のひとつである．DAP は LPC の拡張であるとみなすことができるた

2.3 複合自己回帰モデル
複合自己回帰モデルは，図 3 に示す通り，多重音のスペ

め，その確率モデルは式 (8) と同じ形をとるが，一部の周クトログラムを I 個の微細構造（ソース）と J 個のスペク

波数ビン Ω 上にのみ定義されている点が異なる．

トル包絡（フィルタ）とに分解することができるソース・Xm ∼ Exponential(σ 2 Fm ) m ∈ Ωフィルタ NMF の一種である [5]．いま，観測パワースペク(9)

ここで，Ω = {1, · · · , M } とすると，DAP は LPC と等価になる．調波構造に対してスペクトル包絡を推定するためには，離散的な倍音周波数の集合 Ω を定義すればよい．こうすることで，図 2 に示すように，調波構造のピークの近トログラムを X ∈ RM ×N する．ここで，M は周波数ビンの数，N はフレーム数である．非負値行列 X を，三つの因子 S ，F ，H に分解することを考える．

Xmn ≈

J
I 


Sim Fjm Hnij = Ymn
def

(16)

i=1 j=1

くを正確に通り，倍音周波数に影響されにくい極を持つス

こ こ で ，{Sim }M
m=1 は ソ ー ス i の パ ワ ー ス ペ ク ト ル ，

ペクトル包絡を推定することができる．式 (9) で示される{Fjm }M
m=1 はフィルタ j のパワースペクトル，Hnij は2

尤度を最大化するため，a および σ を交互に反復最適化

フレーム n におけるソース i・フィルタ j の組み合わせのする手法が提案されている [17]．この手法は，乗法更新ア

音量を表す．これらすべての変数は，X から推定すること

ルゴリズムの一種とみなすことができる [2, 9]．

ができる．

c 2014 Information Processing Society of Japan

3

Vol.2014-MUS-104 No.9
2014/8/26

情報処理学会研究報告
IPSJ SIG Technical Report
線形予測分析
ソースとフィルタの各周波数ビンの積調波構造⊗

こ こ で ，E[Xmn ] = Ymn が 成 立 し ，{Xmn }M
m=1 お よ び
全極型スペクトル包絡

幅スペクトルである．ただし，KL ダイバージェンスを用いる場合には，「パワー」ではなく「振幅」スペクトログラ

基底⊗ムに対して定式化することが一般的である [14, 18].

音量

線形周波数 [Hz]線形周波数 [Hz]

…

…

⊕

ソース (音高){Ymn }M
m=1 はフレーム n における観測および再構成振CAR の拡張として，ノンパラメトリックベイズ理論を用いることで，観測データ X に合わせて適切な個数のソースとフィルタを推定することができる無限複合自己回帰モフィルタ (音色)デル (iCAR) が提案されている [6]．

各フレームで観測される多重音スペクトル

図 3

線形周波数領域における複合自己回帰モデル (CAR)

Xmn ≈

I→∞
→∞
 J
i=1

θi φj Sim Fjm Hnij

(22)

j=1

ここで，θi および φj はそれぞれ，はソース i およびフィル

2.3.1 基本的な定式化
CAR の確率モデルを定式化するためには，音源信号の

タ j の大域的な重みを表す．ガンマ過程 NMF [10] と同様

スペクトルが従う確率分布を仮定する必要がある．LPC でに，θ および φ に対してガンマ過程事前分布を仮定するこは式 (6) で示すように音源信号がガウス性白色雑音であるとにより，理論的には無限個存在するソース・フィルタのと仮定していたのに対し，CAR では各周波数ビンごとにうちで，観測データに合わせて高々有限個だけが実質的に異なるパラメータを持つ独立な確率分布を許容する．

アクティベートされる機構を実現できる．

Si (eiωm ) ∼ Nc (0, Sim )(17)

ここで，{Si (eiωm )}M
m=1 はソース i の複素スペクトルである．式 (5) および式 (17) を用いると，

Xijmn (eiωm ) ∼ Nc (0, Sim Fjm Hnij )CAR の別の拡張として，ソーススペクトル {Sim }M
m=1が調波構造をもつような制約を導入する手法も提案されている [6, 19]. もし，音源信号が周期的インパルスの系列であるとすると（声道の理想的なモデル），{Sim }M
m=1 は同じ強度で線形周波数軸に等間隔に並んだ倍音成分からなる調(18)を得る．ここで，{Xijmn (eiωm )}M
m=1 はフレーム n におけ波構造をもつことになる．F0 は，式 (21) で定義される尤
度を最大化するように最尤推定で求めることができる．こるソース i・フィルタ j の組み合わせから生成されたガウのモデルにおいては，スペクトル包絡推定と同時に F0 推

ス分布の再生成を考慮すると，定を行うことができるため，DAP の問題点を解決するこ

Xmn (eiωm ) ∼ Nc (0, Ymn )とができる可能性を秘めている．

(19)を得る．{Xmn (eiωm )}M
m=1 はフレーム n における観測複素スペクトルである．式 (19) は次式と等価である．

3. 無限重畳離散全極型モデル
本章では，ウェーブレット変換で得られる対数周波数ス
ペクトログラムに対してソース・フィルタ分解を行うため

Xmn ∼ Exponential (Ymn )(20)

ここで，E[Xmn ] = Ymn が成立しており，{Xmn }M
m=1 および {Ymn }M
m=1 はフレーム n における観測および再構成パ
ワースペクトルである．これは，式 (16) において Xmn とYmn の近似誤差を評価するコスト関数として，理論的には
板倉-齋藤 (IS) ダイバージェンスが適切であることを示している [5]．ただし，IS ダイバージェンスは Ymn に対して
凸関数ではないため，Ymn の最適化は局所解に陥りやすい傾向がある．

Xmn と Ymn との近似誤差を評価するコスト関数として，
IS ダイバージェンスの代わりにカルバック・ライブラ (KL)ダイバージェンスを用いると異なる確率モデルが得られる [6]．すなわち，式 (20) の代わりに次式を仮定する．

c 2014 Information Processing Society of Japan

デルでは，与えられた音楽音響信号に対し，各フレームにおける複数の F0 を推定すると同時に，複数の典型的なス
ペクトル包絡（楽器の音色）を発見することができる．これを実現するため，スペクトル包絡推定のための離散全極型モデル (DAP) [17] を，F0 とスペクトル包絡の同時推定のための複合自己回帰モデル [5, 6] の枠組みに統合する．
本研究の主な貢献は以下の通りである．

( 1 ) 重畳離散全極型モデル：多重音に含まれる複数の調波
構造それぞれに対してスペクトル包絡を推定すること

2.3.2 拡張Xmn ∼ Poisson (Ymn )の無限重畳離散全極型モデル (iSDAP) を提案する．本モ(21)ができる．通常の DAP では，単独音のスペクトルに対するスペクトル包絡推定を目的としていた．

( 2 ) 基本周波数モデリング：各フレームごとに異なる F0（ソーススペクトル）の集合が含まれることを許容することにより，ビブラートやポルタメントなどの基本周
波数の微細変動を精緻にモデル化することができる．

4

Vol.2014-MUS-104 No.9
2014/8/26

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

すべてのフレームで共有されているため，各ソースス

振幅

通常の CAR においては，ソーススペクトルの集合が全極型スペクトル包絡が倍音成分の重みを決定

ペクトルは半音レベルの F0 に対応していることが期
待されている．
・・・( 3 ) 対数周波数領域での定式化：本モデルは，人間の聴覚特性にあった対数周波数スペクトログラムを分解する

図 4

・・・

対数周波数 [cents]

対数周波数領域におけるフレーム n におけるソース i とフィことができる初めてのソース・フィルタ NMF である．

ルタ j との組み合わせによる楽器音スペクトルの生成

この実現のためには，スペクトル包絡推定の際に離散
1
2
Smnir = exp − 2 (fm − (μni + 1200 log2 r)) (25)2σ的な周波数のみを考慮することができる DAP が不可欠である．従来は，対数周波数領域における単独発話
スペクトルのモデル化に利用されていた [20].
本研究では，観測スペクトログラムに合わせてソースと

ここで，μni はフレーム n におけるソース i の対数基本周
波数 F0 [cents]，fm はスペクトログラムにおける m 番目

フィルタの個数を適切に調節するため，ノンパラメトリッの周波数ビンに対応する対数周波数，σ 2 は各倍音を中心とクベイズ理論に基づく確率モデルを定式化する．本モデル

したメインローブの広がりである．

は，観測データが無限にあれば，理論的には無限個のソースとフィルタが存在することを仮定している．一方，現実的には有限の観測データが与えられると，そこに含まれる
高々有限個のソースとフィルタを推定することができる．
ノンパラメトリックベイズ理論を用いると，無限次元の空間内でスパースな学習を行うことができる．

数周波数領域における振幅スペクトログラムを X ∈ RM ×Nとする．ここで，M は周波数ビンの個数であり，N はフレームの個数である．非負値行列 X を，二つの因子 W お
よび H に分解することを考える．
⎞
⎛
→∞
I→∞
 J
Xmn ∼ Poisson ⎝θni φj Wnijm Hnij ⎠ここで，aj ≡ [aj0 , · · · , ajP ]T であり，ωnir はフレーム n

(23)

j=1

み，φj は全フレームにおけるフィルタ j の大域的な重みを表す．Hnij はフレーム n におけるソース i・フィルタ j の組み合わせの音量を表す．{Wnijm }M
m=1 は，フレーム n におけるソース i・フィルタ j から生成されたある楽器音の振幅スペクトルである．もし，θni および Wnijm が時間変化しないとすると，すなわち，ソーススペクトルの集合がすべてのフレームで共通であるとすると，式 (16) で示される無限複合自己回帰モデルと同じ形式をとる．

重畳離散全極型モデル (SDAP) の重要な特徴は，ソースとフィルタを組み合わせる際に，ソースの倍音成分の重みのみが全極型スペクトル包絡によって制御される点である．
示す通り，全周波数帯域にわたって周波数ビンごとにソースとフィルタとの要素積を計算する．したがって，スペク
トル包絡を推定する際に，SDAP では倍音成分（調波構造のピーク）のみが参照されるのに対して，CAR では調波
構造スペクトル全体が参照される．これは，SDAP が離散
全極型モデル (DAP) の多重音拡張であるとみなせるのに対して，CAR は線形予測分析 (LPC) の多重音拡張とみなせることを意味する．

3.1.2 事前分布の設計
無限次元のベクトル θ n = [θn1 , · · · , θnI ]T および φ =[φ1 , · · · , φJ ]T に対してスパースな学習を行うため，ガンマ過程事前分布を仮定する [6, 10]．これを近似的に実現する

3.1.1 ソース・フィルタの組み合わせ
図 4 で示す通り，振幅スペクトル {Wnijm }M
m=1 は対数
周波数領域において調波構造を持つことを仮定する．

Fnijr Smnir

要素は [U (ω)]pq = cos(ω(p − q)) で与えられる．Fnijr は
パワーではなく，振幅を表すことに注意する．

一方，複合自己回帰モデル (CAR) においては，式 (16) で

ここで，θni はフレーム n におけるソース i の局所的な重Wnijm =

1−1
 = aTj U (ωnir )aj 2 (26)
Fnijr = P p=0 ajp e−ωnir pi [rad] であり，U (ω) は (P + 1) × (P + 1) の行列であり，各本節では iSDAP の確率モデルの定式化を行う．いま，対R

において全極型伝達関数を用いて表現する．

におけるソース i の第 r 次倍音に対応する正規化角周波数

3.1 確率モデルの定式化

i=1

ここで，倍音成分の重み {Fnijr }R
r=1 は，対数周波数領域

(24)

r=1

ひとつの方法として，θn および φ の各要素に対して独立なガンマ事前分布を仮定する．

αθ
, αθ
θni ∼GammaIαφφj ∼Gamma
, αφ
J

(27)
(28)

ここで，R は倍音の個数であり，{Smnir }M
m=1 はフレーム

ここで，αθ および αφ は超パラメータであり，ガンマ過程

n におけるソース i の第 r 次倍音成分の単峰的なスペクトの集中度と呼ばれる．打ち切りレベル I を無限に大きくしルであり，次式で表すものとする．

ていけば，ベクトル θn は集中度 αθ をもつガンマ過程から

c 2014 Information Processing Society of Japan

5

Vol.2014-MUS-104 No.9
2014/8/26

情報処理学会研究報告
IPSJ SIG Technical Reportのランダムサンプルとみなすことができる．ここで，任意の正の実数 に対して θni > を満たす実効的な要素数 I +はほとんど確実に有限であることが証明されている．現実的には，I を αθ に比べて十分大きく設定すれば，θn の I
個の要素のいくつかだけがゼロよりある程度大きな値をとることが期待できる．一方，音量 Hnij に関しては事後分
布計算の容易さから，ガンマ事前分布を仮定する．

ここで，以降の説明のしやすさのため，Xmnijr およびYmnijr を次式で定義しておく．
Xmnijr =λmnijr Xmn(34)

Ymnijr =E[θni φj Smnir Fnijr Hnij ](35)

3.3 確率変数 θ, φ, H に対する変分ベイズ推論
一般に，式 (30) で分解されたそれぞれの確率変数に関す

Hnij ∼ Gamma(aH , bH )(29)

ここで，aH および bH は超パラメータである．

3.2 確率モデルの変分推論
我々の目標は，ベイズの定理に従って確率変数の事後分布,H ;μ,a)p(θ ,φ,H )p(θ, φ, H|X; μ, a) = p(X |θ ,φp(を求め，パX ;μ,a)る変分事後分布は次式で求められる．

q(θ) ∝ exp(Eq(φ,H ) [log p(X, θ, φ, H; μ, a)])(36)q(φ) ∝ exp(Eq(θ ,H ) [log p(X, θ, φ, H; μ, a)])(37)q(H) ∝ exp(Eq(θ ,φ) [log p(X, θ, φ, H; μ, a)])(38)

対数周辺尤度の変分下限 L は各確率変数に関する線形項おラメータ μ および a については周辺尤度 p(X; μ, a) を最よび対数項を含んでいることから，各確率変数に対応する

大化するよう最尤推定（経験ベイズ法による最適化）を行

変分事後分布は，事前分布と同じ形であるガンマ分布をとうことである．しかし，周辺尤度 p(X; μ, a) を解析的に計ることが分かる．したがって，算することはできないため，変分ベイズ法を用いて，事後

q(θni )=Gamma(aθni , bθni )(39)

q(φj )=Gamma(aφj , bφj )(40)

Hq(Hnij )=Gamma(aHnij , bnij )(41)

分布を因子分解できる形

q(θ, φ, H) =q(θni )nijq(φj )

q(Hnij )

(30)

nij

であると仮定し，対数周辺尤度 log p(X; μ, a) の変分下限

L を最大化するように最適化を行う．まず，対数周辺尤度
log p(X; μ, a) に対する第一の変分下限 L0 は
log p(X; μ, a) ≥ E[log p(X|θ, φ, H; μ, a)](31)

ないが，対数関数に対する Jensen の不等式を利用することで，さらなる変分下限を求めることができる．



≥

mnijr

−



(32)

す．等号が成立する，すなわち L0 が最大化される条件は，(33)で与えられる．式 (32) を用いることで，最大化すべき目的

c 2014 Information Processing Society of Japan

Xmnijr(46)

E[θni φj ]Wnijm

(47)

mr


mr

が最小となるような x を求める問題を考える．C の x に関
ここで，λmnijr は補助変数であり， ijr λmnijr = 1 を満た

関数 L は，L0 の下限として求めることができる．



カラパラメータ x に依存するコスト関数 C が与えられ，Cmnijr

λmnijr ∝ exp(E[log(θni φj Smnir Fnijr Hnij )])aH
nij =aH +ルゴリズム [2, 9] を用いることができる．一般に，あるス

E[θni φj Smnir Fnijr Hnij ]λmnijrE[θni φj Smnir Fnijr Hnij ] + const.

(45)

パラメータ μ および a を推定するためには，乗法更新アE[θni φj Smnir Fnijr Hnij ] + const.

λmnijr Xmn log(44)

3.4 パラメータ μ, a に対する乗法更新

ijr

mnijr



(43)

mnir

bH
nij =bH +E[log p(X|θ, φ, H; μ, a)]
Xmn logE[θni φj Smnir Fnijr Hnij ]
=
−

αφ 
+
XmnijrJ mnirbφj =αφ +E[θni Hnij ]Wnijm

aφj =で与えられる．右辺第一項は依然として解析的に計算でき

mn

(42)

mjr

+ E[log p(θ)] + E[log p(φ)] + E[log p(H)]− E[log q(θ)] − E[log q(φ)] − E[log q(H)] ≡ L0とすると，変分パラメータは次式で求められる．
αθ +aθni =
XmnijrI
mjr

E[φj Hnij ]Wnijmbθni =αθ +する偏微分が ∂C∂x = R − R という二つの正の項の差で表現
可能であるとき，x の乗法更新則は x ← RR x で与えられる．この更新によって，x が正であるとき，正値性は自動的に保持される．本研究では，周辺尤度の変分下限 L の符号を反転したものをコスト関数とみなせば，乗法更新アルゴリズムを適用できる．
まず，−L を基本周波数 μni に関して偏微分すると−∂L= Rni − Rni
∂μni

(48)

6

Vol.2014-MUS-104 No.9
2014/8/26

情報処理学会研究報告
IPSJ SIG Technical Report[cents]
10,000
9,000
8,000
7,000
6,000
5,000
4,000
3,000
2,000
1,000

40
観測スペクトログラム

30

[dB]

20
10
0
-10

最も優勢なフィルタ

2番目に優勢なフィルタ

3番目に優勢なフィルタ

-20
-30
5

[cents]
10,000
9,000
8,000
7,000
6,000
5,000
4,000
3,000
2,000
1,000

10

15

20

25

30 [s]（アタック音に対応）-40
1,000 2,000 3,000 4,000 5,000 6,000

7,000 8,000 9,000 10,000 [cents]

再構成スペクトログラム

図 6 “MUS-mz 570 1 ENSTDkCl” から推定された全極型スペク
トル包絡（フィルタ）

4. 評価実験
5

10

15

20

25

30 [s]

本章では，無限重畳離散全極型モデル (iSDAP) を多重
図 5 “MUS-mz 570 1 ENSTDkCl” の観測ウェーブレットスペク

音解析に適用した結果について報告する．

トログラムと無限重畳離散全極型モデルの推定結果から生成した再構成スペクトログラム

4.1 実験条件
実験には，MAP ピアノデータベース [21] に含まれるピを得る．ここで，Rni および Rniは正の項であり，

Rni =

アノ演奏 “MUS-mz 570 1 ENSTDkCl” から，文献 [8] と(μni + 1200 log2 r)Xmnijr + fm Ymnijr (49)

mjr

=
Rni



fm Xmnijr + (μni + 1200 log2 r)Ymnijr (50)mjr

同様に冒頭 30 秒を切り出した 44.1kHz・モノラルの音響信号を用いた．対数周波数領域の振幅スペクトログラムを得るため，ガボール関数を用いた連続ウェーブレット変換を行った．その際，時間方向に 10 [ms] ずつシフトさせながで求まる．したがって，スカラパラメータである μni の乗ら，900 [cents] から 10400 [cents] まで 10 [cents] 間隔で解

法更新則は次式で与えられる．

．超パラメータは析を行った（M = 960 および N = 3000）−1 
μni ← Rni
Rni μni(51)

I = 88, J = 3, αθ = αφ = 1, R = 20, P = 12, σ 2 = 100,aH = 1, bH = Eemp [Xmn ]−1 と設定した．観測スペクトロ

この更新によって μni の正値性は保証されるが，μni は対グラム X にはピアノ音しか含まれていないため，フィル

数基本周波数 [cents] であるから理論的には負の値をとっタ数は J = 3 とした．各フレームにおける対数基本周波数てもよい．しかし，本研究は過去の研究を参考に，対数周{μni }Ii=1 は，ピアノの 88 鍵の音高に対応するように初期波数軸の基準となる 0 [cents] をピアノの最低音よりもさら化し，他のパラメータはランダムに初期化した．

に低い 440 × 23/12−5 = 16.3516 [Hz] に設定しているため，μni は正の値であると考えて実用上は問題ない．
次に，関連研究 [2, 9] で行われているのと同様に，−L を各 フ レ ー ム に 含 ま れ る 複 数 の F0 は 閾 値 処 理 に よ って 求 め た ．具 体 的 に は ，ソ ー ス i の 基 本 周 波 数 成 分j θni φj Fnij1 Smni1 Hnij がある閾値より大きい場合，フaj に関して微分すると

レーム n には μni で示される F0 が含まれているとした．

−∂L= (Rj − Rj )aj∂aj

閾値は，フレームレベルの適合率と再現率のバランスを(52)を得る．ここで，Rj および Rj は正定値行列であり，

Rj =



2
Xmnijr Fnijr
U (ωnir )(53)

3
Ymnijr Fnijr
U (ωnir )

(54)

mnir

Rj =



とって F 値が最大化するように設定した．

4.2 実験結果
図 5 に示す通り，変分ベイズ法と乗法更新アルゴリズムを用いることにより，観測スペクトログラムをよく近似すmnir

るパラメータを学習することができた．しかし，主に発音
時刻付近に見られる周波数方向に広く分布しているノイズで求まる．したがって，ベクトルパラメータである aj に成分（例：ハンマーがピアノ弦をたたく際の打撃音）に対関するベクトル型の乗法更新則は次式で与えられる．

しては，必要以上に多数の調波構造の重ね合わせによって
aj ← R−1
j R j aj(55)

できるだけ正確に近似しようとする振る舞いが見られた．
したがって，現状では最新の手法 [7, 8] に比べて F0 の推定

これは，DAP を用いたスペクトル包絡推定における乗法

精度は遠く及ばないが，音響理論に裏付けられた統計的モ更新則式 (15) と同様の形式をしていることが分かる．

デリングアプローチは有望であると考えている．図 6 に示

c 2014 Information Processing Society of Japan

7

Vol.2014-MUS-104 No.9
2014/8/26

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

すように，推定されたスペクトル包絡は，調波構造のピークに影響を受けず，なだらかな形状であった．このことは，[6]

観測データのみから，倍音成分の重みは指数的に減衰する傾向があることが学習できたことを示す．

iSDAP の問題点として，ソースやフィルタの実効的な個[7]

数は自動調節されるものの，F0 推定のためには最終的に閾値処理が必要になることが挙げられる．この問題を解決するには，音符の存在を表現するバイナリ潜在変数を導入する方法が考えられる．F0 推定精度を改善するには，文[8]

献 [6, 11] と同様に，音高をもつ楽器音と音高を持たない
打楽器音に対してそれぞれ異なるモデルを定式化し，それらの重ね合わせで音楽音響信号を表現する方法が考えられる．ただし，ベイズモデルにおいては，事前分布の影響で[9]異なる性質を持つモデルを適切な重みで安定的に組み合わせることは簡単ではなく，さらなる研究が必要である．

[10]

5. おわりに
本稿では，無限重畳離散全極型モデル (iSDAP) と呼ぶ
ノンパラメトリックベイズソース・フィルタ NMF につい[11]て述べた．本モデルは，ウェーブレット変換で得られる対数周波数スペクトログラムを，音源の調波構造スペクトル（音高），フィルタの全極型スペクトル包絡（音色）
，および

[12]

それらの組み合わせの音量に分解することができる．事後
分布の推定およびパラメータの推定のために変分ベイズ法と乗法更新アルゴリズムを統合する手法を提案し，実験で[13]動作を確認した．今後は，多重音基本周波数推定と自動採譜の間にあるギャップを埋めるため，音符配置に関する事
前分布をベイズ的に統合することを検討している．
[14]
謝辞: 本研究の一部は，JSPS 科研費 26700020, 24220006，

24700168 および JST CREST「OngaCREST プロジェクト」[15]の支援を受けた．

参考文献
[1]

[2]

[3]

[4]

[5]

Smaragdis, P., Févotte, C., Mysore, G., Mohammadiha,N. and Hoﬀman, M.: Dynamic Source Separation Using Nonnegative Factorizations: A Uniﬁed View, IEEE
Signal Proc. Magazine, Vol. 31, No. 3, pp. 66–75 (2014).
Hennequin, R., Badeau, R. and David, B.: NMF With
Time-Frequency Activations to Model Nonstationary
Audio Events, IEEE Trans. on ASLP, Vol. 19, No. 4,
pp. 744–753 (2011).
Durrieu, J.-L., Richard, G., David, B. and Févotte, C.:
Source/Filter Model for Unsupervised Main Melody Extraction From Polyphonic Audio Signals, IEEE Trans.
on ASLP, Vol. 18, No. 3, pp. 564–575 (2010).
Virtanen, T. and Klapuri, A.: Analysis of Polyphonic
Audio Using Source-Filter Model and Non-Negative Matrix Factorization, NIPS Workshop on Advances inModels for Acoustic Proc. (2009).
Kameoka, H. and Kashino, K.: Composite Autoregressive System for Sparse Source-Filter Representation of
Speech, IEEE International Symposium on Circuits

c 2014 Information Processing Society of Japan


[16]

[17]

[18]

[19]

[20]

[21]

and Systems (ISCAS), pp. 2477–2480 (2009).
Yoshii, K. and Goto, M.: Inﬁnite Composite Autoregressive Models for Music Signal Analysis, International
Conference on Music Information Retrieval (ISMIR),
pp. 79–84 (2012).
Carabias-Orti, J. J., Virtanen, T., Vera-Candeas, P.,
Ruiz-Reyes, N. and Cañadas-Quesada, F. J.: Musical Instrument Sound Multi-Excitation Model for NonNegative Spectrogram Factorization, IEEE Journal of
Selected Topics in Signal Proc., Vol. 5, No. 6, pp. 1144–
1158 (2011).
Bertin, N., Badeau, R. and Vincent, E.: Enforcing Harmonicity and Smoothness in Bayesian Non-Negative Matrix Factorization Applied to Polyphonic Music Transcription, IEEE Trans. on ASLP, Vol. 18, No. 3, pp.
538–549 (2010).
Badeau, R. and Ozerov, A.: Multiplicative Updates
for Modeling Mixtures of Non-stationary Signals in theTime-Frequency Domain, European Signal Proc. Conference (EUSIPCO) (2013).
Hoﬀman, M., Blei, D. and Cook, P.: Bayesian Nonparametric Matrix Factorization for Recorded Music, International Conference on Machine Learning (ICML), pp.
439–446 (2010).
Benetos, E., Ewert, S. and Weyde, T.: Automatic Transcription of Pitched and Unpitched Sounds from Polyphonic Music, International Conference on Acoustics,
Speech, and Signal Proc. (ICASSP), pp. 3131—3135(2014).
Vincent, E., Bertin, N. and Badeau, R.: Adaptive Harmonic Spectral Decomposition for Multiple Pitch Estimation, IEEE Trans. on ASLP, Vol. 18, No. 3, pp. 528–
537 (2010).
Bryan, N. J., Mysore, G. and Wang, G.: Source Separation of Polyphonic Music with Interactive User-Feedback
on a Piano Roll Display, International Conference onMusic Information Retrieval (ISMIR), pp. 119–124(2013).
Lee, D. and Seung, H.: Algorithms for Non-Negative Matrix Factorization, Neural Information Processing Systems (NIPS), pp. 556–562 (2000).
Cemgil, A. T.: Bayesian Inference for Nonnegative Matrix Factorisation Models, Computational Intelligence
and Neuroscience, Vol. 2009, Article ID 785152 (2009).
Itakura, F. and Saito, S.: Analysis Synthesis Telephony based on the Maximum Likelihood Method, International Congress on Acoustics (ICA), pp. C17–C20(1968).
El-Jaroudi, A. and Makhoul, J.: Discrete All-Pole Modeling, IEEE Trans. on Signal Proc., Vol. 39, No. 2, pp.
411–423 (1991).
Févotte, C., Bertin, N. and Durrieu, J.-L.: Nonnegative
Matrix Factorization with the Itakura-Saito Divergence:With Application to Music Analysis, Neural Computation, Vol. 21, No. 3, pp. 793–830 (2009).
安良岡直希, 奥乃博: 調波・非調波・音色構造因子分解による音響信号分析と音源分離インターフェースへの応用, 情報処理学会研究報告, Vol. 2012-MUS-94, pp. 1–8 (2012).
亀岡弘和: 全極型声道モデルと F0 パターン生成過程モデ
ルを内部にもつ統一的音声生成モデル, 日本音響学会秋季
研究発表会講演論文集, pp. 211—214 (2010).
Emiya, V., Badeau, R. and David, B.: Multipitch Estimation of Piano Sounds Using a New Probabilistic
Spectral Smoothness Principle, IEEE Trans. on ASLP,
Vol. 18, No. 6, pp. 1643–1654 (2010).

8