卒業研究論文
論文題目
多重解像度時間周波数表現に基づく独立低ランク行列分析

提出年月日

令和 4 年 2 月 25 日

学

科

氏

名

電気情報工学科
細谷 泰稚

印

指導教員（主査）

北村 大地 講師

印副査

重田 和弘 教授印長

辻 正敏 教授

印学科

香川高等専門学校

Independent low-rank matrix analysis based on
multi-resolution time-frequency representations
Taichi Hosotani
Department of Electrical and Computer Engineering
National Institute of Technology, Kagawa College

Abstract
Blind source separation (BSS) is a technique for estimating each original audio sourcefrom an observed signal, and independent low-rank matrix analysis (ILRMA) can achieve
high performance. Optimization in ILRMA consists of estimating spatial and source models (training of spatial demixing ﬁlters and low-rank approximation of a sourcewise timefrequency structure). In conventional ILRMA, the same time-frequency resolutions areused for both models. However, it is reported that the separation performance of ILRMA
strongly depends on the resolution of time-frequency representation. Thus, the performance of ILRMA should be improved by utilizing diﬀerent resolutions of time-frequency
representation for the both models. In this thesis, I propose a new ILRMA algorithm that
introduces multiple time-frequency representations into spatial and source model. Also, I
indicate experimentally that the same resolution of time-frequency representations in theboth models does not always provide the best performance.
Keywords: blind source separation, independent low-rank matrix analysis, spectrogram
resolutions, window function（和訳）
ブラインド音源分離（BSS）は，複数の音源が混合した観測信号から混合前の個々の音源を推定する技術であり，独立低ランク行列分析（ILRMA）が高い分離性能を示している．ILRMA
における最適化は，空間モデル推定（空間分離フィルタの学習）及び音源モデル推定（各音源の時間周波数構造の低ランク近似）からなる．従来の ILRMA では，各モデルの推定において同一の解像度の時間周波数表現を用いている．しかしながら，ILRMA の分離性能は時間周波数領域の解像度に強く依存することが報告されている．このことから，モデル毎に最適な解像度の時間周波数表現を用いた方が良い分離をもたらすと推測される．そこで，本論文では，

ILRMA の空間モデル及び音源モデルの推定に異なる解像度の時間周波数表現を導入した手法を提案する．そして，各モデルの推定に同一の解像度の時間周波数表現を用いることが，必ずしも最大の分離性能を与えるとは限らないことを実験的に示す．

目次
第1章

緒言

1

1.1

本論文の背景 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

本論文の目的 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.3

本論文の構成 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

BSS の基礎技術と従来手法

4

2.1

まえがき . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.2

STFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2.3

時間領域及び周波数領域における BSS の定式化 . . . . . . . . . . . . . . .

6

2.4

ICA 及び FDICA の概要. . . . . . . . . . . . . . . . . . . . . . . . . . .

92.5

NMF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

2.5.1

NMF の概要 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

2.5.2

ISNMF における最適化問題及び反復更新式の定式化 . . . . . . . . .

12

2.6

ILRMA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

2.7

Consistent ILRMA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

2.7.1

スペクトログラム無矛盾性 . . . . . . . . . . . . . . . . . . . . . . .

16

2.7.2

スペクトログラム無矛盾性に基づく ILRMA . . . . . . . . . . . . .

18

本章のまとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

提案手法

21

3.1

まえがき . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

3.2

動機 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

3.3

異なる時間周波数解像度を扱う上での問題点 . . . . . . . . . . . . . . . . .

23

3.4

多重解像度時間周波数表現に基づく ILRMA . . . . . . . . . . . . . . . . .

25

3.5

Chebyshev 窓に基づく 2 種類の窓関数の設計 . . . . . . . . . . . . . . . . .

26

3.6

最適化アルゴリズム . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

3.7

本章のまとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

実験

32

4.1

まえがき . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

4.2

実験条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

第2章

2.8
第3章

第4章

4.3

実験結果 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

4.4

本章のまとめ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

結言

38

第5章
謝辞

39

参考文献

39

付録 A

4 章の実験に対する全結果

45

1

第1章

緒言
1.1 本論文の背景
ブラインド音源分離（blind source separation: BSS）[1] は，混合系や音源情報が未知の条件下で，複数の音源が混合した観測信号から混合前の各音源信号を推定する技術である．

Fig. 1.1 に BSS の概略図を示す．BSS は補聴器や自動採譜，音声認識等の様々な技術に応用されている．
チャネル数（観測に利用したマイクロホン数）が混合している音源数以上となる観測信号を扱う優決定条件 BSS は，音源間の統計的な独立性に基づいた独立成分分析（independent

component analysis: ICA）[2, 3, 4] の登場以降，盛んに研究されている．ICA を周波数毎に適用することで耐残響性を高めた周波数領域独立成分分析（frequency-domain independent

component analysis: FDICA）[5] が提案され，FDICA で推定される周波数毎の分離信号の順番を適切に並び替えるパーミュテーション問題の解決法が検討された [6, 7, 8]．2006 年には，FDICA に対して音源の時間周波数構造仮定を導入することで，パーミュテーション問題を回避しながら分離信号を推定する独立ベクトル分析（independent vector analysis: IVA）

[9, 10, 11] が登場した．その後，補助関数法 [12, 13, 14] 及び反復射影法（iterative projection:IP）[15] に基づく安定かつ高速な IVA（auxiliary-function-based IVA: AuxIVA）[16] が提案されている．

2016 年には，パーミュテーション問題を回避するために何らかの音源モデルを導入するという IVA の画期的なアイデアを拡張して，非負値行列因子分解（nonnegative matrix

factorization: NMF）[12, 17] に基づく低ランク時間周波数構造を ICA の音源モデルに取り入れた独立低ランク行列分析（independent low-rank matrix analysis: ILRMA）[1, 18, 19]が提案された．また，これらの音源モデルを plug-and-play で変更可能な最適化アルゴリズムを採用した BSS [20, 21, 22] も提案されている．さらに，時間周波数領域におけるスペクトロ
グラム無矛盾性 [23, 24] と呼ばれる性質を FDICA 及び IVA に導入した BSS [25] や ILRMAに導入した BSS（consistent ILRMA）[26, 27] も提案されている．

ILRMA の最適化は，時間周波数領域における空間モデル（周波数毎の分離行列）の更新と音源モデル（NMF による低ランク時間周波数構造）の更新からなる．信号の時間周波数表現

2

第 1 章 緒言Fig. 1.1. Overview of BSS.

Fig. 1.2. Comparison between conventional and proposed ILRMA in terms of window
length used in STFT. Conventional ILRMA uses same window length for both
spatial and source models, whereas proposed ILRMA can set window length to
optimal values in each of models.

は短時間 Fourier 変換（short-time Fourier transform: STFT）によって得ることができる．
過去の実験的な調査 [28] により，ILRMA の分離性能は信号の時間周波数表現を求める際のSTFT の窓長に，強く依存していることが分かっている．具体的には，STFT の窓長が極端に短い場合，周波数領域での瞬時混合仮定が成立せず性能が劣化 [29] し，逆に窓長が極端に長い
場合，時間フレーム数が減少することによる統計的推定の不安定さから性能が劣化する．従って，STFT の窓長にはトレードオフが存在する．

1.2 本論文の目的
本論文では，1.1 節で述べた STFT の窓長に対する BSS 性能の傾向に基づき，ILRMA の性能向上を目的として，ILRMA における空間モデル及び音源モデルの最適化に，異なる解像度の時間周波数表現を導入したアルゴリズムを提案する．Fig. 1.2 に提案手法と従来手法の違いを表す図を示す．従来の ILRMA では，空間モデル及び音源モデルの最適化において，同一の窓長で STFT して得られる観測信号（すなわち，1 種類の解像度の時間周波数表現）のみ

1.3 本論文の構成

3を用いているため，Fig. 1.2 に示すように，両モデルの窓長が互いに一致する組しか選べないという制約がある．本論文では，空間モデル及び音源モデルにおける最適な窓長をそれぞれ設
定できるアルゴリズムを提案する．これはすなわち，異なる複数の解像度（多重解像度）の時
間周波数表現に基づく BSS であり，従来の ILRMA を時間周波数解像度という観点から一般化したアルゴリズムに対応している．また提案手法では，STFT による時間周波数解析時にパラメトリックな窓関数である Chebyshev 窓を用いる．これにより，時間周波数表現の（行列としての）サイズを変えることなく見かけ上の窓長を変化させることができ，シンプルな最適化アルゴリズムを得ることができる．実験では，各モデルの最適化における見かけ上の窓長を様々に変化させ，それぞれの窓長に対する音源分離性能の比較を行う．得られた実験結果か
ら，従来手法の consistent ILRMA と提案手法の分離性能を比較し，各モデルの時間周波数表
現の違いが分離性能に与える影響について調査する．

1.3 本論文の構成
2 章では，提案手法を理解する上で重要な，STFT や NMF 等の BSS の基礎技術及びILRMA を含む従来手法について述べる．3 章では，ILRMA の空間モデルと音源モデルの最適化にそれぞれ異なる解像度の時間周波数表現を導入するという提案手法の動機及び詳細について述べる．具体的には，ILRMA の各モデルに異なる解像度の時間周波数表現を導入することがもたらす分離性能への影響，導入の際の問題点，及び問題に対する解決策について説明する．4 章では，提案手法の空間モデル及び音源モデルにおける窓長を変化させることで，従来手法である consistent ILRMA と提案手法の比較実験を行い，得られた結果の傾向について述べる．最後に，5 章で本論文の総括を行い，今後の課題を述べる．

4

第2章

BSS の基礎技術と従来手法
2.1 まえがき
本章では，BSS の基礎技術と従来手法について説明する．時間領域の信号を捉える上で，その信号を時間的に変化するスペクトルとして表現すること，すなわち時間周波数領域で表現することは非常に有効な手段である．特に，音源信号の混合は，時間領域で捉えるよりも時間周波数領域で捉える方が適している．このことから，2.2 節では，時間領域の信号から，その信号の時間周波数領域の表現を得る手法である STFT について取り上げる．また，2.3 節では，時間領域及び周波数領域における BSS の定式化を行う．さらに，2.4 節では，時間領域における音源分離手法である ICA 及び ICA を周波数領域に拡張した FDICA の概要を述べる．
本論文で提案する手法の核となる ILRMA は，FDICA に対して，音源の時間周波数構造における低ランク性を導入した手法である．音源の時間周波数構造を低ランク近似する場合には，NMF という手法が有効である．2.5 節では，この NMF の概要を述べ，ILRMA で用いられる NMF の一種である ISNMF の定式化を行う．そして，2.6 節では，ILRMA の原理について述べる．提案手法との比較を行う consistent ILRMA は，ILRMA に対して，自然なスペク
トログラムの持つ性質であるスペクトログラム無矛盾性を導入した手法である．そこで，2.7
節では，この consistent ILRMA の原理について述べる．最後に，2.8 節で本章の総括を行う．

