一般社団法人 電子情報通信学会
THE
INSTITUTE
OF ELECTRONICS,
社団法人
電子情報通信学会
INFORMATIONAND OFCOMMUNICATIONENGINEERS
THE INSTITUTE
ELECTRONICS,

信学技報
IEICE Technical Report
信学技報
EA2017-56 (2017-10)
TECHNICALREPORT OF IEICE.

INFORMATION AND COMMUNICATION ENGINEERS［招待講演］独立低ランク行列分析に基づくブラインド音源分離
北村

大地†

小野

順貴††

宏††† 亀岡澤田弘和††† 猿渡洋†

† 東京大学 〒 113–8656 東京都文京区本郷 7–3–1†† 国立情報学研究所 〒 101–8430 東京都千代田区一ツ橋 2–1–2††† 日本電信電話株式会社 〒 619–0237 京都府相楽郡精華町光台 2–4
あらまし 本稿では，新しいブラインド音源分離手法である独立低ランク行列分析（ILRMA）を提案する．ILRMA
は，従来の有名なブラインド音源分離手法の独立ベクトル分析（IVA）と，非負行列を対象とした低ランク近似理論の非負値行列因子分解（NMF）の二つを融合したアルゴリズムであり，独立成分分析を起源とする統計的独立性に基づく信号分離理論の正当な拡張手法である．また，本手法は NMF を多次元観測信号用に拡張した多チャネル NMF（MNMF）の特殊形とも解釈できることから，近年の最先端音源分離アルゴリズム（IVA 及び MNMF）の関係性を，

ILRMA を通して理論的に明らかにすることができる．
キーワード ブラインド音源分離，独立成分分析，非負値行列因子分解

Blind Source Separation Based on Independent Low-Rank Matrix Analysis
Daichi KITAMURA† , Nobutaka ONO†† , Hiroshi SAWADA††† , Hirokazu KAMEOKA††† , andHiroshi SARUWATARI†† The University of Tokyo 7–3–1 Hongo, Bunkyo-ku, Tokyo 113–8656, Japan†† National Institute of Informatics 2-1-2 Hitotsubashi, Chiyoda-ku, Tokyo 101–8430, Japan††† Nippon Telegraph and Telephone Corporation 2–4 Hikaridai, Seika-cho, Kyoto 619-0237, Japan
Abstract In this paper, we propose a new eﬃcient algorithm for blind source separation problem (BSS) called
independent low-rank matrix analysis (ILRMA). ILRMA is a uniﬁed method of independent vector analysis (IVA)
and nonnegative matrix factorization (NMF) and can be interpreted as a natural extension of the source separation
theory based on statistical independence. Also, ILRMA is identical to a special case of multichannel NMF (MNMF),
and this fact reveals the relationship between the state-of-the-art BSS techniques including IVA and MNMF.
Key words blind source separation, independent component analysis, nonnegative matrix factorization
ICA（frequency-domain ICA: FDICA）[2] が提案された．し1. は じ め にかしながら，ICA は分離信号の順序に任意性があるため，周ブラインド音源分離（blind source separation: BSS）とは，波数ビン毎に推定される分離信号の順序を適切に並び替えるポ音源位置や混合系が未知の条件で観測された信号のみから混合

スト処理が新たに必要となる．このポスト処理はパーミュテー

前の音源信号を推定する音響信号処理技術である．優決定条
件（音源数 <
= 観測チャネル数）における BSS では，独立成分ション問題と呼ばれ，これまでに様々な解決法が提案されている [3]– [7]．やがて，パーミュテーション問題の解決と周波数毎分析（independent component analysis: ICA）[1] に基づく手の分離行列の推定を同時に達成する手法として，独立ベクトル法が歴史的に主流である．ICA は，音源間の統計的独立性と分析（independent vector analysis: IVA）[8]– [10] が提案され音源信号の非ガウスな生成モデルを仮定することで，混合行列

た．IVA は ICA を多変量に拡張した理論であり，各音源の周の逆系である分離行列を推定する理論である．音響信号では混波数成分を一つにまとめた周波数ベクトルの生成モデル（多変合系が部屋の残響に起因する畳み込み混合となることから，短量確率分布）を仮定する．このとき，球対称な性質を持つ非ガ

時間 Fourier 変換（short-time Fourier transform: STFT）でウスな多変量分布を仮定することで，同一音源の周波数成分間得られる複素スペクトログラムの各周波数ビンに対して ICAの共変関係（高次相関）[10] を考慮でき，パーミュテーション問を適用することで周波数毎の分離行列を推定する周波数領域

題を避けながら分離行列を推定できる．また，ICA や IVA の—1—- 73 -

This article is a technical report without peer review, and its polished and/or extended version may be published elsewhere.

Copyright ©2017

by

IEICE

最適化に対して補助関数法 [11] を用いた高速かつ安定なアルゴ
リズムも提案されている [12], [13]．

この混合モデルは，時不変混合系の残響時間が STFT の窓長よりも十分短い場合に成立する．このとき，M = N かつ Ai が正一方，モノラル信号等を対象とした劣決定条件（音源数 > 観測
チャネル数）下では，非負行列の低ランク近似理論である非負値行列因子分解（nonnegative matrix factorization: NMF）[14]則であれば，分離ベクトル wi,n = (wi,n1 · · · wi,nM )T で構成

H
N ×N
される分離行列 A−1が存在
i ≈ Wi = (wi,1 · · · wi,N ) ∈ Cし，分離信号は次式で与えられる．

を応用した音源分離手法 [15]– [17] が注目を集めている．NMF

yij = Wi xijに基づく手法では，混合信号のパワースペクトログラム等を非(5)負行列とみなして少数の非負パーツに分解することで，音源間ここで，H はエルミート転置を示す．優決定条件 BSS では，式のスペクトルの違いを用いた音源分離とが可能となる．しか(5) 中の分離行列 Wi を全周波数において推定することが最終しながら，分解された非負パーツを音源毎に選別することは的な目標となる．本稿では，以後常に決定的な系（M = N ）を

