Vol.2021-MUS-132 No.24
2021/9/17

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

非負値行列因子分解を用いた被り音の抑圧
溝渕 悠朔1

北村 大地1

中村 友彦2

猿渡 洋2

高橋 祐3

近藤 多伸3

概要：音楽の生演奏を録音する場合，各音源に近接するようにマイクロホンを配置することが一般的である．これは，各音源のみの音響信号を得ることを目的としているが，実際には他の音源からの音も少なからず混入してしまう．これは一般に被り音と呼ばれる．本研究では，音楽信号の録音時に生じる被り音の抑圧を目的として，非負値行列因子分解を用いた新しい手法を提案する．本手法では，周波数毎の被り音成分の漏れゲインに対してガンマ事前分布を導入し，混合行列を教師なし推定する．推定された混合行列からWiener フィルタを構成することで，被り音を抑圧した信号を得ることができる．シミュレーションによる音楽信号の被り音抑圧実験では，提案手法が従来手法と比べて高い精度で被り音を抑圧できることを示す．

1. はじめに音楽の生演奏を録音する場合，音源数以上のマイクロホンを配置することが一般的である．特に，Fig. 1 に示すよ
うに，楽器本体やボーカル，アンプ等の各音源に近接させてマイクロホンを配置することが多い．これらのマイクロホンは，近接している音源（以後，目的音源と呼ぶ）からの音響信号のみを観測することを意図して配置されている．
しかし，通常はこれらの近接マイクロホンに，目的音源以外の音源（以後，非目的音源と呼ぶ）の音が漏れて混入してしまう．この問題は，一般にクロストークや被り音等と呼ばれている．
ライブ演奏のステージ上でのミキシングでは，サウンド

Fig. 1 Spatial arrangement of sources and close microphones,

エンジニアが各音源の音量バランスや音質を調節し，会場

where M = N = 4. Target sound is contaminated with

用スピーカ及びモニタースピーカを通して，それぞれ観客

bleeding sound from other non-target sources.

及び演奏者に提供する．このようなサウンドエンジニアによる音量及び音質の調整は sound reinforcement（SR）と音の反射を抑えることで，被り音の低減を試みることもあ呼ばれる．近接マイクロホンから得られる各観測信号に被るが，同じ環境で複数の音源が同時に鳴る状況では，被りり音が含まれている場合は SR が困難となり，会場の音質音を完全に防ぐことはほとんど不可能である．

や演奏の質を低下をさせてしまう可能性がある．また，録複数のマイクロホンで観測された信号に対する被り音の音時の被り音は，ライブ演奏後のリミックスにおいても悪

抑圧は，多チャネル音源分離（multichannel audio source

影響を及ぼす問題がある．以上の理由から，サウンドエン

separation: MASS）[1], [2], [3], [4] と呼ばれる問題と類似

ジニアは録音時の被り音を極力回避するために，近接マイしている．但し，次に示す点において MASS とは異なる特クロホンの配置を慎重に行う必要がある．音楽スタジオで徴がある．

の収録では，音源間に遮音壁を設けることや，吸音材等で(a) 観測信号の信号対ノイズ（signal-to-noise: SN）比は比較的高い．ここで，SN 比における信号は目的音源

1

2

3

香川高等専門学校
National Institute of Technology, Kagawa College
東京大学
The University of Tokyo
ヤマハ株式会社
Yamaha Corporationc 2021 Information Processing Society of Japan⃝の成分を指し，ノイズは非目的音源からの混入成分である．

(b) 各マイクロホンにおける目的音源は既知である（観測信号の各チャネルには音源のラベルが付与される）．

1Vol.2021-MUS-132 No.24
2021/9/17

情報処理学会研究報告
IPSJ SIG Technical Report(c) 一般的なマイクロホンアレイと異なり，マイクロホン
間隔が空間的に離れており（例えば 2 m 以上），観測
信号には空間エイリアシングが発生する．

(d) 音楽信号を対象するため，その芸術的価値を低下させないよう高品質な被り音抑圧が要求される．
条件 (a) と (b) は，MASS と比較した際に被り音抑圧の方が容易と考えられる性質である．一方，条件 (c) と (d) はFig. 2 Instantaneous mixture model for bleeding-sound reduc-

被り音抑圧における困難な点であり，特に条件 (c) の空

tion, where M = N = 4. Color brightness in mix-

間エイリアシングはより深刻である．ビームフォーミン

ing matrix Ã shows amplitude level of each element

グ [1], [2] やブラインド音源分離（blind source separation:(brighter is larger). Due to close miking setup, diag-

BSS）[5], [6], [7], [8], [9], [10], [11] 等の多くの高品質な

onal elements in Ã have larger amplitudes compared

MASS はマイクロホン間の位相差を利用しており，空間エwith oﬀ-diagonal elements.

イリアシングが生じる場合は音源を抑圧・分離することが入する．この事前分布の導入については，Cemgil [24] に極めて困難になる．

よって提案された最大事後確率（maximum a posteriori:

観測信号の位相を無視し振幅やパワーのみを用いる

MAP）推定 NMF に基づいている．

MASS もいくつか提案されている [12], [13], [14], [15]．
Togami ら [12] は，多チャネル観測信号の周波数毎の時間