2.2 STFT
STFT は，Fig. 2.1 に示すように，時間領域の信号から時間的に変化する音色（スペクトル）
としての表現を得る手法である．STFT において，時間領域から周波数領域への変換（解析）
時の窓関数の長さ（すなわち，短時間区間信号の長さ）及びシフト長をそれぞれ Q 及び τ とする．このとき，時間領域の信号 z = [ z[1], z[2], . . . , z[l], . . . , z[L] ]T ∈ RL の j 番目の短時間区間（時間フレーム）信号は次式で表される．

Tz[j] = z[(j − 1)τ + 1], z[(j − 1)τ + 2], . . . , z[(j − 1)τ + Q]T
= z [j] [1], z [j] [2], . . . , z [j] [q], . . . , z [j] [Q] ∈ RQ

(2.1)

2.2 STFT

5

Fig. 2.1. Mechanism of STFT.

ここで，l = 1, 2, . . . , L，j = 1, 2, . . . , J ，及び q = 1, 2, . . . , Q はそれぞれ離散時間のインデクス，時間フレーム，及び時間フレーム内のサンプルのインデクスであり，·T は転置を表す．
なお，本論文では，信号の時間領域及び時間周波数領域における二つの表現の混同を避けるた
め，信号の時間領域のベクトルや行列をローマン体，信号の時間周波数領域におけるベクトル
や行列をイタリック体でそれぞれ示す．また，時間フレーム数 J は次式で与えられる．

J=

L
τ

(2.2)

ただし，J は正の整数となる必要があるため，信号長 L がシフト長 τ で割り切れるように信号をゼロ埋めする．そして，信号 z の STFT によって得られた複素スペクトログラム Z を次式で表記する．

Z = STFTω (z) ∈ CI×J

(2.3)

ここで，ω = [ ω[1], ω[2], . . . , ω[q], . . . , ω[Q] ]T ∈ RQ は STFT で用いる解析時の窓関数を表す．このとき，複素スペクトログラム Z の (i, j) 成分 zij は次式で表される．
QX



2πı(q − 1)(i − 1)zij =z [q]ω[q] exp −
F
q=1
[j]


(2.4)

ここで，F は窓長 Q 以下の正の整数を，i = 1, 2, . . . , I は周波数ビンのインデクスを，ı は虚数単位を示している．また，周波数ビン数 I は次式で与えられる．



F
I=
2


+1

(2.5)

なお，⌊·⌋ は床関数を表す．本稿では，F = Q として扱う．このとき，周波数ビン数 I は窓長

Q に依存する．

6

第 2 章 BSS の基礎技術と従来手法

e とおくとき，逆 STFT を
周波数領域から時間領域への変換（合成）時の窓関数を ωe のペアが次式の完全再構成条件を満たすことを
ISTFTωe (·) と表記する．本論文では，ω と ω仮定する．

z = ISTFTωe (STFTω (z))∀z ∈ RL(2.6)

ここで，ISTFTω
e (Z) は複素スペクトログラム Z を時間信号 z に戻す逆 STFT を表す．なe のペアが式 (2.6) を満たすとき，τ の上限は窓関数 ω の情報（窓関数の各点でのお，ω と ω値と窓長）によって定まる*1 ．このことから，周波数ビン数 J は窓長 Q によって制限される．
式 (2.4) に示すように，STFT は，一定時間毎に信号を切り出し，それぞれの区間信号に解析窓関数 ω を乗じて，離散 Fourier 変換（discrete Fourier transform: DFT）を施すという処理からなる．STFT を適用して得られた複素スペクトログラム Z は，行が周波数，列が時間の複素行列として表すことができる．また，音響信号処理では各時間周波数成分の大きさのみを取り扱うことも多い．その場合は，複素スペクトログラム Z の各要素に関して絶対値を取った振幅スペクトログラム |Z| ∈ RI×J
≥0 や，絶対値の 2 乗をとったパワースペクトログラム|Z|.2 ∈ RI×J
≥0 を処理の対象とする．ここで，行列に対する絶対値記号及びドット付き指数乗はそれぞれ要素毎の絶対値及び要素毎の指数乗を表す．例として，Figs. 2.2(a) 及び (b) にそれぞれ音楽信号及び音声信号のパワースペクトログラムを示す．図中の色の変化は，青色に近づくほどパワーが小さく，黄色に近づくほどパワーが大きいことを表している．

2.3 時間領域及び周波数領域における BSS の定式化
N 個の音源信号が M 個のマイクロホンで観測される状況を考える．この状況での信号の残
響長を L′ とする．n 番目のチャネルの音源信号，m 番目のチャネルの観測信号，及び n 番目のチャネルの分離信号をそれぞれ次式で表す．

T
sn = sn [1], . . . , sn [l], . . . , sn [L] ∈ RLT
xm = xm [1], . . . , xm [l], . . . , xm [L] ∈ RLT
yn = yn [1], . . . , yn [l], . . . , yn [L] ∈ RL

(2.7)
(2.8)
(2.9)

また，時間領域の多チャネルの音源信号，観測信号，及び分離信号をそれぞれ次式で表す．

Ts[l] = s1 [l], . . . , sn [l], . . . , sN [l] ∈ RN
T
x[l] = x1 [l], . . . , xm [l], . . . , xM [l] ∈ RMTy[l] = y1 [l], . . . , yn [l], . . . , yN [l] ∈ RN

(2.10)(2.11)(2.12)*1 STFT を適用する際に情報の欠落が生じていない場合に，時間信号の再構成が可能となる．情報の欠落が生じないようにシフト長 τ を定める必要がある．

2.3 時間領域及び周波数領域における BSS の定式化

Fig. 2.2. Power spectrogram of (a) music and (b) speech signals.

7

8

第 2 章 BSS の基礎技術と従来手法

ここで，n = 1, 2, . . . , N 及び m = 1, 2, . . . , M はそれぞれ音源及びチャネルのインデクスを表す．さらに，各音源信号に対する多チャネルの観測信号及び混合係数をそれぞれ次式で表す．

Tcn [l] = cn1 [l], . . . , cnm [l], . . . , cnM [l] ∈ RMTan [l′ ] = an1 [l′ ], . . . , anm [l′ ], . . . , anM [l′ ] ∈ RM

(2.13)(2.14)

ただし，l′ = 0, 1, . . . , L′ − 1 は残響時間のインデクスである．
各チャネルに STFT を適用して得られる音源信号，観測信号，分離信号，及び各音源信号に対する観測信号のスペクトログラムの (i, j) 番目の要素をそれぞれ次式で表す．

sij = [ sij1 , . . . , sijn , . . . , sijN ]T ∈ CN
xij = [ xij1 , . . . , xijm , . . . , xijM ] ∈ C
T(2.15)
M(2.16)

yij = [ yij1 , . . . , yijn , . . . , yijN ]T ∈ CN
cijn = [ cijn1 , . . . , cijnm , . . . , cijnM ] ∈ C
T(2.17)
M(2.18)

また，音源信号，観測信号及び分離信号に対して，各チャネルの時間周波数行列（スペクトログラム）の表記を Sn ∈ CI×J ，Xm ∈ CI×J 及び Yn ∈ CI×J と定義する．
無響室等の特殊な環境を除き，実際の音響信号の収録環境では音源の混合において残響が発生する．このような音源の混合は，次式で表される畳み込み混合となる．

x[l] =NX′

cn [l] =n=1

N min(LX
X,l)−1an [l′ ]sn [l − l′ ](2.19)

l′ =0

n=1

なお，残響が生じない環境では，式 (2.19) で L′ = 1 とした混合となる．
混合系の残響時間 L′ が STFT における窓長 Q よりも十分短い場合，(2.19) は STFT によって，近似的に周波数領域での瞬時混合へ変換される．ILRMA を含む周波数領域 BSS では，この条件を仮定することで，時間領域での畳み込み混合を次式に示す周波数領域での瞬時混合として扱う．

xij =

N
X

cijn ≈n=1

NX

ain sijn(2.20)

n=1

ここで，ain = [ ain1 , . . . , ainm , . . . , cinM ]T ∈ CM は各音源信号についての混合係数に対して，各チャネルに Q 点 DFT を適用して得られるベクトル*2 であり，次式で与えられる．
QX2πı(q − 1)(i − 1)ainm =
anm [q] exp −F
q=1

(2.21)

なお，上式において，anm [q] が範囲外の参照となる場合，すなわち q ≥ L′ となる場合については，anm [q] = 0 として扱う．いま，周波数毎の混合行列を Ai = [ ai1 , . . . , ain , . . . , aiN ] ∈*2 正確には，式 (2.4) と同様に，Q 点 DFT を適用して得られたベクトルから対称な成分を取り除いたものである．

2.4 ICA 及び FDICA の概要

9

CM ×N とすると，式 (2.20) は次式のように書ける．
xij = Ai sij(2.22)

優決定条件 BSS では M = N を仮定でき，BSS は Ai の逆行列を推定する問題となる．この逆行列を Wi ≈ A−1
i とすると，分離信号は次式となる．

yij = Wi xij(2.23)

ここで，Wi = [ wi1 , . . . , win , . . . , wiN ]H ∈ CN ×M は分離行列と呼ばれ，·H はエルミート転置を表す．

2.4 ICA 及び FDICA の概要
統計数理アルゴリズムである ICA は，時間領域での瞬時混合及び信号の統計的性質を仮定した BSS である．ICA では，前述の仮定を用いて，瞬時混合系の逆系を推定することで，音源分離を行っている．

ICA が仮定する統計的性質とは，音源信号間の統計的独立性及び音源信号が従う生成モデルの非ガウス性である．信号のチャネルの順序及びスケール（大きさ）の違いはこれらの性質に影響を与えない．従って，ICA によって推定される分離信号には，以下の任意性が存在する．

1. 分離信号のチャネルの順序には任意性がある
2. 分離信号のスケールには任意性がある
これらの任意性は分離信号に対して Fig. 2.3 のように現れる．上記の任意性 1 より，元々の信号源の順序が入れ替わる可能性がある．また，任意性 2 より，分離信号のスケールが混合前の音源信号のスケールから変化してしまう可能性がある．なお，信号のスケールの任意性に関しては，プロジェクションバック（projection back: PB）法 [30] と呼ばれる補正方法が提案されている．さらに，ICA には上記の問題に加えて，残響の生じた混合信号に対する音源分離性能が著しく悪化するという問題がある．これは，2.3 節で述べたように，残響が生じた音源の混合は瞬時混合ではなく畳み込み混合となり，ICA における時間領域での瞬時混合仮定が成り立たないことに起因する．
残響を含む信号に対する音源分離は，畳み込み混合系の逆系を推定することで達成できる．
しかし一般に，時間領域における畳み込み混合系の逆系の推定は困難である．この問題は，2.2節で述べた STFT を用いることで解決できる．畳み込みは STFT によって積和に変換される
ため，時間領域での畳み込み混合を式 (2.22) で表される時間周波数領域での瞬時混合として扱うことが可能となる．このことを利用して，観測信号を STFT して得られたスペクトログ
ラムの各周波数ビンの複素時系列信号 xi1 , xi2 , . . . , xij , . . . , xiJ に対して，独立な ICA を適用し，分離信号のスペクトログラムを推定する手法が提案されている．この手法は FDICA と呼ばれる．前述した通り，ICA の分離信号にはスケールや順序の任意性がある．FDICA では
周波数毎に独立な ICA が適用されるため，各周波数ビンで推定した分離信号のスペクトログ