一般に容易ではない．多チャネルの観測信号に対しては，音源考える．観測チャネル数が音源数より多い場合は，事前に xij

毎の空間的な伝達系（空間相関行列 [18]）の違いを用いて非負に対して主成分分析を適用し M = N とする方法が用いられる．

パーツを音源毎に選別することで分離を実現する多チャネル

2. 2 FDICA と IVA

NMF（multichannel NMF: MNMF）[19]– [21] が提案されていFigure 1 に，FDICA 及び IVA の音源モデルの違いを示す．

る．MNMF は，劣決定条件における高精度な音源分離を実現

FDICA [2] は残響による畳み込み混合を分離するために，周波しているが，最適化変数の初期値に対する頑健性や反復計算に

数ビン i 毎の複素時系列信号 xi1,m · · · xiJ,m に対して ICA を適用する．このとき，音源信号の生成モデル（統計的な音源モおける膨大な計算コストが解決すべき問題となっている．
本稿では，優決定条件における BSS を対象として，IVA とNMF を融合した新しい音源分離手法の独立低ランク行列分析

デル）として音源間で独立かつ非ガウスな分布 p(s) ≈ p(y)（複素ラプラス分布等）を仮定する．各周波数での分離行列 Wi は（independent low-rank matrix analysis: ILRMA）[22], [23] を他の周波数に関わりなく推定されるため，正しい分離信号を得．ILRMA は，FDICA や IVA の自然な拡張であるためには，FDICA の後に周波数毎の分離信号の順序をそろ

（注 1）

提案する

り，周波数毎の分離行列をより高精度に推定するアルゴリズム

えるパーミュテーション問題を解く必要がある [3]– [7]．

である．また，MNMF に空間的な近似仮定を導入したモデル

IVA [8]– [10] は，パーミュテーション問題に対するエレガン

が ILRMA と等価であることを理論的に示し，近年の最先端音トな解法である．FDICA のように周波数毎に独立な ICA を適源分離アルゴリズム（IVA 及び MNMF）の関係性を俯瞰する．

用するのではなく，全周波数成分を一つにまとめた周波数ベク

音声及び音楽信号を用いた実験により，ILRMA は従来手法よトル x̄j,m = (xij,m · · · xIj,m )T ∈ CI×1 を考え，このベクトルの時系列信号 x̄1,m · · · x̄J,m に対して ICA を適用する．このりも高速かつ高精度な音源分離を達成することを示す．

とき，音源周波数ベクトル s̄j,n = (sij,n · · · sIj,n )T の生成モ

2. 従 来 手 法

デルには，I 次元の非ガウスかつ球対称（確率値 p(s̄) がノルム

2. 1 定 式 化
音源数と観測チャネル数をそれぞれ N 及び M とし，各時間∥s̄∥ にのみ依存する性質）な多変量分布 p(s̄) ≈ p(ȳ) を仮定する．例えば，次式の球対称ラプラス分布がよく用いられる [10]．

 �

�
�
�
�∑ � yij,n �2�� p(ȳj,n ) = ∏exp −�� σi,n �πσi,n周波数における音源信号，観測信号，分離信号をそれぞれ

sij = (sij,1 · · · sij,N )T ∈ CN ×1

xij = (xij,1 · · · xij,M )T ∈ CM ×1
T

yij = (yij,1 · · · yij,N ) ∈ C

N ×1

1

(1)

i

(2)

(6)

i

ここで，σi,n はスケールパラメタである．この分布は球対称性(3)を持つため，同一ベクトル内の成分（同一音源の全周波数成分）と表す．ここで，i = 1, · · · , I は周波数インデクス，j = 1, · · · , J

は時間インデクス，n = 1, · · · , N は音源インデクス，m =が高次相関を持つ．従って IVA は，基本周波数とその倍音のよ
うに同時に生起する周波数成分を一つの音源にまとめる傾向が

1, · · · , M はチャネルインデクスを示し，T は転置を表す．また，あり，パーミュテーション問題を回避しつつ分離行列を推定で及び Yn ∈ CI×J で表す．これらの行列の要素はそれぞれ sij,n ，∏各信号の複素スペクトログラム行列を Sn ∈ CI×J ，Xm ∈ CI×J ，xij,m ，及び yij,n に一致する．混合系が線形時不変であり，時間周波数領域での複素瞬時混合で表現できると仮定すると，周きる．以後，式（6）に基づく IVA をラプラス IVA と表記する．
音 源 周 波 数 ベ ク ト ル 間 の 独 立 性 p(ȳj,1 , · · · , ȳj,N ) =L = −2J（ai,n = (ai,n1 · · · ai,nM )T は各音源のステアリングベクトル）

xij = Ai sijp(ȳi,n ) を仮定すると，IVA の観測信号に対する負対数

尤度関数は次式のように得られる（注 2）．

波数毎の時不変な複素混合行列 Ai = (ai,1 · · · ai,N ) ∈ CM×Nが定義でき，観測信号を次式で表現できる．

n∑ilog | det Wi | +∑G(ȳj,n )

(7)

j,n

（注 2）： [10] を含む多くの文献で，式（7）の第一項の係数 2 が抜けているが，(4)

この係数は複素確率変数の変数変換（観測信号 xij から分離信号 yij への変換）
おけるヤコビアンに由来するものである（複素変数のヤコビアンは実数変数の場
合の二乗となる）．しかしながら，最適化の過程で分離信号が互いに無相関とな（注 1）：本手法は文献 [22], [23] では rank-1 MNMF と呼ばれている．

ると第一項は 0 になるため，この係数の間違いはあまり重要ではなくなる．

- 74 -—2—

