1-7-3
スパース時系列正則化に基づく独立低ランク行列分析を用いた
ブラインド音源分離 ∗☆三井祥幹（東大），北村大地（総研大），高道慎之介（東大），猿渡洋 (東大)1

はじめに音源信号，観測信号および分離信号をそれぞれブラインド音源分離とは，複数の音源からなる混合信号を観測した際に，マイクや音源の位置等の事前情報を利用することなく，混合する前の音源信号を推定する技術である．観測チャネル数が音源数以上となる優決定条件では，音源信号間の独立性を利用して音源分離を行う独立成分分析（independent componentanalysis: ICA） [1] に基づく手法が盛んに研究されている．ICA に基づく音源分離では，時間周波数領域での線形時不変な瞬時混合モデルを仮定する．これは，チャネル間空間相関行列のランクが 1 になることから，ランク 1 空間モデル [2] と呼ばれる．このような仮定の下で，周波数ビンごとに独立な ICA を適用する
周波数領域 ICA（frequency-domain ICA: FDICA）
が提案されているが，分離後にパーミュテーション問題を解く必要がある．分離とパーミュテーション問題を同時に解く手法として，周波数ビンをまとめて 1 つのベクトルとして取り扱う独立ベクトル分析（independent vector analysis: IVA） [3] が提案されている．
一方で，観測信号のスペクトログラムを非負行列の積に分解する非負値行列因子分解（nonnegative matrix factorization: NMF） [4] に基づく音源分離も盛んに研究されており，多チャネルに拡張した MNMF
（multichannel NMF: MNMF） [5] や，IVA と NMFを統合したブラインド音源分離手法である独立低ランク行列分析（independent low-rank matrix analysis:ILRMA） [6, 7] などが提案されている．特に ILRMAは，分離信号の時間周波数構造（NMF モデル）と線形分離行列を同時に推定する手法であり，音楽信号において IVA や MNMF よりも頑健かつ高精度な音源分離を実現している．
本稿では，音源間の独立性のみを用いる ILRMA に対して，時間周波数構造に関する事前分布を新たに導入する．特に，音声信号や音楽信号を対象とし，時間周波数領域での時系列がスパースとなる正則化を新たに付与することで，従来の ILRMA よりも高精度な音源分離を目指す．スパース正則化には，ベイズ型ポストフィルタの一種である MOSIE（MMSE estimation
with Optimizable Speech model and InhomogeneousError criterion） [8] を利用し，その統計モデルパラメータを自己推定する枠組みも併せて提案する．

2

従来の音源分離手法

2.1 表記の定義
観測チャネル数および音源数を，それぞれ M, N とおく．また，時間周波数領域における複素多チャネル∗sij = (sij,1 , sij,2 , . . . , sij,N )T(1)

T

xij = (xij,1 , xij,2 , . . . , xij,M )

(2)

T

(3)

yij = (yij,1 , yij,2 , . . . , yij,N )とおく．ここで，i = 1, . . . , I は周波数ビン，j =1, . . . , J は時間フレーム，n = 1, . . . , N は音源，m =1, . . . , M は観測チャネルのインデクスをそれぞれ表し，T は行列またはベクトルの転置を表す．
短時間フーリエ変換の窓長が各音源から各マイクまでのインパルス応答長より十分長い場合，ランク 1
空間モデルが成立し，観測信号を次式で表せる．

xij = Ai sij(4)

ここで，Ai = (ai,1 , . . . , ai,N ) は混合行列であり，ai,nは各音源のステアリングベクトルである．M = N となる場合は，分離行列 Wi = (wi,1 , . . . , wi,M )H が存在し，分離信号は

yij = Wi xij(5)と書ける．ここで H は行列またはベクトルの共役転置を表し，Wi = A−1
i の関係にある．

2.2 NMF
NMF は，非負値の要素のみで構成される行列を，非負の 2 つの行列の積として分解する手法であり，行列の低ランク構造を抽出できる．NMF では，非負行列DをD ≈ TV

(6)として近似的に分解する．ここで，T および V は非負行列である．音響信号に対しては，D をパワースペクトログラムとして分解する．このとき，T の各列にはパワースペクトログラム中の頻出スペクトルパタンを表す基底ベクトルが得られ，V の各行には各基底ベクトルの時間的な強度変化を表すアクティベーションが得られる．
NMF の目的関数は，観測 D とモデル T V のダイバージェンスで定義される．