10

第 2 章 BSS の基礎技術と従来手法

Fig. 2.3. Uncertainty in ICA. ICA cannot determine order and scales of estimated signals.

ラムの並びとスケールがバラバラになるという問題が生じる．周波数毎のスケールの任意性については，前述の PB 法により解析的に復元可能である．それに対し，Fig. 2.4 に示すような周波数毎の順序の任意性を解決することは困難である．このような問題をパーミュテーション問題と呼び，これを解決することが FDICA における大きな課題である．

2.5 NMF
2.5.1 NMF の概要
NMF とは，1999 年に D. D. Lee と H. S. Seung によって提案された非負行列に対する分
解アルゴリズムである．行列の分解は，線形方程式を解くための LU 分解や QR 分解，ベクトル空間の概念に基づく固有値分解や特異値分解等が代表的である．それに対して，NMF は非負行列を分解対象としている点において，これらの行列分解手法と大きく異なっている．2.2
節で説明した通り，時間領域の信号に対して，STFT を経て得られる振幅（あるいはパワー）スペクトログラムは非負行列である．従って，これらの行列は NMF の分解対象となる．

NMF は，次式に示すように，非負行列を別の 2 つの非負行列の行列積に分解する数理アル2.5 NMF

11

Fig. 2.4. Permutation problem in FDICA.

Fig. 2.5. NMF for audio signals, where K = 2.

ゴリズムである．

Z ≈ TV(2.24)

I×K
ここで，Z ∈ RI×J
≥0 は分解の対象となる非負行列であり，T = [ t1 t2 · · · tK ] ∈ R≥0 及び

V = [ v1 v2 · · · vK ]T ∈ RK×J
はそれぞれ基底行列及びアクティベーション行列と呼ばれ≥0る非負行列である．T の列ベクトルである tk は基底ベクトルと呼ばれ，その本数 K は通常K ≪ min(I, J) となるように設定される．ただし，k = 1, 2, . . . , K は基底ベクトルのインデクスを示す．行列 T V のランクは基底ベクトルの本数 K が最大であるため，行列 Z のラン
クが K を上回る場合，Z を T V によって完全に再構成することができない．従って，この場
合において，NMF は低ランク近似分解となる．これは，Z 中に頻出する少数（K 個）の潜在的なパターンを基底ベクトルとして抽出できる教師無し学習であると解釈できる．

NMF を音響信号に適用する場合，振幅（あるいはパワー）スペクトログラムを非負観測行
列 Z とすることが一般的である．この場合，Fig. 2.5 に示すように，音響信号中の頻出スペク

12

第 2 章 BSS の基礎技術と従来手法

トルが tk として得られ，さらに各スペクトルの時間的強度変化が vk として現れる．このように，NMF は音響信号中のスペクトルパターンを学習できるため，音楽信号解析 [31] や音源分
離 [32, 33] 等に頻繁に適用される．なお，NMF では，T 及び V の推定問題を解析的に解くことができないため，次節に示す最適化問題を解く必要がある．

2.5.2 ISNMF における最適化問題及び反復更新式の定式化
NMF では，次式の最適化問題を解くことで非負の変数行列 T 及び V を推定する．
Minimize D(Z|T V ) s.t. tik , vkj ≥ 0 ∀i, j, k
T ,V

(2.25)

ここで，tik 及び vkj はそれぞれ T 及び V の要素である．また，D(·|·) は 2 つの行列間の乖離度を測る関数である．そのような関数としては，二乗 Euclid 距離，一般化 Kullback–Leibler ダイバージェンス，Itakura–Saito ダイバージェンスが用いられる．本論文では，Itakura–Saito
ダイバージェンスに基づく NMF（ISNMF）を取り扱う．

A 及び B を I × J 型の非負行列とするとき，これらの非負行列に対する Itakura–Saito ダイバージェンスは次式で定義される．

DIS (A|B) =

I X
J 
X
aij
i=1 j=1


aij
− log
−1
bij
bij

(2.26)

I×J
ただし，aij 及び bij はそれぞれ A ∈ RI×J
≥0 及び B ∈ R≥0 の要素を表す．式 (2.25) にItakura–Saito ダイバージェンスを導入することで，ISNMF の最適化問題は次式となる．
!
J
I X
X
zij
zij
Minimize DIS (Z|T V ) =− log PK−1PK
T ,V
k=1 tik vkj
k=1 tik vkj
i=1 j=1
s.t. tik , vkj ≥ 0 ∀i, j, k

(2.27)式 (2.27) で表される最適化問題に補助関数法を適用して得られる反復更新式は次式のよう
になる．

v
P
−2
uP
[t] [t]
[t]
u J
K
z
vkj
t
v
u
′
′
′
j=1 ijk =1 ik k j
[t+1]
[t] u
tik = tik t P−1[t] PK[t] [t]Jv
tv
j=1 kjk′ =1 ik′ k′ j
v
P
−2
uP
[t] [t]
[t]
u I
K
z
t
v
tik
u i=1 ij
k′ =1 ik′ k′ j
[t+1]
[t] u
vkj = vkj t P−1[t] PK[t] [t]Ii=1 tikk′ =1 tik′ vk′ j(2.28)(2.29)

ここで，各変数の上付き文字 ·[t] は変数更新の反復回数である．式 (2.28) 及び (2.29) は，次の

2.6 ILRMA

13

ように行列形式で表現することもできる．

(T ←T ⊙Z ⊙ (T V )

.−2 

VT

). 12
(2.30).−1

V ←V ⊙(T V ) V T()1.−2  . 2
T T Z ⊙ (T V )T T (T V )(2.31).−1

ここで，行列間の演算 ⊙ 及び行列間の分数はそれぞれ要素毎の積及び商を表す．なお，式(2.30) 及び (2.31) では反復回数 t の表記は省略し，変数更新を表す演算子 ← を用いている．

2.6 ILRMA
2.4 節では，FDICA にはパーミュテーション問題という大きな問題が存在することを述べた．それに対して，2006 年に提案された IVA では，FDICA の音源モデルに「同一音源の全周波数成分は連動して生起する傾向にある」という仮定を導入することで，パーミュテーショ
ン問題を回避している．このように，FDICA の音源モデルに他の仮定を導入する方法はパーミュテーション問題の回避策として有効である．その一方で，IVA が仮定する全周波数成分の共変性は，周波数成分の部分的な強度変化を表現することができない．つまり，IVA は楽器音信号のような，基本周波数とその倍音のみが強い共変性を持つ音響信号，すなわち，共変性の程度が各時間周波数に依存して変化する音響信号の分離には適していない．この問題を克服した BSS として，ILRMA が提案されている．ILRMA は，FDICA の音源モデルに「音源の時間周波数構造は低ランクで近似できる」という仮定を導入した BSS である．この仮定によって，ILRMA では，各音源の時間周波数領域における共変性をより詳細に取り扱うことができ，明確な倍音構造を有する音響信号に良く適合した分離が可能となる．

ILRMA では，分離信号の各チャネルにおける複素スペクトログラム Yn ∈ CI×J に対して，
次式の複素ガウス分布を生成モデルとして仮定する．

p(Yn ) =Yp(yijn )

i,j



1
|yijn |2
=
exp −
πrijn
rijn
i,j
Y

(2.32)

ただし，rijn は各音源の時間周波数毎の分散であり，パワーの期待値として rijn = E[|yijn |2 ]と書ける．なお，E[·] は期待値を表す．式 (2.32) の仮定の下での，分離行列 Wi に関する観測
信号の負対数尤度関数は次式で与えられる．

L = −2J

X
ilog | det Wi | +X  |yijn |2
i,j,n

rijn


+ log rijn

(2.33)

ここで，本節の冒頭で述べた ILRMA の音源モデルに対する仮定を導入すると，分散 rijn は

14

第 2 章 BSS の基礎技術と従来手法

NMF によって次式のように低ランク近似される．
X
rijn =
tikn vkjn(2.34)k
I×J
また，rijn ，tikn ，及び vkjn は非負の値であり，それぞれを要素として持つ分散行列 Rn ∈ R>0，
K×J
基底行列 Tn ∈ RI×K
によって，次式のように行≥0 ，及びアクティベーション行列 Vn ∈ R≥0

列形式で表現できる．

Rn = T n V n(2.35)式 (2.34) 及び (2.35) は，各音源の時間周波数構造 Rn がランク K の非負行列で近似されることを示している．ここで，式 (2.33) に式 (2.34) を代入することで次式が得られる．

L = −2J

X
ilog | det Wi | +

X
i,j,n

X
|y |2
P ijn
+ log
tikn vkjn
k tikn vkjn!
(2.36)kN
ILRMA では，分離行列 {W }li=1 ，基底行列 {T }N
n=1 ，及びアクティベーション行列 {V }n=1に関して，式 (2.36) の最小化を行っている．式 (2.36) の第 2 項（変数 i，j ，及び k に関する総和）の最小化は，ISNMF の最適化問題 (2.27) において，Z ，T ，及び V をそれぞれ |Yn |.2 ，Tn ，及び Vn とした式に対応する*3 ．これは，ILRMA における分散の最尤推定が，分離信号の各チャネルにおけるパワースペクトログラム |Yn |.2 を低ランク近似しながら分離行列 Wiの推定を行うこと，すなわち分離信号のスペクトログラムを低ランクに誘導することに等価であることを示している．また，式 (2.36) の第 1 項（分離行列 Wi の行列式の負対数を取った項）及び第 2 項内部の yijn に関する部分は IVA における負対数尤度関数に対応する．このことから，ILRMA は IVA と同様に，周波数成分の共変性，及び，FDICA の性質である分離信号間の独立性を考慮しているともいえる*4 ．つまり，ILRMA は，音源の時間周波数構造 Rn
を低ランク行列として近似しながら，時間周波数構造の共変性と分離信号間の独立性を加味した分離行列 Wi の推定を行っていると解釈できる．本論文では，分離行列 Wi を空間モデル（あるいは空間分離フィルタ），音源の時間周波数構造 Rn を音源モデルと呼ぶ．ここまでで説明した，ILRMA における音源分離の原理を Fig. 2.6 に示す．
先に述べたことから，ILRMA における音源モデルの推定には式 (2.30) 及び (2.31) で表される ISNMF の更新式を，空間モデルの推定には AuxIVA で提案された反復最適化手法である IP を適用することができる．ここで用いる IP は，初期の IVA での反復最適化手法である
自然勾配法 [34, 35] よりも，高速かつ安定な分離行列の最適化手法であることが実験的に示されている．空間モデルの推定では，分離行列 Wi は行毎に更新され，その更新式は IP を用い*3 最小化を考える上で，定数項の差異は問題とならない．
*4 既に述べたように，IVA と ILRMA の間で共変性の構造は異なる．IVA では全周波数で一様な強度変化を持