Fig. 1 Source models assumed in (a) FDICA and (b) IVA, where x̃m and ỹn shows
time-domain signals of xij,m and yij,n , respectively.

こ こ で ，G(ȳj,n ) は コ ン ト ラ ス ト 関 数 と 呼 ば れ G(ȳj,n ) =− log p(ȳj,n ) で定義される．ラプラス IVA の場合は σi,n = 1 とおくと G(ȳj,n ) = log π+∥ȳj,n ∥2 となる（∥ · ∥2 は L2 ノルム）．

K は基底数であり，低ランク近似するために K ≪ min(I, J) と設定される．分解行列 T 及び V は，次式の最小化問題の解として推定される．

H
分離信号は yij,n = wi,nxij として最適化変数 Wi を含む点に注

min D (|B|.p ∥T V ) s.t. tik , vkj >= 0 ∀i, j, k

意する．式（7）の最適化には，補助関数法 [11] を用いた高速かつ安定な手法が提案されている [12], [13]．なお，FDICA やIVA のような独立性基準の音源分離では，分離信号 yij,n のスケールが推定できない．最終的な分離信号のスケールは，プロ
ジェクションバック法 [24] で復元できるため，ラプラス IVA のスケールパラメタは常に σi,n = 1 とおいてかまわない．
式（6）の球対称ラプラス分布以外の優ガウス分布を音源モデルに用いた IVA として，次式のような時変分散 rj,n を持つ原点対称複素ガウス分布に基づく手法が提案されている [25]．

p(ȳ1,n , · · · , ȳJ,n ) =
=

∏
j

∏
j

(9)

T ,V

ここで，D(·∥·) は 2 つの行列引数の要素間のダイバージェンスの総和，tik 及び vkj は行列 T 及び V の非負要素，k = 1, · · · , K

は基底インデクスを示す．行列 T V が低ランクである制約から，基底行列 T の列ベクトルは，|B|.p 中に頻出するスペクト

ルパタンを表し，アクティベーション行列 V の行ベクトルは各スペクトルパタンの時間的な強度変化を表す．これまでに様々なダイバージェンスの活用が提案されているが，板倉斎藤ダイ
バージェンスに基づく NMF（Itakura–Saito NMF: ISNMF）

p(ȳj,n )によるパワースペクトログラムの低ランク近似には，次に示す(∥ȳj,n ∥22
1
exp −
πrj,n
rj,n

)

特殊な生成モデルによる解釈が存在する [27]．

ISNMF のコスト関数は次式で定義される．

(8)

LISNMF =時変分散 rj,n は周波数に関しては共通の値である．式（6）と同様に，式（8）も ȳj,n のノルムのみに依存するため，球対称な分布である．以後，式（8）に基づく IVA を時変ガウス IVA と表記する．式（8）はガウス分布で構成されているが，分散 rj,nが時間的に変動することから，時間フレーム j に関する周辺分

i,j)|bij |p|b |p− log ∑
−1t v
t v
k ik kj
k ik kj

∑ ij

(10)

ここで，bij は B の複素要素である．今，K 個の複素スペクト
ログラム要素 cij1 , · · · , cijK の生成モデルを原点対称複素ガウス分布と仮定する．