2. 従来手法

チャネル領域に対して非負値行列因子分解（nonnegative

2.1 混合モデル

matrix factorization: NMF）[17], [18] を適用する時間チャ

M 及び N をそれぞれマイクロホン数及び音源数とし，

ネル NMF（time-channel NMF: TCNMF）を提案してい音源信号，観測信号，及び推定信号をそれぞれ次式で表す．

る．TCNMF では，周波数ビン毎の非負混合行列と各音源のアクティベーションを推定しており，条件 (c) や非同期

s̃(t) = [s̃1 (t), · · · , s̃n (t), · · · , s̃N (t)]T ∈ RN

(1)

録音多チャネル信号に対しても音源分離できることが確認

x̃(t) = [x̃1 (t), · · · , x̃m (t), · · · , x̃M (t)]T ∈ RM

(2)

されている [13]．また，TCNMF と同じく位相情報を無視

ỹ(t) = [ỹ1 (t), · · · , ỹn (t), · · · , ỹN (t)] ∈ R(3)した BSS として，線形分離領域多チャネル NMF（linearTN

demixed domain multichannel NMF: DMNMF）[15] があこ こ で ，t = 1, 2, · · · , T ，n = 1, 2, · · · , N ，及 び m =る．DMNMF も TCNMF と同じく周波数ビン毎の非負混

1, 2, · · · , M はそれぞれ離散時間，音源，及びマイクロ合行列を推定する BSS であるが，いずれの手法も音楽信

ホンのインデクスである．被り音抑圧では，1 章で述べた

号の被り音抑圧に対する有効性は検討されていない．音楽

録音条件より，混合系は優決定系（M ≥ N ）となる．本稿

信号の被り音抑圧を目的とした手法では，周波数毎の非負では，以後 M = N の条件のみ取り扱う．

混合行列の対角成分（近接マイクロホン毎の目的音源のゲイン）と非対角成分（近接マイクロホン毎の非目的音源の漏れゲイン）を事前計測する教師あり手法が提案されている [16]．しかし，録音現場での SR のコストを極力抑える
ためには，このような事前計測はあまり望ましくない．ま複数の音源信号が瞬時混合する場合，観測信号及び推定
信号はそれぞれ次式でモデル化される．

x̃(t) = Ãs̃(t)(4)ỹ(t) = W̃ x̃(t)(5)

た，事前計測された混合行列と実際の状態に大きなミス

ここで，Ã ∈ RM ×N 及び W̃ ∈ RN ×M は，それぞれ時不マッチがある場合，被り音抑圧の性能は著しく低下してし変な混合行列及び分離行列である．この時不変瞬時混合モ

まう．

デル (4) を Fig. 2 に示す．観測信号 x̃(t) は，1 章の条件

本研究では，音源やマイクロホンの空間的な配置が不明(b) より各チャネルに目的音源のラベルが付与されている（ブラインド）の条件下で被り音を抑圧することを目的とことから，x̃m (t) は Figs. 1 及び 2 に示すように，m 番目する．但し，訓練データとテストデータのミスマッチを避の音源 s̃m (t)（n = m）の近接マイクロホンの観測信号とけるために，ソロ演奏のデータセット等の学習データや混定義する．従って，x̃m (t) には目的音源 s̃m (t) の成分が主合系の事前計測等は用いない．従って，教師あり学習に基に含まれており，同時に非目的音源 s̃m′ (t) の被り音成分づいた手法 [19], [20], [21], [22], [23] は本研究の対象外であが含まれている．ここで，m′ ̸= m である．このため，Ãる．本稿では，位相に依存しないブラインド被り音抑圧手の対角成分の絶対値（目的音源のゲイン）は大きく，非対法として，TCNMF の改良手法を新たに提案する．被り音

角成分の絶対値（非目的音源の漏れゲイン）は小さくなるの相対的な漏れゲインをモデル化するために，新たに周波（1 章の条件 (a) に対応）
．被り音抑圧は，混合行列 Ã 又は数毎の非負混合行列の対角及び非対角成分に事前分布を導W̃ = Ã−1 なる分離行列 W̃ を，観測信号 {x̃(t)}Tt=1 のみ

c 2021 Information Processing Society of Japan⃝2

Vol.2021-MUS-132 No.24
2021/9/17

情報処理学会研究報告
IPSJ SIG Technical Reportから推定する問題となる．
実際の録音では，音が各マイクロホンに到達するまでの時間差や録音環境の残響により，複数音源の混合は時不変
畳み込み混合となる．本稿では，時不変畳み込み混合をシンプルにモデル化するため，各マイクロホンと各音源間のインパルス応答長（即ち残響長）が，短時間フーリエ変換（short-time Fourier transform: STFT）で用いる窓長より

Fig. 3 Decomposition model of TCNMF, where M = N = 4

も短いと仮定する．この仮定により，複数の音源が時不変

and I = 6. Note that abs(Xm ) is channel-wise time-

畳み込み混合された観測信号は，時間周波数領域では周波

frequency matrix, but Ai and Si are frequency-wise

数毎の時不変複素瞬時混合として表すことができ，観測信

source-channel and time-source matrices, respectively.

