2-1-9
混合 DOA モデルに基づく多チャンネル複素 NMF による劣決定 BSS ∗☆武田和馬 (筑波大), 亀岡弘和 (東大/NTT CS 研),
澤田宏, 荒木章子 (NTT CS 研), 宮部滋樹, 山田武志, 牧野昭二 (筑波大)

1序 論

ブ ラ イ ン ド 音 源 分 離 (Blind Source Separation;BSS) とは，音源の成分と音源からマイクロホンまでの伝達特性がともに未知のもとで，マイクロホン入力信号から音源成分を復元する技術である。音声信号を対象とした BSS は，ハンズフリーテレビ会議システムなど，多くの応用が期待されている。
以下，まず BSS の問題を定式化する。M 個の信号源から到来する音源信号を J 個のマイクロホンで 観 測 す る 場 合 を 考 え ，j 番 目 の マ イ ク ロ ホ ン で観測される信号の時間周波数成分を yj (ω, t)，m 番目 の 音 源 信 号 の 時 間 周 波 数 成 分 を sm (ω, t) と し ，y(ω, t) = (y1 (ω, t), . . . , yJ (ω, t))T ∈ CJ ，s(ω, t) =(s1 (ω, t), . . . , sM (ω, t))T ∈ CM と す る 。た だ し ，1 ≤ ω ≤ Ω, 1 ≤ t ≤ T は，それぞれ周波数および時刻に対応するインデックスである。ここで，時間周波数分析 (短時間 Fourier 変換など) の分析窓長に比べて，信号源から観測点までのインパルス応答長が十分短いとすると，観測信号モデルは近似的にy(ω, t) = A(ω)s(ω, t) + n(ω, t)(1)のように瞬時混合の形で表すことができる。この近似は，各狭帯域信号において畳み込み混合モデルを仮定していることに相当するため，狭帯域近似という。m 番目の信号源から j 番目のマイクロホンまでの伝達周波数特性 aj,m (ω) を要素にした行列 A(ω) =(aj,m (ω))J×M = (a1 (ω), . . . , aM (ω)) ∈ CJ×M を混合行列と呼び，以下ではこれを時不変と仮定する。
n(ω, t) は，多数の方向から到来する背景雑音や，フレーム長を超える残響成分など，時不変な伝達特性として表現できない成分を表す。このモデルに基づく
BSS は周波数領域 BSS と呼ばれ，時間領域の畳み込み混合モデルに基づく BSS に比べて演算量の少ないアルゴリズムを構成できる点が特長である。
BSS では観測信号だけから未知の混合行列と音源成分を推定する必要があるため，通常は音源に関して何らかの仮定を置き，これにより立てられる規準をもとに両未知変数を最適推定する問題として定式化される。例えば，観測信号数が信号源数以上の場合には，独立成分分析が有効な手法として知られ [1]，音源信号間の独立性を最大化するように分離行列 (混合行列の逆行列) を推定することが目的となる。しかし，信号源の数が観測信号の数よりも多い劣決定の条件下では，たとえ混合過程が既知であったとしても式 (1) を満たす音源成分の解は無数に存在し，独立性の規準ではこの中から一意に解が決められない。このため，劣決定の条件では，音源に関してスパース性などの独立性よりさらに強い仮定が必要となる [2–6]。また一方
で，仮に ω ごとに分離が達成され，M 個の分離信号∗s1 (ω, t), . . . , sM (ω, t) が得られたとしても，これら M個の成分がそれぞれどの信号源のものであるかは一意に決定できず，その対応づけのためには別の手がかりが必要となる。この問題はパーミュテーション整合問
題と呼ばれ，周波数領域 BSS では通常 BSS の後段処理としてパーミュテーションの整合が行われ，分離信号が最終的に得られる。パーミュテーション整合の方法として，これまで音源到来方向 [8]，帯域間の振幅相
関 [9]，それらの組み合わせと調波性 [10]，デルタ振幅成分の帯域間の同期性 [11] などを利用した様々な解法が検討されている。
以上のように，劣決定条件の周波数領域 BSS は，周波数ごとの分離問題の劣決定性とパーミュテーション問題の任意性を同時に有することになる。これまで
我々は，複素 NMF と呼ぶモノラル信号分離手法 [12]において仮定した信号モデルを多チャンネルに拡張し，周波数ごとの信号分離とパーミュテーション整合を同時解決することを目的とした，多チャンネル複素NMF と呼ぶ枠組を提案した。そして，このコンセプトに基づきいくつかの手法を検討してきた [13–16]。一方，近年，非負値行列因子分解 (Non-negative MatrixFactorization; NMF) を多チャンネル化した BSS の手法が提案され，周波数ごとの信号分離とパーミュテーション整合を同時解決できる手法として注目されている [17, 18]。多チャンネル複素 NMF や多チャンネルNMF では，同一音源の信号の周波数成分はコヒーレントに時間変化する，という仮定に基づき観測信号がモデル化されるが，本稿では，この仮定に加え，同一音源の信号は単一方向から平面波として到来する，という仮定を組み込んだ統合モデルを提案する。

