3-3-12
非負値周期成分分析手法による音楽音響信号の音源分離∗◎林亜紀 (NTT EV 研), 亀岡宏和 (NTT CS 研), △松林達史, 澤田宏 (NTT EV 研)
はじめに「サステイン」「リリース」などの楽音の状態変化を基底スペクトルに仮定したモデル [5]，ビート（拍）に音楽音響信号の多くは，複数の楽器などの音源か応じたリズムをアクティベーションに仮定した制約付ら構成され，混合信号の分離は音楽構造理解や自動NMF [6] が提案されている．
採譜，楽器やピッチの部分的な編集，雑音除去，欠本稿では，音楽音響信号に含まれる，各音符やドラ損値補間など多くの場面で必要になる．近年，音響
ム音の発音タイミングにおける，リズムの周期性を信号処理分野において，自動採譜や音源分離などへ仮定した制約付 NMF を提案する．特にベースやドラの有効性から非負値行列分解（Non-negative Matrix
ムなどの楽器で周期性が想定される．単純かつ単一Factorization, NMF) [1] への注目が高まっている [2]． な周期成分を抽出する手法 [7] が提案されているが，NMF では，混合信号（観測信号）のスペクトログラ
周期性は単一であるとは限らず，複数の音源で異なるムを非負値行列 Y で表現した上で，2 つの非負値行複数の周期性が存在する可能性がある．そこで，周期列 H ，U の積の形に分解する．これは，各時間窓で性を想定した一部の基底で，どの音源がいつ発音す観測されたスペクトルを，時間変化する係数行列 Uるかを示すアクティベーションに基底毎に異なる周期（アクティベーション）の値でスケールされた基底ス性を仮定し，目的関数に制約として加えることによペクトル H の線形和で近似することに該当し，観測り，周期的な発音パターンの抽出を促進する．周期性
行列はランク 1 のスペクトログラム Hω,k Uk,t の和のを仮定する制約を加えた目的関数において，最適解形で表現される．NMF の重要な特徴として，非負値を求めるアルゴリズムを定式化した上で，教師あり・制約により分解行列がスパースな値を持つことが挙
教師なし両方の枠組みにおいて音源分離の性能を評げられる．例えば，U には多くの 0 要素が含まれる． 価する．
この特徴により，各時間窓の観測スペクトルは，多く2 提案モデルの場合少数の基底スペクトルの和で表現される．従って，各基底スペクトルには観測信号全体において複数
本稿では，提案する周期を考慮したパターン抽出方回出現する典型的なスペクトル（多くは頻出する楽法を非負値周期成分分析手法（Non-negative Periodic音）が抽出され，効率の良い分解表現が可能になる． Component Analysis, NPCA) と呼ぶ．NPCA の目
しかしながら，音楽スペクトログラムに NMF を的関数は以下のように定義される．
適用する際には，2 つの課題がある．第 1 に，精度とI(H, U ) + w1 J(U , α, β) + w2 L(α, β)(1)
分離信号の有意性との間にトレードオフが発生する．
NMF では X = HU と観測行列 Y の間の近似誤差w1 , w2 > 0 は正則化パラメータである．I(H, U ) はを目的関数として定義し，その目的関数を最小化す
通常の NMF の目的関数で考慮されるものと同じ，観る．基底数を増やすことにより近似誤差は減少する
測行列 Y と，NMF モデルによる分解行列の積 X =が，基底数を増やしすぎると分解行列の解釈性が下HU との近似誤差である．以下のような NMF モデがる．例えば，極端な例として Y = HI や Y = IU
ルを考える．などの無意味な分解結果が出力されてしまう場合がKある．ここで I は単位行列を示す．したがって，多∑Y≃X=Hω,k Uk,t(2)ω,tω,tくの基底を用いて精度（少ない近似誤差）を保ちながk=1
ら，優位性の高い分離を行うためには，非負値性以外の有用な制約を HU に課したモデルが必要になる．
ここで，k は基底インデックス，K は基底数を示し，NMF の第 2 の課題として，局所的最適解問題が挙ω と t はそれぞれ周波数と時間窓を示す．Y と X とげられる．既存の NMF の最適化アルゴリズムの多くの間の距離尺度としてＩダイバージェンスを用いる
は収束が保証されているものの，ある限られた領域と，I(H, U ) は以下のように書ける．
においてのみ最適な局所的最適解へ収束してしまう)∑(Yω,t場合が多く，必ずしも全ての領域において最適な大Yω,t log
I(H, U ) =− Yω,t + Xω,t(3)Xω,t域的最適解が求まるとは限らない．音楽音響信号へω,tの適用において，この局所的最適解は，例えば分離J(U , α, β) はアクティベーション Uk,1 , . . . , Uk,T に信号が聴覚で捉えられる音楽事象に対応しないなど，周期性を仮定するための制約を表す項である．T は音楽的に解釈不可能な解となり得る．
時間-周波数表現における時間窓数を示す．提案法で上記の問題を解決するには，データの特徴を反映は，U の非負性を保つために，Uk,t の対数を取ったした制約をモデル HU や目的関数に組み込むことがものが，倍音構造を考慮した正弦曲線の重み付き和考えられる．例えば，音色の類似性を考慮して基底で近似されること，すなわちをクラスタリングする枠組みを導入した手法 [3] が提∑案されている. 音楽の時間変化特徴を考慮した制約付
kk(αm,ncos nPm t + βm,nsin nPm t) (4)log Uk,t ≃NMF も提案されており，発音が一定時間持続するこm,nとを仮定したモデル [4] や，「アタック」「ディケイ」

1∗

Non-negative Periodic Component Analysis for Music Source Separation. by HAYASHI, Aki (NTT ServiceEvolution Laboratories), KAMEOKA, Hirokazu (NTT Communication Science Laboratories), MATSUBAYASHI Tatsushi, SAWADA Hiroshi (NTT Service Evolution Laboratories)

日本音響学会講演論文集- 639 -

2016年3月

めの項L(α, β) =∑ {( ∑kkと表現できることを仮定する．αm,n, βm,n∈ R が，周期 1/Pm の繰り返しパターンの形状を決定する．
ここで，各正弦曲線の周期は基本周波数 Pm の整数倍に限定される形を取っており，n は各周波数の n倍音を表現するための変数である．この制約により，∑kkn (αm,n cos nPm t + βm,n sin nPm t) がちょうど周波数 Pm の周期構造を持つことが保証される．以上により，log Uk,t は M 個の周期関数の和で表現される．
なお，式 (4) は U を正弦曲線の重み付き和の指数の形で近似していると捉えられる．
周期 {Pm } について，例えば Pm = logTM log m +1.0, などと定義することが考えられる．また，倍音を表す定数の集合 {n} について，{1, 2, ..., 9, 10} などの定義が考えられる．上述の例では，想定できる周期の数 M は 20 となる．
提案モデルでは，自動的に Uk,t を最適に表現する周期を推定することができる．図 1 に，ある短い音
楽データにおける Uk,t と，Uk,t の太枠で囲まれた時kk
間窓内に対応する αm,n, βm,nの推定例を示す．例においては，(A)-(D) で示されるいくつかの m と n のkk組み合わせに該当する αm,nと βm,nの値が他の m やn に比べて，大きい値を持つと推定された．下の周期関数の図は，抽出された m, n に対応する cos nPm tまたは sin nPm t の U の太枠内の時間窓における離散曲線を示す．Sum と書かれた曲線は，(A)-(D) に該当する曲線の重み付き和を示し，Exp と書かれた尖った曲線は，Sum の指数を取ったものである．Exp の曲線は実際に推定された Uk,t に近い形をしていることが分かる．
提案手では，アクティベーション U に周期性を仮定した目的関数 J(U , α, β) を，式 (4) の両辺間のユークリッド距離を用いて，以下のように定義する．
∑2J(U , α, β) =log Uk,t − fk,t(5)k,t

ここで，以下のように置いた．
)∑(kkαm,ncos nPm t + βm,nsin nPm tfk,t =+(∑kβm,n) }2 p/2n(7)を定義する．L(α, β) は，高々数個の α と β のグループのみが非ゼロの値を取り，他の全てのグループでは値が 0 になっている場合，小さい値を取る．従って，L(α, β) を目的関数に含めて最小化することにより，Uk,t を最小限の周期で表現できる．図 1 は group
sparsity を考慮した α と β の推定結果を示している．
数種類の m に関連する α, β のグループでのみ，α とβ の値が大きくなっている．グループ m ごとに，倍音 n 方向の二乗和を最小化することにより，n 方向の二乗和が 0 に近くなった場合，全ての n で該当するα, β の値が強制的に 0 に変換される．スパース性を保ちつつ，n 方向では複数の係数が非ゼロになることを許容している．
提案手法では，以上のように各基底スペクトルのアクティベーションに対して，十分な数の周期候補の中から適切に選択された，最小限の周期関数の和で表現される任意の周期を持つ周期的イベントを検出・推定した上で，その周期性を行列分解の枠組みで考慮することが可能になる．

最適化アルゴリズム3

目的関数である式 (1) を最小化する H, U , α, β の値を求めるアルゴリズムについて説明する．繰り返し計算により目的関数を最小化することができる．各ループでは，それぞれ H, U , α と β について最適化を行う 3 段階のステップを実施する．最適化を解析的に行うことは不可能であるが，補助関数法を用いて，目的関数の各項に対して設定した上界関数 I + ，J + ，L+ を最小化することにより，間接的に H, U , α, βの値を最適化する更新式を導出することができる．上界関数を H, U , α, β に関して最小化するステップと，目的関数と上界関数の値が等しくなるように補助変数を更新するステップを繰り返すことで，目的関数を単調に降下させることができる．
まず，I(H, U ) の上界関数 I + については，通常のNMF と同様に，Jensen の不等式を用いて以下のように設定することができる．
∑(I + (H, U ) =Yω,t log Yω,t− Yω,t∑kω,tλω,k,t log)
Hω,k Uk,t− Yω,t + Xω,t(8)λω,k,t∑補助変数 λ は k λω,k,t = 1 を満たす正の重み定数でHUk,tあり，等号成立条件は λω,k,t = ∑ ω,kHω,k Uk,t である．
k(6)m,n式 (4) で U に周期性を仮定すると，m, n の組み合わせによっては，想定する正弦曲線が冗長になる場合がある．Uk,t を最小限の周期関数の重み付き和で表現するために，本発明では α と β に以下に示す groupsparsity を考慮する．ここで，係数 α と β を M 個の異なるグループに分け，group sparsity を誘導するた

日本音響学会講演論文集)
2 p/2

nk,mkFig. 1 An example of the estimates of Uk,t , αm,nkand βm,n .

kαm,n次に，J(U , α, β) の上界関数 J + は以下のように設計できる．
∑{ 1
J + (U , α, β) =+ b(ξ) Uk,t + c(ξ)
Uk,tk,t}(9)− 2 log Uk,t fk,t + (fk,t )2 .
ξ122
ここで，b(ξ) = 2 log
ξ + ξ 2 ，c(ξ) = (log ξ) −2 log ξ− ξである．補助変数 ξ の等号成立条件は ξ = Uk,t であ- 640 -2016年3月

る．なお，上記の上界関数を設計するにあたって使用した不等式は [3] で用いられている．
最後に，f (x) = xp/2 に接する放物線が x > 0 のとき常に f (x) 以上の値を取ることを利用し，L(α, β)に関する上界関数 L+ を以下のように設計する．
)∑ { p p/2−1 ( ∑p/2k2αm,n− ϕk,m + ϕk,mL+ (α, β) =ϕk,m2nk,m)}
p p/2−1 ( ∑ k 2p/2+ νk,mβm,n − νk,m + νk,m2n(10)ϕ0，νk,m > ∑0 であり，等号成立条件は ϕk,m =∑k,m >k2k 2α，
ν=k,mn m,nn βm,n である．
続いて，上述した上界関数を用いて H, U , α, β の更新式を導出する．更新式は，上界関数を偏微分したものを 0 と置いて，各変数について解くことにより得られる．H ，U の公式は以下のように導出される．
∑Yω,t λω,k,tHω,k ← t∑,(11)t Uk,t√dk,t + dk,t 2 + 4w1 ak,t,(12)Uk,t ←
2 dk,t∑こ=ω Hω,k + w1 b(ξ)，dk,t∑ こ で ak,t =ω Yω,t λω,k,t + w1 2fk,t である．α，β は，以下のよkk
うな αm,nと βm,nに関する連立方程式を解くことにより更新できる．
∑ (∑kw1cos n′ Pm′ t cos nPm t · αm′ ,n′m′ ,n′+∑tksin n′ Pm′ t cos nPm t · βm′ ,n′tp/2−1 kαm,n =+w2 p ϕk,mw1∑ (∑m′ ,n′+∑∑log Uk,t cos nPm t (13)tk
cos n′ Pm′ t sin nPm t · αm′ ,n′tksin n′ Pm′ t sin nPm t · βm′ ,n′tp/2−1 kβm,n =+w2 p νk,m)∑)log Uk,t sin nPm t (14)t

連 立 方 程 式 AX = B の 解 は ，例 え ば Matlab で は X = A\B ，Python 言 語 の NumPy
では X =numpy.linalg.solve(A, B)，R 言語ではsolve(A, B) とすることにより解くことができる．

4

評価実験

本章では，提案手法を使用した教師なし音源分
離，教師あり音源分離の評価実験の結果をそれぞれ示す．音楽データセットとしては，SiSEC20151 で“Professionally-produced music recordings (MUS)”として公開されている “The Mixing Secret Dataset
100 (MSD100)” を用いた．データセットは異なるジャンルを含む 100 のヴォーカル付き楽曲で構成されている．混合信号に加え，歌手または楽器毎に分離された分離信号も公開されている．
実験データとして，下記の 3 楽曲のデータを用いた．
1 https://sisec.inria.fr/professionally-produced-musicrecordings/)

日本音響学会講演論文集

Fig. 2 Estimates of Uk,t for musicA1 for unsupervised audio source separation.
• musicA: AM Contra - Heart Peripheral (Electronic Dance Pop)• musicB: Actions - Devil’s Words (Power Pop)• musicC: Actions - South of the Water (PowerPop)
3 楽曲の分離信号は，全てベース，ドラム，その他の楽器，ヴォーカルの 4 トラックであった．
本実験では，musicA1，musicA2，musicB，musicCとして，それぞれ上記楽曲の一部を抽出した．抽出区間は，順番に 1:50-2:00, 0:50-0:60, 1:47-1:57, 1:071:17 である．各抽出区間は，大きな旋律の変化を避けて，人手で決定した．音楽データはモノラルでサンプ
リング周波数は 22.050kHz とした．短時間フーリエ変換におけるフレーム長は 32ms として，フレームシフト長は 16ms とした．NPCA における H, U の初期値には，通常の NMF による分解結果を使用した．
NMF における H, U の初期値には乱数を使用した．
正則化パラメータ w1 ，w2 はともに 0.2 とし，p = 1.0とした．
まず，教師なし音源分離において，提案手法 NPCAが分離信号を使用することなく周期的に発生する楽音を抽出できることを検証する．図 2 は musicA1 に対するアクティベーション Uk,t の推定結果を示す．基底数 K は 24 とした．図中下部の太枠内の 12 個の基底にのみ，周期性制約を適用した．以下，この 12 基底を周期基底と呼ぶ．α，β の推定は，周期基底のみで行った．図中で各基底のアクティベーションの左に付けられた楽器ラベルは，分離信号を人手で聴いて付加した．NPCA による分離結果は右側のラベルに示されている．周期的に発音される音と，非周期的に発音される音を分離して，主に周期基底において周期的に発音される音が抽出されていることが確認できる．周期的なリズムが想定されるドラムやベースの音は主に周期基底に，ランダムなリズムが想定されるボーカルは主にそれ以外の基底に割り当てられている．その他の楽器についても周期成分と非周期成分がそれぞれ周期基底，それ以外の基底に分類されていることが分かる．
図 2 の左側の楽器ラベルは，同じデータに対して通常の NMF を適用した場合の結果を示す．多くの基底のアクティベーションが，ドラムやベースの周期的なリズムに影響され，周期的でないその他の楽器など- 641 -

2016年3月

Table 1 SNR improvement (dB) under each sectionof music applying NMF and NPCA.
datamodel
bass drum
othervocal averagemusicA1 NMF
1.59
9.66 11.27 15.02
9.39
NPCA
2.35
9.61
11.14
14.55
9.41
musicA2 NMF
2.40
10.48 10.89 12.91
9.17
NPCA
3.19
10.44
10.71
12.79
9.28
musicB NMF
10.41 15.80
4.45
9.19
9.96
NPCA 10.56
15.87
4.40
9.17
10.00
musicC NMF
9.80
18.11
7.02
8.96
10.97
NPCA
9.77
18.16
6.51
9.20
10.91
average NMF
6.05
13.51
8.41 11.52
9.87
NPCA
6.47 13.52
8.19
11.43
9.90

においても，ドラムやベースのリズムに該当するアクティベーションの微小な増加がみられる．アクティ
ベーションの左側に示された楽器ラベルにおいても，多くの基底がベース（11 基底）とドラム（8 基底）に割り当てられており，ヴォーカルは部分的に 1 基底に抽出されているのみである．これに対して，提案手法 NPCA では，周期成分と非周期成分が分離され，4基底でヴォーカルの大部分を表現できており，分離性能も NMF と比べて向上していることが確認できた．
実験を通して，提案モデルが周期的にアクティベートされる基底を限られた基底数の中にまとめて扱うことができることが示され，この特長は精度と分離信号の解釈性とのトレードオフを解決するのに有効であることが示唆される．加えて，提案モデルを使用すれば，例えば図 1 のような形で α と β の値を参照することにより，各基底にどのような周期と位相を持つ周期性が推定されたのかを把握することもできる．
このようなリズム構造の把握は，自動採譜や楽曲理解のための分析において重要である．
続いて，教師あり音源分離において NPCA の性能を通常の NMF と比較する．教師あり音源分離の学習フェーズでは，学習データとして分離信号を用いて各基底スペクトルを事前に学習する．分離フェーズでは，テストデータとして混合信号を用いて事前に学習した基底スペクトルを用いて音源分離を行う．この時，基底行列 H の更新は行わず，U , α, β のみ最適化する．NPCA・NMF ともに全ての楽器に対して基底数を 6 とした．したがって，4 楽器による混合信号の分離時の基底数 K は 24 となる．定量的に分離性能の比較を行うために，signal-to-noise ratio (SNR)
を評価指標として用いる．SNR は∑2ω,t si,ω,tSN R(dB) = 10 log10 ∑(15)2
ω,t (si,ω,t − ŝi,ω,t )で定義される．ここで，si,ω,t は入力された楽器 iの分離信号の複素スペクトログラムを，ŝi,ω,t =∑HUω,k k,tYω,t ∑k∈iHω,kUk,t は Winner filter により得られる分k離信号の複素スペクトログラムを示す．SNR が高いほど，分解性能が高いことを示す．以下，本稿では，
SNR improvement，すなわち SNR と SNR-all の差を評価のために用いる．SNR-all は ŝi,ω,t の代わりに混合信号の複素スペクトログラムを用いて算出した
SNR の値である．表 1 は NMF と NPCA による音源分離結果から算出した SNR improvement の値の一覧である．結果から，特にベースやドラムなどの周期性が想定される楽器における分離性能の向上が確認できる．musicC においては，ヴォーカルパートにおいても SNR improvement の値が向上している．実

日本音響学会講演論文集

際，musicC のヴォーカルの旋律は抽出した 10 秒間で 2 回繰り返されており，提案法が多くの周期的なリズムを奏でる楽器に対して分離性能を向上できることが分かる．楽器によっては SNR improvement の値が低下している場合があるものの，周期性という強い制約を与えているにも限らず SNR improvementの平均値が NMF よりも僅かながら増加していることが確認でき，教師あり音源分離の枠組みにおいても提案手法が音源分離性能を向上できることが示せた．

5まとめ

本稿では，NMF の目的関数において，どの音源がいつ発音するかを示すアクティベーションに任意の周期性を仮定する項を付加し，その周期性形状を推定しながら最適化を行うアルゴリズムを導出した．教師あり・教師なし音源分離の両方で提案手法の有効性を評
価し，周期的に発生するイベントを抽出することにより音源分離性能が向上することを示した．今後は，特に教師あり音源分離において，各音源における周期性の強さに応じて w1 の値を自動的に調整するよう手法を改良したい．また，提案手法を可変基底 NMF[5] と融合させることにより，周期的に基底行列を変化させる手法を考案したい．

参考文献[1] Daniel D Lee and H Sebastian Seung, “Learning the parts of objects by non-negative matrix
factorization,” Nature, vol. 401, no. 6755, pp.
788–791, 1999.
[2] Paris Smaragdis and Judith C Brown, “Nonnegative matrix factorization for polyphonicmusic transcription,” in IEEE Workshop on
Applications of Signal Processing to Audio andAcoustics, pp. 177–180, 2003.
[3] Hirokazu Kameoka et al., “Constrained andregularized variants of non-negative matrix
factorization incorporating music-specific constraints.,” in ICASSP 2012, pp. 5365–5368,
2012.
[4] Tuomas Virtanen, “Monaural sound source
separation by nonnegative matrix factorization
with temporal continuity and sparseness criteria,” IEEE Transactions on Audio, Speech, andLanguage Processing, vol. 15, no. 3, pp. 1066–1074, 2007.
[5] Masahiro Nakano et al., “Nonnegative matrix factorization with markov-chained bases
for modeling time-varying patterns in music
spectrograms,” in Latent Variable Analysis andSignal Separation, pp. 149–156. Springer, 2010.
[6] Kazuki Ochiai et al., “Explicit beat structure
modeling for non-negative matrix factorizationbased multipitch analysis,” in ICASSP 2012,pp. 133–136, 2012.
[7] Rafii Zafar and Bryan Pardo. “Repeating pattern extraction technique (REPET): A simple
method for music/voice separation,” in IEEE
Transactions on Audio, Speech, and Language
Processing, Vol. 21, no.1, pp. 73–84, 2013.

- 642 -

2016年3月