号と推定信号はそれぞれ次式となる*1 ．

対値演算 abs(·) はそれぞれ要素毎の q 乗及び要素毎の絶対

xij = Ai sij(6)

yij = Wi xij(7)

ここで，sij ，xij ，及び yij はそれぞれ次式で定義される．

sij = [sij1 , · · · , sijn , · · · , sijN ]T ∈ CN

(8)

xij = [xij1 , · · · , xijm , · · · , xijM ]T ∈ CM

(9)yij = [yij1 , · · · , yijn , · · · , yijN ] ∈ C

(10)

T

N

また，i = 1, 2, · · · , I 及び j = 1, 2, · · · , J はそれぞれ周波数ビン及び時間フレームのインデクスであり，Ai ∈ CM ×N
は周波数毎の複素混合行列である．さらに，sijn ，xijm ，及び yijn はそれぞれ，複素数の音源，観測，及び推定スペク
トログラム Sn ∈ CI×J ，Xm ∈ CI×J ，及び Yn ∈ CI×J の
時間周波数要素である．
一般的なビームフォーミング [1], [2] や BSS [5], [6], [7],[8], [9], [10], [11] は，各マイクロホンへ間の位相差を用いて複素分離行列 Wi を推定する．マイクロホン間隔が大きく.2
M値を表す．従って，x.2
ij 及び sij はそれぞれ {Xm }m=1 と{Sn }Nn=1 のパワースペクトログラムの時間周波数要素となる．DMNMF では，位相情報を無視したパワースペクト
ログラム領域において，式 (6) を周波数毎の非負混合行列

Ai で式 (11) として近似する．さらに，各音源のパワース
ペクトログラムは，NMF を用いて低ランク行列でモデル.2化される．x.2
ij から Ai 及び sij を推定した後に，Wiener

フィルタを用いて推定信号 yij を復元する．

2.3 TCNMF
TCNMF [12] は位相を無視し振幅のみ用いる BSS であり，音声強調に適用されている．一般的な NMF が時間周波数行列を低ランク分解するのに対し，TCNMF は周波数毎の時間チャネル行列を次のように分解する．

Xi ≈ Ai Si ∀i

(15)

M ×J
Xi = [xi1 · · · xij · · · xiJ ] ∈ R≥0(16)

離れている場合は空間エイリアシングが生じるため，ビー

Fig. 3 に示すように，Xi は振幅領域の周波数ビン毎の時

ムフォーミングや BSS では Wi を正確に推定できなくな

×J
間チャネル観測信号，Si ∈ RN
は各音源のゲインの時間≥0る．従って，位相差に基づくビームフォーミングや BSS で被り音を抑圧することは困難となる．

変化を行ベクトルとして含む行列である．Xi から Ai 及
び Si を推定した後に，Wiener フィルタを用いて推定信号を復元する．

2.2 DMNMF

変数 Ai 及び Si の推定は，次の最適化問題となる [18]．

空間エイリアシングに対処するために，位相を無視しパ
ワーのみ用いる BSS の DMNMF が提案されている [15]．

DMNMF は，独立低ランク行列分析（independent low-rankmatrix analysis: ILRMA）[4], [10], [11] の位相非依存な手法である．DMNMF では，観測信号を次式でモデル化する．
.2x.2
ij ≈ Ai sij ∀i, j

(11)

×N
Ai = abs(Ai ) ∈ RM
≥0

(12)

xij = abs(xij ) ∈ RM
≥0

(13)

sij = abs(sij ) ∈ RN
≥0

(14)

ここで，ベクトルや行列のドット付き指数演算 ·.q 及び絶*1

本稿では，ローマン体の信号は複素数を表し，イタリック体の信号は実数を表す．

c 2021 Information Processing Society of Japan
⃝

min
A,S

X

DKL (Xi |Ai Si )is.t. aimn ,sinj ≥ 0 ∀i, m, n, j

(17)

ここで，DKL (·|·) は次式のように定義される．
Xximj
DKL (Xi |Ai Si ) =ximj log P
n aimn sinjm,j!
X
− ximj +aimn sinj(18)n

式 (18) は，Xi 及び Ai Si の類似性を測る一般化 Kullback–Leibler（KL）ダイバージェンスである．また，A 及び S
は，それぞれ {Ai }Ii=1 及び {Si }Ii=1 の集合であり，ximj ，3

Vol.2021-MUS-132 No.24
2021/9/17

情報処理学会研究報告
IPSJ SIG Technical ReportPaimn ，及び sinj は，それぞれ Xi ，Ai ，及び Si の要素であ

ここで，zimnj ∈ N は ximj = e +る．決定系では Ai が M × N の正方行列であるため，最変数，P(z; λ) は確率変数 z ∈ N とパラメータ λ > 0 を持小化問題である式 (17) は，全ての i に対して Ai = I（I はつポアソン分布，Γ(z + 1) = z! はガンマ関数，e は [0, 1)

単位行列）という自明解を持つ．この自明解を回避するたの範囲で一様分布に従う確率変数である．また，zimnj はめに，L0.5 ノルムに基づくスパース正則化が各時間フレー

i, m, n, j に関して互いに独立と仮定する．ポアソン分布

