Vol.2015-MUS-108 No.2
2015/8/31

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

非ガウス性モノラル音響信号に対する音源分離のための非負値行列分解と半正定値テンソル分解
吉井 和佳1,a)

糸山 克寿1,b)

後藤 真孝2,c)

概要：本稿では，非ガウス性モノラル音響信号に対して音源分離を行うのに適した，非負値行列分解 (NMF)と半正定値テンソル分解 (PSDTF) の新しい確率モデルを提案する．従来，複素スペクトルは平均を 0 とする複素ガウス分布に従うことを仮定するのが一般的であった．このとき，各フレームにおける複素ガ
ウス分布の分散行列（半正定値行列）を少数の基底分散行列（半正定値行列）の錐結合で表現するのがLog-Determinant ダイバージェンスに基づく PSDTF (LD-PSDTF) であり，半正定値行列を対角行列（対角成分はパワースペクトル密度なので非負ベクトル）に限定すると板倉・斎藤ダイバージェンスに基づく
NMF (IS-NMF) に帰着する．しかし，実際の音響信号は非ガウス性（多くの場合優ガウス性）を持つた
め，ガウス分布より裾の重い確率分布を尤度関数に用いることが望ましい．本研究では，複素ガウス分布
を極限形式として含む複素ガウススケール混合 (GSM) 分布を尤度関数とした PSDTF および NMF を提案する．具体的には，複素 GSM 分布のうちで複素 t 分布（特殊形として複素コーシー分布を含む）を尤度関数とした t-PSDTF を定式化し，基底とアクティベーションを同時に最尤推定するための効率的な乗法更新アルゴリズムを導出する．実験の結果，高品質な音源分離ができることを確かめた．

1. はじめに

我々は以前，NMF の自然な数学的拡張となる，半正定値テンソル分解 (positive semideﬁnite tensor factorization:

音楽情報処理分野において，多重音に対する音源分離に

PSDTF) [7] と呼ぶ新しい因子分解法を提案した．NMF とおける有用性から，非負値行列分解 (nonnegative matrix

同様，PSDTF もコスト関数の違いにより複数の変種が考

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

えられるが，Log-Determinant ダイバージェンスに基づく

NMF では，多重音の振幅あるいはパワースペクトログラ

PSDTF (LD-PSDTF) が音響信号解析に有効である．まム （非負値行列）を，周波数方向の基底スペクトルの集合

ず，連続時間上の入力音響信号は局所的に定常なガウス過（非負値行列）と各基底スペクトルに対応する時間方向の程に従うことを仮定する．これは，フーリエ変換は線形変音量変化の集合（非負値行列）との積に分解を行う．この換であるから，各フレームの複素スペクトルは連続周波数

とき，観測スペクトログラムと，得られた 2 つの非負値行上の複素ガウス過程に従うことを仮定することと等価であ

列の積である再構成スペクトログラムとの誤差を表すコスる．複素ガウス過程の任意の周辺分布は複素ガウス分布にト関数を最小化するため，効率的な乗法更新アルゴリズム

なることから，離散周波数上の複素スペクトルの尤度は複が導出されている [5]．また，コスト関数の最小化は，あ素ガウス分布で与えられる．このとき，各フレームにおける確率モデルの最尤推定（対数尤度の最大化）としても解る複素ガウス分布の分散行列（半正定値行列）を少数の基釈可能である [6]．NMF の結果を用いれば，多重音の複素底分散行列（半正定値行列）の錐結合で表現し，それら基底

スペクトログラムに対して，時間周波数ビンごとに独立に行列と非負の重みを同時に最尤推定するのが LD-PSDTF

ウィナーフィルタを適用することで，音源スペクトログラ

である．ただし，半正定値行列を対角行列（対角成分はパムを求めることができる（位相は同一となる）．

ワースペクトル密度なので非負ベクトル）に限定すると板倉・斎藤ダイバージェンスに基づく NMF (IS-NMF) に帰

1
2
a)
b)
c)

京都大学 大学院情報学研究科 知能情報学専攻
Yoshida-honmachi, Sakyo, Kyoto, 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 2015 Information Processing Society of Japan

着する．LD-PSDTF の結果を用いれば，多重音の複素ス
ペクトログラムに対して，周波数ビン間の相関を考慮した
ウィナーフィルタを適用することで，適切な位相をもつ音
源スペクトログラムを推定可能となった．

1Vol.2015-MUS-108 No.2
2015/8/31

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

IS-NMF や LD-PSDTF において，多重音である入力音
響信号がガウス過程（離散時間上ではガウス分布）に従う

2. 複素ガウススケール混合分布という仮定にはある程度妥当性が認められる．各音源信号

本章では，NMF と PSDTF に対する新しい確率モデル

はガウス過程に従わない場合でも，中心極限定理により，を導出するうえで必要なガウススケール混合 (Gaussian

十分に多くの音源信号が重畳することで得られる音響信号

scale mixutre: GSM) 分布について説明する．GSM 分布

はガウス過程に従うとみなせる．実際，単独発話の音声信とは，ガウス分布の分散行列のスケールに対してある事前

号はガウス分布よりもラプラス分布などの優ガウス分布で分布を考え，スケールを積分消去することで得られる複合

よく近似できることが知られており，独立成分分析では，

分布 (compound distribution) のことをさす [13]．もとにそのような各音源信号の非ガウス性が，ガウス分布に従う

なる分布を多変量複素ガウス分布とすれば，得られる複合

観測信号を分離するときの手がかりとなる [8]．

分布は多変量複素 GSM 分布となる．一般に，GSM 分布

ただし，音源信号の性質によっては，通常の中心極限定はガウス分布より裾が重く (heavy-tailed) ，外れ値に頑健理が成立せず，多重音である入力音響信号がガウス分布にな性質を持ち，具体例として t 分布や K 分布などが知られ従わない場合が存在する．例えば，音源信号がコーシー分ている．コーシー分布は t 分布の，ラプラス分布は K 分布

布（裾が非常に重く，平均も分散も定義されない）に従うの特殊な場合である．ガウス分布は，GSM 分布（あるい場合は，そのような音源信号が多数重畳しても依然としては t 分布や K 分布）の極限形式としての特殊形である．

