2-R1-5
局所時間周波数構造に基づく深層パーミュテーション解決法の実験的評価 ∗☆山地修平, 北村大地 (香川高専)

1

はじめにブラインド音源分離（blind source separation:BSS）とは，複数の音源が混合した観測信号から，混合前の音源を推定する技術である．特に独立成分分析（independent component analysis: ICA） [1] の登場以降，音源信号数 ≥ マイクロホン数である優決定条件下の音源分離問題に広く適用されている．
一般的には，周波数領域での瞬時混合を仮定し ICAを周波数毎に適用する周波数領域 ICA（frequencydomain ICA: FDICA） [2] に基づく手法がよく用いられる．しかし，ICA は分離信号の順番が不定であり，FDICA は周波数毎に分離を行うため，周波数毎の分離信号の順番がバラバラになってしまう．これを正しい順番に並び替える問題は一般にパーミュテーション問題 [3] と呼ばれている．近年では FDICA に対して音源の時間周波数成分の共起関係を新たに仮定して，パーミュテーション問題を起こさず分離信号を推定する BSS が登場している．例えば，独立ベクトル分析（independent vector analysis: IVA） [4] や独立低ランク行列分析（independent low-rank matrixanalysis: ILRMA） [5, 6] が有名である．さらに，深層学習（deep neural networks: DNN）を用いて音源の時間周波数構造の仮定を学習し，FDICA に適用する独立深層学習行列分析 [7] も提案されている．しかしながら，上記いずれの手法を用いても完璧にパーミュテーション問題を解くことは難しく，とくに複数音声の混合信号では IVA や ILRMA 等でも分離に失敗してしまうことが多い．
文献 [8] では，局所的な時間周波数構造から分離信号の正しいパーミュテーションを予測する DNN を新たに提案し，FDICA が周波数毎に理想的な分離を行った状況では高精度な BSS が達成できることを示している．本稿では，上記の DNN に基づくパーミュテーション解決法を拡張し，FDICA の分離結果に対してより頑健な予測を達成できるアルゴリズムを新たに提案する．また，提案手法を実際に FDICA のポスト処理として用いて，従来の BSS と音源分離性能を比較する．提案するパーミュテーション解決法は，モノラル音源の推定値のみを必要とするため，劣決定条件の場合にも一般性を損なうことなく適用できる．

22.1 定式化
音源数と観測チャネル数（マイクロホン数）をそれぞれ N 及び M とし，各時間周波数における音声信号，混合信号，及び分離信号をそれぞれT

xij = (xij,1 , · · · , xij,m , · · · , xij,M )zij = (zij,1 , · · · , zij,n , · · · , zij,N )T∈ CNT

∈C∈CN

(1)
M(2)
(3)と表す．ここで，i = 1, 2, · · · , I ，j = 1, 2, · · · , J ，
n = 1, 2, · · · , N ，及び m = 1, 2, · · · , M はそれぞれ∗

xij = Ai sij(4)

以後，決定的な系（M = N ）を仮定する．Ai が正則であれば，分離行列 Wi = A−1∈ CN ×M を用いて，
i
分離信号を次式で表せる．

zij = Wi xij(5)

従って，音源分離は分離行列 Wi の推定問題となる．

2.2 パーミュテーション問題
FDICA は音源間の統計的独立性のみに基づいて分離行列を推定するため，分離音源のスケール及び順番に関しては任意である．従って，FDICA の推定分離行列を Ŵi とすると，次式のような不定性が残る．
Ŵi = Di Pi Wi(6)

ここで，Pi は任意のパーミュテーション行列，Di は任意の対角行列である．また，式 (8) の音源毎の複素スペクトログラム行列を Yn ∈ CI×J と定義する．すなわち，FDICA で推定される分離信号

yij = Ŵi xij(7)
T= (yij,1 , · · · , yij,n , · · · , yij,N )

∈ CN

(8)

は，推定音源の順番やスケールが周波数毎にばらばらになっている状態である．このうち，Di によって生じるスケールの任意性は，プロジェクションバック
法 [9] で復元可能である．一方で，Pi によって生じる分離信号の順番の任意性（パーミュテーション）を復元することは容易ではない．
パーミュテーション問題を解決して得られる分離信号は次式となる．

zij = Pi−1 Di−1 yij(9)

本稿では，Pi−1 を推定することが課題となる．

FDICA とパーミュテーション問題