ムに対して導入されている [12]．
XX
DKL (Xi |Ai Si ) + µ
min∥sij ∥0.5に従う確率変数は再生性を持つ．即ち，zn ∼ P(zn ; λn ) かPP
つ x = n zn のとき，周辺分布は p(x) = P(x; n λn ) でA,Si与えられる．従って，Xi の周辺対数尤度は次式となる．

i,js.t. aimn , sinj ≥ 0 ∀i, m, n, j

n zimnj を満たす確率(19)

ここで，µ は正則化のための重み係数，sij は Si の時間軸

log p(Xi ; Ai , Si )
Y X
= log
p(xinm ; zimnj )p(zimnj ; aimn sinj )方向のベクトル（Si = [si1 · · · sij · · · siJ ]）である．

m,j zimnj

Y

= log

3. 提案手法

P(ximj ;P

n aimn sinj )m,j3.1 動機前述の通り，被り音抑圧では観測信号の位相情報が信=−ローチと考えられる．しかしながら，実験で確認する通り精度な推定は難しい問題である．実際に，文献 [15] では，
ステアリングベクトル（Ai の列ベクトル）の事前情報や位相依存の BSS による事前推定を用いて推定精度を安定させている．一方で TCNMF は，位相情報を用いずにある程度の精度で音声強調ができる [13] が，音楽の BSS や被り
音抑圧に対する性能は調査されていない．特に，式 (19) のP
スパース正則化項 i,j ∥sij ∥0.5 は，時間周波数領域における W-disjoint-orthogonality [25]（1 つの時間周波数グリッドに生起する音源数は混合信号であっても高々 1 つ）という仮定に基づいている．これは音声の混合信号には妥当だ"

Si のスパース正則化は，音楽信号の場合には音質を低下させる可能性がある．
提案手法では，TCNMF における Ai の自明解を避けるために，Si ではなく非負混合行列 Ai の対角成分と非対角成分をそれぞれ正則化する．本手法は，被り音の相対的な漏れゲインがガンマ事前分布から生成されると仮定した

MAP 推定と解釈できる．

X

aimn sinjnX#aimn sinj − log Γ(ximj + 1)

aimn と sinj に関する式 (22) の最大化（即ち最尤推定）は，
式 (18) の最小化と等価な問題である [24]．

3.3 被り音の相対漏れゲインの事前分布生成モデル
提案手法では，Ai の自明解を避けるために，Ai の対角成分及び非対角成分に次の事前分布生成モデルを導入する．
δ(aimn − 1)
(m = n)aimn ∼(23)
G(a(m ̸= n)imn ; k, θ)G(a; k, θ) =1
ak−1 e−a/θΓ(k)θk率変数 a ≥ 0，形状母数 k > 0，及び尺度母数 θ > 0 から成るガンマ分布である．また，aimn は i，m，n に関して互
いに独立と仮定する．式 (23) より，Ai の事前分布は次式となる．

p(Ai ; k, θ) =

Y
m,n=m

Y

m,n=m

p(aimn )Yにしている [24]．KLNMF の最小化問題は，ポアソン生成モデルを仮定した最尤（maximum likelihood: ML）推定と等価である．式 (17) は，次の生成モデルを仮定している．

zimnj ∼ P(zimnj ; aimn sinj )P(z; λ) =1
e−λ λzΓ(z + 1)c 2021 Information Processing Society of Japan
⃝

(20)
(21)

p(aimn ; k, θ)m,n̸=mYδ(aimn − 1)G(aimn ; k, θ)m,n̸=m

3.2 KL ダイバージェンスに基づく NMF の生成モデル
divergence-based NMF: KLNMF）の生成モデルを明らか(24)

ここで，δ(a) は Dirac のデルタ関数であり，G(a; k, θ) は確=

Cemgil は，KL ダイバージェンスに基づく NMF（KL-

(22)

n

が，時間周波数領域で多分に重なり合う音楽信号では不適切と考えられる．従って，従来の TCNMF で用いられる

ximj logm,j用できため，位相に依存しない DMNMF は合理的なアプ

DMNMF の最適化は不安定であり，非負混合行列 Ai の高X(25)
式 (25) では，Ai の対角成分は全て 1 となるよう制限され，非対角成分（被り音の相対漏れゲイン）の生起確率は kや θ によって制御される．Fig. 4 に示すように形状母数を

k > 1 とすることで，Ai の全非対角成分において aimn = 0
を避けることができ，Ai = I なる自明解を回避できる．本
稿では，以下 k > 1 のみを考える．
一方，アクティベーション行列 Si には明示的な構造を仮定せず，次式の非負制約事前分布のみを導入する．

4

Vol.2021-MUS-132 No.24
2021/9/17

情報処理学会研究報告
IPSJ SIG Technical Reportsinj ∼ lim

1

β→∞ β

I[0 ≤ sinj ≤ β]

∝ I[0 ≤ sinj ](26)

ここで，β は正規化係数，I[·] は引数が真ならば 1，偽なら0 となる 2 値の関数である．Ai と同様に，sinj も i，n，jに関して互いに独立と仮定されるため，Si の事前分布は次式となる．

p(Si ) =Yp(sinj )Fig. 4 Probabilistic density function of gamma distribution.