コーシー分布に従う [9]．ガウス分布もコーシー分布も安複素 GSM 分布は，複素楕円対称 (complex elliptically

定分布 (stable distribution) の一種であり，それぞれ安定

symmetric: CES) 分布 [13] と呼ばれる，実用上重要なクラ分布の特性指数が 1 および 2 のときの特殊な場合である．

スの分布のうちで代表的なものである．そのため，CES 分

安定分布に従う確率変数の和は，同じ特性指数をもつ安定

布がもつ好ましい性質を多く引き継いでおり，周辺分布や分布に従うという性質があり，音源信号の重畳過程に対す

条件付き分布が解析的に計算可能である点で都合がよい．

る生成モデルを定式化するうえで都合がよい．

GSM 分布を用いると，非ガウス分布およびガウス分布の本研究では，入力音響信号がガウス分布より裾の重い

統一的な取り扱いが可能になる．

分布であるガウススケール混合 (Gaussian scale mixture:GSM) 分布に従うと仮定した場合の PSDTF（およびその特殊形としての NMF）の定式化に取り組む．多変量ガウス分布と同様に，多変量 GSM 分布は任意の周辺分布が多変
量 GSM 分布となるため，連続時間上の音響信号に対する確率過程として，ガウス過程の代わりに GSM 過程を考え

2.1 分布の導出
多変量 CGSM 分布の導出について説明する．まず，多変量複素ガウス分布に従う確率変数 x ∈ Cd を考える．

x | Σ ∼ Nc (x|0, Σ)(1)

ることができる．すなわち，連続周波数上では，複素スペク

ここで，Σ ∈ Cd×d は半正定値行列である．次に，x のアトルが複素 GSM 過程に従うことを仮定するのと等価であフィン変換で得られる新たな確率変数 z ∈ Cd×d を考える．

る．本稿では，まず，GSM 分布の中でもとりわけ重要性の高い t 分布に着目する．ガウス分布は自由度が無限の t 分

z =μ+

√

τx

(2)

布であり，コーシー分布は自由度が 1 の t 分布に相当する．

ここで，μ ∈ Cd×d は移動パラメータであり，τ ≥ 0 はス

t 分布に基づく確率モデルの最尤推定には，LD-PSDTF に

ケールパラメータである．このとき，ガウス分布は確率変対する乗法更新則を内包する，t-PSDTF のための乗法更数のアフィン変換に対して閉じていることから，z も多変新則を導出することができる．

量複素ガウス分布に従う．

提案手法では，入力音響信号中に含まれる外れ値に対して頑健な低ランク分解が可能になる．NMF や PSDTF のz | μ, Σ, τ ∼ Nc (z|μ, τ Σ)

(3)

基本的な考え方は，対象となるデータをできるだけ少数の

ここで，分散行列のスケールである τ に対して，任意の事意味のある基底の組み合わせで表現しようとするものであ

前分布（ガンマ分布，逆ガンマ分布，逆ガウス分布，一般

るから，入力データ中に突発的な雑音が含まれていたとし

化逆ガウス分布など）を考える．

ても，それを表現するためにわざわざ新たな基底を準備すτ ∼ p(τ |θ)ることは得策ではない．これまでに，回帰問題において，ガウス過程の代わりに t 過程を用いたり [10]，行列分解問
題において，雑音がガウス分布ではなく t 分布やラプラス分布に従うと仮定することで [11, 12]，外れ値に頑健なモデルが学習できることが報告されている．音楽音響信号には突発的な外れ値が多く，頑健性向上が期待できる．

c 2015 Information Processing Society of Japan(4)

ここで，θ は事前分布のパラメータである．ただし，スケー
ルの任意性を解消するため，E[τ ] = 1 とする．この事前分
布のもとで，式 (3) における τ を積分消去を行う．
 ∞p(z|μ, Σ, θ) =N (z|μ, τ Σ)p(τ |θ)dτ
(5)
0

2

Vol.2015-MUS-108 No.2
2015/8/31

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

ここで，p(z|μ, Σ, θ) が多変量複素 GSM 分布の確率密度関数となる．ただし，スケール τ の事前分布がデルタ分布である（τ の値を一意に定める）場合には，式 (3) で与え

ガウス分布
コーシー分布られる多変量複素ガウス分布が得られることから，多変量
複素 GSM 分布は多変量複素ガウス分布を極限形式として含むことが分かる．
多変量複素 GSM 分布の具体的な確率密度関数を数式で表す際には，量 Q(z) を定義しておくと便利である．
H−1

Q(z) = (z − μ) Σ(z − μ)(6)

これは，平均ベクトルが μ で，共分散行列が Σ であるようなベクトル z で表される一群の値に対するマハラノビス距
離の二乗となっている．

2.2 アフィン変換・周辺分布・条件付き分布
多変量複素 GSM 分布は多変量複素ガウス分布と同じくいくつかの好ましい性質を持つ．まず，多変量複素 GSM
分布はアフィン変換に関して閉じている．すなわち，ある多変量複素 GSM 分布に従う確率変数 z ∈ Cd をz ∼ GSMc (z|μ, Σ, θ)(7)とする．ここで，μ は位置ベクトル，Σ はスケール行列，θ

図 1 実数軸でスライスした（虚部を周辺化した）一変量複素 t 分布

2.3 多変量複素 t 分布
式 (4) の事前分布に逆ガンマ分布を与えた場合に，式 (5)で得られる多変量複素 GSM 分布は多変量複素 t 分布となる．参考のため標準 t 分布の確率密度関数を図 1 に示す．
具体的には，λ = τ −1 にガンマ事前分布を仮定する．
 ν νλ ∼ G λ ,(15)
2 2
ここで，G(a, b) は形状パラメータ a，逆スケールパラメー
タ b をもつガンマ分布を表す．このとき，式 (5) を具体的に計算すると以下を得る．
 ν ν
dλp(z|μ, Σ, ν) = Nc (z|μ, λ−1 Σ)G λ ,2 2
− 2d+νd
2
Γ( 2d+ν
2
−1
2 ) 2
=
1
+
|Σ|
(16)
Q(z)
ν
d
Γ( 2 ) (πν)ν