QNMF = D(D | T V )(7)

ここで，D(· | ·) はダイバージェンスであり，行列の要素ごとの二乗ユークリッド距離，一般化カルバック・ライブラー情報量，板倉斎藤擬距離などが用いられ，その最適化手法として乗算更新型の反復アルゴリズムが提案されている [4, 9]．

2.3 ILRMAIVA は FDICA の独立なスカラー変数をベクトル変数に拡張したモデルであるが，ILRMA は IVA のベクトル変数をさらに低ランク行列変数（NMF 変数）に拡張したアルゴリズムである．また，MNMF [5] にBlind source separation using independent low-rank matrix analysis with sparse regularization for timeseries activations, by Yoshiki Mitsui (The University of Tokyo), Daichi Kitamura (SOKENDAI), Shinnosuke Takamichi (The University of Tokyo), Hiroshi Saruwatari (The University of Tokyo).

日本音響学会講演論文集- 325 -2016年9月

対してランク 1 空間モデルを仮定した理論と等価なため，Rank-1 MNMF とも呼ばれる [7]．
ILRMA の目的関数は次式で表される [7]．ただし，
M = N およびランク 1 空間モデルを仮定している．
[∑ ∑|y |2∑ ij,nQILRMA =− 2 log | det Wi |k znk tik vkj
ni,j]
∑∑
+log
znk tik vkj(8)nk

ここで，znk は n 番目の音源信号に対する k 番目の基底の寄与を表す隠れ変数であり，0 ≤ znk ≤ 1 かつ∑n znk = 1 が成り立つ．tik および vkj は，NMF の基底行列 T およびアクティベーション行列 V の非負要素である．式 (8) の第一項と第二項は IVA の目的関数に一致し，第一項と第三項は板倉斎藤擬距離に基づく NMF（Itakura-Saito NMF: ISNMF）の目的関数に一致する．文献 [7] では式 (8) の効率的な最適化手法が提案されている．

Fig. 1

Process flow of proposed method.の最適化を図る．

3.1 動機ILRMA は分離信号の音源モデル（時間周波数モデル）を T と V で捉えながら分離行列 Wi を推定する．
しかし，より複雑な時間周波数構造を持つ音源信号に対しては，NMF の音源モデル推定に失敗し，分離精度が上がらないという問題がある．そこで「音源信号の時間周波数構造の時系列（スペクトログラム行列上の任意周波数ビンにおける行ベクトル）はスパースになる」という事前知識を踏まえ，ILRMA の音源モデル T V の時系列がスパースになる正則化を付与することで，音声や音楽等の複雑な構造を持つ音源信号の分離を促進する．ここで，一般の音響信号に関しては，各周波数ビンにおける時系列のスパース度合いはその周波数位置ごとに大きく異なる．よって，アクティベーション行列 V を直接スパースに誘導するのではなく，各周波数ビンのスパース度合いをMOSIE の事前分布の形状母数として推定する手法も提案する．

3.3 処理の流れ
Figure 1 は提案手法の処理の流れを示している．まず，従来の ILRMA を数回反復更新し，分離行列 Wi ，基底行列 T およびアクティベーション行列 V を得る．この基底行列 T のスペクトル形状を補償した上で，n 番目の音源に寄与する基底ベクトルのみで構成される部分基底行列 Tn を構成する．次に，観測信号
xij と ILRMA の暫定的な分離信号 yij を用いて周波数ビンごとの適切な形状母数推定および MOSIE を適用し，スパース性を強調した n 番目の分離信号 Ỹn を得る．このとき，反復初期においては，推定される形状母数値よりもさらに強い度合いで強制的にスパース性を誘導するためのハイパーパラメータ η を用いることで，音源分離を促進する．ここで得た観測信号 Ỹn に対して，基底行列 Tn を固定した ISNMF を適用し，スパース正則化された分離信号に対応するアクティベーション行列 V̂n を得る．これを ILRMAの内部モデル V に戻し，再び ILRMA を数回反復更新する．この手順を繰り返すことで，ILRMA の音源モデル T V の時系列がスパースになるように誘導し，音源分離を促進する．以降では，各ステップにおける処理の詳細を説明する．

3.2 定式化
提案手法は，以下に示す二段階最適化として定式化できる．ここでは，正則化を表す制約式 (10) も最小化問題を含むのが特徴である．