sij = (sij,1 , · · · , sij,n , · · · , sij,N )周波数，時間，音源，チャネルのインデクスを示し，·Tは転置を表す．また，各信号の複素スペクトログラム行列を Sn ∈ CI×J ，Xm ∈ CI×J ，及び Zn ∈ CI×Jで表す．周波数毎の時不変な混合行列 Ai ∈ CM ×Nを定義し，混合信号が次式で表現できると仮定する．

3提案手法

3.1 動機文献 [10] では，BSS の短時間フーリエ変換（shorttime Fourier transform: STFT）おける最適な窓長を実験的に検討している．文献 [10] 中の Fig. 6 (b) より，IVA 及び ILRMA では，残響状態 T60 = 470 msの条件では音声を分離できないことがわかる．一方で，FDICA に対して，音源信号 si,j を用いる理想的なパーミュテーション解決法（ideal permutationsolver: IPS）を適用した結果では 10 dB 以上の信号対歪み比（signal-to-distortion ratio:SDR）[11] の改Experimental evaluation of deep permutation solver based on local time-frequency structure. By Shuhei
YAMAJI and Daichi KITAMURA (NIT Kagawa).

日本音響学会講演論文集- 265 -2020年9月

T T

r̃i,γ = (ri,γ,1 , ri,γ,2 )

(10)

×1
∈ R2τ
≥0(11)

DNN
outputs

Fig. 2

1 : Diff.

DNN

1 : Diff.

0 : Same

1 : Diff.

DNN

Frequency
Time

DNN

Input vector…

T

Input vector of DNN.

…

×1T
T
di,ω,γ = (r̃i,γ, g̃i,ω,γ)T ∈ R4τ
≥0Fig. 1…

3.2 DNN の入出力提案する DNN モデルの入力ベクトルを Fig. 1 に示す．観測された混合信号 Xn に FDICA を適用すると，パーミュテーション問題が生じた分離信号 Yn が得られる．|Yn |.2 から，2 つの周波数（i,i + ω ）の短時間時系列パワー（長さ τ ）を以下のように集める．

Time…

• 分離信号 Y1 及び Y2 から同じ 2 つの周波数の時系列パワーをそれぞれ抽出し DNN に入力する• DNN は入力された 2 つの時系列パワーが同一音源か否かを予測し 0 または 1 として出力する• Y1 及び Y2 の全時間及び全周波数に対して DNNが適用される• 最終的な推定値 Pi−1 は，予測値の時間方向及び周波数方向の多数決結果から決定される

Frequency

善を達成している．この事実は，高残響下での音声混合信号であっても，Ŵi は FDICA で正確に推定でき，Pi−1 の推定のみ失敗していることを示している．
そこで，文献 [8] において，著者らは DNN に基づく教師付きパーミュテーション解決法を提案し，FDICAが周波数毎に理想的な分離を行った状況では高精度な BSS が達成できることを示した．本稿では，音源数 N = 2 及びチャネル数 M = 2 を仮定したうえで上記のパーミュテーション解決法を拡張し，FDICAの分離結果に対してより頑健な予測を達成できるアルゴリズムを新たに提案する．この手法の概要は以下の通りである．

0 : Same

DNN predictions in subband frequency bins.

3.3 DNN の構造提案する DNN の構造は，入力層，隠れ層 6 層，及び出力層の計 8 層からなる全結合構成となっており，1〜5 番目の隠れ層には ReLU 関数，最終隠れ層にはsigmoid 関数を適用している．1 6 番目の隠れ層の次元数はそれぞれ 128, 128, 128, 64, 64, 1 である．予測結果 qi,ω,γ と正解ラベルとの誤差関数には，平均二乗誤差（mean squared error: MSE）を使用している．