はその他のパラメータとする．このとき，任意の A ∈ Cm×d

ここで，μ ∈ Cd は位置ベクトル，Σ ∈ Cd×d はスケール行

および b ∈ Cm に対して以下が成立する．

列，ν > 0 は自由度である．ただし，ν < 1 では平均と分

Az + b ∼ GSMc (Az + b|Aμ + b, AΣAH , θ)(8)すなわち，確率分布のパラメータは μ → Aμ + b および

Σ → AΣAH のように変換される．

散が，ν < 2 では分散が定義されない（無限大に発散する）
ため，ν < 2 の多変量複素 t 分布に従う確率変数に対しては通常の中心極限定理は成立しないことに注意する．
多変量複素 t 分布は，自由度が ν = 1 のとき多変量複素

また，多変量複素 GSM 分布の周辺分布・条件付き分布

コーシー分布，ν → ∞ のとき多変量複素ガウス分布となはいずれも多変量複素 GSM 分布となる．まず，式 (7) にり，いずれも再生性を有する．すなわち，確率変数 z 1 お

おいて，確率変数 z ∈ Cd を 2 つの部分 z 1 ∈ Cd1 および

よび z 2 が以下の多変量複素ガウス分布z 2 ∈ Cd2 に分割し（ただし d = d1 + d2 ），位置ベクトルや
スケール行列も対応するように分割しておく．
     
 
 

z1
Σ11 Σ12z 1  μ1
∼ GSMc,, θ (9)z2z 2  μ2
Σ21 Σ22
このとき，z 1 および z 2 の周辺分布は以下で与えられる．

z 1 ∼ Nc (z 1 |μ1 , Σ1 )(17)

z 2 ∼ Nc (z 2 |μ2 , Σ2 )(18)に従うとすると，確率変数の和も多変量複素ガウス分布

z 1 + z 2 ∼ Nc (z 1 + z 2 |μ1 + μ2 , Σ1 + Σ2 )(19)

z 1 ∼ GSMc (z 1 |μ1 , Σ11 , θ)(10)に従う．同様のことが多変量複素コーシー分布についても

z 2 ∼ GSMc (z 2 |μ2 , Σ22 , θ)

(11)成立する．すなわち，確率変数 z 1 および z 2 が以下の多変

一方，z 2 が与えられた場合の z 1 の条件付き分布は以下で与えられる．

z 1 | z 2 ∼ GSMc (z 1 |μ1|2 , Σ1|2 , θ 1|2 )(12)

ここで，θ 1|2 は式 (4) で与えられる事前分布に依存して決まる量であるが，位置ベクトル μ1|2 およびスケール行列

Σ1|2 は常に以下で与えられる．
μ1|2 = μ1 + Σ12 Σ−1
22 (z 2 − μ2 )(13)

Σ1|2 = Σ11 − Σ12 Σ−1
22 Σ21(14)

c 2015 Information Processing Society of Japan

量複素コーシー分布z 1 ∼ Cc (z 1 |μ1 , Σ1 )(20)

z 2 ∼ Cc (z 2 |μ2 , Σ2 )(21)に従うとすると，確率変数の和も多変量複素コーシー分布

z 1 + z 2 ∼ Cc (z 1 + z 2 |μ1 + μ2 , Σ1 + Σ2 )(22)に従う．t 分布のうちで再生性を有するのはガウス分布とコーシー分布のみであり，一般の自由度 ν に対しては再生性は成立しないことに注意する．

3

Vol.2015-MUS-108 No.2
2015/8/31

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

観測 と再構成との距離を最小化するように と を反復最適化

ガウス分布
ラプラス分布

周波数
時間

図 2 実数軸でスライスした（虚部を周辺化した）一変量複素 K 分布

図 3

IS-NMF による混合音スペクトログラムの低ランク分解

2.4 多変量複素 K 分布
X ≈ W ×2 H = W H T = Y
def

式 (4) の事前分布にガンマ分布を与えた場合に，式 (5)で得られる多変量複素 GSM 分布は多変量複素 K 分布となる．参考のため標準 K 分布の確率密度関数を図 2 に示す．具体的には，τ にガンマ事前分布を仮定する．

τ ∼ G (τ |ν, ν)

ここで，×2 は 2 モード積 [15] であり，wk ∈ RM
+ および

hk ∈ RN+ はそれぞれ基底ベクトルおよび対応するアクティ
ベーションベクトルである．式 (25) は次式で書き直せる．

(23)

K

xn ≈
このとき，式 (5) を具体的に計算すると以下を得る．

Kc (x|μ, Σ, ν) = Nc (z|μ, τ Σ)G (τ |ν, ν) dτ
(24)



ν−d
2ν d
=
|Σ|−1
νQ(x)
Kν−d 2 νQ(x)d
Γ(ν)π
ここで，μ ∈ C は位置ベクトル，Σ ∈ C
d

d×d

(25)

はスケール

行列，ν > 0 は自由度である．多変量複素 t 分布と異なり，クトルの錐結合で近似することである．
観測ベクトル xn と再構成ベクトル y n との間の誤差D(x|y) を評価する尺度として，Kullback-Leibler (KL) ダ
イバージェンス [5] や板倉・斉藤 (IS) ダイバージェンス [16]がよく用いられる．本稿では，特に後者に着目する．

M

DIS (xn |y n ) =m=1

分布，ν → ∞ のとき多変量複素ガウス分布となるが，再本章では，モノラル音響信号の音源分離の観点から，半(26)

すなわち，NMF の目標は，非負ベクトルを少数の非負ベ複素 K 分布は，自由度が ν = 1 のとき多変量複素ラプラス

3. 複素ガウス分布に基づく非負値行列分解と半正定値テンソル分解

def

k=1

任意の ν について平均と分散が定義できる．また，多変量生性を持つのは ν → ∞ のときのみである．

hkn wk = y n

xnm
xnm
− log
−1
ynm
ynm


(27)