2多チャンネル複素 NMF

まず，多チャンネル複素 NMF の基本モデルについて述べる。観測信号に関し，以下を仮定する。
仮定 1. 観測信号は，振幅成分が周波数間でコヒーレントに時間変化する，I 個の要素信号 si (ω, t) からなる。
ここで，要素信号とは，各音源信号をより細かく分解した成分を指す。よって，各要素信号はいずれかの音源に属することになる。i 番目の要素信号の周波数成分比を Hi (ω) とし，これに対応する時刻 t におけるゲインおよび位相スペクトルをそれぞれ Ui (t), ϕi (ω, t)とすると，仮定 1 より，要素信号は

si (ω, t) = Hi (ω)Ui (t)eȷϕi (ω,t)(2)と表される。ただし，∀i :∑Hi (ω) = 1(3)ωUnderdetermined blind source separation with multichannel complex NMF incorporating DOA mixture
model. by Kazuma TAKEDA (University of Tsukuba), Hirokazu KAMEOKA (University of Tokyo/NTTCorporation), Hiroshi SAWADA, Shoko ARAKI (NTT Corporation), Shigeki Miyabe, Takeshi YAMADA,
Shoji MAKINO (University of Tsukuba)

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

2012年3月とする。周波数成分比 Hi (ω) が時不変であることから，si (ω, t) の振幅成分が周波数間でコヒーレントに時間変化するような拘束を与えていることになる。
以上より，i 番目の要素信号と j 番目のマイクロホンとの間の伝達周波数特性を aj,i (ω) とし，ai (ω) =(a1,i (ω), . . . , aJ,i (ω))T とすると，観測信号 y(ω, t) は，y(ω, t) =I∑

ai (ω)si (ω, t) + n(ω, t)θ i = Θz i

ただし，π = (π1 , . . . , πK ) である。zi ∈ {1, . . . , K}
は i 番目の要素信号にどの DOA 候補値が割り当てられるかを表すインジケータ変数であり，上式はこれが離散分布 (各確率値が π1 , . . . , πK ) から生成されることを意味している。このプロセスにより要素信号のDOA が決定され，伝達周波数特性 ai (ω) は，(4)aj,i (ω)|zi ∼ NC ([ℓθi (ω)]j , ν 2 (ω)), (j ≥ 2)i=1と表される。ここで，ai (ω) と si (ω, t) のスケールの任意性を除くため，a1,i (ω) = 1 とする。雑音成分n(ω, t) は，平均が 0，共分散が σ 2 (ω)I の複素正規分
布 NC (0, σ 2 (ω)I) に従うと仮定する。

n(ω, t) ∼ NC (0, σ 2 (ω)I)(5)

また，[12] と同様，ゲイン関数 Ui (t) をスパースにする目的で事前分布として一般化正規分布を仮定する。

Ui (t) ∼ GN (0, υ, γ)(6)