ri,γ,n = (|yi,(γ−1)η+1,n | , |yi,(γ−1)η+2,n | ,2

23.4 サブバンド周波数帯域における DNN 推定
DNN の 予 測 例 を Fig. 2 に 示 す．こ こ で ，f,2τ ×1
TT T
1 f2 , · · · , f5 は サ ブ バ ン ド 内 の 周 波 数 で あ り，
g̃i,ω,γ = (gi,ω,γ,1 , gi,ω,γ,2 ) ∈ R≥0(13)簡 単 の た め に γ は 省 略 し て い る ．こ の 図 で は ，gi,ω,γ,n = (|yi+ω,(γ−1)η+1,n |2 , |yi+ω,(γ−1)η+2,n |2 ,参照周波数を i = f3 とし，その近傍周波数を
τ ×12 Ti + ω = f1 , f2 , · · · , f5 と定義している．Y1 の参照周· · · , |yi+ω,(γ−1)η+τ,n | ) ∈ R≥0(14)波数 f3 の成分は赤色の音源成分であり，Y1 の f2 ，ここで，行列の |·|.2 は，要素ごとの二乗絶対値を返す． f4 ，及び f5 は青色の音源成分である．この内，2 本
また，ω = −Ω, −Ω+1, · · · , −1, 0, 1, · · · , Ω は，ri,γ,nの短時間時系列パワー (r̃i , g̃i,ω ) の全組み合わせがと gi,ω,γ,n の周波数の差であり，η は，短時間セグメ
DNN に入力される．入力された短時間時系列パワー
ントの時間軸に沿ったストライド幅，γ = 1, 2, · · · , Γ（i と i + ω ）のパーミュテーションが正しい場合，は，短時間セグメントのインデクスである．なお，Γ
DNN は理想的には「0」を出力する．逆に，パーミュは，短時間のアクティベーションの長さ τ とストラテーションが正しくない場合は，DNN は理想的にはイド幅 η によって決まる．ベクトル ri,γ,n は，参照「1」を出力する．その結果 Fig. 2 の右側に示すよう周波数 i の短時間時系列パワーに対応し，ベクトルに，参照周波数 i に基づくパーミュテーション問題
gi,ω,γ,n は，Fig. 1 に示すように，隣接又は局所周波の発生個所の推定が可能となる．以後このベクトル数 i + ω の短時間時系列パワーに対応する．DNN のは，サブバンドベクトルと呼ぶ．実際には，DNN の入力ベクトルは，(10) を正規化したものである．
出力は [0, 1] の範囲内の値となるので，サブバンドベクトルの生成時は以下の閾値処理を行う．
di,ω,γ×1
d˜i,ω,γ =∈ R4τ(15)≥0∥di,ω,γ ∥2q̃i,ω,γ = round(qi,ω,γ )∈ {0, 1} (17)
ここで，∥ · ∥2 は，L2 ノルムを表す．
q̃i,γ = (q̃i,−Ω,γ , q̃i,−Ω+1,γ , · · · , q̃i,−1,γ , q̃i,0,γ ,提案 DNN モデルは，0 または 1 を出力する 2 値
Tq̃i,1,γ , · · · , q̃i,Ω,γ ) ∈ {0, 1}2Ω+1 (18)
分類器である．推定結果が「0」の場合は，ri,γ,1 とgi,ω,γ,1 が同一音源であることを意味し，同様に ri,γ,2
ここで，round(·) は，丸め演算子である．と gi,ω,γ,2 も同一音源である，一方，推定結果が「1」
音声信号は本来，無音区間が多く存在することかの場合は ri,γ,1 と gi,ω,γ,1（同様に ri,γ,2 と gi,ω,γ,2 ）がら，長さ τ の短時間時系列パワー r̃i,γ や g̃i,ω,γ はほ異なる音源成分であることを意味している．これら
ぼ零ベクトルになる可能性があり，その場合 DNN のの推定処理を Fig. 2 に示す．実際には，DNN の予測予測は不安定になる．この問題に対処するために，提結果は 2 値ではなく次式のような値である．
案手法では，Fig. 3 に示すように，長さ τ の入力ベ()˜
クトルをストライド幅 η でシフトさせて，全時間フqi,ω,γ = DNN di,ω,γ ∈ [0, 1](16)
×1· · · , |yi,(γ−1)η+τ,n |2 )T ∈ Rτ≥0

日本音響学会講演論文集(12)- 266 -

2020年9月

Subbandpermutation
vectorDNN
outputs

11111100111000…

01

Majority
decision

Frequency

Stride1010TimeFig. 3 DNN predictions for all short-time subbands
and their majority decision.

(a)Stride

DNN andMajority decision

Frequency

Subbandpermutation
vectors
11011010011001011

100(b)

TimeFig. 4 Estimation of subband permutation vectors
in all frequency bins.
レームに DNN 予測を走査する．そして，DNN の予測結果を時間軸に関して多数決することで，より信頼性の高いサブバンドベクトル ṽi を得る．この処理は，次のように示される．

1∑q̃i,γ ∈ {0, 1}2Ω+1Γ γ(19)ṽi = round(vi ) ∈ {0, 1}2Ω+1(20)vi =(c)

実際に，音源の順番は時間に依存しないため，多数決によって予測誤差の悪影響を大幅に軽減できる．