つ音源構造を仮定しているが，ILRMA ではより詳細な音源構造を仮定している．

2.6 ILRMA

15

Fig. 2.6. Principle of BSS based on ILRMA.

て次式で与えられる．

Uin ←

1X
1
xij xHij
J j [Tn Vn ]i,jwin ← (Wi Uin )−1en(2.37)(2.38)

Hwin ← win win
Uin win− 12(2.39)
N

ここで，[·]r,c は行列の (r, c) 番目の要素を表す．なお，en ∈ {0, 1}

は n 番目の要素が 1，他

要素が 0 の単位ベクトルである．また，win は分離行列 Wi の行ベクトルであり，分離ベクトルと呼ばれる．各分離ベクトルを更新した後は，分離信号 Yn を次式で更新する．
H
yijn ← win
xij(2.40)

音源モデルの推定では，基底行列 Tn 及びアクティベーション行列 Vn が式 (2.30) 及び (2.31)
によって更新され，式中の変数を ILRMA の変数で書き換えた更新式は次式で与えられる．

(Tn ← Tn ⊙(Vn ← Vn ⊙|Yn |.2 ⊙ (Tn Vn )(Tn Vn )

.−1

.−2 

VnT

VnT

). 12

)1

.−2  . 2
TnT |Yn |.2 ⊙ (Tn Vn )TnT (Tn Vn ).−1

(2.41)(2.42)

Tn 及び Vn の更新後に，分散行列 Rn の更新を明示的に行うのであれば，それは次式で表される．

Rn ← T n V n(2.43)

以上の空間モデルと音源モデルの更新式を交互に反復することで，式 (2.33) を最小化できる．
また，ILRMA によって推定された分離信号には，FDICA や IVA と同様に，周波数毎の順序とスケールの任意性がある．従って，反復最適化後は，次式の PB 法を適用することで，分離信号のスケールを補正する．

ỹijn = Wi−1 (en ⊙ yij ) = yijn λin ,(2.44)

16

第 2 章 BSS の基礎技術と従来手法

ここで，ỹijn = [ ỹijn1 , ỹijn2 , . . . , ỹijnm , . . . , ỹijnM ]T ∈ CM はスケール補正後の分離信号の(i, j) 番目の成分，λin = [ λin1 , λin2 , . . . , λinm , . . . , λinM ]T ∈ CM はスケール補正係数を表す．また，ベクトル間の演算 ⊙ は，行列に対する演算と同様に，要素毎の積を表す．

2.7 Consistent ILRMA
2.7.1 スペクトログラム無矛盾性
STFT を適用する際の窓関数を乗じるという操作は，その窓関数のスペクトルが時間周波数領域において，（時間領域での各短時間区間信号のスペクトルを重み係数として）周波数方向に畳み込まれることに相当する．従って，時間周波数領域での大きなパワー値は周波数方向に
滲み，共起性が生まれる．また，STFT を適用した際の時間フレーム間のオーバーラップにより，時間周波数領域において，時間方向の冗長性が生じる．ここで述べた冗長性とは，互いに交わる 2 つの時間フレームが共有している区間の情報は，両方のフレームに含まれるという性質のことである．これによって，時間周波数領域での大きなパワー値は時間方向にも滲み，共起性が生まれる．結果的に，時間周波数領域では，各時間周波数成分の近傍で一貫した共起性が生じていることが自然である．この共起性はスペクトログラム無矛盾性と呼ばれる．また，
時間周波数領域における何らかの信号処理によって，共起性が崩された状態はスペクトログラ
ム矛盾と呼ばれる．Figs. 2.7(a) 及び (b) に，共起性が一貫していない矛盾したスペクトログ
ラム，及び，このスペクトログラムに対応する無矛盾なスペクトログラムをそれぞれ示す．ここで，Sart は，人工的に作られたスペクトログラム，すなわち一貫した共起性が考慮されていない矛盾したスペクトログラムを表す．Fig. 2.9(a) は，中心部分の時間周波数グリッドのみが大きなパワー値を持つ矛盾したスペクトログラムである．また，Fig. 2.9(b) は，Fig. 2.9(a)のスペクトログラムに対して，後述する方法により無矛盾性を担保したスペクトログラムである．
任意のスペクトログラムに対するスペクトログラム無矛盾性の担保は，逆 STFT 及び STFT
を続けて適用することで実現できる．この操作を視覚的に表現した図を，Fig. 2.8 に示す．

Fig. 2.8 内の時間領域の信号 s は，STFT によって時間周波数領域の信号 S に写される．このスペクトログラム S に対し，何らかの信号処理を加えることで，新たに得られたスペクトログ
ラムを S ′ とする．S ′ が共起性の崩れた矛盾したスペクトログラムとなった場合，S ′ に直接対応する時間波形は時間領域に存在しない．この矛盾したスペクトログラム S ′ に対して，逆STFT を適用すると，時間周波数領域における無矛盾なスペクトログラム S ′′ に射影されたう
えで時間領域に変換される．その結果，新たな時間信号 s′′ が得られる．従って，時間周波数
領域のいかなるスペクトログラムも，一度逆 STFT を適用して時間領域に戻し，再び STFTを適用して時間周波数領域の信号に変換することで，無矛盾なスペクトログラムに変換することができる．

Figs. 2.9(a) 及び (b) に，共起性に一貫性のない矛盾した音楽信号のスペクトログラム，及び，このスペクトログラムに対応する無矛盾なスペクトログラムをそれぞれ示す．Fig. 2.9(a)2.7 Consistent ILRMA

Fig. 2.7. (a) inconsistent power spectrograms |Sart |2 and (b) their consistent version obtained by applying inverse STFT and STFT. Spectrogram of (a) is artiﬁcially
produced with random phase.

Fig. 2.8. Spectrogram (in)consistency with STFT and inverse STFT.

17

18

第 2 章 BSS の基礎技術と従来手法

Fig. 2.9. (a) inconsistent power spectrograms |Sart |2 and (b) their consistent version obtained by applying inverse STFT and STFT. Spectrogram of (a) is music andspeech signals with random dropout.

は，音楽信号の無矛盾なスペクトログラムに対して，時間周波数グリッドをランダムに選び，それらのパワーを 0 にしたスペクトログラムである．Figs. 2.7 及び 2.9 から，無矛盾なスペク
トログラムにおいては，時間と周波数の両方向に信号のパワーの滲みが生じていることが確認できる．これにより，スペクトログラム無矛盾性の担保は，スペクトログラムの時間方向及び周波数方向に対して，スムージングを施す処理であると解釈できる．このことから，スペクト
ログラム無矛盾性は，スペクトログラムにおける各時間周波数近傍の連動性とも解釈できる．

Fig. 2.8 の通り，逆 STFT は，矛盾したスペクトログラムから無矛盾なスペクトログラムを復元する操作である．すなわち，スペクトログラム Z の無矛盾性は次式によって特徴付けることができる．

E(Z) = Z − STFTω (ISTFTωe (Z))(2.45)

このとき，(2.45) のノルム ∥E(Z)∥ が 0 となるスペクトログラム Z を無矛盾と呼ぶ．

2.7.2 スペクトログラム無矛盾性に基づく ILRMA
前節で述べたスペクトログラム無矛盾性は，BSS の性能向上に寄与することが明らかにされている [25]．これは，スペクトログラム無矛盾性の担保によって与えられる時間周波数領域の共起性が，パーミュテーション問題を緩和する働きを持っていることに起因する．
パーミュテーション問題が生じたスペクトログラムに対して，スペクトログラム無矛盾性の担保が与える影響を Fig. 2.10 に示す．Fig. 2.10(a) は，音楽信号のパワースペクトログ
ラム |S|.2 である．Fig. 2.10(b) は，スペクトログラム S に対して，人工的にパーミュテー
ション問題を引き起こしたもののパワースペクトログラム |S (perm) |.2 である．Fig. 2.10(c)は，S (perm) に対してスペクトログラム無矛盾性を担保したもののパワースペクトログラム

2.7 Consistent ILRMA

19

Fig. 2.10. Smoothing eﬀect of spectrogram consistency applied to permutation misaligned
music signal. (a) is power spectrogram of original source signal |S|.2 . (b) is
randomly permuted version of (a), which simulates the permutation problem(perm)and is denoted as Sn. (c) is the consistent version of (b).

|STFTω (ISTFTωe (S (perm) ))|.2 である．これらの図から分かるように，スペクトログラムにおける時間周波数領域の共起性は，周波数方向の連続性を強調する．その結果，周波数毎の音源成分が正しく整列された状態に誘導され，パーミュテーション問題を緩和することができる．
スペクトログラム無矛盾性を ILRMA の最適アルゴリズムの毎反復に導入した BSS として，consistent ILRMA が提案されている [26, 27]．ILRMA の空間モデルの反復最適化更新
式 (2.37)–(2.40) 及び NMF 音源モデルの反復最適化更新式 (2.41)，(2.42) において，以下の演算を挿入することで，毎回の反復においてスペクトログラム無矛盾性を担保する．

Yn ← STFTω (ISTFTωe (Yn ))(2.46)式 (2.46) は，Fig. 2.11 に示すように，分離信号のスペクトログラム Yn を無矛盾なスペクト
ログラムの集合へと射影していることに対応する．ここで，赤色，青色，橙色，及び紫色の矢
印はそれぞれ STFT の適用，ILRMA における分離行列 Wi の反復更新，Consistent ILRMA
におけるスペクトログラム無矛盾性の担保，及び ISTFT の適用（すなわち，無矛盾なスペク
トログラム集合への射影及び時間領域への変換）を表す．また，x，X ，Y ，及び S はそれぞれ時間領域での観測信号，観測信号 x のスペクトログラム，分離信号のスペクトログラム，及び音源信号のスペクトログラムを表す．仮に Yn が無矛盾であれば，式 (2.46) は Yn に影響を与えない．また，Yn に矛盾があれば，Figs. 2.7，2.9，及び 2.10 に示すように，式 (2.46) は

Yn の時間方向と周波数方向の両方にスムージングをかける形で作用する．
上記の新しい処理の導入により，最適化の毎反復において，矛盾したスペクトログラムが無矛盾なスペクトログラムに射影される．つまり，Fig. 2.11 に示すように，Consistent ILRMA
では，ILRMA よりも真の音源信号に近づきながら音源分離を進めることができる．

2.6 節で述べた通り，ILRMA で推定した分離信号には周波数毎のスケールの任意性がある．
この任意性は，スペクトログラム無矛盾性を崩す原因となる．従って，Consistent ILRMA では，周波数毎のスケールの任意性に起因する矛盾の影響を最小限に抑える必要がある．これを目的として，最適化の毎反復で式 (2.46) を行う直前に，式 (2.44) の PB 法を適用する．PB 法

20

第 2 章 BSS の基礎技術と従来手法

Fig. 2.11. Comparison between ILRMA and consistent ILRMA, where arrows in red show
STFT, arrows in dark blue show iterative update of parameters in ILRMA,
arrows in orange show ensuring process of spectrogram consistency, and arrows
in purple show projection onto set of consistent spectrograms applied in inverse
STFT. Separated signal estimated by consistent ILRMA tends to approach tooracle source signal S.