3.3.1 ILRMA の更新
ILRMA の乗算更新式 [7] に基づき，分離行列の Wi ，基底行列 T およびアクティベーション行列 V を数回更新する．

3提案手法

min

QILRMA()∑s.t. min
DIS |Ỹn |2 | Tn Vn
Wi ,T ,V ,ZV(9)(10)n

ここで，DIS (· | ·) は板倉斎藤擬距離であり，行列 Xに対する演算 |X|2 は行列の要素ごとの絶対値の 2 乗を表す．また，Y˜n は，時系列のスパース性が強調された n 番目の分離信号のスペクトログラムであり，3.3.3
節で詳しく述べる．Tn および Vn は，それぞれ n 番目の音源信号に対する寄与が大きい基底とアクティべーションのみで構成される行列であり，3.3.4 節で詳しく述べる．この二段階最適化を直接解くのは困難であるが，本稿では文献 [10] に倣い，目的関数の最適化更新を数回行い，制約式の最適化更新を数回行う，といった流れを繰り返し行うことにより，全体

日本音響学会講演論文集

3.3.2 基底ベクトルの補償
ILRMA は ICA と同様に分離信号のスケールを決められないため，分離行列 Wi および基底 tik は周波数ごとにスケールがばらばらになっている．そこで，観測信号の m 番目のチャネルのスケールに補償した(m)
基底 tik を次式のように計算する．ただし，本稿では簡略化のため，n = m の場合のみを考える．
(n)tik = |ai,nn |2 znk tik(11)

ここで，ai,nn は ILRMA で推定した分離行列 Wi の逆行列から得られる．

3.3.3 MOSIE による時系列スパース性の強調
MOSIE はベイズ型ポストフィルタの一種であり，目的音源信号および妨害音信号に対してそれぞれ事- 326 -

2016年9月

前分布を仮定し，妨害音の分散の推定値を利用して目的音源信号の強調および妨害音の抑圧を行うフィル
タである [8]．MOSIE では，各周波数ビンにおいて，目的音の振幅スペクトルの事前分布としてカイ分布を仮定する．ここで，i 番目の周波数ビンのある時系列に対し，χ ∈ R≥0 を目的音の振幅スペクトルを表す確率変数とすると，カイ分布の確率密度関数は次式で与えられる．
()ρi()2ρiρi
2ρi −12p(χ) =χexp −χΓ(ρi ) E[χ2 ]E[χ2 ](12)
ここで，ρi は i 番目の周波数ビンにおけるカイ分布の形状母数，E[·] は期待値を表す．
まず，式 (12) で表される時系列統計モデルに関して，そのスパース性を制御する周波数別形状母数 ρiの決定方法について述べる. 3.3.2 節と同様に m = nの場合を考えるとき，n 番目のチャネルにおける観測信号 xij,n のうち，ILRMA の出力から推定される妨害音成分を qij,n とおく．本研究においては，ρi の推定値を以下の式で算出する．

ρi = η · {µ22 (|xij,n |) + µ22 (|qij,n |)− 2µ2 (|xij,n |)µ2 (|qij,n |)}· {µ4 (|xij,n |) − µ4 (|qij,n |) − µ22 (|xij,n |)− 2µ2 (|xij,n |)µ2 (|qij,n |) + 3µ22 (|qij,n |)}−1 (13)
ここで，µd (·) は時間フレーム
d 次のモー∑ j に関するdメントを表し，µd (αij ) = ( j αij)/J である．また，
η は後述するハイパーパラメータを表す．式 (13) の右辺から η を除いた因子は，現状の分離音から得られる形状母数の empirical 推定値であり，モーメントマッチングによって各周波数別に最尤推定されたものである [11]．ただし，ILRMA の反復初期においては音源分離が不十分な状態であるため，中心極限定理により，この推定値は過大推定されていることに留意する（これは，音源のスパース性を過小評価していることと等価である）．このため，ハイパーパラメータ
0 < η ≤ 1 を乗ずることによって，上述の empirical推定値を補正し，スパース性を強調する．この η は，反復処理が進行するにつれて音源分離の精度が向上することを考慮し，次第に 1 へと収束する簡単な数列を適当に定めればよい．本稿ではこれを η = l/Lと与える．ここで，l は提案手法の現在の反復回数を，L は総反復回数を表す．
前述の ρi が与えられるとき，目的音複素スペクト
ル ỹij,n の最適な推定値は以下のように書ける [8]．
√ξij
µ2 (|qij,n |)·ỹij,n = xij,n ·|xij,n |2
ρi + ξij)
 (ξij|xij,n |2Φ 1 − β2 − ρi , 1; − µ2 (|q|)ρ+ξij,niij)(·|xij,n |2ξijΦ 1 − ρi , 1; − µ2 (|qij,n |) ρi +ξij
] β1
Γ( β2 + ρi )·(14)
Γ(ρi )
ここで，β はドメインパラメータ，Γ(·) はガンマ関数，
Φ(e, f ; g) = 1 F1 (e, f ; g) は第 1 種合流型超幾何関数である．ξij は事前 SN 比であり，以下のように定義さ日本音響学会講演論文集れる．