この値は常に非負であり，xn = y n のときのみ 0 をとる．
通常の距離尺度と異なり，対称性が成り立たない，すなわ
ち DIS (xn |y n ) = DIS (y n |nn ) であることに注意する．

3.2 コスト関数最小化としての LD-PSDTF

正定値テンソル分解 (PSDTF) [7, 14] のための新しい確

PSDTF は，半正定値行列 *1 の集合（テンソル）X =率モデルを提案する．まず，コスト関数最小化の立場か[X 1 , · · · , X N ] ∈ CM ×M ×N に 対 し て ，少 数 の 半 正 定ら，NMF と PSDTF の一般的な定式化を説明し，PSDTF

値 行 列 の 集 合（ テ ン ソ ル ）W = [W 1 , · · · , W K ] ∈が NMF の自然な数学的拡張になっていることを示す．

CM ×M ×K お よ び H = [h1 , · · · , hK ]T ∈ RK×Nの積+

NMF と同様に PSDTF にも様々なコスト関数が考えう

Y = [Y 1 , · · · , Y N ] ∈ CM ×M ×N で近似する．

るが，Log-Determinant (LD) ダイバージェンスに基づく

X ≈ W ×3 H = Y
def

PSDTF (LD-PSDTF) は，板倉・斉藤 (IS) ダイバージェ(28)

ンスに基づく NMF (IS-NMF) の数学的拡張となっている．

ここで，×3 は 3 モード積 [15] であり，W k ∈ CM ×M

ここで，コスト関数の最小化は，多変量複素ガウス分布でおよび hk ∈ RN
+ はそれぞれ基底行列およびアクティベー

与えられる尤度関数の最大化（確率モデルの最尤推定）とションベクトルである．式 (28) は次式で書き直せる．

等価であることを示す．

KXn ≈3.1 コスト関数最小化としての IS-NMF
×N
NMF は，非負値行列 X = [x1 , · · · , xN ] ∈ RMに対+

hkn W k = Y ndef

すなわち，PSDTF の目標は，半正定値行列を少数の半正定値ベクトルの錐結合で近似することである．

×Kの積 Y = [y 1 , · · · , y N ] ∈よび H = [h1 , · · · , hK ] ∈ RN+

*1

予め与えるものとし，X を低ランク行列 Y で近似する．
c 2015 Information Processing Society of Japan


(29)

k=1

×K
お
して，二つの非負値行列 W = [w 1 , · · · , wK ] ∈ RM+
×N
RM
で近似する（図 3）．ただし，K  min(M, N ) は+0

半正定値行列と非負値行列とは異なる概念であることに注意．あるエルミート行列 A = AH ∈ RM ×M が半正定値行列であるとは，A = V V H となる行列 V が存在すること，あるいは，Aの M 個の固有値が全て非負であることと同値である．

4

Vol.2015-MUS-108 No.2
2015/8/31

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

観測行列 X n と再構成行列 Y n との間の誤差 D(X n |Y n )合に: 離散フーリエ変換行列: フレーム における局所信号: フレーム における複素スペクトル

結の凸を評価する尺度として，von-Neumann ダイバージェン
ス [17] や Log-Determinant (LD) ダイバージェンス [18] が

分解

PSDTF
列
底行

基

NMFは半正定値行列の対角要素のみに着目(本来存在するはずの周波数成分間の相関は無視)考えられる．本稿では，特に後者に着目する．

 − M (30)
− log X n Y −1
DLD (X n |Y n ) = tr X n Y −1
n
n

解

分

に

合

結

凸

の

ル

[dB]

この値は常に非負であり，xn = y n のときのみ 0 をとる．

ト

ク

ベ

NMF

底

テンソルデータ: 局所的な共分散行列の集合 行列データ: パワースペクトルの集合 基(半正定値行列の集合)(非負値ベクトルの集合)

式 (25) と式 (28)，式 (26) と式 (29)，式 (27) と式 (30) を

図 4

IS-NMF と LD-PSDTF の比較

比較すると，PSDTF は NMF の自然な数学的拡張となって．ここで，半正定値行列を対角行いることが分かる（図 4）

このとき，音源信号が重畳して得られる観測信号も局所列に限定すると，対角成分には固有値が並び，それらはすべ

的に定常なガウス過程に従う．式 (31) に着目すると，複素て非負であることから，X n = diag(xn )，Y n = diag(y n )，ガウス分布の再生性から以下を得る．

W k = diag(wk) とすれば *2 ，PSDTF は NMF に帰着すx̃n | Y n ∼ Nc (x̃n |0, Y n )る．すなわち，NMF では PSDTF のように，M 個の要素
ただし，Y n =間の相関を考慮することができない．



(33)

k Y kn である．これが LD-PSDTF の確率モデルであり，対数尤度関数は以下で与えられる．

3.3 尤度最大化としての LD-PSDTF
モノラル音響信号の音源分離を題材として，確率モデル

log p(xn |Y n ) = −M log(π) − log |Y n | − xn Y −1
n xnの最尤推定の観点から LD-PSDTF を捉え直す．まず，観(34)

測される音響信号（混合音）の複素スペクトログラムをX̃ = [x̃1 , · · · , x̃N ] ∈ CM ×N ，k 番目の音源信号の複素ス
ペクトログラムを X̃ k = [x̃k1 , · · · , x̃kN ] ∈ CM ×N とする．
ここで，M は周波数ビン数，N はフレーム数である．観= − log |Y n | − tr X n Y −1nc

ここで，X n = xn xH
n0 とした．式 (34) を式 (30) と比較すると，X n は観測データで一定であるから，Y n に関する対数尤度の最大化は，LD ダイバージェンスの最小化と等価であることが分かる．

測した混合音が K 個の音源信号の瞬時混合であると仮定
すると，以下が成り立つ．
K

3.4 乗法更新アルゴリズム

k=1

LD-PSDTF のパラメータ W および H を求めるため，

K

x̃kn

X̃ k ⇔ x̃n =

X̃ =

(31)

k=1