の適用によって，目的関数 (2.33) の値が変動する．この変化分を補正するために，他の変数を次のように更新する．

win ← λinmref win(2.47)

H
yijn ← win
xij(2.48)

tikn ← |λinmref |2 tikn

(2.49)

ここで，mref は PB 法で用いるリファレンスチャネルのインデクスである．

2.8 本章のまとめ
本章では，BSS の基礎技術と従来手法について説明した．次章では，2.6 節で説明した

ILRMA の空間モデル及び音源モデルの更新に，それぞれ異なる解像度の時間周波数表現を導入した方法を提案する．また，提案手法が持つ利点や，提案手法のアイデアを実現する上で生じる問題点とそれに対する解決策について詳しく述べる．

21

第3章

提案手法
3.1 まえがき
本章では，提案手法である多重解像度時間周波数表現に基づく ILRMA について説明する．

3.2 節では，ILRMA の空間モデルと音源モデルの最適化にそれぞれ適した窓長が存在するという予想に基づき，提案手法の動機について詳しく述べる．3.3 節では，提案手法のアイデアを実現する上で生じる問題について述べる．3.4 節では，3.3 節で述べた問題点を回避するために，窓関数の有効幅を変化させるという方法を導入し，その具体的な内容について説明する．3.5 節では，窓関数の有効幅を設定できる Chebyshev 窓の説明を行う．3.6 節では，提案手法における最適化アルゴリズムを示す．最後に，3.7 節で本章の総括を行う．

3.2 動機
2.3 節で述べた通り，STFT における窓長が観測信号の残響時間より短い場合，時間領域の残響による畳み込み混合を時間周波数領域の瞬時混合に変換できず，式 (2.22) の仮定が成り立たない．このことから，式 (2.23) での高精度な BSS は原理的に困難となる．そのため，STFTの窓長は残響時間を十分超える長さ（例えば 256 ms 以上）に設定されることが一般的である[28]．よって，ILRMA では，空間モデルの最適化を行う際に適切な窓長が存在すると予想される．一方で，音声や楽器音等の音響信号を NMF で時間周波数解析する場合は，16〜128 ms
程度の窓長で STFT を適用する場合が多い（例えば，文献 [31]，[33]，[36] 等）．ILRMA においても，各音源のスペクトログラムを NMF でモデル化する関係上，低ランク行列として良く近似できる窓長が存在するはずである．例として，Figs. 3.1(a) 及び (b) に，ボーカルの音響信号を異なる窓長（32 ms 及び 256 ms）で STFT して得られるパワースペクトログラムを示す．この図から，NMF によるスペクトログラムの近似精度は STFT の窓長の長さに大きく依存することが予想される．従って，ILRMA においては音源モデルの最適化の観点からも適切な窓長が存在し，それは前述の空間モデルの窓長の最適値とは異なる可能性がある．

22

第3章

提案手法

Fig. 3.1. Vocal spectrograms using (a) narrow (32 ms) and (b) wide (256 ms) windows.

3.3 異なる時間周波数解像度を扱う上での問題点

23

3.3 異なる時間周波数解像度を扱う上での問題点
前節で述べたことを基に，ILRMA の空間モデルと音源モデルの最適化に異なる窓長の窓関数を導入することを考える．いま，空間モデル及び音源モデルの最適化における STFT の窓関
数を，それぞれ ω (Ex) 及び ω (In) とする．なお，本節に限り，二つの窓関数 ω (Ex) と ω (In) の
長さは互いに異なっているものとして扱う．観測信号 xm を窓関数 ω (Ex) 及び ω (In) で STFT
(Ex)

して得られるスペクトログラムを，それぞれ Xm
を窓関数 ω

(Ex)

及び ω

(In)

(In)

及び Xm

と表す．また，分離信号 yn(Ex)

で STFT して得られるスペクトログラムを，それぞれ Yn

及び

(In)
(Ex)
(Ex)
Yn と表す．分離信号の更新式 (2.40) にスペクトログラム Xm 及び Ynを導入した式を以下に示す．
(Ex)

(Ex)

H
yijn ← win
xij

(3.1)

(Ex)

ただし，yijn は，分離信号の n 番目のチャネルを窓関数 ω (Ex) で STFT して得られる観測信(Ex)号のスペクトログラムの (i, j) 番目の要素である．また，xij

は，観測信号の各チャネルを

窓関数 ω (Ex) で STFT して得られる観測信号のスペクトログラムの (i, j) 番目の要素であり，次式で表される．
(Ex)

xij

(Ex)

(Ex)

(Ex)

= [ xij1 , . . . , xijm , . . . , xijM ]T ∈ CM

(3.2)

(In)

を導入した式を以下に

音源モデルの反復更新式 (2.41) 及び (2.42) にスペクトログラム Yn示す．

(Tn ← Tn ⊙(Vn ← Vn ⊙(In) .2|Yn| ⊙ (Tn Vn )(Tn Vn )

.−1

.−2 

VnT

). 21
(3.3)

VnT

)1
 (In).−2  . 2
TnT |Yn |.2 ⊙ (Tn Vn )TnT (Tn Vn )

(3.4)

.−1

(Ex)

また，空間モデルの反復更新式の 1 つである (2.37) にスペクトログラム Xmを導入した式を以下に示す．

Uin ←

1X
1
(Ex) (Ex)Hxxij
J j [Tn Vn ]i,j ij(3.5)

2.2 節で述べた通り，スペクトログラムの周波数ビン数は STFT の窓長に依存する*1 ．よっ(Ex)て，ω (Ex) と ω (In) の長さが異なる状況において，二つのスペクトログラム Xn

(In)

及び Yn

のサイズ（周波数ビン数）は互いに異なる．つまり，式 (3.1) で得られたスペクトログラム*1 スペクトログラムの時間フレーム数は，STFT の窓関数の情報（窓関数の各点での値と窓長）によって制限される．STFT の窓長が時間フレーム数に直接影響することはないため，ここでは周波数ビン数に対する影響のみを考える．

第3章

24

提案手法

Fig. 3.2. Process ﬂow for converting time-frequency resolution via time-domain waveform.

(Ex)

Yn

(In)

を式 (3.2) の Ynとして用いることはできない．これは，ILRMA の空間モデル及び(Ex)

音源モデルの最適化更新に，異なる窓長で得られたスペクトログラム Yn

(In)

及び Yn

をそ

のまま用いることはできないことを意味する．従って，ILRMA における各モデルの最適化に(Ex)

異なる窓長の窓関数を導入するためには，二つのスペクトログラム Yn

(In)

と Yn

とでサイ

ズが異なっているという問題を何らかの方法で解決する必要がある．
スペクトログラムのサイズを一致させるには，様々な方法が考えられる．本論文では，スペ(Ex)

クトログラム Yn

(In)

及び Yn

のサイズを一致させる方法の一つとして，Fig. 3.2 に示す方法を提案する．これは，スペクトログラムに逆 STFT 及び STFT を施すことで，時間領域を介してサイズを一致させる方法である．

Fig. 3.2 の処理の詳しい説明は次節に回し，ここからは残った問題について説明する．前(Ex)

述の通り，二つのスペクトログラム Xm

(In)

及び Yn

の周波数ビン数は互いに異なる．ここ

(In)
(Ex)
の周波数ビン数をそれぞれ I (Ex) 及び I (In) とする．
で，スペクトログラム Xm 及び Yn

式 (3.3) 及び (3.4) から，音源モデルの最適化によって得られる基底行列 Tn は I (In) × K 型の非負行列，アクティベーション行列 Vn は K × J 型の非負行列となる．これらの行列を式(3.5) に用いると，式中の Tn Vn のサイズは I (In) × J となる．しかし，式 (3.5) で使っている(Ex)

スペクトログラム Xmのサイズは I (Ex) × J である．I (Ex) ̸= I (In) であるから，式 (3.5)(Ex)

において，[Tn Vn ]i,j と xijの間で周波数ビンのインデクス i のとる範囲が異なる．つまり，

ILRMA の空間モデル及び音源モデルの最適化に用いる STFT の窓長が互いに異なっている
場合，アルゴリズム内部の計算に不都合が生じ，各モデルの反復更新を進めることができない．

3.4 多重解像度時間周波数表現に基づく ILRMA

25

前述の問題を解決するために，Fig. 3.2 と同様の方法を適用し，音源モデルの更新式により(Ex)得られた分散行列 Rn = Tn Vn とスペクトログラム Xmのサイズを一致させることを考える．分散行列 Rn は音源モデルの窓関数 ω (In) から得られた表現であるため，ω (In) による時
間周波数表現から ω (Ex) による時間周波数表現への変換，すなわち Fig. 3.2 の逆変換を施す形となる．なお，Fig. 3.2 において，空間モデルと音源モデルの立ち位置を入れ替えることで，逆変換が得られる．ここで，分散行列 Rn はパワースペクトログラムであるため，これに対して逆 STFT を適用する際には，位相を付与しなければならない*2 ．このとき，分散行列 Rnに付与する位相としては，様々なものが考えられる．つまり，Fig. 3.2 の方法をそのまま適用
すると，分散行列 Rn に対する位相付与の任意性という問題に直面する．スペクトログラムの周波数ビン数は STFT の窓長に依存するため，ω (Ex) と ω (In) の長さを揃えることで，Tn Vn
(Ex)

と Xm

のサイズが一致し，位相付与の任意性を回避できる．しかし，この解決法は提案手法における「ILRMA の空間モデル及び音源モデルの最適化に異なる窓長の時間周波数表現を用いる」という戦略と矛盾する．そこで，提案手法では，空間モデルと音源モデルの最適化更新に異なる時間周波数解像度のスペクトログラムを用いるとともに，Fig. 3.2 内の二つの窓関数の実際の長さを一致させるという方法によって，この二律背反の状態を解決する．その方法の具体的な内容については，次節で説明する．なお，位相付与の任意性に関しては，本論文の趣旨から外れるため，これ以上は踏み込まないこととする．

3.4 多重解像度時間周波数表現に基づく ILRMA
提案手法では，ILRMA の空間モデルと音源モデルにそれぞれ異なる窓長の窓関数を導入
し，各変数の最適化を行う．その際，前節で述べた位相付与の任意性を回避するため，各モデルの変数更新において，Fig. 3.3 に示すように，STFT における窓関数の有効幅（細さ）のみを変化させることにより，見かけ上の窓長を変化させるという方法を提案する．この方法で(Ex)は，実際の窓長は変化しないため，二つのスペクトログラム Xm

(In)

と Yn

のサイズが一致(Ex)する．これによって，二つの行列 Tn Vn と Xm のサイズも一致し，式 (3.5) における問題が解消される．このとき，提案手法のアルゴリズムでは，Fig. 3.2 の逆変換を適用する必要がなくなる．つまり，常に Fig. 3.2 に示す一方向の変換のみを反復して適用することになる．なお，窓関数の有効幅を変化させる方法については，次節で説明する．
改めて，空間モデル及び音源モデルの最適化に用いる窓関数をそれぞれ ω (Ex) ∈ RQ 及び