(|cijk |2
1
exp −p(cijk ) =
πrijk
rijk

布は優ガウスとなり [26]，IVA の音源モデルとして活用できる．

2. 3 単一チャネル信号を対象とした NMF
単一チャネルの音響信号を対象とした NMF では，複素スペ

∑(

)

(11)

ここで，式（11）は時間周波数の各グリッド i, j に対して互いに独を別の二つの立と仮定している．式（11）は原点対称であるため，確率値は複非負行列 T ∈ R>= 0 （基底行列）及び V ∈ R>= 0 （アクティベー

素成分 cijk の位相によらず振幅 |cijk | あるいはパワー |cijk |2 にクトログラム B ∈ C

I×J

I×K

を非負化した行列 |B|

.p

K×J

ション行列）の行列積に分解する．ここで，行列に対する絶対値
記号は要素毎に絶対値を取った行列を表し，ドット付きの指数
は要素毎の累乗を表す．従って |B|.p は，p = 1 が振幅スペクト

ログラム，p = 2 がパワースペクトログラムに対応する．また，のみ依存する．従って，分散はパワーの期待値 rijk = E[|cijk |2 ]に対応する．また，行列 B の複素要素が bij =∑kcijk と成分分解できるとすると，ガウス分布の再生性より，B に関する生成モデルも原点対称複素ガウス分布で表される．

- 75 -

—3—

=

i,j

π

(

|bij |2
exp − ∑
r
r
k ijk
k ijk

∑1

)

(12)

Vector-activated timefrequency source model(a)
Arbitrary number of bases

この事実は，一般に複素数のパワー値の加法性は成り立たな
いが，上記の生成モデルが成立する場合はパワースペクトル

Frequency|cijk |2 の加法性が期待値の意味で保存されることを意味している（注 3）．従って，式（12）において rijk = tik vkj とおくことで，
2

パワースペクトル |bij | を NMF で∑

t v に加法分解する
k ik kj− log p(B) =i,j(log π + log∑
k|bij |2
tik vkj + ∑t v
k ik kj

)

Time

Time

Low-rank timefrequency source model
(b)

Basis

ことの妥当性が正当化される．このモデルの負対数尤度関数は

∑

Time

Frequencyi,j∏Frequency

Uniform variance over frequency binsp(bij )

Basis∏

Frequencyp(B) =

Time

Fig. 2 Structure of source models in (a) Time-varying Gaussian
IVA and (b) ILRMA, where grayscale depicts intensity ofvariance.

(13)

この事実に基づいて，本稿では，ISNMF を IVA の音源モデとなり，式（13）の最小化問題は式（10）のそれと等価である．
なお，式（12）に示す ISNMF の生成モデルは式（8）と同様に，

ルとして導入することを提案する．本手法では，音源信号の生成モデルとして式（12）を仮定し，その分散 rij,n は，Fig. 2∑

時間周波数に依存して分散が変動することから，時間フレーム（b）のように，ISNMF で rij,n =j に関する周辺分布は優ガウスになる．

な形で分離行列 Wi と同時にブラインドに推定される．ここで，til,n 及び vlj,n は各音源に対応する L 本の基底を持つ基底

3. 独立低ランク行列分析
3. 1 動t v
として低ランク
l il,n lj,nL
×J
行列 Tn ∈ RI×
及びそのアクティベーション行列 Vn ∈ RL
の
>
>
=0
=0

機

音声の BSS においては，ラプラス IVA や時変ガウス IVA は非負要素であり，l = 1, · · · , L は新たな基底のインデクスである．従って，本手法は時変ガウス IVA の自然な拡張であり，1

ブラインドにパーミュテーション問題を解きつつ分離行列 Wi

本の基底であった時変ガウス IVA の音源モデルを，L 本の基底を比較的高い精度で推定できる．しかしながら，式（6）や式で構成される低ランク行列に拡張していると解釈できることか（8）で定義される周波数成分間の高次相関は，Fig. 2（a）に示ら，独立低ランク行列分析（ILRMA）と呼ぶ．

3. 2 尤 度 関 数

すように，全周波数成分に対して一様な強度の共変関係しか表

ILRMA は ISNMF と同様に，次式の複素ガウス分布を音源現できない．実際の音響信号では，基本周波数とその倍音は強い高次相関を持つが，それ以外の成分や周波数差の大きい成分信号の生成モデルとして仮定する．

間の高次相関は弱くなる等，共変関係の度合いは時間周波数のp(ȳj,1 , · · · , ȳj,N ) =グリッド i, j に依存して大きく変動する．このような理由から，

IVA は特に楽器音信号等のように明確な時間周波数構造を持つ=

音響信号の分離には適していない．一方で NMF は，Fig. 2（b）に示すように，限られた数のスペクトルパタンから成る低ラン
クな時間周波数構造の表現に良く適合しており，離散的な音階や音符の構造を持つ音楽信号等を効率的にモデル化できる．

の時間周波数の分散 rij,n の推定値として用いるモデルベース
散値（時間周波数成分の共変関係の強度に関する構造）を考慮しながら分離行列 Wi を推定できる．さらに，正確な時間周波

L = −2J∑ilog | det Wi | +rij,n =

いが，式（11）のモデルでは E[|bij | ] =kE[|cijk | ] =∑kる．即ち，E[|c1 +c2 |2 ] = E[|c1 |2 ]+E[|c2 |2 ] が成立している．

rijk が成立す∑ ( |yij,n |2

i,j,n

rij,n

+ log rij,n

)

(15)

∑

t v
とおくと，ISNMF のコスト関数（9）におl il,n lj,n

ける bij = yij,n の場合に等しいことがわかる．但し，分散 rij,nの低ランク分解については，次式の 2 通りを考える．

rij,n =∑

til,n vlj,n(16)

znk tik vkj

(17)

l

rij,n =

（注 3）：複素数 c1 と c2 に対して |c1 +c2 |2 = |c1 |2 +|c2 |2 は一般に成立しな
2(14)で与えられる．式（15）の第一項及び第二項はモデルベー値）が分離行列の推定精度を向上させうることを示唆している．

∑)

ス IVA のコスト関数に対応し，第二項及び第三項は分散を数毎の分散 rij,n （即ち音源毎のパワースペクトログラム推定

2(|yij,n |21
exp −πrij,nrij,n式（14）に基づく観測信号の負対数尤度関数は式（7）よりの事前分布として用いている．また，時間と周波数両軸の分散

IVA [28] も提案されている．これらの手法では，与えられた分∏

ここで，rij,n は音源毎の分散であり，rij,n = E[|yij,n |2 ] である．

れた分離したい音源の時間的なアクティビティを時変分散 rj,nク法等何らかの手法で一時的に推定した分離信号を，式（12）p(ȳj,n )nn,i時変ガウス IVA を用いた手法 [25] では，ユーザから与えら変動を用いる手法として，スペクトル減算法やバイナリマス∏∑k

ここで，znk ∈ [0, 1] は基底行列 T 中の全音源を表現する K 本- 76 -—4—

STFT

Update demixing matrix so that estimated signals
are mutually independent and have low-rank
structure in time-frequency domain

Frequency
Frequency

Frequency-wisedemixing matrix

til,n 及び vlj,n の更新後は，推定分散を次式で更新する．

Low-rank approximation

Estimated signal

Observed signal

rij,n ←
Time

n

�
(∑
)−2
�∑
� i,j |yij,n |2 tik vkj
z ′t ′v ′k′ nk ik k j�(25)
znk ← znk
(∑
)−1
∑
t v
i,j ik kj

znk = 1 を満たす．式（16）を用いる場合，各音源の端に異なる場合は最適な基底数が音源毎に異なる．そこで，式（17）のように分割関数 znk を導入することで，K 本の基底を

ILRMA における音源分離の原理を Fig. 3 に示す．Wi 及び
NMF 音源モデル最適化の過程では，分離信号のパワースペクの時間周波数構造を共変関係として加味した分離行列 Wi を推定する．混合前の各音源のパワースペクトログラム |Sn |.2 が低

ここで，式（26）はランクに誘導することでパーミュテーション問題を避けつつ．

ICA や IVA における分離行列 Wi の最適化では，補助関数法に基づく手法が提案されており，従来の自然勾配法よりも高速かつ安定に収束することが示されている [12], [13]．式（15）同じであることから，分離ベクトル wi,n の更新式は反復射影法（iterative projection: IP）[13] を用いて次式で得られる．

wi,n ← wi,n

)