n,j

∝

Y

I[0 ≤ sinj ]

(27)

n,j

3.4 MAP 推定のコスト関数

布 (24) に対応する正則化項であり，次式となる．


1
(32)
R(aimn ; k, θ) = −(k − 1) log aimn + aimn
θ

前節の事前分布に基づき，変数 Ai と Si を MAP 推定で求める．事後分布は次のように得られる．

Y

p(Ai , Si ; Xi ) ∝

Y

i

i

3.5 最適化アルゴリズムの導出

p(Xi ; Ai , Si ) p(Ai ; k, θ)p(Si )
|
{z
}|
{z
}
PriorsLikelihood(28)
式 (28) の負の対数を取ると，右辺を以下のように分解で
きる．

J =−

式 (31) の最小化問題は，NMF の最適化で頻繁に用いられる majorization-minimization（MM）アルゴリズム [18], [26]で解くことができる．目的関数 DKL (Xi |Ai Si ) の上限関数
は Jensen’s の不等式を用いて次のように設計される．

DKL (Xi |Ai Si )X

[log p(Xi ; Ai , Si ) + log p(Ai ; k, θ) + log p(Si )]

c

=

i

式 (22)，(25)，(27) を式 (29) に代入すると，コスト関数 Jが次式として得られる．
"
X
X
J =
−ximj log
aimn sinj
n

i,m,j

+

X

+

i,m,n=m+X
i,m,n̸=m+X

aimn sinj X+
aimn sinj−ximj log
ξimnj
ξimnj
n
n

X

X

i,m,j

X

−ximj

X

!

aimn sinj Xξimnj log+
aimn sinj
ξimnj
n
n

!

(33)

ての i，m，j ，n に対する ξimnj の集合を示す．式 (33) の1−(k − 1) log aimn + aimnθI[0 ≤ simn ]

等号成立条件は次式となる．

ξimnj = P(30)

aimn sinj∀i, m, j, nn′ aimn′ sin′ j(34)

式 (33) より，上限関数の最小化問題は次のようになる．

を取る指示関数である．Ai と Si の MAP 推定は式 (30) の最小化問題であり，Ai と Si に関する最小化は次の問題と等価である．

DKL (Xi |Ai Si ) +aimn sinjnc

ここで，I[·] = − log I[·] は引数が真のとき 0，偽のとき ∞i!

ここで，= は定数項の違いを除いて等しいことを表す．まPた，ξimnj > 0 は n ξimnj = 1 を満たす補助変数，Ξ は全

i,n,j

X

aimn sinj +X≡ D+ (Ai , Si , Ξ)#I[aimn = 1]


X
n

i,m,j

aimn sinj + log Γ(ximj + 1)

X

=
≤

n

A,S

−ximj log

i,m,j

(29)

min

X

X

min

A,S,Ξ

X
i

D+ (Ai , Si , Ξ) +X

R(aimn ; k, θ)

i,m,n̸=ms.t. aimn , sinj ≥ 0 ∀i, m, n, j, ξimnj > 0 ∀i, m, n, j,
and aimn = 1 ∀m = n

R(aimn ; k, θ)

i,m,n̸=ms.t. aimn , sinj ≥ 0 ∀i, m, n, j and aimn = 1 ∀m = n(31)
ここで，R(aimn ; k, θ) は Ai の非対角成分のガンマ事前分

c 2021 Information Processing Society of Japan⃝(35)

従って，式 (35) の上限関数を aimn 及び sinj で偏微分し 0 と置き，さらに式 (34) の統合成立条件を ξimnj に代入すること
で，各変数の更新式を導出できる．正則化項 R(aimn ; k, θ)は sinj での偏微分に影響しないため，sinj の更新式は次式のように単純な KLNMF [18] の更新式と等しくなる．

5Vol.2021-MUS-132 No.24
2021/9/17

情報処理学会研究報告
IPSJ SIG Technical ReportPsinj ← sinj

m

∑

ximj

n′ aimn′ sinj ′P

タ化する．具体的には，前処理として観測信号 x̃(t) に対しaimn

m aimn(36)て次のような処理を施す．

一方，非対角成分 aimn （m ̸= n）については，式 (35) のx̃(t) ←

αx̃(t) ∀tv

(44)


上限関数の偏微分から次式が得られる．

X
1
ξimnj
1
−ximj
+ sinj − (k − 1)
+ =0
aimn
aimn
θ

ここで，max(·) は入力集合の最大値を返す．即ち，式 (44)(37)

は {x̃(t)}Tt=1 のダイナミックレンジが ±α となるようにゲj

式 (37) を aimn について整理すると，次式となる．
P(k − 1) + j ximj ξimnj
P
aimn =(38)1
j sinjθ +v = max {abs(x̃(t))}Tt=1(45)

インを調整している．α は式 (19) における µ と同様に，データ近似項と正則化項のバランスを制御するハイパーパ
ラメータとなる．例えば，α を小さな値に設定すると，最適化に対してより強い正則化が課せられる．

式 (38) に統合成立条件 (34) を代入することで，非対角成分 aimn （m ̸= n）の更新式は次式として導出される．
Px(k − 1) + aimn j ∑ ′ a imj′ s ′ sinj
imn in jnPaimn ←
(39)
1
+
s
j inj
θ