ズムが提案されている [7, 14]．本稿では結果のみ記すと，我々の目標は，観測変数 X̃ を潜在変数 X̃ k に分解することである．しかし，これは不良設定問題であるので，何らかの制約を導入して最適解を同定する必要がある．
我々は，各音源信号は局所的に定常なガウス過程に従うと仮定する．すなわち，各フレーム n における局所的な連続信号は定常であり，時不変な共分散関数（カーネル）をもつガウス過程に従うと仮定する．このとき，ガウス過程の任意の周辺分布は多変量ガウス分布になることから，あるサンプリング周波数に従ってサンプルされた離散信号は多変量ガウス分布に従う．したがって，フーリエ変換はアフィン変換であることから，各フレームにおける複素スペ
クトル x̃kn は多変量複素ガウス分布に従う．

x̃kn | Y kn ∼ Nc (x̃kn |0, Y kn )補助関数法に基づく収束性の保証された乗法更新アルゴリ
乗法更新則は以下で与えられる．

W k ← W k Lk (LTk W k P k W k Lk )− 2 LTk W k
 12

−1
tr X n Y −1
n W kY n
hkn ← hkn
tr W k Y −1
n
1

(35)
(36)

ここで，Lk はコレスキー分解 Qk = Lk LTk で求まる下三角行列であり，P k および Qk は次式で求まる．
N

Pk =
n=1

hkn Y −1
Qk =
n

N

−1
hkn Y −1n X nY n

(37)

n=1

hkn の非負性と W k の半正定値性は自然に保たれているが，tr(W k ) = 1 を満たすよう，反復のたびに W k および(32)

hk のスケールを調整しておく．式 (35) および式 (36) において，半正定値行列を対角行列に限定すると，よく知られ

ここで，共分散行列を Y kn = hkn W k とした．すなわち，

た IS-NMF の乗法更新アルゴリズムに帰着する [19, 20]．

各音源信号において，W k で定まる基本的なダイナミクス
は定常であるが，そのスケール hkn のみが時変化すると考えることで，音源信号の非定常性を制限する．
*2

diag(a) はベクトル a を対角成分に持つ対角行列を意味する．

c 2015 Information Processing Society of Japan


3.5 LD-PSDTF に基づくウィナーフィルタリング
最終的に，式 (32) および式 (33) から，x̃n が与えられた
ときの x̃kn の事後分布は多変量複素ガウス分布になること

5Vol.2015-MUS-108 No.2
2015/8/31

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

が分かり，その平均と分散は次式で求めることができる．

E[x̃kn |x̃n ] = Y kn Y −1
n x̃n

(38)

V[x̃kn |x̃n ] = Y kn − Y kn Y −1
n Y kn(39)

適用するうえでは，振幅スペクトルに加法性が成立するという強い仮定を置く必要があった．しかし，実験的には，

KL-NMF の方が IS-NMF より優れた音源分離結果を示す
ことが多い [9]．このことは，音源信号の混合過程を数学

このウィナーフィルタリングでは，X̃ k の位相は X̃ の位相的に正しくモデル化することが，必ずしもよい分離結果にとは異なる点に注意する．各フレーム n ごとに周波数ビン

つながるとは限らないことを示している．

間の相関を考慮しながら一挙に分離を行うことで，高品質

実際の音響信号に対して適切な尤度関数を決定する問題

な分離が可能となる．一方，IS-NMF では式 (38) や式 (39)は，モデル選択の問題として定式化できる．将来的に，パにおける半正定値行列はすべて対角行列に制限されているラメータに事前分布を導入すれば，エビデンスに基づくベ

ため，各周波数ビン n, m ごとに独立に行われることにな

イズ的なモデル選択も可能になると考えられる．

り，X̃ k の位相は X̃ の位相と同一となる．

4. 複素ガウススケール混合分布に基づく非負値行列分解と半正定値テンソル分解
我々は，多変量複素ガウス分布の代わりに，より裾の重い

4.2 尤度最大化としての定式化
PSDTF を頑健化するため，式 (33) で与えられる尤度関数において，多変量複素ガウス分布の代わりに，裾の重い多変量複素 GSM 分布を用いることを考える．

多変量複素 GSM 分布を尤度関数に用いることで，外れ値x̃n | Y n ∼ GSMc (x̃n |0, Y n )に対して頑健な PSDTF（およびその特殊形である NMF）

(40)

K

を提案する．本研究では，多変量複素 GSM 分布のなかで

ここで，Y n =も，多変量複素 t 分布（多変量複素コーシー分布を特殊形

タ X 全体に渡る対数尤度は次式で与えられる．

k=1 hkn W k である．このとき，観測デー

として含む）を尤度関数に用いる場合に着目する．多変量N

複素コーシー分布は多変量複素ガウス分布と同様に再生

log p(X |Y) =log p(x̃n |Y n )(41)

性を持つため，多変量複素コーシー分布に基づく PSDTF

n=1

は，混合音の生成過程として妥当な解釈を与えることがで

ここで，p(x̃n |Y n ) は式 (40) で与えられる各フレームごと

きる．誌面の都合上，多変量複素 K 分布（多変量複素ラプの尤度である．この対数尤度の一般形を最大化する乗法更ラス分布を特殊形に含む）を尤度関数に用いた PSDTF に新則を導出することは困難であるため，以降，複素 t 分布

ついては項を改めて議論したい．

である場合について議論する．

4.1 複素ガウススケール混合分布の妥当性

4.3 多変量複素 t 分布に基づく t-PSDTF

一般に，どのような尤度関数が適切であるかは，音源信

多変量複素 t 分布を尤度関数に用いる場合，最大化すべ号やそれらの和である観測信号の性質によって異なる．例

き対数尤度（式 (41) 参照）は具体的に次式で与えられる．

えば，音源信号がガウス過程に従うならば，観測信号もガウス過程に従う．また，音源信号がコーシー過程に従うな
らば，観測信号もコーシー過程に従う．ガウス過程もコー
シー過程も t 過程の特別な場合であるが，音源信号が一般の t 過程に従う場合は，観測信号は t 過程に従うとはいえない．一般に，GSM 過程に従う音源信号の和は必ずしもlog p(X |Y)
(42)


N 

2
ν
c
log 1 + tr X n Y −1
− log |Y n | − M +
=
n
2
ν
n=1
c

ここで，= は定数項を除いて等号が成立することを示す．