ω (In) ∈ RQ とする．これらの窓関数は，Fig. 3.3 に示すように，有効幅は異なっているが，実(Ex)際の窓長は同一である．提案手法では，前者の Yn

が空間モデルの最適化に，後者の ω (In)が音源モデルの最適化に用いられる．これらの解像度が異なる二つの時間周波数表現の変換に*2 逆 STFT の対象は複素スペクトログラムであるため，位相情報を失ったパワースペクトログラムに対しては，

逆 STFT をそのまま適用することができない．

26

第3章

提案手法

は，次式で表される，時間領域を経由した方法を使用する．

Yn(In) ← STFTω(In) (ISTFTωe (Ex) (Yn(Ex) ))(3.6)

この処理の概念は，Fig. 3.2 に示す方法と一致する．式 (3.6) の時間周波数表現の変換は，従来の ILRMA における分離行列の反復更新則 (2.37)–(2.40) と，NMF 音源モデルの反復更新
則 (2.41) 及び (2.42) の間に，PB 法と共に挿入される．すなわち，空間モデルの最適化に用い(Ex)た時間周波数表現 Yn(In)から，音源モデルの最適化に用いる分離信号の時間周波数表現 Yn

へと，時間領域を介して解像度の変換を行っている．Fig. 3.3 の方法を基にした上記の処理により，ILRMA の空間モデル及び音源モデルのそれぞれで，独立した時間周波数表現を利用することが可能となる．
式 (3.6) により，提案手法ではスペクトログラム無矛盾性の担保も行っている．これは，分離信号のスペクトログラム Yn を異なる解像度の無矛盾なスペクトログラムの集合に射影する操作であるともいえる．特に，空間モデル最適化用の窓関数 ω (Ex) と音源モデル最適化用の窓
関数 ω (In) が一致する場合，式 (3.6) は consistent ILRMA におけるスペクトログラム無矛盾性を担保する式 (2.46) に一致する．つまり，この場合において，提案手法のアルゴリズムは

consistent ILRMA のアルゴリズムに一致する．従って，提案手法は，consistent ILRMA を複数の解像度の時間周波数表現に一般化した手法とも解釈できる．提案手法を視覚的に表現した図を Fig. 3.4 に示す．ここで，図中の赤色，青色，及び紫色の矢印はそれぞれ STFT の適用，ILRMA における分離行列 Wi の反復更新，及び ISTFT の適用を表す．この図から，提案手法では，空間モデル及び音源モデルの最適化に用いる時間周波数表現が二層構造をなしていることが見て取れる．なお，図中の青色の矢印は，音源モデルの変数更新で得られた分散行
列 Rn を空間モデルの変数更新に使っていることから生じる表現であり，Fig. 3.2 のような時
間周波数表現の変換は行っていない．

3.5 Chebyshev 窓に基づく 2 種類の窓関数の設計
一般的な窓関数の例として Hann 窓を挙げ，その時間波形及び周波数特性を Fig. 3.5 に示す．Fig. 3.5 のように，窓関数の周波数特性にはメインローブ及びサイドローブと呼ばれるピークが存在する．窓関数には，メインローブの幅が狭いこと，そしてサイドローブの大きさが小さいことが求められる．しかし，この二つの特徴の間にはトレードオフが存在する．このことから，トレードオフの限界となる窓関数は何かという問題，すなわちメインローブ幅が与えられた下でのサイドローブピークの最大値が最も小さい窓関数は何かという問題が提起される．窓関数，窓関数の周波数特性，及びメインローブ幅をそれぞれ

ω = [ ω[1], ω[2], . . . , ω[q], . . . , ω[Q] ]T ∈ RQ ，Ω = [ Ω[1], Ω[2], . . . , Ω[i], . . . , Ω[I] ]T ∈ RI ，及び 2ic とすると，この問題は次式で表される．

Minimize ∥sidelobes(Ω)∥∞ s.t.
ω

Q
X
q=1

ω[q] = 1(3.7)

3.5 Chebyshev 窓に基づく 2 種類の窓関数の設計

27

Fig. 3.3. Narrow and wide window functions used in proposed algorithm.

ここで，∥ · ∥∞ 及び sidelobes(·) はそれぞれ Chebyshev ノルム及び窓関数の周波数特性におけるサイドローブの集合であり，次式で定義される．

sidelobes(Ω) = {Ω[i] | i > ic }(3.8)∥A∥∞ = max |x|

(3.9)

x∈A

ただし，A は複素数値からなる集合である．よって，式 (3.7) は次式のようにも書き表せる．

Minimize max |Ω[i]| s.t.
ω

i>ic

Q
X

ω[q] = 1(3.10)

q=1

なお，式 (3.7) を満たす窓関数には定数倍の任意性があるため，Pq ω[q] = 1 を制約条件として定めている．また，式 (3.7) はサイドローブピークの最大値 ∥sidelobes(Ω)∥∞ が定められた
下での，メインローブ幅が最も狭い窓関数を求める問題と等価であるから，次式のように表すこともできる．

Minimize ic s.t. Ω[0] = 1, |Ω[i]| ≤ α ∀|i| ≥ ic
ω(3.11)

ただし，α は与えられたサイドローブピークの内，最大のピーク値である．ここで，式 (3.11)における条件 Ω[0] = 1 は式 (3.7) における条件

Pq ω[q] = 1 に対応している．式 (3.11) を満たす窓関数は Dolph-Chebyshev 窓 [37] と呼ばれる．本論文では，Dolph-Chebyshev 窓を単に Chebyshev 窓と呼ぶ．長さ Q の Chebyshev 窓に DFT を施した関数は次式で与えられる．

Ω[i] =


TQ−1 β cos π(i−1)
Q
TQ−1 (β)(3.12)

28

第3章

提案手法

Fig. 3.4. Overview of proposed ILRMA. Red arrows indicate STFT using ω (In) . Purplee (Ex) . Blue arrows indicate iterative update ofarrows indicate ISTFT using ωparameters in ILRMA.

ここで，β は次式で表される．

β = coshα1cosh−1 10 20
Q(3.13)

3.6 最適化アルゴリズム

29

Fig. 3.5. Waveform and frequency response of Hann window.

また，Tb (·) は第一種 Chebyshev 多項式であり，任意の実数 x に対して次式で定義される．



−1

cos b cos x 
Tb (x) = cosh b cosh−1 x(−1)b cosh b cosh−1 (−x)

if − 1 ≤ x ≤ 1
if x > 1
if x < −1(3.14)

ただし，b は任意の非負整数である．Chebyshev 窓は，式 (3.12) に逆 DFT を施した関数を正規化して得られる．本論文では，サイドローブピークの最大値 α をサイドローブレベルという窓関数のパラメタで与える．このとき，サイドローブレベルは窓関数の有効幅を左右するパラ
メタとなる．

Figs. 3.6(a)–(c) に Q = 64，α = 60 [dB], 200 [dB], 500 [dB] の正規化を行っていない
Chebyshev 窓の時間波形及び周波数特性を示す．Figs. 3.6(a) 及び (b) から，メインローブのピークとサイドローブのピークの差がそれぞれ 60 dB と 200 dB であることが確認できる．
なお，Fig. 3.6(c) において，メインローブのピークとサイドローブのピークの差が 500 dB となっていないのは，計算機イプシロンに起因する現象である．また，Chebyshev 窓はサイドローブが等リプルであるという性質を持っていることも確認できる．さらに，サイドローブレ
ベルが大きくなるにつれて，Chebyshev 窓の有効幅は細くなっていくことも見て取れる．

3.6 最適化アルゴリズム
提案手法における最適化アルゴリズムの擬似コード*3 を Algorithm 1 に示す．Algorithm 1
中の 6 及び 7 行目が NMF 音源モデルの更新，8–10 行目が IP に基づく分離行列の更新，

11–14 行目が mref 番目のチャネルへの PB 法，15 行目が異なる解像度の時間周波数表現の変換を表す．

*3 記法は algorithms という TeX パッケージに準拠する．

30

第3章

提案手法

Fig. 3.6. Waveforms and frequency responses of non-normalized Chebyshev window. Side
lobe levels are (a) 60 dB, (b) 200 dB, and (c) 500 dB. As value of α increases,
window function gets narrower.

3.7 本章のまとめ
本章では，多重解像度時間周波数表現に基づく ILRMA を提案し，その具体的な処理について説明した．その中で，提案手法の利点，実現における問題点，及びその問題への解決策についても述べた．次章では，提案手法の空間モデル及び音源モデルの最適化に用いる窓長を変化

3.7 本章のまとめ

31

Algorithm 1 Proposed ILRMA
Input: {x[l]}L
l=1 , maxIter
Output: {y[l]}L
l=1
NI
1: Initialize {Tn }N
n=1 , {Vn }n=1 , {Wi }i=1
(In)
2: Xm = STFTω (In) (x)
(Ex)
3: Xm
= STFTω(Ex) (x)(In)
H (In)
4: yijn = win
xij ∀i, j, n
5: for iter = 1, 2, . . . , maxIter do (In) .2
. 1[|Yn | ⊙(Tn Vn ).−2 ]VnT 2
6:
Tn ← Tn ⊙∀n(Tn Vn ).−1 VnT T (In) .2
. 12
Tn [|Yn | ⊙(Tn Vn ).−2 ]7:Vn ← Vn ⊙∀n
TnT (Tn Vn ).−1P(Ex) (Ex)H8:
Uin ← J1 j [Tn V1n ] xij xij
∀i, n
9:
10:
11:
12:
13:
14:

win ← (Wi Uin )

−1

i,j

en ∀i, n− 1Hwin ← win win Uin win 2 ∀i, n
λin ← Wi−1 m ,n ∀i, n
refwin ← λin win ∀i, n(Ex)
H (Ex)
yijn ← win
xij ∀i, j, n
[Tn ]i,k ← |λin |2 [Tn ]i,k ∀i, k, n
(In)

(Ex)

Yn ← STFTω(In) (ISTFTωe (Ex) (Yn
16: end for(Ex)
17: y = ISTFTω
)e (Ex) (Yn
15:

)) ∀n

させ，音源分離の実験を行う．そして，得られた結果における傾向の評価，及び，従来手法である consistent ILRMA と提案手法の分離性能の比較を行う．

32

第4章

実験
4.1 まえがき
本章では，提案手法の空間モデル及び音源モデルの最適化に用いる窓長を，サイドローブレ
ベルというパラメタによって，様々に変化させることで，従来手法である consistent ILRMAと提案手法の比較実験を行う．4.2 節では，実験に使用する音源や提案手法におけるパラメタ
設定等の実験条件について述べる．4.3 節では，全実験結果の内，傾向がよく確認できる代表的なデータの実験結果について述べる．4.3 節に掲載しなかった実験結果については付録 A に示す．最後に，4.4 節で本章の総括を行う．