3.7 推定信号の再構成提案手法では，従来の TCNMF と同様に，複素数の観測
信号 xijm に対して次式の Wiener フィルタを適用すること
で，複素数の推定信号 Yn を復元する．

式 (36) 及び (39) より，aimn 及び sinj の非負制約は全変数(aimm simj )2xyijn = P
2 ijm
n (aimn sinj )の初期値を非負値とすることで常に満たされる．また，Aiの対角成分の値は制約条件を満たすように 1 で初期化し，他の変数の更新中も常に 1 に固定する．
式 (36) 及び (39) は，以下のように行列形式で実装すると，多くの言語で効率的に計算できる．
(k − 1) + Ai ⊙ AXi Si i SiT
∀i
Ai ←
1
T
θ + 1Sidiag(Ai ) ← [1, 1, · · · , 1]T ∀i
Si ← Si ⊙Xi
ATi Ai S i

AT
i 1

∀i

(46)

なお，aimm = 1 より，式 (46) は次式で実装できる．

Si.2
yijn =
xijm
(47)
.2
A.2
i Si
m,j
ここで，[·]m,j は入力行列の (m, j) 要素を表す．その後，推(40)

定信号 Yn に逆 STFT を適用することで，時間領域の信号ỹn (t) が得られる．推定信号のゲインは次式で復元する．
(41)
ỹ(t) ←
(42)

ここで，⊙ と行列の除算はそれぞれ要素毎の乗算と除算を表す．また，1 は 1 だけを含む M × J の行列，diag(·) はvỹ(t) ∀t
α

(48)

4. 実験
4.1 実験条件

入力された正方行列の対角成分をベクトルとして表す．な

提案手法の性能を評価するために，音楽信号の被り音抑お，式 (40) は Ai の対角成分の値を変更するが，式 (41) に圧の実験を行った．観測信号となる音楽の混合信号は，人よって対角成分は直ちに 1 に上書きされる．また，式 (40)

工的な音楽データセットである songKitamura [27], [28] を及び (42) の反復計算は，コスト関数 (31) を単調減少させ用いた．混合前の音源信号 Sn には，クラリネット（Cl.），

ることが理論的に保証されている．

オーボエ（Ob.），ピアノ（Pf.），トロンボーン（Tb.）の 4 種類の楽器音を用い，M = N = 4 となるように 4 チャネルの観測信号 xij を作成した．このとき，被り音を含む観測信号

3.6 データ近似項と正則化項のバランス提案手法では，Ai の対角成分を 1 に制約しているため，
非対角要素は被り音の相対漏れゲインに対応する．また，KL ダイバージェンスは次式のスケール依存性を持つ．
DKL (αXi |αAi Si ) = αDKL (Xi |Ai Si )

(43)

×N
を模擬するために，周波数毎の非負混合行列 Ai ∈ RM
≥0を用いて，楽器音 sij を次式のように混合した．

xij = Ai sij(49)

但し，Ai の対角要素は 1 とし，非対角要素は (0, 0.2) の範

ここで，α ≥ 0 は任意の係数である．これらの事実より，観囲の一様分布から生成される乱数に設定した．また，異な測信号 {Xi }Ii=1 のゲインは，式 (31) におけるデータ近似項PPi DKL (Xi |Ai Si ) 及び正則化項
i,m,n̸=m R(aimn ; k, θ)

る擬似乱数シードを用いて，上記の方法で 10 種類の観測信間のバランスに影響を与えることが分かる．そこで，提案た. 信号のサンプリング周波数は 44.1 kHz とし，STFT

手法では，観測信号を正規化したうえでゲインをパラメー

では，窓長 4096 点（約 92.9 ms）の Hamming 窓をハーフ

c 2021 Information Processing Society of Japan⃝

号（即ち 10 種類の非負混合行列 Ai ）を用いて実験を行っ6

Vol.2021-MUS-132 No.24
2021/9/17

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

4.2 実験結果
Fig. 5 は 5 つの比較手法の平均 SDR 改善量を示している．ここで，Fig. 5 における従来の TCNMF と提案手法の性能は，ハイパーパラメータを実験的に決定し，本実験において最も性能が高かった µ = 0.56，k = 1.25，θ = 0.6，及び α = 0.006 に設定した場合の結果を掲載している．位相を考慮した BSS である IVA と ILRMA では，被り音を全く抑圧できないことが確認できる．これは，本実験で観測された混合信号が，非負乱数混合行列 Ai で生成されたものであることに起因しており，観測信号の位相が分離行
Fig. 5 Comparison of SDR improvements, where each bar is

列の推定に全く役に立たないためと考えられる．結果的

average over 10 diﬀerent observed mixtures and 4 in-に，IVA や ILRMA の推定信号には多くの人工的な歪みがstrumental sources.

生じ，SDR が大幅に劣化している．DMNMF は被り音を抑圧できる可能性があるが，その平均性能は 0 dB を超え

オーバーラップシフトさせた．周波数ビン数及び時間フなかった．これは，DMNMF における最適化の難しさに起

レーム数はそれぞれ I = 2049 及び J = 109 であった．