同じ形の GSM 過程に従うとは限らない．中心極限定理に

式 (42) を Y に関して直接最大化することは困難である

より，理論的には，音源数が増加するほど観測信号はガウ

ため，補助関数法（A.1 節）を用いて式 (42) の下限関数をス過程に漸近する場合も考えられる．

設計し，その下限関数を最大化することで，式 (42) を間接

多変量複素ガウス分布と多変量複素コーシー分布に基づ
く PSDTF のみが「数学的に正しい」モデルであったとしても，より広いクラスの裾の重い分布に基づく PSDTF を的に最大化することを考える．

4.3.1 下限関数導出のための不等式
下限関数導出のために必要な不等式を説明する．まず，定式化しておくことは応用上重要である．これまで，音源

f (Z) = log |Z| が凹関数であることに着目すると，f (Z) に信号がガウス過程に従うという仮定のもとでは，音源信号

対して 1 次のテイラー展開を行うことで，次式を得る．

のパワースペクトル密度の加法性が成立するため，混合
音のパワースペクトログラムに対して IS-NMF を適用す

− log |Z| ≥ − log |Ω| − tr(Ω−1 Z) + M(43)

ることには一定の妥当性があることが知られていた．一

ここで，Ω は任意の半正定値行列（展開点）であり，M は Z

方，混合音の振幅スペクトログラムに対して KL-NMF をのサイズである．等号成立条件は，Ω = Z で与えられる．

c 2015 Information Processing Society of Japan

6

Vol.2015-MUS-108 No.2
2015/8/31

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

次に，任意の半正定値行列 A に対して g(Z) = tr(Z −1 A)

したがって，hkn の非負性と V k の半正定値性は自然に保は凸関数であることに着目すると，澤田らの提案する不等

たれているが，tr(W k ) = 1 を満たすよう，反復ごとに W k

式 [21] を適用可能である．

および hk をスケーリングしておく．

−tr




K
k=1 Z k

−1



A

K

≥−

式 (49) および式 (50) は，式 (35) および式 (36) の自然
Ttr Z −1k Φk AΦk(44)

な拡張となっていることが分かる．実際，多変量複素 t 分
布の自由度を ν → ∞ とすると，πn → 1 となり，多変量複k=1

ここで，{Z k }K
{Φk }K
k=1 は任意の半正定値行列であり，k=1 は素ガウス分布に基づく PSDTF (LD-PSDTF) に帰着する．

足すと単位行列になるような補助行列である ( k Φk = I)．

等号成立条件は，Φk = Z k ( k Z k )−1 で与えられる．

PSDTF が得られる．すなわち，提案する多変量複素 t 分



4.3.2 下限関数の導出
これら不等式を用いると，式 (42) に対する補助関数は以4.3.4 t-NMF の導出
多変量複素 t 分布に基づく t-PSDTF の特殊形として，複− log |Ωn | − tr(Ω−1n Y n) + M

n=1





2
−1
−1
−1
ψn + ψn 1 + tr X n Y n
ν



c
−1
−1
−1 2M
+ 1 tr X n Y n−tr(Ωn Y n ) − ψn
=
ν
n=1

N


ν
− M+
2

N

−tr(Ω−1n Y n)≥n=1− ψn−1

であるとして LD-PSDTF を行っていると解釈できる．た
反復ごとに更新されることに注意する．

log p(X |Y)
N≥

布に基づく t-PSDTF では，πn X n を疑似的な観測データ
だし，πn は現在の再構成データ Y n に依存しているため，

下の通り導出できる．

c

一方，ν = 1 とすると，多変量複素コーシー分布に基づく

2M
+1
ν

K


−1
T
tr h−1
kn W k Φkn X n Φkn(45)

素 t 分布に基づく NMF を導出しておく．3.2 節で説明した通り，半正定値行列を対角行列に限定する場合を考えばよい
ので，X n = diag(xn )，Y n = diag(y n )，W k = diag(wk )とすると，乗法更新則は以下で与えられる．

 12
