1-R-14
非負値テンソル因子分解に基づく分散マイクロホンアレイを用いたスポットフォーミング ∗☆綾野翔馬 (香川高専), 李莉, 関翔悟 (サイバーエージェント), 北村大地 (香川高専)

1

はじめに複数の音源の音が混ざりあった信号から特定の音源の信号を抽出する音源分離技術は広く用いられている．例えば，会議において特定の話者の音声のみを取り出すことができれば議事録の作成が容易になる．
音源分離の手法の一つに，複数の同期マイクロホンで構成されるマイクロホンアレイを用いたビームフォーミング (beamforming: BF) があり，遅延和 BF や最小分散無歪応答 (minimum variance distortionlessresponse: MVDR) BF 等が活用される．BF はマイクロホンアレイから見て特定の方位に存在する音源を強調するため，同一方位に複数の音源がある場合はそれらを分離することができない．そこで，Fig. 1のように複数のマイクロホンアレイを用いて，空間上の特定の領域に存在する目的音源を強調するスポットフォーミングと呼ばれる技術が提案されている [1]．
スポットフォーミングはいくつかの手法が提案されており，マイクロホンアレイの配置を最適化する手法[2] や，分散して置かれたマイクロホンアレイ内の全マイクで空間フィルタを設計する手法 [1] がある．しかしながら，前者はマイクロホンアレイが容易に移動可能でなければならない点，後者はマイクロホンアレイ間の完全な同期が必要な点において大がかりなシステムが要求される．そこで，非負値行列因子分解 (nonnegative matrix factorization: NMF) [3] を用いた，マイクロホンアレイ同士の同期の条件が緩和されたスポットフォーミング [4] が提案された．本論文では，分散マイクロホンアレイを用いたスポットフォーミングとして，非負値テンソル因子分解 (nonnegativetensor factorization: NTF) を用いた手法を提案する．

2

分散マイクロホンアレイを用いたスポットフォーミング

本稿では Fig. 1 のように目的音源と干渉音源が存在し，また複数のマイクロホンアレイが分散して置かれている状況を考える．ここで，各マイクロホンアレイ内では同期録音しているが，マイクロホンアレイ間では非同期である状況を想定している．Fig. 1の場合，各マイクロホンアレイは正面方向に BF を行うことで，目的音源を強調できるが，同時に同一方位上に存在する干渉音源も強調されてしまう．Fig. 1のような状況では，この干渉音源がマイクロホンアレイ毎に異なるため，各 BF の出力信号間で共通する∗音源成分を抽出できれば，目的音源のみの信号が得られると考えられる．

2.2

信号モデル

本稿では，文字 D について，D を 3 階テンソル，D を行列，D をベクトルとし，(a, b, c) 要素を Da,b,cで表す．また信号については，小文字と大文字でそれぞれ時間領域と時間周波数領域を表す．マイクロホンアレイのインデクスを a = 0, 1, · · · , A − 1，マイクロホンアレイ内のマイクロホンのインデクスを
m = 0, 1, · · · , M − 1，離散時間インデクスを t =0, 1, · · · , T − 1 とする．a 番目の各マイクロホンアレ
イの M 個のマイクの観測信号 x(a) ∈ RM ×T に対して，方位 θa を強調する任意の BF fθa : RM ×T → RTを適用する処理を次式で表す．

y(a) = fθa (x(a) )(1)

BF の出力 y(a) ∈ RT に対して短時間 Fourier 変換(short-time Fourier transform: STFT) を行い，スペクトログラム Y (a) ∈ CI×J を得る．ここで，i =0, 1, · · · , I − 1 は周波数ビン，j = 0, 1, · · · , J − 1 は時間フレームである．
2.3

従来法(a)

想定する状況

2.1

Fig. 1 Situations and signals estimated by twobeamformers.

文献 [4] の手法では，BF の出力 Yi,j のマイクロホンアレイ方向と時間フレーム方向を結合して 1 次元とし，時間周波数マスクを作成することでスポットフォーミングを行う．結合した振幅スペクトログラムを C (conv) ∈ RI×N
≥0 とすると次式で与えられる．
(conv)

Ci,n(a)(conv):= Ci,aJ+j= Yi,j(2)

ここで，N = AJ であり，n = 0, 1, · · · , N − 1 はC (conv) の列インデクスである．さらに，この行列C (conv) に対して次式のように NMF を適用する．
P(conv)Ci,n≃ k Ti,k Ṽn,k (C (conv) ≃ T Ṽ T )(3)