ξij = κ|ỹi(j−1),n |2µ2 (|qij,n |)[]|xij,n |2+ (1 − κ) max− 1, 0µ2 (|qij,n |)(15)

また，κ は忘却係数である．
以上より，全ての周波数 i に関して MOSIE を通すことで，時系列のスパース性が強調された分離信号の振幅スペクトログラム |ỹij,n | が得られる．これを要素に持つ n 番目の分離信号の振幅スペクトログラ
ム行列を Ỹn と定義する．

3.3.4 アクティベーション行列の再推定
基底インデクスの集合 Kn を次のように定義する．
}
{
1(16)
Kn = k | znk >N
Kn は n 番目の分離信号に寄与する基底インデクスを(n)表す集合である．T (n) （要素に tik を持つ行列）のk ∈ Kn 列目の基底ベクトルから構成される部分行列を Tn とすると，これは n 番目の音源のみを表現する部分基底行列となる．同様に，n 番目の音源のみを表現する部分アクティベーション行列 Vn も得られる．
時系列のスパース性が強調された分離信号 Ỹn に対して式 (10) のように Tn を固定した ISNMF を適用することで，Ỹn のアクティベーション行列 V̂n を再推定できる．この NMF の更新式は，文献 [9] より次式で与えられる．
√∑∑−22′ tik ′ ,n v̂k ′ j,n )i tik,n |ỹij,n | (v̂kj,n ← v̂kj,n∑∑ k−2′′
i tik,n (k′ tik ,n v̂k j,n )(17)
ここで tik,n および v̂kj,n は，それぞれ行列 Tn およびV̂n の各要素である．

3.3.5 スパース正則化の反映
全音源に対する V̂n を統合して得られるアクティ
ベーション行列 V̂ を，ILRMA の音源モデル中の Vに代入することで，MOSIE によるスパース正則化をILRMA に反映できる．ただし，基底インデクスの順番が V と一致するように V̂n を統合する必要がある．

4評価実験

4.1 実験条件
提案手法の有効性を確認するため，音声と音楽の混合信号の音源分離実験を行った．比較対象法は，
IVA および従来の ILRMA である．観測信号として，
SiSEC [12] に収録されている音声および音楽信号を，
RWCP database [13] に収録されている E2A インパルス応答（残響時間 300 ms）と畳み込むことで，擬似的な 2 チャネル 2 音源の信号を作成した．音源信号の到来方向として，正面を 90◦ として，70◦ と 110◦ ，50◦ と 130◦ ，50◦ と 70◦ の三種の組み合わせについて検証した．分離性能の客観評価尺度には signal-todistortion ratio（SDR） [14] の改善量を用いた．提案
手法は，ILRMA を 20 回反復するごとに一度 MOSIEによる正則化を行うように設定した．また，式 (17)の ISNMF の更新回数は 100 回とした．その他の条件- 327 -

2016年9月

Table 1

謝辞 本研究は，総合科学技術・イノベーション会
議による革新的研究開発推進プログラム（ImPACT）およびセコム科学技術振興財団の助成を受けた．

Experimental conditions

Sampling frequency
FFT length
Window shift
Number of bases K
Number of ILRMA iterations L
Domain parameter βForgetting factor κ

16000 Hz
256 ms (4096 points)
64 ms (1024 points)
20200
1.00.98

参考文献[1] P. Comon, “Independent component analysis, A
new concept?” Signal processing, vol. 36, no. 3, pp.
287–314, 1994.
[2] N. Q. K. Duong, E. Vincent and R. Gribonval,“Under-determined reverberant audio source separation using a full-rank spatial covariance model,”IEEE Trans. ASLP, vol. 18, no. 7, pp. 1830–1840,
2010.
[3] T. Kim, H. T. Attias, S.-Y. Lee and T.-W. Lee,“Blind source separation exploiting higher-orderfrequency dependencies,”
IEEE Trans. ASLP,
vol. 15, no. 1, pp. 70–79, 2007.
[4] D. D. Lee and H. S. Seung, “Algorithms for nonnegative matrix factorization,” Proc. NIPS 2000,pp. 556–562, 2001.
[5] 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.
[6] H. Kameoka, T. Yoshioka, M. Hamamura, J. Le
Roux and K. Kashino, “Statistical model of speech
signals based on composite autoregressive systemwith application to blind source separation,” Proc.
LVA/ICA 2010, pp. 245–253, 2010.

Fig. 2 Comparison of signal separation performance among IVA, conventional ILRMA and proposed method, extracting (a) instrumental signaland (b) speech signal.
およびパラメータ設定を Table 1 に示す．
文献 [7] では，znk を 0 または 1 に固定する基底数固定型 ILRMA と，znk を同時に推定する基底数適応法型 ILRMA の二手法が提案されている．本稿の実験では，前者の基底数固定型 ILRMA を用いている．

4.2 実験結果
Figure 2 は，異なる初期値に対して 20 回の試行を行った際の平均 SDR 改善量を示している．いずれの音源配置に関しても，提案法は従来 IVA と比較して3 dB から 13 dB 程度，従来 ILRMA と比較して 1 dBから 2 dB 程度の SDR の改善を達成しており，より高精度な音源分離が実現されている．以上より，提案するスパース時系列正則化が音源分離に関して有効に機能していることが分かる．

5

おわりに

本稿では，より高精度なブラインド音源分離手法を目的とし，音源の独立性だけを用いる ILRMA に対して，分離信号の時間周波数の時系列がスパースになるような正則化を施す手法を提案した．正則化には，ベイズ型ポストフィルタの MOSIE を用いており，周波数ビンごとのスパース性を同時に推定する．
音声と音楽の混合信号を用いた実験では，提案手法の分離精度の向上を確認した．

日本音響学会講演論文集[7] D. Kitamura, N. Ono, H. Sawada, H. Kameoka andH. Saruwatari, “Determined blind source separation unifying independent vector analysis and nonnegative matrix factorization,” IEEE/ACM Trans.
ASLP, vol. 24, no. 9, pp. 1622–1637, 2016.
[8] C. Breithaupt and R. Martin, “Analysis of thedecision-directed SNR estimator for speech enhancement with respect to low-SNR and transientconditions,” IEEE Trans. ASLP, vol. 19, no. 2, pp.
277–289, 2011.
[9] M. Nakano, H. Kameoka, J. Le Roux, Y. Kitano, N. Ono and S. Sagayama, “Convergenceguaranteed multiplicative algorithms for nonnegative matrix factorization with β-divergence,” Proc.
IEEE MLSP 2010, pp. 283–288, 2010.
[10] H. Nakajima, D. Kitamura, N. Takamune,S. Koyama, H. Saruwatari, N. Ono, Y. Takahashiand K. Kondo, “Music signal separation using supervised NMF with all-pole-model-based discriminative basis deformation,” Proc. EUSIPCO, 2016.
[11] Y. Murota, D. Kitamura, S. Nakai, H. Saruwatari,
S. Nakamura, Y. Takahashi and K. Kondo, “Musicsignal separation based on bayesian spectral amplitude estimator with automatic target prior adaptation,” Proc. ICASSP 2014, pp. 7490–7494, 2014.
[12] S. Araki, F. Nesta, E. Vincent, Z. Koldovský,
G. Nolte, A. Ziehe and A. Benichoux, “The2011 signal separation evaluation campaign(SiSEC2011): – Audio Source Separation –,” Proc.
LVA/ICA 2012, pp. 414–422, 2012.
[13] S. Nakamura, K. Hiyane, F. Asano, T. Nishiura
and T. Yamada, “Acoustical sound database in realenvironments for sound scene understanding andhands-free speech recognition,” Proc. LREC 2000,pp. 965–968, 2000.
[14] E. Vincent, R. Gribonval and C. Févotte, “Performance measurement in blind audio source separation,” IEEE Trans. ASLP, vol. 14, no. 4, pp.
1462–1469, 2006.

- 328 -2016年9月