4.2 実験条件
RWCP データベース [38] 収録のインパルス応答 E2A（T60 = 300 ms）による 2 音源の畳み込み混合を行い，10 曲分の 2 チャネル観測信号を生成した．ここで，用いたインパルス応答の収録条件は Fig. 4.1 に示す通りである．そして，これらの観測信号に対する BSS 性能を従来手法の consistent ILRMA 及び提案手法の 2 手法で比較した．Table 4.1 に実験で用いた音
源信号（ドライソース）の詳細を示す．評価指標は音源対歪み比（source-to-distortion ratio:

SDR）[39] の改善量（improvement of SDR: SDRi）を用いた．STFT における窓関数 ω (Ex)及び ω (In) には長さ 256 ms（4096 点）の Chebyshev 窓を用い，サイドローブレベルパラメタ

α を変化させることにより，窓関数の有効幅を変更した．空間モデル及び音源モデルの最適化
における両方の Chebyshev 窓のサイドローブレベルには，Table 4.2 に示す値を用い，これらを総当たり的に変化させた．その他の実験条件は Table 4.2 に示す通りである．なお，STFTの実装には DGTtool [40, 41] を用いた．

4.3 実験結果
各サイドローブレベルに対して，異なる乱数シードを用いて 10 回実験を行った際の平均

SDRi を付録 A に示す．ここで，これらの図における色の濃淡は，暗いほど平均 SDRi が低4.3 実験結果

33

Fig. 4.1. Recording conditions of impulse response E2A.
Table 4.1. Music sources obtained from SiSEC2011 dataset used as dry sourcesSong number

Data nameSource (1/2)

Length [s]1another dreamer-the ones we love

drums/guitar

25.6

2

another dreamer-the ones we love

guitar/vocals

25.6

3

bearlin-roads

acoustic guit main/vocals

14.6

4

bearlin-roads

drums/bass

14.6

5

bearlin-roads

piano/acoustic guit main

14.6

6

fort minor-remember the nameviolins synth/vocals

24.6

7

fort minor-remember the namevocals/drums

24.6

8

tamy-que pena tanto faz

guitar/vocals

13.6

9

ultimate nz tour

drums/vocals

18.6

10

ultimate nz tour

guitar/synth

18.6

く，明るいほど平均 SDRi が高いことを示している．なお，Figs. 4.3–4.6 の太線で囲んだ対角部分は，空間モデルと音源モデルの両方で同一の窓関数を使用していることから，従来のconsistent ILRMA の結果に対応する．また，それ以外の部分（非対角部分）は異なる解像度の時間周波数表現を用いる提案手法の結果に対応する．
まず，実験結果の全体的な傾向を述べる．行方向及び列方向を音源モデル及び空間モデルのサイドローブレベルとして実験結果を表した図は，平均 SDRi が比較的高い場所の分布によって，Figs. 4.2(a) 及び (b) の 2 種類に大別される．これらの図は，Figs. A.1–A.11 に対応し，図中の赤色の部分は，比較的高い平均 SDRi が観測されたパラメタに相当する場所を示してい

34

第4章

実験

Table 4.2. Experimental conditions
Window shift length

32 ms

Side lobe level of{20, 30, 40, 50, 60, 70, 80, 90,
100, 120, 150, 200, 300, 500,

Chebyshev window [dB]

700, 1000, 1500, 2000, 3000}
Number of NMF bases K

10
Wi : identity matrix

Initialization ofparameters

Tn and Vn : uniformly distributed
random values in the range (0, 1)

Number of iterations

200

for parameter updates
Number of trials

10 with diﬀerent random seeds

Reference channel mref

1

for projection back techniqueる．Fig. 4.2(a) は，比較的高い平均 SDRi を示す部分が対角部分の下側に広がっているという性質を表している．また，音源モデルのサイドローブレベルがおよそ 50〜60 dB であると
きに，空間モデルのサイドローブレベルが大きくなっても，比較的高い分離性能を示し続けるという特徴も有している．Figs. A.2 及び A.3 より，Fig. 4.2(a) の傾向は Song 1 及び Song 9の結果に見られる．Fig. 4.2(b) は，対角部分のやや下側から上側にかけて，比較的高い平均

SDRi を示す部分が広がっているという性質を表している．Figs. A.4，A.7，A.8，及び A.11
より，Fig. 4.2(b) の傾向は Song 3，Song 6，Song 7，及び Song 10 の結果に見られる．また，Fig. A.3 では，対角部分の下側と上側両方に平均 SDRi が高い部分が広がっている．よって，Song 2 の結果は Figs. 4.2(a) 及び (b) の傾向を部分的に含んでいるといえる．さらに，Fig. A.6 では，音源モデルのサイドローブレベルが 50 dB，空間モデルのサイドローブレベルが 70〜500 dB という横に細長い領域で，平均 SDRi が高い部分が広がっている．このことから，Song 5 の結果は Fig. 4.2(a) の傾向を部分的に含んでいるといえる．なお，残りの Song

4 及び Song 8 の結果である Figs. A.5 及び A.9 では，平均 SDRi の分布は Figs. 4.2(a) 及び(b) のどちらにも該当しないことが見て取れる．
次に，10 曲の観測信号全ての平均 SDRi，Fig. 4.2(a) の傾向を示す Song 9 の平均 SDRi，及び Fig. 4.2(b) の傾向を示す Song 3 と Song 6 の平均 SDRi を取り上げて，分離性能の評価を行う．10 曲の観測信号全て，Song 3，Song 6，及び Song 9 の平均 SDRi のうち，サイ
ドローブレベルが 40〜1500 [dB] のデータをそれぞれ Figs. 4.3–4.6 に示す．Figs. 4.4–4.6 において，最も高い平均 SDRi を示しているのは，空間モデル用の窓関数 ω (Ex) 及び音源モデ
ル用の窓関数 ω (In) のサイドローブレベルがそれぞれ 50 dB と 60 dB，60 dB と 90 dB，及び 70 dB と 120 dB の場合であり，対角成分（consistent ILRMA の結果）から外れた場所に位置していることが分かる．また，Fig. 4.3 より，10 曲の観測信号全ての平均 SDRi においても，空間モデル用の窓関数 ω (Ex) のサイドローブレベルが 80 dB や 90 dB の場合は，音源

4.3 実験結果

35

Fig. 4.2. Rough classiﬁcation of experimental results by distribution of better SDRigroup. Red parts indicate group of relatively better SDRis. (a) is trend ofSong 1 and Song 9. (b) is trend of Song 3, Song 6, Song 7, and Song 10.

Fig. 4.3. Average SDRi using 10 diﬀerent initializations of all 10 music. Range of side
lobe level is from 40 to 1500.

モデル用の窓関数 ω (In) のサイドローブレベルを 100 dB や 120 dB とした方が，consistent

ILRMA よりも高い性能を示すことが読み取れる．以上から，音源モデル及び空間モデルの最適化に同一の窓関数を用いることが常に最良の結果を与えるとは限らないことが実験的に示さ
れた．

36

第4章

実験

Fig. 4.4. Average SDRi using 10 diﬀerent initializations of Song 3. Range of side lobelevel is from 40 to 1500.

Fig. 4.5. Average SDRi using 10 diﬀerent initializations of Song 6. Range of side lobelevel is from 40 to 1500.

4.4 本章のまとめ
本章では，提案手法の空間モデル及び音源モデルの最適化における窓長を様々に変化させ，音源分離の実験を行った．そして，得られた結果の傾向を評価し，従来手法である consistent

ILRMA との比較を行った．実験結果から，ILRMA の空間モデル及び音源モデルの最適化に

4.4 本章のまとめ

37

Fig. 4.6. Average SDRi using 10 diﬀerent initializations of Song 9. Range of side lobelevel is from 40 to 1500.

同一の窓関数を用いることが必ずしも最良の結果を与えるとは限らないことが示された．また，音源分離の対象とする観測信号に応じて，分離性能の傾向が変化し，そこには少なくとも2 つ以上のパターンが存在することを確認した．

38

第5章

結言
本論文では，ILRMA に対して，空間モデル及び音源モデルの最適化にそれぞれ異なる窓長の窓関数を導入する手法を提案した．これは，ILRMA における空間モデルと音源モデルのそれぞれに適切な窓長が存在するのではないかという考えに基づくものである．提案手法では，
窓関数の有効幅を変化させるという方法により，スペクトログラムのサイズを変えることな
く，見かけ上の窓長を変化させ，空間モデルと音源モデルに異なる時間周波数解像度のスペク
トログラムを用いることで生じる問題を解消した．窓関数のサイドローブレベル（サイドロー
ブピークの最大値）とメインローブ幅（有効幅）がトレードオフの関係にあることに注目し，窓関数の有効幅を変化させるために，サイドローブレベルを設定可能な Chebyshev 窓を用いた．実験結果から，提案手法は consistent ILRMA と比較して，分離精度が向上する場合があることが確認された．
最後に，今後の展望を述べる．4.3 節で述べたように，実験結果には，観測信号毎に少なく
とも 2 つ以上のパターンが存在していた．このことから，実験結果の傾向と観測信号の構成
要素（音源信号の楽器の種類や時間周波数構造等）の関係性を調べることで，各々の観測信号が持つ固有の特徴が提案手法の分離性能に及ぼす影響を明確にできると予想される．これを基に，音源分離に用いる信号の特徴に応じて，空間モデルと音源モデルに最適な窓関数を設定できると考える．また，式 (2.36) で表される ILRMA のコスト関数に，複数の解像度の時間周波数表現を導入することで，提案手法のコスト関数を設計することがおそらく可能である．このコスト関数の設計により，提案手法をより理論的な面から捉えることができると思われる．
以上の，観測信号に応じた窓関数の最適化及び提案手法のより理論的な解析によって，提案手法の更なる分離性能の向上が期待できる．

39

謝辞
本論文は，香川高等専門学校電気情報工学科北村研究室にて行われた研究に基づくものです．
まず，本研究を進めるにあたり，ご多忙のところ熱心にご指導くださいました指導教員の北村大地講師に心より感謝申し上げます．北村大地講師には，論文執筆や研究に関する議論など，細部にわたるまで丁寧にご指導いただきました．研究室に配属されるまで，研究というものに触れたことがなかった私にとっては，大変貴重な経験となりました．また，実験を行う際の計算機サーバや，Notion 上での充実したナレッジ，学生が相談しやすい環境など，研究を行う上で快適な環境作りを徹底して行っていただけたことに対しては，感謝の念に堪えません．
本論の副査である重田和弘教授には，論文の構成や記述に関して大変有益な助言を頂き，大変お世話になりました．ここに厚く御礼申し上げます．
早稲田大学の矢田部浩平講師には，共同研究を通じ，研究に関する議論から，国内学会の原稿の記述に関することまで，多数のご助言をいただきました．心より感謝申し上げます．
北村研究室の先輩方が執筆された論文は本論文を執筆する上で参考にさせていただきました．また，北村研究室同期の川口翔也氏・蓮池郁也氏・溝渕悠朔氏・村田佳斗氏には，ゼミや日頃のディスカッションのほか，1 年に亘る研究室生活を様々な面で支えていただきました．
ここに感謝申し上げます．
最後になりますが，現在に至るまで私の学生生活を金銭的に支え，暖かく見守って下さった
両親には感謝の念に堪えません．これまで本当にありがとうございました．

40