因していると思われる．従来の TCNMF と提案手法では，

本実験では，独立ベクトル解析（independent vector

平均 SDR 改善量が 0 dB を超えていることが確認できる．

analysis：IVA）[9]，ILRMA [11]，DMNMF [15]，従来の特に提案手法は，従来の TCNMF よりも 2.5 dB 以上改善

TCNMF [12]，及び提案手法の 5 つの手法を比較した．IVAされている．この改善は音楽演奏の高品質な後処理や SR

及び ILRMA は複素数の分離行列 Wi を推定する手法で位を実現する上で重要といえる．

相を考慮した BSS であり，その他の手法は振幅又はパワー
スペクトログラムのみを用いる位相非依存な BSS である．

5. まとめ

IVA 及び ILRMA の Wi の初期値は，DMNMF や従来の本稿では，音楽の生演奏を近接マイクロホンで録音する

TCNMF 及び提案手法で用いた初期混合行列の逆行列と際に生じる被り音の除去を目指し，TCNMF を改良した新した．また IVA と ILRMA では，反復音源ステアリング

しい手法を提案した．提案手法では，被り音の相対漏れゲ法（iterative source steering：ISS）[30] と呼ばれる数値的インを正則化しており，これはガンマ事前分布を仮定したに安定な分離行列の更新式を用いた．分離行列 Wi 推定後

MAP 推定に基づく TCNMF と解釈できる．模擬的な観測は，式 (7) を用いて推定信号を復元し，さらにプロジェク

信号を用いた実験では，提案手法が最も被り音を抑圧でき

ションバック [31] を適用して周波数毎のスケールを復元しることが実験的に示された．提案手法は 3 つのハイパーパ

た．一方，DMNMF，従来の TCNMF，及び提案手法につラメータ（k ，θ，及び α）を持っているため，ハイパーパ

いては，非負混合行列 Ai の初期値に対角要素が 1，非対角ラメータに対する提案手法の性能変化の解析と，効率的な要素が (0, 0.1) の範囲の一様分布から生成される乱数を用

パラメータのチューニング法が今後の課題として挙げらいた．その他のパラメータについては，(0, 1) の範囲で一れる．

様分布から生成される乱数を用いて初期化した．非負混合行列 Ai 推定後は，Wiener フィルタ (46) を用いて推定信号を得た．ILRMA と DMNMF における NMF 音源モデルの基底ベクトル数 L は 10，30，及び 80 と設定した．各手

謝辞
本研究の一部は，JSPS 科研費 19K20306 及び 19H01116の助成を受けた．

法の最適化アルゴリズムは更新式を 200 回反復計算し，いずれもコスト関数値が収束していることを確認している．
客 観 評 価 尺 度 と し て ，source-to-interference ratio と

参考文献
[1]

sources-to-artiﬁcial ratio の 2 つを加味した総合的な分離品質を示す source-to-distortion ratio（SDR）[29] を用いた．

[2]

1 章の条件 (a) で述べたように，被り音抑圧における観測
信号の SDR（入力 SDR）は高い値を示した．今回の実験[3]では，10 種類の観測信号の Cl.，Ob.，Pf.，及び Tb. の平均入力 SDR がそれぞれ 18.8，15.0，14.7，及び 8.6 dB となった．これらの入力 SDR からの改善量を楽器音毎に計算し，各手法での平均性能を評価した．

c 2021 Information Processing Society of Japan⃝[4]

M. Brandstein and D. Ward, Microphone Arrays: Signal Processing Techniques and Applications, SpringerVerlag Berlin Heidelberg, 2001.
H. L. Van Trees, Optimum Array Processing, John Wiley and Sons, New York, 2002.
X. Yu, D. Hu, J. Xu, Blind Source Separation: Theory and Applications, John Wiley and Sons, New York,2014.
H. Sawada, N. Ono, H. Kameoka, D. Kitamura, and
H. Saruwatari, “A review of blind source separation
methods: Two converging routes to ILRMA originating

7

Vol.2021-MUS-132 No.24
2021/9/17

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

[5]

[6]

[7]

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

[19]