3.5 フルバンドベクトルの推定提案する DNN パーミュテーション解決法は（a）サブバンドのストライドによる全周波数のサブバンドベクトルの推定（3.5.1 項及び Fig. 4）及び（b）類似度比較と多数決に基づくフルバンドベクトルの構築（3.5.2 項及び Fig. 5）で構成される．
3.5.1 全周波数におけるサブバンドベクトルの推定サブバンドベクトル ṽi は，基準周波数 i をシフトすることにより全周波数で推定する．この処理をFig. 4 に示す．全部で I 個のサブバンドベクトルṽ1 , ṽ2 , · · · , ṽI を推定している．ここで，各サブバン
ドベクトル ṽ1 , ṽ2 , · · · , ṽI 内の 2 値（「0」及び「1」）は，同じ意味を持たない．これはサブバンド内の周波数成分が，参照周波数 i の成分と同一音源か否かを示しているに過ぎず，参照周波数 i の変化（サブバンドのシフト）とともに対応音源が変化するためである．
例えば，Fig. 4 の ṽf3 の 0 と 1 は，それぞれ赤色と青色の音源成分を示している．一方で，ṽf4 の 0 と 1はそれぞれ青色と赤色の音源成分を示している．これ
らのサブバンドベクトルの整列は次項で処理される．
3.5.2 フルバンドベクトルの構成推定されたサブバンドベクトル ṽ1 , ṽ2 , · · · , ṽI から，次式で定義されたフルバンドベクトル u を構成する．
u = (u1 , u2 , · · · , uI )T ∈ {0, 1}I

日本音響学会講演論文集(21)

Fig. 5 Reconstruction of fullband permutation label: (a) initialization, (b) second, and (c) last steps.
u の構成処理を Fig. 5 に示す．前項で述べた通り，サブバンドベクトル ṽ1 , ṽ2 , · · · , ṽI の 2 値は同じ意味を持たない. そのため，「0」と「1」の値がそれぞれ赤色の音源と青色の音源を示すように，全てのサブバンドベクトルを統一する必要がある．
Fig. 5 (a) は，フルバンドベクトル u の構成における最初のステップを示している．図に示すように，最も低い周波数のサブバンドベクトル ṽis が，フルバンドベクトル u の対応する周波数に挿入される．
ここで，is は，最も低い参照周波数のインデクスを表す．Fig. 5 (a) では，is = f3 及び Ω = 2 であり，u1 , u2 , · · · , u5 は ṽf3 により決定される．
Fig. 5 (b) は Fig. 5 (a) の次のステップを示している．このステップでは，最も低い周波数のサブバンドに隣接するサブバンドベクトルを算出する．推定されたサブバンドベクトル ṽis +1 及びその論理反転ベクトル ṽis +1 （Fig. 5 (b) 中の ṽf4 と ṽf4 ）が用意される．
次に，u の一部を
ǔi = (ui−Ω , ui−Ω+1 , · · · , ui+Ω−1 )T ∈ {0, 1}2Ω (22)の 様 に 定 義 し ，MSE(ṽis +1 , ǔis +1 )及 びMSE(ṽis +1 , ǔis +1 ) を比較する．ここで MSE(·, ·) は，2 つの入力ベクトル間の MSE 値である．その結果より，MSE が小さくなるベクトルを ṽis +1 又は ṽis +1から選択してメモリに格納する．フルバンドベクト
ル u は，Fig. 5 (b) に示すように，メモリに格納されたベクトルを基に周波数毎に多数決処理を行って更新される．

- 267 -

2020年9月

Fullband
permutation
vector

1111111100011000110

Replace0Frequency

001

Majority decision00

SDR improvement

112103.6 フルバンドベクトルに基づく成分の交換
フルバンドベクトル u は，Pi−1 の推定値そのものである．従って，Fig. 6 の様に，u に基づいて周波数毎の分離信号成分を入れ替えることで，パーミュテーション問題が解決できる．

実験