参考文献
[1] H. Sawada, N. Ono, H. Kameoka, D. Kitamura, and H. Saruwatari, “A review ofblind source separation methods: Two converging routes to ILRMA originating from
ICA and NMF,” APSIPA Transactions on Signal and Information Processing, vol.
8, no. e12, pp. 1–14, 2019.
[2] P. Comon, “Independent component analysis, a new concept?,” Signal Processing,vol. 36, no. 3, pp. 287–314, 1994.
[3] A. Hyvärinen, “Independent component analysis: recent advances,” PhilosophicalTransactions of the Royal Society A, vol. 371, 2013.
[4] 浅野 太, “音のアレイ信号処理-音源の定位・追跡と分離-,” 日本音響学会編 音響テクノロ
ジーシリーズ, コロナ社, 2018

[5] P. Smaragdis, “Blind separation of convolved mixtures in the frequency domain,”
Neurocomputing, vol. 22, pp. 21–34, 1998.
[6] N. Murata, S. Ikeda, and A. Ziehe, “An approach to blind source separation basedon temporal structure of speech signals,” Neurocomputing, vol. 41, no. 1-4, pp. 1–24,
2001.
[7] H. Sawada, R. Mukai, S. Araki, and S. Makino, “A robust and precise method forsolving the permutation problem of frequency-domain blind source separation,” IEEE
Transactions on Speech and Audio Processing, vol. 12, no. 5, pp. 530–538, Sep. 2004.
[8] H. Saruwatari, T. Kawamura, T. Nishikawa, A. Lee, and K. Shikano, “Blind source
separation based on a fast-convergence algorithm combining ICA and beamforming,”
IEEE Transactions on Audio, Speech, and Language Processing, vol. 14, no. 2, pp.
666–678, 2006.
[9] A. Hiroe, “Solution of permutation problem in frequency domain ICA using multivariate probability density functions,” Proc. International Conference on Latent
Variable Analysis and Signal Separation, pp.601–608, 2006.
[10] T. Kim, T. Eltoft, and T.-W. Lee, “Independent vector analysis: an extension of
ICA to multivariate components,” Proc. International Conference on Latent Variable
Analysis and Signal Separation, pp. 165–172, 2006.
[11] T. Kim, H. T. Attias, S.-Y. Lee, and T.-W. Lee, “Blind source separation exploiting

参考文献

41

higher-order frequency dependencies,” IEEE Transactions on Audio, Speech, andLanguage Processing, vol. 15, no. 1, pp. 70–79, 2007.
[12] D. D. Lee and H. S. Seung, “Algorithms for non-negative matrix factorization” Proc.
International Conference on Neural Information Processing Systems, pp. 556–562,
2000.
[13] D. R. Hunter and K. Lange, “A tutorial on MM algorithms,” The American Statistician, vol. 58, no. 1, pp. 30–37, 2004.
[14] M. Nakano, H. Kameoka, J. Le Roux, Y. Kitano, N. Ono, and S. Sagayama,“Convergence-guaranteed multiplicative algorithms for nonnegative matrix factorization with β-divergence,” Proc. IEEE International Workshop on Machine Learning
for Signal Processing, 2010, pp. 283–288.
[15] N. Ono and S. Miyabe, “Auxiliary-function-based independent component analysis
for super-Gaussian sources,” Proc. Latent Variable Analysis and Signal Separation,
pp. 165–172, 2010.
[16] N. Ono, “Stable and fast update rules for independent vector analysis based on auxiliary function technique,” Proc. IEEE Workshop on Applications of Signal Processing
to Audio and Acoustics, pp. 189–192, 2011.
[17] 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.
[18] D. Kitamura, N. Ono, H. Sawada, H. Kameoka, and H. Saruwatari, “Determinedblind source separation unifying independent vector analysis and nonnegative matrix
factorization,” IEEE/ACM Transactions on Audio, Speech, and Language Processing,vol. 24, no. 9, pp. 1626–1641, 2016.
[19] D. Kitamura, N. Ono, H. Sawada, H. Kameoka, and H. Saruwatari, “Determinedblind source separation with independent low-rank matrix analysis,” in Audio Source
Separation, S. Makino, Ed., pp. 125–155. Springer, Cham, 2018.
[20] K. Yatabe and D. Kitamura, “Determined blind source separation via proximal splitting algorithm,” Proc. IEEE International Conference on Acoustics, Speech and Signal Processing, pp. 776–780, 2018.
[21] K. Yatabe and D. Kitamura, “Time-frequency-masking-based determined BSS withapplication to sparse IVA,” Proc. IEEE International Conference on Acoustics,
Speech and Signal Processing, pp. 715–719, 2019.
[22] K. Yatabe and D. Kitamura, “Determined BSS based on time-frequency masking
and its application to harmonic vector analysis,” IEEE/ACM Transactions on Audio,
Speech, and Language Processing, vol. 29, pp. 1609–1625, 2021.
[23] J. L. Roux, H. Kameoka, N. Ono, and S. Sagayama, “Fast signal reconstruction frommagnitude STFT spectrogram based on spectrogram consistency,” Proc. DAFx, 2010.

42

参考文献

[24] J. L. Roux and E. Vincent, “Consistent Wiener ﬁltering for audio source separation,”
IEEE Signal Processing Letters, vol. 20, no. 3, pp. 217–220, 2013.
[25] K. Yatabe, “Consistent ICA: Determined BSS meets spectrogram consistency,” IEEE
Signal Processing Letters, vol. 27, pp. 870–874, 2020.
[26] 豊島直，北村大地，矢田部浩平, “スペクトログラム無矛盾性を用いた独立低ランク行列
分析,” 日本音響学会 2020 年秋季研究発表会講演論文集, pp. 291–294, 2020.

[27] D. Kitamura and K. Yatabe, “Consistent independent low-rank matrix analysis fordetermined blind source separation,” EURASIP Journal on Advances in Signal Processing, vol. 2020, no. 46, 35 pages, 2020.
[28] D. Kitamura, N. Ono, and H. Saruwatari, “Experimental analysis of optimal window
length for independent low-rank matrix analysis,” Proc. European Signal Processing
Conference, pp. 1210–1214, 2017.
[29] S. Araki, R. Mukai, S. Makino, T. Nishikawa, and H. Saruwatari, “The fundamental
limitation of frequency domain blind source separation for convolutive mixtures ofspeech,” IEEE Transactions on Speech and Audio Processing, vol. 11, no. 2, 2003.
[30] K. Matsuoka and S. Nakashima, “Minimal distortion principle for blind source separation,” Proc. International Conference on Latent Variable Analysis and Signal Separation, pp. 722–727, 2001.
[31] C. Févotte, N. Bertin, and J.-L. Durrieu, “Nonnegative matrix factorization with the
Itakura-Saito divergence. With application to music analysis,” Neural Computation,
vol. 21, no. 3, pp. 793–830, 2009.
[32] D. FitzGerald, M. Cranitch, and E. Coyle, “On the use of the beta divergence formusical source separation,” Proc. IET Irish Signals and Systems Conference, 2009.
[33] D. Kitamura, H. Saruwatari, K. Yagi, K. Shikano, Y. Takahashi, and K. Kondo,“Music signal separation based on supervised nonnegative matrix factorization withorthogonality and maximum-divergence penalties,” IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences, vol. E97-A, no. 5,
pp. 1113–1118, 2014.
[34] A. J. Bell and T. J. Sejnowski, “An information-maximization approach to blindseparation and blind deconvolution,” Neural Computation, no. 7, pp. 1129–1159,
1995.
[35] S. Douglas and S. Amari, “Natural-gradient adaptation,” in Unsupervised adaptiveﬁltering, Ed. S. Haykin, vol. I, pp. 13–61, Wiley, 2000.
[36] Y. Iwase and D. Kitamura, “Supervised audio source separation based on nonnegative matrix factorization with cosine similarity penalty,” IEICE Transactions on
Fundamentals of Electronics, Communications and Computer Sciences, vol. E105-A,
no. 6, 2022.

参考文献

43

[37] C. L. Dolph, “A Current Distribution for Broadside Arrays Which Optimizes theRelationship between Beam Width and Side-Lobe Level,” Proceedings of the IRE,
vol. 34, pp. 335–348, 1946.
[38] S. Nakamura, K. Hiyane, F. Asano, T. Nishiura, and T. Yamada, “Acoustical sounddatabase in real environments for sound scene understanding and hands-free speech
recognition,” Proc. International Conference on Language Resources and Evaluation,
pp. 965–968, 2000.
[39] E. Vincent, R. Gribonval, and C. Févotte, “Performance measurement in blind audio
source separation,” IEEE Transactions on Audio, Speech, and Language Processing,vol. 14, no. 4, pp. 1462–1469, 2006.
[40] K. Yatabe, DGTtool. Zenodo, 2021, doi:10.5281/ZENODO.5010751.
[41] 矢田部浩平, “短時間フーリエ変換および離散ガボール変換の MATLAB 実装について,”
日本音響学会 2021 年秋季研究発表会講演論文集, pp. 253–256, 2021.

44

発表文献一覧
国内学会
1. 細谷泰稚, 北村大地, 矢田部浩平, “解像度の異なる複数の時間周波数表現を用いた独立
低ランク行列分析,” 日本音響学会 2022 年春季研究発表会講演論文集, 1-1P-2, 2022.

45

付録 A

4 章の実験に対する全結果
Figs. A.1–A.11 に 4 章の実験により得られた結果を示す．ここで，Fig. A.1 は 10 曲の観測信号全ての平均 SDRi，Figs. A.2–A.11 は Song 1〜Song 10 の平均 SDRi である．なお，

4.3 節に掲載した Figs. 4.3–4.6 は，サイドローブレベルが 40〜1500 [dB] のデータである．
それに対して，ここで示す Figs. A.1–A.11 は，サイドローブレベルが 20〜3000 [dB] のデー
タという点において異なっている．つまり，Figs. A.1，A.4，A.6，及び A.10 は，それぞれFigs. 4.3–4.6 を含んだデータとなっている．

Fig. A.1. Average SDRi using 10 diﬀerent initializations of all 10 music.

46

付録 A 4 章の実験に対する全結果

Fig. A.2. Average SDRi using 10 diﬀerent initializations of Song 1.

Fig. A.3. Average SDRi using 10 diﬀerent initializations of Song 2.

47

Fig. A.4. Average SDRi using 10 diﬀerent initializations of Song 3.

Fig. A.5. Average SDRi using 10 diﬀerent initializations of Song 4.

48

付録 A 4 章の実験に対する全結果

Fig. A.6. Average SDRi using 10 diﬀerent initializations of Song 5.

Fig. A.7. Average SDRi using 10 diﬀerent initializations of Song 6.

49

Fig. A.8. Average SDRi using 10 diﬀerent initializations of Song 7.

Fig. A.9. Average SDRi using 10 diﬀerent initializations of Song 8.

50

付録 A 4 章の実験に対する全結果

Fig. A.10. Average SDRi using 10 diﬀerent initializations of Song 9.

Fig. A.11. Average SDRi using 10 diﬀerent initializations of Song 10.