from ICA and NMF,” APSIPA Trans. Signal and Info.
Process., vol. 8, no. e12, pp. 1–14, 2019.
P. Smaragdis, “Blind separation of convolved mixtures in
the frequency domain,” Neurocomputing, vol. 22, pp. 21–
34, 1998.
H. Saruwatari, T. Kawamura, T. Nishikawa, A. Lee, andK. Shikano, “Blind source separation based on a fastconvergence algorithm combining ICA and beamforming,” IEEE Trans. Audio, Speech, and Lang. Process.,
vol. 14, no. 2, pp. 666–678, 2006.
A. Hiroe, “Solution of permutation problem in frequency
domain ICA using multivariate probability density functions,” Proc. Int. Conf. Independent Compon. Anal.
Blind Source Separation, pp. 601–608, 2006.
T. Kim, H. T. Attias, S.-Y. Lee, and T.-W. Lee, “Blindsource separation exploiting higher-order frequency dependencies,” IEEE Trans. Audio, Speech, and Lang.
Process., vol. 15, no. 1, pp. 70–79, 2007.
N. Ono, “Stable and fast update rules for independent vector analysis based on auxiliary function technique,” Proc. IEEE Workshop Appl. Signal Process.
Audio Acoust., pp. 189–192, 2011.
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. Audio, Speech,
and Lang. Process., vol. 24, no. 9, pp. 1626–1641, 2016.
D. Kitamura, N. Ono, H. Sawada, H. Kameoka, and
H. Saruwatari, “Determined blind source separation withindependent low-rank matrix analysis,” in Audio Source
Separation, S. Makino, Ed., pp. 125–155. Springer,
Cham, 2018.
M. Togami, Y. Kawaguch, H. Kokubo, and Y. Obuchi,“Acoustic echo suppressor with multichannel semi-blind
non-negative matrix factorization,” Proc. Asia-Paciﬁc
Signal Info. Process. Assoc. Annu. Summit Conf.,
pp. 522–525, 2010.
H. Chiba, N. Ono, S. Miyabe, Y. Takahashi, T. Yamada,and S. Makino, “Amplitude-based speech enhancementwith nonnegative matrix factorization for asynchronous
distributed recording,” Proc. Int. Workshop Acoustic
Signal Enhancement, pp. 203–207, 2014.
Y. Murase, H. Chiba, N. Ono, S. Miyabe, Y. Takahashi, T. Yamada, and S. Makino, “On microphone arrangement for multichannel speech enhancement basedon nonnegative matrix factorization in time-channel domain,” Proc. Asia-Paciﬁc Signal Info. Process. Assoc.
Annu. Summit Conf., 2014.
T. Taniguchi and T. Masuda, “Linear demixed domain
multichannel nonnegative matrix factorization for speechenhancement,” Proc. IEEE Int. Conf. Acoust., Speech
Signal Process., pp. 476–480, 2017.
O. Das, J. O. Smith, and J. S. Abel, “Microphone crosstalk cancellation in ensemble recordings with maximum
likelihood estimation,” Proc. Audio Eng. Soc. Convention, 2021.
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.
D. D. Lee and H. S. Seung, “Algorithms for non-negative
matrix factorization” Proc. Neural Info. Process. Syst.,
pp. 556–562, 2000.
A. A. Nugraha, A. Liutkus, and E. Vincent, “Multichannel audio source separation with deep neural networks,”
IEEE/ACM Trans. Audio, Speech, and Lang. Process.,c 2021 Information Processing Society of Japan
⃝

[20]

[21]

[22]

[23]

[24]

[25]

[26]

[27]

[28]

[29]

[30]

[31]

vol. 24, no. 9, pp. 1652–1664, 2016.
N. Makishima, S. Mogami, N. Takamune, D. Kitamura,
H. Sumino, S. Takamichi, H. Saruwatari, and N. Ono,
“Independent deeply learned matrix analysis for determined audio source separation,” IEEE/ACM Trans.
Audio, Speech, and Lang. Process., vol. 27, no. 10,
pp. 1601–1615, 2019.
H. Kameoka, L. Li, S. Inoue, and S. Makino, “Supervised determined source separation with multichannel
variational autoencoder,” Neural Comput., vol. 31, no. 9,
pp. 1891–1914, 2019.
N. Makishima, Y. Mitsui, N. Takamune, D. Kitamura,
H. Saruwatari, Y. Takahashi, and K. Kondo, “Independent deeply learned matrix analysis with automatic
selection of stable microphone-wise update and fastsourcewise update of demixing matrix,” Signal Process.,
vol. 178, 107753, 2021.
T. Nakamura, S. Kozuka, and H. Saruwatari, “Timedomain audio source separation with neural networks
based on multiresolution analysis,” IEEE/ACM Trans.
Audio, Speech, and Lang. Process., vol. 29, pp. 1687–
1701, 2021.
A. T. Cemgil, “Bayesian inference for nonnegative matrix factorisation models,” Computational Intelligence
and Neuroscience, vol. 2009, no. 785152, 2009.
O. Yılmaz and S. Rickard, “Blind separation of speech
mixtures via time-frequency masking,” IEEE Trans.
Signal Process., vol. 52, no. 7, pp. 1830–1847, 2004.
Y. Sun, P. Babu, and D. P. Palomar, “Majorizationminimization algorithms in signal processing, communications, and machine learning,” IEEE Trans. Signal
Process., vol. 65, no. 3, 2017.
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. Audio, Speech, and Lang. Process.,
vol. 23, no. 4, pp. 654–669, 2015.
D. Kitamura, “Open dataset: songKitamura,” http://d-kitamura.net/dataset_en.html. Accessed 17 August 2021.
E. Vincent, R. Gribonval, and C. Févotte, “Performancemeasurement in blind audio source separation,” IEEE
Trans. Audio, Speech, and Lang. Process., vol. 14, no. 4,
pp. 1462–1469, 2006.
S. Robin and N. Ono, “Fast and stable blind source separation with rank-1 updates,” Proc. IEEE Int. Conf.
Acoust., Speech Signal Process., pp.236–240, 2020.
N. Murata, S. Ikeda, and A. Ziehe, “An approach toblind source separation based on temporal structure ofspeech signals,” Neurocomputing, vol. 41, no. 1–4, pp. 1–
24, 2001.

8