4.1 実験条件
高残響下の音声信号の混合を用いて提案パーミュテーション解決法の性能を評価した．実験では IPSを用いた FDICA [10]，ILRMA [6]，及び提案手法を用いた FDICA の 3 手法を比較した．IPS は，sij を利用しており，FDICA に基づく BSS の上限性能を示している．本実験では JVS コーパス [12] の音声信号（nonpara30）を使用し，これらの音声信号に RWCP
データベース [13] の JR2 インパルス応答（詳細は文献 [10]）を畳み込んで，1 ファイル当たり 20 s，残響
長 470 ms の音声信号を男性 46 名 95 ファイル及び女性 48 名 95 ファイル分作成した．STFT は，窓長512 ms 及びシフト長 128 ms に設定した．
学習データとして，上記の畳み込み音声信号から 2話者選び，観測信号を作成した．これらに FDICA を適用し分離した後，IPS でパーミュテーション問題を解決した Z1 及び Z2 を作成し，もう一度周波数成分をランダムにシャッフルすることで，パーミュテーション問題が残る信号 Y1 及び Y2 と正解ラベルを作成した．
上記の処理で計 40 万ファイルの学習データを作成し，学習用及び検証用に 2 等分した．DNN 学習時の条件は文献 [8] と同様である．短時間区間長は τ = 40 とし，短時間区間のストライド幅は η = 4 とした（Γ = 37となった）．テストデータは，SiSEC2011 [14] の音
声信号 (男性 4 名及び女性 4 名) を使用した．これらの音声信号を JR2 インパルス応答で畳み込み処理を行い，28 種類の観測信号を生成した．評価指標には，SDR [11] の改善量を用いた．
4.2 実験結果
各手法の平均 SDR 改善量を Fig. 7 に示す．文献[10] で報告されているように，ILRMA の分離性能は4 dB 以下であるのに対し，IPS を用いた FDICA では 10 dB 以上の改善を達成している．提案手法を用いた FDICA は，ILRMA を上回る 8 dB 以上の改善を達成しており，上限性能に比較的近いものである．
これは，DNN に基づくパーミュテーション解決法が，正しくを予測を行っている結果である．

日本音響学会講演論文集

642FDICAwith IPS

Solving permutation problem.

以降，前述のステップを繰り返し，完全なフルバンドベクトル u が，Fig. 5 (c) のように構成される．
提案手法では，u の構成過程の反復的多数決により，DNN 推定誤差の悪影響を軽減している．

480TimeFig. 610Fig. 75ILRMA(2 bases)

ILRMA(3 bases)

ILRMA(4 bases)

ProposedmethodAverage SDR improvements.

まとめ

本論文では，DNN に基づく頑健なパーミュテーション解決法を提案し，その性能を評価した．提案手法
は，DNN の予測値の示す音源成分を全周波数で統一するアルゴリズムであり，さらに DNN の予測誤差の影響を軽減するため時間軸と周波数軸に沿った多数決を行う．実験結果より，提案手法の有効性が確認された．今後の課題として，N ≥ 3 の音源数に対応できるアルゴリズム構築が挙げられる．
謝辞 本研究の一部は JSPS 科研費 19K20306 及び19H01116 の助成を受けたものである．

参考文献[1] P. Comon, “Independent component analysis, a new concept?,” Signal Process., 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] H. Sawada, R. Mukai, S. Araki, and S. Makino, “A robustand precise method for solving the permutation problem offrequency-domain blind source separation,” IEEE Trans.
SAP, vol. 12, no. 5, pp. 530–538, Sep. 2004.
[4] 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.
[5] 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. 1626–1641, 2016.
[6] D. Kitamura, N. Ono, H. Sawada, H. Kameoka, andH. Saruwatari, “Determined blind source separation withindependent low-rank matrix analysis,” in Audio Source
Separation, S. Makino, Ed., pp. 125–155. Springer, Cham,2018.
[7] N. Makishima , S. Mogami , N. Takamune, D. Kitamura,
H. Sumino, S. Takamichi, H. Saruwatari, and N. Ono, “Independent deeply learned matrix analysis for determinedaudio source separation,” IEEE/ACM Trans. ASLP, vol.
27, no. 10, pp. 1601–1615, 2019.
[8] 山地修平, 北村大地, “局所時間周波数構造に基づく深層パーミュ
テーション解決法,” 日本音響学会 2020 年春季研究発表会講演
論文集, pp. 317–320, 2020.
[9] K. Matsuoka and S. Nakashima, “Minimal distortion principle for blind source separation,” Proc. ICA, pp. 722–727,
2001.
[10] D. Kitamura, N. Ono, and H. Saruwatari, “Experimentalanalysis of optimal window length for independent low-rankmatrix analysis,” Proc. EUSIPCO, pp. 1210–1214, 2017.
[11] E. Vincent, R. Gribonval, and C. F, “Performance measurement in blind audio source separation,” IEEE Trans.
ASLP, vol. 14, no. 4, pp. 1462–1469, 2006.
[12] S. Takamichi, K. Mitsui, Y. Saito, T. Koriyama, N. Tanji,and H. Saruwatari, “JVS corpus: free Japanese multispeaker voice corpus,” arXiv preprint, 1908.06248, 2019.
[13] 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.
[14] S. Araki, F. Nesta, E. Vincent, Z. Koldovsky, G. Nolte,
A. Ziehe and A. Benichoux, “The 2011 signal separationevaluation campaign (SiSEC2011): -Audio source separation,” Proc. LVA/ICA, pp. 414-422, 2012.

- 268 -2020年9月