−1
−2
wi,n ← wi,n λ−1
n , yij,n ← yij,n λn , rij,n ← rij,n λn(20)

及び vlj,n に関する偏微分は ISNMF のコスト関数（9）の偏微[(IJ)−1t ′ v′
l′ il ,n l j,nt
i il,nt ′ v′
l′ il ,n l j,n∑znk λ−2n(32)nznk λ−2nznk ← ∑z ′ λ−2n′ n k n′

(33)

∑

i,j

|yij,n |2 ]1/2 等が適当である．これらの正規化はコスト関数（15）の値を変えない点に注意する．全変数の推定
後は次式のプロジェクションバック法 [24] を適用する．

ŷij,n = Wi−1 (en ◦ yij )分と同じである．従って，til,n 及び vlj,n は次の ISNMF にお

vlj,n

(31)

で あ る ．こ こ で ，λn は 任 意 の 正 規 化 係 数 で あ る が ，λn =(21)

ける乗算型反復更新式で最適化できる．

(30)

音源モデルの正規化は，分割関数がない場合

tik ← tik

される．分割関数の無い式（16）を用いた場合，式（15）の til,n

j

(29)

であり，分割関数を含む場合

音源モデルの更新式は，分割関数の有無に応じて二種類導出

�
(∑
)−2
�∑
� i |yij,n |2 til,nt ′ v′
l′ il ,n l j,n�
vlj,n ← vlj,n
(∑
)−1
∑

znk tik vkj

k規化を反復毎に行うことで，この問題を避ける．

ルである．分離ベクトル更新後は，分離信号を次式で更新する．

�
(∑
)−2
�∑
� j |yij,n |2 vlj,nt ′ v′
l′ il ,n l j,n�
til,n ← til,n
(∑
)−1
∑

∑

任意性が存在し，数値計算が不安定になる．そこで，次式の正1
ここで，en ∈ RN×{0, 1} は n 番目の要素が 1，他要素が 0 のベクトHxij
yij,n ← wi,nznk = 1 を保証するために計算される．

til,n ← til,n λ−2
n

(19)

−1
H
wi,n
Ui,n wi,n 2nと rij,n はいずれも分離信号 yij,n のスケールを決定できるため

(18)

j