2
nm xnm )hkn /ynmn (πwkm ← wkm
n hkn /ynm

 12
2
nm xnm )wkm /ynmm (πhkn ← hkn
m wkm /ynm

条件（下限関数の最大化条件）は以下で与えられる．

Ωn = Y n(46)

2
ψn = 1 + tr X n Y −1
n
ν

(47)

Φkn = Y kn Y −1
n

(48)

(54)

ここで，πnm は以下で与えられる．

k=1

ここで，Ωn ，ψn および Φkn は補助変数であり，等号成立

(53)

πnm =

2+ν
−1
2xnm ynm+ν(55)

このとき，疑似的な観測データ πnm xnm は次式となる．
−1

2
ν
−1
πnm xnm =
+
(56)
ynm
x−1
2+ν
2 + ν nm
これは，観測データ xnm と再構成データ ynm に対して，重み ν : 2 で調和平均をとることを意味し，自由度 ν が大きくなるほど，xnm を重視するようになる．反対に，ν が小

4.3.3 乗法更新則の導出
この下限関数を W k および hkn に関してそれぞれ偏微
分しゼロとおくことで，以下の更新式を得る．

W k ← W k Lk (LTk W k P k W k Lk )− 2 LTk W k
 12

−1
tr πn X n Y −1
n W kY n
hkn ← hkn
tr W k Y −1n1

さくなるほど，ynm を重視するようになる．

t-PSDTF では，疑似的な観測データ πn X n に対して，低ランクな再構成データ Y n をフィットさせることが目的(49)

であるが，すでに πn X n には Y n が取り込まれており，ある種の平滑化が行われていると見ることができる．この結(50)

果，観測データに対して過剰にフィットすることが抑制され，外れ値に対して頑健な分解ができると考えられる．

ここで，Lk はコレスキー分解 Qk = Lk LTk で求まる下三角行列であり，P k および Qk は次式で求まる．
NPk =hkn Y −1
Qk =
n

n=1

N

4.4 多変量複素 K 分布に基づく K-PSDTF−1
hkn Y −1n (πn X n )Y n (51)n=1

2M + ν
2tr(X n Y −1n )+νc 2015 Information Processing Society of Japan

べき対数尤度（式 (41) 参照）は，式 (24) を用いて表現できる．4.3 節同様に，補助関数法（A.1 節）を用いて対数尤

また，πn は次式で求められる．

πn =

多変量複素 K 分布を尤度関数に用いる場合，最大化す度の下限関数を設計し，その下限関数を最大化することを(52)

考える．このとき，第二種変形ベッセル関数 Kν (x) が x に関して対数凸であることを利用する．

7

Vol.2015-MUS-108 No.2
2015/8/31

情報処理学会研究報告
IPSJ SIG Technical Report
与えられた混合音

5. 評価実験
多変量複素 t 分布に基づく t-PSDTF を用いた予備的な

C4

E4

0

音源分離実験について報告する．

G4

2

C4+E4

4

観測行列

C4+G4

6

E4+G4

8

観測行列

C4+E4+G4

10

12

観測行列

14[s]

観測行列

5.1 実験条件
実験には，RWC 研究用音楽データベース：楽器音 [22]に収録されているピアノ (011PFNOM)，エレキギター
対数尤度関数を最大化するような基底行列の凸結合を見つける(131EGLPM) およびクラリネット (311CLNOM) の単独音

基底行列

基底行列

基底行列

を用いた．各楽器ごとに，異なる 3 つの音高 (C4, E4, G4)をもつ 2 秒間の音響信号を準備し，それらを 7 つの異なる組み合わせで重畳したもの (C4, E4, G4, C4+E4, C4+G4,

E4+G4, C4+E4+G4) を連結することで 14 秒の音響信号を合成した．サンプリング周波数は 16[kHz] とした．

128

128

128

256

256

256

384

384

C4の周期
C4の周期

512

128

256

384

512

音量変化

0.10

384

E4の周期
E4の周期

512

128

256

384

512

音量変化

0.14

G4の周期
G4の周期

512

128

256

384

512

音量変化

0.35

次に，t-PSDTF を用いて，与えられた混合音を C4, E4,

G4 に対応する音源信号に分離することを試みた．まず，ガT
ウス窓を用いて局所信号 {on }N
n=1 を切り出し，X n = on onとすることで観測データ {X n }N
n=1 を得た．窓幅は 512 点，

0

0

2

4 6

図 5

8 10 12 14

0

0

2

4 6

8 10 12 14

0

0

2

4 6

2

8 10 12 14 x10

時間領域 t-PSDTF を用いたピアノ信号の分解結果

シフト長は 160 点とした (M = 512, N = 1400)．ここ

的に避けられない．PSDTF はこの影響を考慮することでで，t 分布などの GSM 分布はアフィン変換に対して閉じ優れた音源分離性能を達成している．また，t-PSDTF にている (2.2 節参照）ことから，周波数領域だけではなく，対する乗法更新アルゴリズムが，ν の値に関わらず数値的

時間領域においても等価な分解が可能である [7]．このとにも安定して収束することを確認した．

き，ウィナーフィルタリングを行うと，各音源信号の複

PSDTF の主な課題は，計算コストが O(KN M 3 ) であ

素スペクトルではなく，時間信号が直接求まる．比較のり，NMF の O(KN M ) よりもはるかに大きいことである．

ため，パワースペクトログラムに対する t-NMF も評価し計算時間を短縮し，局所解を回避するため，実際には NMF

た．各手法に対して，基底数は K = 3，反復回数は 100

で PSDTF を初期値する，すなわちある程度収束が進むま回とし，乗法更新アルゴリズムを用いた．音源分離結果

で基底行列を対角行列に制限して反復更新を行う方法が推は，BSS Eval Toolbox [23] を用いて，source-to-distortion

奨される．

ratio (SDR)，source-to-interferences ratio (SIR) および

6. おわりにsources-to-artifacts ratio (SAR) で評価した．

本稿では，非ガウス性モノラル音響信号に対して音源分

5.2 実験結果

離を行うのに適した，NMF と PSDTF の新しい確率モデ

音源分離タスクにおいて t-PSDTF は t-NMF に対する明ルを提案した．実際の音響信号は非ガウス性（多くの場合

確な優位性を示した．SDR, SIR, SAR の平均は，IS-NMF

優ガウス性）を持つため，ガウス分布より裾の重い確率分（ν → ∞ とした t-NMF）では 19.1 dB, 24.0 dB, 21.0 dB で布である複素ガウススケール混合 (GSM) 分布を尤度関数あったのに対し，LD-PSDTF（ν → ∞ とした t-PSDTF）に用いることが望ましい．そこで、本稿では，複素 GSM

では 23.0 dB, 27.7 dB, 25.1 dB であった．また，コーシー

分布のうちで特に重要な分布である複素 t 分布（特殊形と尤度に基づく NMF と PSDTF（ν = 1 とした t-NMF として複素コーシー分布を含む）を尤度関数として，基底とt-PSDTF）はそれぞれ，IS-NMF と LD-PSDTF とほぼ同

アクティベーションを同時に最尤推定するための効率的な等の性能を示した．この予備実験で用いたデータは単純で乗法更新アルゴリズムを導出した．実験の結果，高品質なあったが，突発性の雑音が多く含まれるより現実的な音楽

音源分離ができることを確かめた．

音響信号に対する評価を今後実施予定である．

今後は，多変量複素安定分布を尤度関数とした PSDTF

図 5 に示す通り，PSDTF を用いると，減衰音と持続音の導出に取り組みたい．安定分布はガウス分布やコーシーのいずれに対しても基底行列 W および音量変化 H を適

分布を特殊な場合として含み，音源信号の重畳をモデル化

切に推定できた．ここで，各基底行列 W k における斜め縞

するうえで数学的に妥当性がある [9]．

の間隔は周期を表しており，W k の中心付近は巡回行列に謝辞: 本研究の一部は，JSPS 科研費 24220006, 26700020,

近くなっている．しかし，窓関数の影響で周辺部はそうは26280089，JST CREST OngaCREST，および栢森情報科学

ならないので，周波数成分間に相関が発生することは原理振興財団の支援を受けた．

c 2015 Information Processing Society of Japan

8

Vol.2015-MUS-108 No.2
2015/8/31

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

[5]

付

録

[6]

A.1 補助関数法
目的関数 F (θ) を変数 θ に関して最大化する問題を考え[7]

る．補助関数法 [4, 25] を用いると，F (θ) を「間接的に」最大化することができる．まず，F (θ) に対する補助関数
として，以下の下限関数 F + (θ, φ) を考える．

[8]

F (θ) ≥ F + (θ, φ)

[9]

(A.1)

ここで，φ は補助変数である．このとき，以下の反復更

[10]

新則
[11]

φnew ← argmaxφ F + (θ old , φ)(A.2)

θ new ← argmaxθ F + (θ, φnew )

(A.3)

[12]

を用いると，F (θ) は単調非減少となる．このアルゴリズ[13]

ムの収束性は保証されている.

A.2 確率密度関数
確率分布の確率密度関数を以下に整理しておく．

[14]

ガンマ分布

G(x|a, b) =[15]

ba a−1 −bx
x
e
Γ(a)

[16]

多変量複素ガウス分布

Nc (x|μ, Σ) =1
|Σ|−1 exp (−Q(x))
πd

[17]

多変量複素 t 分布[18]

tc (x|μ, Σ, ν)=
− 2d+νd
2
Γ( 2d+ν
2
−1
2 ) 2
1
+
|Σ|
Q(x)
ν
d
Γ( 2 ) (πν)ν

[20]

多変量複素 K 分布

[21]

Kc (x|μ, Σ, ν)
=


2ν d
−1
|Σ|
Γ(ν)π d

[19]

νQ(x)

ν−d


Kν−d 2

νQ(x)
[22]

参考文献
[1]

[2]

[3]

[4]

P. Smaragdis, C. Févotte, G. Mysore, N. Mohammadiha,and M. Hoﬀman. Dynamic source separation using nonnegative factorizations: A uniﬁed view. IEEE Signal
Processing Magazine, 31(3):66–75, 2014.
H. Kameoka and K. Kashino. Composite autoregressive
system for sparse source-ﬁlter representation of speech.
In ISCAS, pages 2477–2480, 2009.
K. Yoshii and M. Goto. Inﬁnite composite autoregressive models for music signal analysis. In ISMIR, pages
79–84, 2012.
M. Hoﬀman, D. Blei, and P. Cook. Bayesian nonparametric matrix factorization for recorded music. In ICML,

c 2015 Information Processing Society of Japan


[23]

[24]

[25]

pages 439–446, 2010.
D. Lee and H. Seung. Algorithms for non-negative matrix factorization. In NIPS, pages 556–562, 2000.
A. T. Cemgil. Bayesian inference for nonnegative matrix
factorisation models. Computational Intelligence and
Neuroscience, 2009:Article ID 785152, 2009.
K. Yoshii, R. Tomioka, D. Mochihashi, and M. Goto. Inﬁnite positive semideﬁnite tensor factorization for source
separation of mixture signals. In ICML, pages 576–584,
2013.
A. Hyvärinen, J. Karhunen, and E. Oja. Independent
Component Analysis. John Wiley & Sons, 2004.
A. Liutkus and R. Badeau. Generalized Wiener ﬁlteringwith fractional power spectrograms. In ICASSP, 2015.
A. Shah and A. G. Wilson Zand . Ghahramani. Student-t
processes as alternatives to gaussian processes. In AISTATS, pages 877–885, 2014.
B. Lakshminarayanan, G. Bouchard, and C. Archambeau. Robust Bayesian matrix factorisation. In AISTATS, pages 425–433, 2011.
N. Wang and D.-Y. Yeung. Bayesian robust matrix factorization for image and video processing. In ICCV,
pages 1785–1792, 2013.
E. Ollila, D. E. Tyler, V. Koivunen, and H. V.
Poor. Complex elliptically symmetric distributions:
Survey, new results and applications. IEEE TASLP,
60(11):5597–5625, 2012.
K. Yoshii, R. Tomioka, D. Mochihashi, and M. Goto. Beyond nmf: Time-domain audio source separation withoutphase reconstruction. In ISMIR, pages 369–374, 2013.
T. G. Kolda and B. W. Bader. Tensor decompositions
and applications. SIAM Review, 51(3):455–500, 2009.
C. Févotte, N. Bertin, and J.-L. Durrieu. Nonnegative
matrix factorization with the Itakura-Saito divergence:With application to music analysis. Neural Computation, 21(3):793–830, 2009.
K. Tsuda, G. Rätsch, and M. K. Warmuth. Matrix exponentiated gradient updates for on-line learning and Bregman projection. JMLR, 6:995–1018, 2005.
B. Kulis, M. Sustik, and I. Dhillon. Low-rank kernel learning with Bregman matrix divergences. JMLR,
10:341–376, 2009.
亀岡弘和. 非負値行列因子分解. 計測と制御, 51(9):835–
844, 2012.
亀岡弘和. 非負値行列因子分解の音響信号処理への応用.
日本音響学会誌, 68(11):559–565, 2012.
H. Sawada, H. Kameoka, S. Araki, and N. Ueda. Eﬃcientalgorithms for multichannel extensions of Itakura-Saito
nonnegative matrix factorization. In ICASSP, pages
261–264, 2012.
M. Goto, H. Hashiguchi, T. Nishimura, and R. Oka.
RWC music database: Music genre database and musicalinstrument sound database. In ISMIR, pages 229–230,
2003.
E. Vincent, R. Gribonval, and C. Févotte. Performancemeasurement in blind audio source separation. IEEETASLP, 14(4):1462–1469, 2006.
J. Le Roux, E. Vincent, Y. Mizuno, H. Kameoka, N. Ono,and S. Sagayama. Consistent Wiener ﬁltering: Generalized time-frequency masking respecting spectrogram
consistency. In LVA/ICA, pages 89–96, 2010.
M. Nakano, H. Kameoka, J. Le Roux, Y. Kitano, N. Ono,and S. Sagayama. Convergence-guaranteed multiplicative algorithms for non-negative matrix factorization
with beta divergence. In MLSP, pages 283–288, 2010.

9