ただし，x ∼ GN (0, υ, γ) ∝ e−|x| /υ である。
観測信号を以上のような要素信号の和で表す考え方
は，複素 NMF [12] における観測信号のモデル化に基づくため，式 (2) の要素信号モデルに基づく BSS 方式
を多チャンネル複素 NMF と呼ぶ。
γ3

次に，以下を仮定する。

ここまで各要素信号の伝達周波数特性を周波数 ω ごとの独立な変数であるかのように扱っていたが，もし上記の仮定が成立するならば，例えばマイクロホン数が 2 の場合，伝達周波数特性の各 ω 間の関係は，到来方向 (Direction-of-Arrival; DOA)θ の関数としてℓθ (ω) =1eȷωd cos θ/czi |π ∼ Discrete(π1 , . . . , πK )(11)を仮定する。ただし，Dirichlet(α0 /K, . . . , α0 /K) はDirichlet 分布であり，0 < α0 < K のときスパースなπ を生成しやすくなる。
以上の，混合モデルに基づく ai (ω) の生成モデルを「混合 DOA モデル」と呼ぶ。

最大事後確率推定

以 上 の 生 成 モ デ ル に 基 づ き ，観 測 信 号 Y :={y(ω, t)}1≤ω≤Ω,1≤t≤T が与えられたもとでの事後確率が最大となる未知パラメータを推定するアルゴリズムを導く。観測信号 Y が与えられた下でのパラメータa := {ai (ω)}1≤i≤I,1≤ω≤Ω
H := {Hi (ω)}1≤i≤I,1≤ω≤ΩU := {Ui (t)}1≤i≤I,1≤t≤T
ϕ： = {ϕi (ω, t)}1≤i≤I,1≤ω≤Ω,1≤t≤T](7)と陽に表される。ただし，0 ≤ θ < 2π ，d をマイクロホンの間隔 (m)，c を音速 (m/s) とする。実際には残響や狭帯域近似の影響で，ai (ω) は上記の理論式からは逸脱する。そこで，ai (ω) は，到来方向 θi が既知のとき，ℓθi (ω) を中心とした複素正規分布より生成されたものと仮定する。到来方向 θi は実際には観測することができないため，これを潜在変数と見なすことにすると，ai (ω) の生成モデルは，DOA を潜在変数とした混合モデルとなる。これを多チャンネル複素 NMF に基づく観測信号の生成モデルに組み込み，観測信号が与えられた下で全体の生成モデルのパラメータ推論を行うことで音源分離とパーミュテーション整合を同時に達成しようというのが提案法の要点である。
まず，DOA を K 個の候補値 Θ1 , . . . , ΘK からなる離散的な変数とする。例えば，180 度を K 等分した角度の集合，Θk = (k − 1)π/K とする。これらの DOAの候補値の中から一つの値が選択され，各要素信号のDOA が決定される，というプロセスをモデル化したい。このプロセスは以下のように記述できる。

日本音響学会講演論文集

π ∼ Dirichlet(α0 /K, . . . , α0 /K)

4

仮定 2. 要素信号は単一方向から平面波到来する(10)により生成される。ただし，[·]j はベクトルの j 番目の要素を表す。
ここで，音源が空間的にスパースに点在するなら，要素信号の到来方向はスパースとなるはずである。従って，zi が従う離散分布をできるだけスパースになるようにしたい。すなわち，1, . . . , K のうち，一部の少数を除くほとんどのものにほぼ 0 の確率値が割り当てられているような離散分布にしたい，ということである。
このためにはスパースな π を生成しやすいような事前分布を設定すれば良い。そこで，

γ