(

∑

z ′t ′v ′k′ nk ik k j

以上より，これらの分離行列と音源モデルの更新式を交互にの Wi に関する偏微分は，ラプラス IVA の補助関数の偏微分と

wi,n ← (Wi Ui,n )−1 en

znk tikz ′t ′v ′k′ nk ik k j

反復することで，式（15）を最小化できる．しかしながら，Wi

3. 3 更新式の導出

1∑ 1
xij xH
ij
J
rij,n

i,n

rij,n ←

ランクは基本的に増加することから，ILRMA は分離信号を低Ui,n =

znk vkj

znk ，tik ，及び vkj の更新後は，推定分散を次式で更新する．

ランクであれば，混合信号のパワースペクトログラム |Xn |.2 の互いに独立となる分離信号を推定することができる．

j,n

�
(∑
)−2
�∑
� i,n |yij,n |2 znk tik
z ′t ′v ′k′ nk ik k j�
vkj ← vkj
(28)
(∑
)−1
∑

N 個の音源に適応的に振り分けることができる．

を低ランク行列としてモデル化しながら，そ

(26)

�
(∑
)−2
�∑
� j,n |yij,n |2 znk vkj′ znk′ tik′ vk′jktik ← tik �
(27)
(∑
)−1
∑

の混合信号等ように，音源毎のスペクトログラムのランクが極

トログラム |Yn |z ′t ′v ′k′ nk ik k jznkznk ← ∑z ′n′ n k基底数は共通の L 本となる．しかしながら，ボーカルとドラム

.2

(24)

び vkj の更新式を補助関数法で導出できる [23]．

の基底を，N 個の音源に振り分ける潜在変数（分割関数）であ∑

til,n vlj,nl

また，分割関数を含む式（17）を用いた場合は，znk ，tik ，及Time

Fig. 3 Principle of source separation based on ILRMA.

り，

∑

(34)

ここで，ŷij,n = (ŷij,n1 · · · ŷij,nM )T は推定された多チャネル(22)

分離信号であり，◦ は要素毎の積を示す．なお，一連の最適化

アルゴリズムの具体的な実装は [29] で公開されている．

3. 4 ILRMA，IVA，及び MNMF の関係性
(23)

MNMF [19]– [21] は NMF を多チャネル観測信号用に拡張した音源分離手法であり，特に Sawada et al. による定式化 [21]- 77 -—5—

Flexibleが最も一般化されている．本手法は，各時間周波数の瞬時空間

X̂ij =∑

Spatial model

相関行列 Xij = xij xH
ij を観測信号とし，次のモデルを仮定する．

(35)

rij,n Ri,n

ここで，Ri,n ∈ C

M×Mは音源毎の空間相関行列 [18] であり，音源からマイクロホンまでの伝達系及び空間的な拡散を表す時不変な特徴量である．また，分散 rij,n は時変な変動成分，即ち音源毎のパワースペクトログラムであり，式（17）の形で低ラIntroduce rank-1
spatial model

Restricted

n

MNMF

IVA

Restricted

Introduce NMF
source model

ILRMA

Source model

FlexibleFig. 4 Relationship between ILRMA, IVA, and MNMF fromviewpoint of flexibility of spatial and source models.

ンク近似される．MNMF では，X̂ij を多変量複素ガウス分布の共分散行列として多チャネル観測信号の生成モデル [18], [21]仮定）を導入したモデルが，本質的に ILRMA と等価であるこを考える．そのときの最尤推定は次式の最小化問題となる．

とを示している．しかしながら，MNMF は混合系を推定する

LMNMF =∑[(

−1log det X̂ij + tr Xij X̂ij

i,j

)]

のに対して ILRMA は分離系を推定する双対な関係となってお(36)

り，この違いが IP に基づく効率的な最適化を実現している．

MNMF では，音源毎の空間相関行列 Ri,n と全音源の低ランクなパワースペクトログラム T V を同時に推定し，T V を分割関数 znk で選別することで音源分離を達成する．空間相関行列

Ri,n はフルランク行列として推定されるが，ここに混合系（4）の近似モデルを導入すると，伝達系は全てステアリングベクト
ル ai,n で記述でき，空間相関行列はランク 1 行列

Ri,n = ai,n aHi,nILRMA，IVA，及び MNMF の関係性を Fig. 4 に示す．3. 1
節で述べたように，IVA のベクトル音源モデルを NMF による
低ランク行列音源モデルに拡張した手法が ILRMA であり，それは MNMF に対して混合系（4）を仮定したモデルと等価である．この事実は，独立に発展してきた IVA 及び MNMF の両理論が，実は密接に関連していることを初めて明らかにしている．

4. 音声及び音楽信号の音源分離実験(37)

4. 1 実 験 条 件に制約される．式（37）を式（35）に代入すると，

X̂ij =∑ILRMA による音源分離の妥当性を示すために，音声及び音楽信号を用いた音源分離実験を示す．比較手法として，補助関

rij,n ai,n aHi,n数法に基づくラプラス IVA [13]，Ozerov et al. の MNMF [20]，n= Ai Dij AH
iSawada et al. の MNMF [21]，分割関数無しの ILRMA，分割(38)

D = diag (rij,1 , · · · , rij,N )関数有りの ILRMA の 5 手法を用いた．音声信号の実験で用い(39)

るデータは，SiSEC2011 [30] の UND タスクで公開されていると表現でき，これをコスト関数（36）に代入すると，次式を得る．

LMNMF =
=

∑[

(

−H
−1
log det Ai Dij AH
i + tr Xij Ai Dij Aii,j∑[log (det Ai ) (det Dij ) (det Ai )∗及び 250 ms である．このデータセットは劣決定条件 BSS を想定しているため，本実験では常に 1 番目及び 2 番目の音声信号のみを用いて M = N = 2 となる混合信号を生成した．音楽信

i,j

(

H
−1
+ tr Wi−1 yij yij
Wij−H WiH Dij
Wi

=

∑[

log | det Ai |2 + log det Dij

i,j

(

−1
H
Dij
+ tr Wi Wi−1 yij yij= − 2J
+

∑

∑
i,j

i

[

log | det Wi |log∏

クで公開されている 4 楽曲を，RWCP データベース [31] に収録されている E2A（T60 = 300 ms）及び JR2（T60 = 470 ms）のインパルス応答で畳み込み混合して生成した．音声信号は音楽信号よりも時間周波数構造が高ランクなため，正確なモデル化にはより多くの基底数が必要であるが，ILRMA での音声信H−1
rij,n + tr yij yij
Dij

n号の実験で用いるデータは，同じく SiSEC2011 の MUS タス

)]

)]

(

∑(

実録音信号（dev1 及び dev2）である．このデータセットは男
性及び女性の 12 個の音声信号を含み，残響時間 T60 は 130 ms)]−1

号分離では，基底数が多いと分離性能が劣化することが確認さ]れている [23]．従って，音声信号の分離では，ILRMA の基底)数を L = 2 及び K = 4 とし，音楽信号の分離では L = 30 及びK = 60 とした．客観評価尺度には信号対歪み比（SDR）[32] の)改善量を用いた．その他の実験条件は [23] に記載されている．

を用いた．式（40）はここで，xij = Wi−1 yij 及び Wi = A−1i

6 に示す．また，Tables 1 及び 2 は，全観測信号の平均分離性

ILRMA のコスト関数（15）と一致することから，MNMF に能をまとめている．なお，これらの結果は変数の乱数初期値を= − 2J∑ilog | det Wi | +i,j,n|yij,n |2+ log rij,nrij,n4. 2 実 験 結 果(40)

音声信号及び音楽信号の典型的な分離結果例を Figs. 5 及び対して空間相関行列のランクが 1 となる制約（混合系（4）の変えて 10 回試行した際の平均及び標準偏差である．この結果- 78 -—6—

SDR improvement [dB]

16
14
12
10
8
6
4
2
0
-2
-4

Source 1

Table 1 Total average SDR improvements (dB) of speech signals

Source 2

Recording conditions(a)(rev. time andOzerov’s
MNMF

Sawada’s
MNMF

ILRMA w/opartitioning
function

ILRMA withpartitioning
function

SDR improvement [dB]SDR improvement [dB]partitioning

partitioning

functionfunction

MNMF

MNMF

130 ms and 1 m

2.98

1.35

0.68

11.91

4.88

130 ms and 5 cm

2.86

2.13

1.13

8.97

3.48

250 ms and 1 m

2.03

0.49

0.48

7.34

2.09

250 ms and 5 cm

2.43

0.91

0.47

6.43

1.91

Impulse(b)

Laplace Ozerov’s Sawada’s

partitioning

functionfunction

10.32

12.29

12.29

6.11

6.62

7.40

IVA

MNMF

MNMF

E2A

5.72

5.73

JR2

1.77

2.37

Laplace IVA

Laplace
IVA

Ozerov’s
MNMF

Sawada’s
MNMF

ILRMA w/opartitioning
function

20
18
16
14
12
10
8
6
4
2
0
-2

ILRMA withpartitioning
function

Source 1

Source 2(a)Laplace
IVA

Ozerov’s
MNMF

Sawada’s
MNMF

ILRMA w/opartitioning
function

ILRMA w/o ILRMA withpartitioning

response

Sawada’s MNMFILRMA withpartitioning function

where reverberation time is (a) 130 ms and (b) 250 ms.

16
14
12
10
8
6
4
2
0
-2
-4

ILRMA w/o ILRMA with

Table 2 Total average SDR improvements (dB) of music signals
Laplace
IVA

Fig. 5 Examples of average SDR improvements for female speech,

16
14
12
10
8
6
4
2
0
-2
-4

Laplace Ozerov’s Sawada’s
IVAmic. spacing)

Ozerov’s MNMFILRMA w/opartitioning function

19.8 sSDR improvement [dB]SDR improvement [dB]

16
14
12
10
8
6
4
2
0
-2
-4

24.3 s8.3 s233.9 s28.1 s

0

50

100
150
Iteration step200

Fig. 7 SDR convergence and examples of actual calculation time

ILRMA withpartitioning
functionwith impulse response E2A, where signal length is 18.6 s.

(b)

CPU を用いている．ラプラス IVA と ILRMA は MNMF よりも高速な収束を示しており，計算時間も比較的短い．Sawada

et al. の MNMF は，1 回の更新において IN 回のサイズ 2M
行列の固有値分解や IJ 回のサイズ M 行列の逆行列演算が必要であり，より多くの計算時間が必要となる．
Laplace
IVA

Ozerov’s
MNMF

Sawada’s
MNMF

ILRMA w/opartitioning
function

5. ま と め

ILRMA withpartitioning
functionFig. 6 Examples of average SDR improvements for music signal
with impulse response (a) E2A and (b) JR2.

本稿では，IVA と NMF を融合した新しい BSS の ILRMAを提案した．ILRMA は，IVA の音源モデルをベクトルから低ランク行列に拡張した手法であり，効率的な最適化とより高精より，ラプラス IVA はシンプルな音源モデルに起因して音楽信度な音源分離を達成している．また，MNMF に対して周波数号の分離には不向きであることが確認できる．2 種の MNMF

領域での複素瞬時混合仮定を導入したモデルと ILRMA が本質

は，低ランクな音源モデルを用いることで，音楽信号を高い精的に等価であることを示し，IVA 及び MNMF の二大音源分離

度で分離しているが，音声信号に対しては十分な分離精度が得理論の関係性を理論的に明らかにした．本稿に記載した内容以られていない．一方で，分割関数無しの ILRMA は，音声信号

外にも，観測チャネル数が十分多い状況における残響のモデルに対して高精度かつ頑健な分離を達成しており，分割関数有り化手法 [33]，適応的なスパース正則化を付与した ILRMA [34]，の ILMRA は音楽信号に対するさらなる性能改善に有効である

STFT の窓長に関する実験的な解析 [35]，生成モデルを複素Student’s t 分布へと一般化した ILRMA [36]，分離行列と音源ことが確認できる．
さらに，Fig. 7 は，各手法の反復更新に対する SDR 改善量の収束例と実際の計算時間例を示している．計算には MAT-

モデルの最適化速度を調整し初期値への頑健性を向上させるパ
ラメトリックな最適化手法 [37] 等が検討されている．

LAB 8.3（64 bit）環境で，Intel Core i7-4790（3.6 GHz）の- 79 -

謝辞 本研究は，セコム科学技術振興財団及び JSPS 科研費—7—

17H06572 の助成を受けた．
文

献[1] P. Comon, “Independent component analysis, a new concept?,” Signal processing, vol.36, no.3, pp.287–314, 1994.
[2] P. Smaragdis, “Blind separation of convolved mixtures in
the frequency domain,” Neurocomputing, vol.22, pp.21–34,
1998.
[3] S. Kurita, H. Saruwatari, S. Kajita K. Takeda and
F. Itakura, “Evaluation of blind signal separation method
using directivity pattern under reverberant conditions,”
Proc. ICASSP, vol.5, pp.3140–3143, 2000.
[4] N. Murata, S. Ikeda and A. Ziehe, “An approach to blindsource separation based on temporal structure of speech signals,” Neurocomputing, vol.41, no.1–4, pp.1–24, 2001.
[5] H. Sawada, R. Mukai, S. Araki and S. Makino, “A robust
and precise method for solving the permutation problem
of frequency-domain blind source separation,” IEEE Trans.
ASLP, vol.12, no.5, pp.530–538, 2004.
[6] H. Sawada, R. Mukai, S. Araki and S. Makino, “Convolutive blind source separation for more than two sources in thefrequency domain,” Proc. ICASSP, pp.III-885–III-888, 2004.
[7] H. Saruwatari, T. Kawamura, T. Nishikawa, A. Lee andK. Shikano, “Blind source separation based on a fastconvergence algorithm combining ICA and beamforming,”
IEEE Trans. ASLP, vol.14, no.2, pp.666–678, 2006.
[8] A. Hiroe, “Solution of permutation problem in frequency domain ICA using multivariate probability density functions,”
Proc. ICA, pp.601–608, 2006.
[9] T. Kim, T. Eltoft, and T.-W. Lee, “Independent vector
analysis: an extension of ICA to multivariate components,”
Proc. ICA, pp. 165–172, 2006.
[10] T. Kim, H. T. Attias, S.-Y. Lee and T.-W. Lee, “Blindsource separation exploiting higher-order frequency dependencies,” IEEE Trans. ASLP, vol.15, no.1, pp.70–79, 2007.
[11] D. R. Hunter and K. Lange, “A tutorial on MM algorithms.”
The American Statistician, vol. 58, no. 1, pp. 30–37, 2004.
[12] N. Ono and S. Miyabe, “Auxiliary-function-based independent component analysis for super-Gaussian sources,” Proc.
LVA/ICA, pp.165–172, 2010.
[13] N. Ono, “Stable and fast update rules for independent vector analysis based on auxiliary function technique,” Proc.
WASPAA, pp.189–192, 2011.
[14] D. D. Lee and H. S. Seung, “Learning the parts of objects
by non-negative matrix factorization,” Nature, vol. 401, no.
6755, pp. 788–791, 1999.
[15] H. Kameoka, M. Nakano, K. Ochiai, Y. Imoto, K. Kashino
and S. Sagayama, “Constrained and regularized variantsof non-negative matrix factorization incorporating musicspeciﬁc constraints,” Proc. ICASSP, pp.5365–5368, 2012.
[16] D. Kitamura, H. Saruwatari, K. Yagi, K. Shikano, Y. Takahashi, and K. Kondo, “Music signal separation based onsupervised nonnegative matrix factorization with orthogonality and maximum-divergence penalties,” IEICE Trans.
Fundamentals, vol. E97-A, no. 5, pp. 1113–1118, 2014.
[17] D. Kitamura, H. Saruwatari, H. Kameoka, Y. Takahashi, K. Kondo, and S. Nakamura, “Multichannel signal separation combining directional clustering and nonnegative matrix factorization with spectrogram restoration,”
IEEE/ACM Trans. ASLP, vol. 23, no. 4, pp. 654–669, 2015.
[18] N. Q. K. Duong, E. Vincent, and R. Gribonval, “Underdetermined reverberant audio source separation using a fullrank spatial covariance model,” IEEE Trans. ASLP, vol. 18,
no. 7, pp. 1830–1840, 2010.
[19] H. Kameoka, T. Yoshioka, M. Hamamura, J. Le Roux, and
K. Kashino, “Statistical model of speech signals based oncomposite autoregressive system with application to blindsource separation,” Proc. LVA/ICA, pp. 245–253, 2010.

[20] A. Ozerov and C. Févotte, “Multichannel nonnegative matrix factorization in convolutive mixtures for audio source
separation,” IEEE Trans. ASLP, vol. 18, no. 3, pp. 550–563,
2010.
[21] H. Sawada, H. Kameoka, S. Araki, and N. Ueda, “Multichannel extensions of non-negative matrix factorization
with complex-valued data,” IEEE Trans. ASLP, vol. 21,
no. 5, pp. 971–982, 2013.
[22] D. Kitamura, N. Ono, H. Sawada, H. Kameoka and
H. Saruwatari, “Eﬃcient multichannel nonnegative matrix factorization exploiting rank-1 spatial model,” Proc.
ICASSP, pp.276–280, 2015.
[23] D. Kitamura, N. Ono, H. Sawada, H. Kameoka, and
H. Saruwatari, “Determined blind source separation unifying independent vector analysis and nonnegative matrix
factorization,” IEEE/ACM Trans. ASLP, vol. 24, no. 9, pp.
1626–1641, 2016.
[24] K. Matsuoka and S. Nakashima, “Minimal distortion principle for blind source separation,” Proc. ICA, pp. 722–727,
2001.
[25] T. Ono, N. Ono, and S. Sagayama, “User-guided independent vector analysis with source activity tuning,” Proc.
ICASSP, pp. 2417–2420, 2012.
[26] N. Ono, “Auxiliary-function-based independent vector analysis with power of vector-norm type weighting functions,”
Proc. APSIPA ASC, 2012.
[27] C. Févotte, N. Bertin and J.-L. Durrieu, “Nonnegative matrix factorization with the Itakura-Saito divergence. Withapplication to music analysis,” Neural Computation, vol.21,
no.3, pp.793–830, 2009.
[28] A. R. López, N. Ono, U. Remes, K. Palom´’aki, and M. Kurimo, “Designing multichannel source separation based on
single-channel source separation,” Proc. ICASSP, pp. 469–
473, 2015.
[29] D. Kitamura, “Algorithms for independent low-rank matrix analysis,” [Online]. http://d-kitamura.net/pdf/misc/
AlgorithmsForIndependentLowRankMatrixAnalysis.pdf[30] S. Araki, F. Nesta, E. Vincent, Z. Koldovsky, G. Nolte,
A. Ziehe and A. Benichoux, “The 2011 signal separation evaluation campaign (SiSEC2011):-audio source separation,” Proc. LVA/ICA, pp.414–422, 2012.
[31] S. Nakamura, K. Hiyane, F. Asano, T. Nishiura and T. Yamada, “Acoustical sound database in real environments forsound scene understanding and hands-free speech recognition,” Proc. LREC, pp.965–968, 2000.
[32] E. Vincent, R. Gribonval and C. Fevotte, “Performancemeasurement in blind audio source separation,” IEEE
Trans. ASLP, vol.14, no.4, pp.1462–1469, 2006.
[33] D. Kitamura, N. Ono, H. Sawada, H. Kameoka, and
H. Saruwatari, “Relaxation of rank-1 spatial constraint in
overdetermined blind source separation,” Proc. EUSIPCO,
pp. 1271–1275, 2015.
[34] Y. Mitsui, D. Kitamura, S. Takamichi, N. Ono, and
H. Saruwatari, “Blind source separation based on independent low-rank matrix analysis with sparse regularization fortime-series activity,” Proc. ICASSP, pp. 21–25, 2017.
[35] D. Kitamura, N. Ono, and H. Saruwatari, “Experimental
analysis of optimal window length for independent low-rankmatrix analysis,” Proc. EUSIPCO, pp. 1210–1214, 2017.
[36] S. Mogami, D. Kitamura, Y. Mitsui, N. Takamune,
H. Saruwatari, and N. Ono, “Independent low-rank matrix
analysis based on complex Student’s t-distribution for blindaudio source separation,” Proc. MLSP, 2017.
[37] Y. Mitsui, D. Kitamura, N. Takamune, H. Saruwatari,
Y. Takahashi, and K. Kondo, “Independent low-rank matrix analysis based on parametric majorization-equalization
algorithm,” Proc. CAMSAP, 2017 (in press).

- 80 -—8—