Spotforming using distributed microphone array based on nonnegative tensor factorization. By Shoma
AYANO (NIT Kagawa), Li LI, Shogo SEKI (CyberAgent, Inc.), and Daichi KITAMURA (NIT Kagawa).

日本音響学会講演論文集- 137 -

2024年3月

P
ここで，T ∈ {T ∈ [0, 1]I×K | i Ti,k = 1} 及び Ṽ ∈×K
RN
は基底行列及びアクティベーション行列であ≥0り，k = 0, 1, · · · , K − 1 は NMF の基底ベクトルのインデクスである．NMF の目的関数は次式となる．

P(conv) P
DCminimizeTṼ
i,kn,ki,ni,nkT ,Ṽs.t. Ti,k , Ṽn,k ≥ 0 ∀i, k, n(4)文献 [4] では，距離関数 D に Euclid 距離を使用している．式 (4) を最小化する基底行列 T 及びアクティ
ベーション行列 Ṽ を求めた後，Ṽ を用いてバイナリ
マスク行列 H̃ ∈ {0, 1}J×K を次式のように作成する．
(1 (if ṼaJ+j,k > µ ∀a)(5)H̃j,k =0 (o/w)
ここで，µ ∈ R≥0 は閾値である．すなわち，このバイナリマスクはすべてのマイクロホンアレイで µ より大きいアクティベーションを持つ成分を共通の目的音源成分とみなし 1 としている．このバイナリマ
スクを用いて目的音源の推定振幅スペクトログラム
E (a)(conv) ∈ RI×J
≥0 を次式で得る．
(a)(conv)

Ei,j=Pk Ti,k H̃j,k ṼaJ+j,k(6)

これは各マイクロホンアレイの信号に対して同じ時間
周波数バイナリマスク H̃ を適用していることに相当する．最後に，各マイクロホンアレイの E (a)(conv) に対して位相を与えて逆 STFT を行い，時間ズレを補正して和をとることで目的音源 e(conv) ∈ RT を得る．
従来法では，NMF における基底ベクトル数 K を観測信号に応じて適切に調整しなければならない．また，式 (5) のバイナリマスクは，マイクロホンアレイ間の時間フレームを跨ぐ（STFT のシフト長を超える）ほどの時間ずれは考慮されていないため，許容できる非同期性（録音開始時刻のずれ及びサンプリング周波数のずれ）の条件がやや厳しい．そこで，次章
では NTF を用いることでこれらの制約を緩和した提案法を説明する．

3提案法

3.1

非負値テンソル表現を用いたモデル化の利点

提案法は，各 BF の出力信号の振幅スペクトログラ
ムで 3 階テンソル C(prop) ∈ RA×I×Jを構成する．
≥0(prop)

Ca,i,j(prop)(a):= Yi,j(7)

この 3 階テンソル C
は式 (2) と異なり，マイクロホンアレイ及び時間フレームの物理的な次元を維持している．提案法は，この C(prop) に対して次節のNTF を適用し，分配行列，基底行列，アクティベーション行列の 3 要素に分解する．K 個の基底ベクトルは分配行列によって自動的に各マイクロホンアレイに振り分けられるため，観測信号の性質に強く依存

日本音響学会講演論文集せず適切に共通成分を抽出できることが期待される．
さらに，分配行列を用いて時間フレームに非依存なバイナリマスクを構成できるため，原理的には完全に非同期な A 個のマイクロホンアレイでもスポットフォーミングが可能となる．

3.2

NTF に基づくスポットフォーミング

提案法の処理の流れを Fig. 2 に示す．提案法では，
式 (7) で定義される観測非負値テンソル C(prop) に対して NTF を適用する．このとき次式のように，分配P行列 Z ∈ {Z ∈ [0, 1]A×K | a Za,k = 1}，基底行列
J×K
T ，及びアクティベーション行列 V ∈ R≥0の積に分解する．
(prop)

Ca,i,j≃Pk Za,k Ti,k Vj,k(8)この NTF における Z は，K 個の基底ベクトルを A個のマイクロホンアレイに分配する役割を持つ．この NTF の目的関数は次式となる．
P(prop) P
minimize
a,i,j D Ca,i,j |k Za,k Ti,k Vj,kZ,T ,Vs.t. Za,k , Ti,k , Vj,k ≥ 0 ∀a, i, j, k(9)次節に示す反復更新式で Z ，T ，V を求めた後，Zを用いてバイナリマスク H ∈ {0, 1}K を構成する．まず，αk = min{Z0,k , Z1,k , · · · , ZA−1,k } を各 k について計算し，α0 , α1 , · · · , αK−1 を降順に並べた際の先頭
γ 個のインデクスの集合 Γ ∈ {Γ ⊆ {k ∈ Z|0 ≤ k <K}|n(Γ) = γ} を求める．ここで，n(·) は．集合 · の要素数を表す．H は Γ を用いて次式で表される．
(1 (if k ∈ Γ)
Hk =(10)0 (o/w)
すなわち，k 番目の基底ベクトルは αk が大きいほど全てのマイクロホンアレイで共通に使用されることを表すため，それらの基底ベクトルが目的音源成分であるとみなすバイナリマスクを構成している．式(10) では αk が大きい順に γ 個の基底ベクトルを目的音源成分として選択している．
このバイナリマスクを用いて目的音源を推定する Wiener フィルタを構成する．この処理は，提案法による目的音源の推定複素スペクトログラムをE (a)(prop) ∈ CI×J とすると次式となる．
P2(a)(prop)(a)k (Hk Za,k Ti,k Vj,k )Ei,j= P
(11)
2 Yi,j(ZTV)a,k i,k j,kk
最後に，各マイクロホンアレイの E (a)(prop) に逆STFT を適用し，時間ズレを補正して和をとること
で目的音源 e(prop) ∈ RT を得る．

3.3

NTF における各変数の反復更新式の導出

本節では，文献 [3] の補助関数法を用いて式 (9) の最適化問題を解く反復更新式を導出する．目的関数 D- 138 -

2024年3月

Fig. 2

Process ﬂow of the proposed spotforming method.には次式で表される定数項を除いた一般化 Kullback–Leibler 擬距離を用いる．

D(y|x) = −y ln x + x
全体の目的関数を新たに G と定義する．
P(prop) P
G = a,i,j D(Ca,i,j | k Za,k Ti,k Vj,k )

(12)

(13)

補助関数法を適用するために補助変数 U ∈ {Ua,i,j,k ∈P[0, 1]A×I×J×K | a,i,j,k Ua,i,j,k = 1} を導入する．
Jensen の不等式を適用して得られる補助関数を G とすると次式となる．

G ≤ G+hP(prop) PZa,k Ti,k Vj,k=−CUloga,i,j
k a,i,j,ka,i,j
Ua,i,j,kP+ k Za,k Ti,k Vj,k ](14)∂G + /∂Za,k = 0 より次式を得る．
oP n(prop)
1i,j −Ca,i,j Ua,i,j,k Za,k + Ti,k Vj,k = 0Za,k =(prop)
i,j Ca,i,j Ua,i,j,kP

i,j Ti,k Vj,k(16)補助変数 U に G + = G となる値（補助変数の等号成立条件）を代入することで，次の反復更新式を得る．
P(prop) P Ti,k Vj,ki,j Ca,i,j
l Za,l Ti,l Vj,lPZa,k ← Za,k(17)i,j Ti,k Vj,k同様にして，Ti,k 及び Vj,k の反復更新式を得る．

P(prop) P Za,k Vj,ka,j Ca,i,j
l Za,l Ti,l Vj,lPTi,k ← Ti,kPVj,k ← Vj,ka,j Za,k Vj,k(prop) P Za,k Ti,ka,i Ca,i,j
l Za,l Ti,l Vj,lPa,i Za,k Ti,k(18)(19)

従って，Z ，T ，V を初期化し，式 (17)–(19) を反復更新することで，Z ，T ，V が推定される．

4実験

4.1

実験条件

本稿では，Pyroomacoustics [5] を用いたシミュレー
ションにより従来法及び提案法のスポットフォーミン

日本音響学会講演論文集

Table 1

Dry sources

File name
84 121123 000008 000002.wav
652 130737 000012 000000.wav
3000 15664 000020 000005.wav
1272 141231 000024 000005.wav

SourceTargetInterf. 0
Interf. 1
Interf. 2

Duration [s]
3.38
4.05
4.07
3.47(15)

これを Za,k について解くと，次式となる．

PFig. 3 Condition of room, source, and microphonearray locations.

グの性能を比較する．部屋の形状と音源及びマイクロホンアレイの配置条件を Fig. 3 に示す．Target は目的音源，各 Interf. は干渉音源を表す．ここでは，壁面反射のみを考慮する 2 次元の鏡像法を用いている．
壁面の反射回数と反射係数を調節し，T60 = 0 ms 及び T60 = 512 ms の 2 種類の部屋を作成した．各音源のドライソースには Table 1 に示す LibriTTS [6] の開発データセット (クリーン) の一部を利用した．サンプリング周波数は 16 kHz とした．
各マイクロホンアレイは 4 個のマイクロホンを間隔 2.83 cm で等間隔に配置して構成した．すべてのマイクロホンアレイは完全同期とし，BF には MVDRを用いた．ただし，理想的な条件として MVDR の目的音源のステアリングベクトルには Target の直接音のインパルス応答を用い，雑音分散共分散行列の計
算には全 Interf. からの信号の和を用いた．STFT の窓長及びシフト長はそれぞれ 32 ms 及び 16 ms とし，窓関数には Hann 窓を用いた．NMF 及び NTF の基
底ベクトル数は K = 30，反復更新回数は 100 回とした．なお本稿では，条件を揃えるために，従来法では NMF の距離関数に式 (12) を用いた反復更新式を採用し，さらにバイナリマスクを作成した後は提案法と同様の Wiener フィルタを適用した．従来法のT 及び Ṽ と提案法の T 及び V の初期値は区間 (0,1) の一様分布乱数を用い，提案法の Z の初期値は全て 0.5 とした．客観評価値には，10 種類の乱数初期- 139 -2024年3月

Fig. 4 Box plots of SDR values with 10 random initializations: (a) conventional and (b) proposed methods
when T60 = 0 ms.

Fig. 5 Box plots of SDR values with 10 random initializations: (a) conventional and (b) proposed methods
when T60 = 512 ms.
値に対する信号対歪み比 (source-to-distortion ratio:SDR) [7] を用いた．

4.2

実験結果

参考文献

SDR の箱ひげ図を Figs. 4 及び 5 に示す．前者はT60 = 0 ms，後者は T60 = 512 ms の結果である．青い箱は四分位範囲，箱の内部の赤線及び青丸はそれぞれ中央値及び平均値，破線は外れ値を除く値の範囲，赤い十字は外れ値を表す．また，黒の水平線は BF 出力の SDR を示す．
両方の残響条件で提案法の SDR の中央値が従来法を上回っていることが分かる．一方で，従来法は µ による SDR のばらつきが小さいのに対し，提案法は同じ γ を用いても SDR のばらつきが大きくなっている．
特に T60 = 512 ms の場合では，提案法のばらつきが顕著である．従って，提案法は NTF の各変数に適切な初期値を与えること，適切な γ を選択することができれば，高い性能が得られると考えられる．

5

おわりに

本稿では，分散マイクロホンアレイと NTF を用いた新しいスポットフォーミング手法を提案した．提案
法は NTF の初期値によって性能のばらつきが大きいものの，中央値では従来法を上回ることを確認した．
今後の課題として NTF に適切な初期値を与える方法の検討，最良の γ を自動決定する手法等が挙げられ

日本音響学会講演論文集

る．また，非同期の実録音信号を用いて，実際の応用における性能についても調査する．

[1] M. Taseska and E. A. P. Habets, “Spotforming: spatial ﬁltering with distributed arrays for positionselective sound acquisition,” IEEE/ACM Trans. Audio, Speech, Lang. Process., vol. 24, no. 7, pp. 1291–1304, 2016.
[2] K. Sekiguchi, Y. Bando, K. Itoyama, and K. Yoshii,“Layout optimization of cooperative distributed microphone arrays based on estimation of source separation performance,” J. Robotics and Mechatronics,
vol. 29, no. 1, pp.83–93, 2017.
[3] D. D. Lee and H. S. Seung, “Algorithms for nonnegative matrix factorization,” Proc. NeurIPS, pp.
556–562, 2000.
[4] Y. Kagimoto, K. Itoyama, K. Nishida, and K.
Nakadai, “Spotforming by NMF using multiple microphone arrays,” Proc. IROS, pp. 9253–9258, 2022.
[5] R. Scheibler, E. Bezzam, and I. Dokmanić, “Pyroomacoustics: a Python package for audio room
simulation and array processing algorithms,” Proc.
ICASSP, pp. 351–355, 2018.
[6] H. Zen, V. Dang, R. Clark, Y. Zhang, R. J. Weiss,
Y. Jia, Z. Chen, and Y. Wu, “LibriTTS: a corpus
derived from librispeech for text-to-speech,” Proc.
Interspeech, pp. 1526–1530, 2019.
[7] E. Vincent, R. Gribonval, and C. Févotte, “Performance measurement in blind audio source separation,” IEEE Trans. Audio, Speech, Lang. Process.,
vol. 14, no. 4, pp. 1462–1469, 2006.

- 140 -

2024年3月