混合 DOA モデル[(9)(12)(13)(14)(15)の事後確率は，

P (a, H, U, ϕ, π|Y ) ∝∑p(a, H, U, ϕ, z, π|Y )z∝ p(Y |a, H, U, ϕ)p(π)p(U )∑p(a|z)p(z|π) (16)zで与えられ，右辺の対数を L(a, H, U, ϕ, π) とすると，

L(a, H, U, ϕπ) = log p(Y |a, H, U, ϕ) + log p(π)∑p(a|z)p(z|π) (17)+ log p(U ) + logzとなる。ただし，z := {zi }1≤i≤I である。ここで，式(4), (5), (6), (8), (10) より，p(Y |a, H, U, ϕ) ∝(8)

- 748 -∏e−∥y(ω,t)−∑

22i ai (ω)si (ω,t)∥2 /σ (ω)ω,tp(a|z) ∝∏∏e−|aj,i (ω)−lj,zi (ω)| /ν (ω)
22j̸=1 i,ω

2012年3月

p(z|π) =∏のとき上不等式の等号は成立する。
最後に，log p(U ) の下界関数を，γ の値に応じて以下のいずれかの不等式

πziip(π) ∝∏πk α0 /K−1kp(U ) ∝∏∏elog p(U ) = −2ρ

(18)log p(Y |a, H, U, ϕ)∑∑2c=−y(ω, t) −ai (ω)si (ω, t) /σ 2 (ω)
ω,t≥−

∑

Vi (t) = Ui (t)L+ (a, H, U, ϕ, π, x, λ, V )∑ 1 ∑1=−2σ (ω) j,i βj,i (ω, t)ω,t+∑xj,i (ω, t) − aj,i (ω)

si (ω, t)のとき上不等式の等号は成立する。
∑次に，log z p(a|z)p(z|π) の下界関数は，負の対数関数が凸関数であることを利用し，Jensen の不等式

p(a|z)p(z|π)

z=∑log∑i≥∑∏∏πk+ (α0 /K − 1)πkλi,k log

e−|aj,i (ω)−lj,k (ω)| /ν (ω)
2λi,k = ∑∏ ∏ −|aj,i (ω)−lj,n (ω)|2 /ν 2 (ω)πnen

日本音響学会講演論文集∑λi,k∏ ∏ −|aj,i (ω)−lj,k (ω)|2 /ν 2 (ω)
ej̸=1 ωfγ [Ui (t)] + Ci,t2∑ Ui,tσ22Hi,ω
+ ω2βj,i,ω,tνωtj̸=1 ωにより設計できる。ここで，λi,k は，i 番目の要素信号がそれぞれの方向からどのくらいの確率で到来し∑ているかを意味する補助変数であり， k λi,k = 1,λi,k > 0 を満たす。また，λi,k がπk∑aj,i,ω =(26)2 ∑∑ xj,i,ω,t Ui,t e−ȷϕi,ω,tσ+ ω2Hi,ωλi,k lj,k,ωβνωj,i,ω,tt∏ ∏ −|aj,i (ω)−lj,k (ω)|2 /ν 2 (ω)
ei,klog πk − 2ρは L(a, H, U, ϕ, π) の補助関数の要件を満たす。ただし，C はパラメータに依らない定数項である。
以下，この補助関数 L+ をもとに 1 < γ ≤ 2 の場合のパラメータの更新則を導く。H に関しては，厳密には式 (3) の条件を考慮に入れる必要があるが，ここでは簡単のため，無制約で L+ を最大化する Hi (ω) を求めたのちに式 (3) を満たすように規格化する方法をとることとする。スペースの都合上，インデックス ω, tを下付き添字で表記することにすると，L+ を最大にするパラメータは(21)2∑kj̸=1 ωk)|aj,i (ω) − lj,k (ω)|2 /ν 2 (ω)

j̸=1 ωilog∑∑−cxj,i (ω, t) = aj,i (ω)1(j̸=1) si (ω, t)(20))(∑+ βj,i (ω, t) yj (ω, t) −aj,i (ω)1(j̸=1) si (ω, t)2(λi,k log πk − log λi,k(19)により設計できる [12]。= は定数項を除いた項の等号を表す。1(·) は括弧内の条件を満たす場合に 1 を返し，満たさない場合に0 を返す指示関数である。ま∑た，βj,i (ω, t) は， i βj,i (ω, t) = 1, βj,i (ω, t) > 0 を満たす範囲で任意に定めても良いパラメータである。
∑xj,i (ω, t) は i xj,i (ω, t) = yj (ω, t) を満たす補助変数であり，∑(25)

xj,i (ω, t) − aj,i (ω)1(j̸=1) si (ω, t)i,k2(24)のとき等号成立する。
以上をまとめると，

1 ∑
12 (ω)σβ(ω,t)j,iω,tj,i1(j̸=1)

fγ [Ui (t)] (23)i,tを用いて設計する。いずれも，補助変数 Vi (t) が2i∑γγ|Vi (t)|γ−2Ui (t)2 + |Vi (t)|γ − γ|Vi2(t)|2
(1 < γ ≤ 2)fγ [Ui (t)] =γ|Vi (t)|γ−1 (Ui (t) − Vi (t)) + |Vi (t)|γ(0 < γ ≤ 1)である。表記の簡略化のため lj,k (ω) := [ℓΘk (ω)]j ,ρ := 1/(2υ γ ) とした。以下，目的関数 L(a, H, U, ϕ)を a, H, U, ϕ に関して最大化する最適化問題を考える。
以上の最適化問題の大域最適解は解析的に得ることは難しいが，[12] と同様，局所最適解を効率的に探索する反復アルゴリズムを補助関数法に基づき導くことができる。補助関数法は，最大化したい目的関数の下界関数 (補助関数と呼ぶ) を設計し，その下界関数を反復的に増加させていくことで間接的に目的関数を増加させていく最適化方法であるが，これを適用するためには目的関数の下界関数を設計する必要がある。
まず，log p(Y |a, H, U, ϕ) の下界関数は以下の不等式|Ui (t)|γ ≥ −2ρi,tti∑c−2ρ|Ui (t)|γHi,ω =j̸=1 ωt∑2Ui,t=1) 2∑ |a1(j̸j,i,ω |j2 ∑ |a1(j̸=1) |2∑ Hi,ω
j,i,ω

σω2

j

βj,i,ω,t(27)βj,i,ω,t=1) ȷϕi,ω,t∑ Hi,ω ∑ Re(x∗j,iω,t a1(j̸)
j,i,ω e2σω jβj,i,ω,tωω- 749 -=1) ȷϕi,ω,t∑ Re(x∗j,i,ω,t a1(j̸)
j,i,ω eβj,i,ω,tj

tUi,t =(22)Ui,tk+ ργ|Vi,t |(28)γ−2

2012年3月

∑∗1(j̸=1)

xj,i,ω,t aj,i,ωjeȷϕi,ω,t = ∑1jλi,k + α0 /K − 1iπk = ∑ ∑λm,n + α0 − K(30)0-1m-2で与えられる。σ 2 (ω) と ν 2 (ω) の更新則も同様に解析的に得ることができるが，以下の実装ではいずれも固定値とした。補助関数法の原理より，以上の更新則によるパラメータの更新と式 (20), (22), (24) による補助変数の更新を交互に繰り返すことで L(a, H, U, ϕ, π)を単調に増加させていくことができる。
上述の反復アルゴリズムにより得られたパラメータと補助変数を用い，以下の手順により分離信号を得る。

-30

1000

2000

3000

4000

5000

6000

7000

8000

5000

6000

7000

8000

fs[Hz]

321degree[rad]n2(29)/βi,j,ω,tdegree[rad]∑∗1(j̸=1)

xj,i,ω,t aj,i,ω

3/βi,j,ω,t0-1

1. (DOA 推定) πk の値より，各音源がどの角度から到来しているらしいかを決定する。

-2

2. (要素信号のグルーピング) 前段の推定結果と λi,kの値から，各要素信号がどの角度から到来していて，どの音源に属しているかを決定する。
3. (分離信号の獲得) 前段のグルーピング結果をもとに，各音源の信号を，補助変数 xj,i (ω, t) を足し合わせることで構成する。

5動作実験

本章では提案法の音源分離性能および DOA 推定性能の検証結果について述べる。
音源信号は，ピアノ，サックスの 2 つの音響信号とし，マイクロホン間の距離を 2.15[cm]，残響時間を130[ms] として測定したインパルス応答を畳込み混合によって人工的に混合した 5[s] の信号を観測信号とした。また，観測信号の時間周波数成分は短時間 Fourier
変換 (標本化周波数 16k[Hz]，フレーム長 128[ms]，オーバーラップ 64[ms]，Hanning 窓) により計算した。提案法の条件については，要素信号数 I = 15，観測雑音の分散 σ 2 (ω) = 1，ステアリングベクトルの分散ν 2 (ω) = 5f /fN ，角度の分割数 K = 90 とした。これらのパラメータは，実験に基づき決定した。補助関数法の繰り返し数は 300 とし，各パラメータ φ の初期値は乱数で設定した。分離信号の精度は，他音声の消し残りによる歪みを表す歪み具合を表す Source toInterference Ratio (SIR) [19] により評価した。単位は dB であり，この数値が高いほど性能が良い。
動作実験の結果，ピアノの SIR は分離前が 2.3[dB]だったのに対し，分離後は 11.4[dB] まで向上し，サックスの SIR は分離前が-2.0[dB] だったのに対し，分離後は 25.8[dB] まで向上した。また，Fig. 1 に混合音の作成に利用したステアリングベクトルの 2 チャンネル間の位相差と，推定ステアリングベクトルのチャンネル間の位相差 arg(a2,1 (ω)) を示す。Fig. 1 より，各音源の到来方向が概ね正しく推定できていることが分かる。

6010002000

30004000fs[Hz]Fig. 1 正解チャンネル間位相差 (上) と推定チャンネ
ル間位相差 arg(a2,i (ω))(下)多チャンネル複素 NMF を提案した。音源分離実験を通して，各音源の伝達周波数特性を概ね正しく推定できたことから，パーミュテーション整合が適切に行えていることを確認した。

参考文献[1] 澤田他, 信学誌, Vol. 91, No. 4, pp. 292–296, 2008.
[2] Yilmaz & Rickard, IEEE Trans. SP, 52(7),pp. 1830–1847, 2004.

[3] Mandel et al., Adv. Neural Inf. Proc. Sys., 2006.
[4] Araki et al., Signal Process., 87, pp. 1833–1847,2007.

[5] Mori et al., Proc. IWAENC’05, pp. 229–232, 2005.
[6] 和泉他, 音講論 (春)’07, 2-1-5, pp. 555–556, 2007.
[7] 小笠原基 音講論 (春)’09,2-9-17, pp. 695–698, 2009[8] Kurita et al., Proc. ICASSP’00, pp. 3140–3143,
2000.

[9] Murata et al., Neurocomputing, 41(1–4), pp. 1–24,2001.

[10] Sawada et al., IEEE Trans. SAP, 12(5), pp. 530–538, 2004.

[11] 小野, 音講論 (秋)’10, 2-10-7, pp. 581–582, 2010.
[12] 亀岡他, 音講論 (秋)’08, 2-8-13, pp. 657–660, 2008.
[13] 北野他, 音講論 (秋)’09, 2-4-15, pp. 645–646, 2009.
[14] Sawada et al., Proc. ICASSP’11, pp. 229–232, 2011.
[15] Sawada et al., Proc. WASPAA’11, pp. 153–156,2011.

[16] 武田他, 音講論 (秋)’11, 1-Q-19, pp. 801–804, 2011.
[17] Ozerov & Févotte, IEEE Trans. ASLP, 18(3),pp. 550–563, 2010.

まとめ

本稿では，劣決定 BSS における，劣決定性およびパーミュテーション問題に起因する解の非一意性を同時に解消する手法として，混合 DOA モデルに基づく

日本音響学会講演論文集-3[18] Sawada et al., Proc. ICASSP’12, to appear.
[19] Vincent et al., Proc. ICA’07, pp. 552–559, 2007.

- 750 -2012年3月