1-Q-19(f)
音源の W-DO 性を仮定した多チャンネル複素 NMF による劣決定 BSS∗
☆武田和馬 (筑波大)，亀岡弘和，澤田宏, 荒木章子 (NTT CS 研)，
山田武志，牧野昭二 (筑波大)る [2–6]。また一方で，仮に ω ごとに分離が達成され，K1

序 論

個の成分 s1 (ω, t), · · · , sK (ω, t) が得られたとしても，こブラインド音源分離 (BSS) とは，音源の成分と音源かれらの成分がそれぞれどの信号源のものであるかは一意

らマイクロホンまでの伝達特性がともに未知のもとで， に決定できず，その対応づけのためには別の手がかりがマイクロホン入力信号から音源成分を復元する技術であ 必要となる。この問題はパーミュテーション問題と呼ばる。音声信号を対象とした BSS は，ハンズフリーテレビ

れ，これまで音源到来方向 [7]，帯域間の振幅相関 [8]，そ会議システムなど，多くの応用が期待されている。

れらの組み合わせと調波性 [9]，デルタ振幅成分の帯域間

以下，K 個の信号源から到来する音源信号を M 個のマの同期性 [10] などを利用した様々な解法が検討されていイクロホンで観測する場合を考え，BSS の問題を定式化ることから分かるとおり，周波数領域 BSS における重要する。m 番目のマイクロホンで観測される信号の時間周課題の一つである。

波数成分を ym (ω, t)，k 番目の音源信号の時間周波数成分を sk (ω, t) とし，y(ω, t) = (y1 (ω, t), · · · , yM (ω, t))以上のように，劣決定の周波数領域 BSS は，周波数ごとの分離問題の劣決定性とパーミュテーション問題の任∈CM ，s(ω, t) = (s1 (ω, t), · · · , sK (ω, t))T ∈ CK とする。 意性を同時に有することになるが，以上の解の非一意性
ただし，1 ≤ ω ≤ Ω, 1 ≤ t ≤ T はそれぞれ周波数および は，(1) 音源の時間周波数成分のスパース性を仮定した観時刻に対応するインデックスである。ここで，時間周波 測 (混合過程) モデル，(2) 振幅が周波数間でコヒーレン
T

数分析 (短時間 Fourier 変換など) の分析窓長に比べて， トに時間変化するような拘束をもつ音源モデル，を統合信号源から観測点までのインパルス応答長が十分短いと 化することで効果的に解消できる可能性がある。これを実現するアプローチとして，これまで我々は多チャンネル

すると，観測モデルは近似的に複素 NMF と呼ぶ劣決定 BSS の解法の枠組を提案し，こy(ω, t) = A(ω)s(ω, t) + n(ω, t)(1)の思想に基づきいくつかの手法を検討してきた [11, 12]。
また，近年，同様の考え方に基づき，NMF(Non-negativeの よ う に 瞬 時 混 合 の 形 で 表 す こ と が で き る 。信 号源 k からマイクロホン m までの伝達周波数特性

am,k (ω) を要素にした行列 A(ω) = (am,k (ω))M ×K =(a1 (ω), · · · , aK (ω)) ∈ CM ×K を混合行列と呼び，以下ではこれを時不変と仮定する。n(ω, t) は，多数の方向かMatrix Factorization) を多チャンネル化した BSS の手法が提案されている [13]。本稿では，この枠組の一環として，音源の混合過程に関し音源成分間の W-DisjointOrthogonality (W-DO 性) [2](スパース性の一種) を仮定した多チャンネル複素 NMF を提案する。

ら到来する背景雑音や，フレーム長を超える残響成分など，時不変な伝達特性として表現できない成分を表す。
このモデルに基づく BSS は周波数領域 BSS と呼ばれ，時間領域の畳み込み混合モデルに基づく BSS に比べて演算量の少ないアルゴリズムを構成できる点が特長である。

BSS では観測信号だけから未知の混合行列と音源成分

2多チャンネル複素 NMF

2.1

音源信号モデル

以下，式 (1) の観測モデルを想定し，n(ω, t) を時間，周波数および空間的に無相関な Gauss 性雑音n(ω, t) ∼ NC (0, σ 2 (ω)I)

を推定する必要があるため，通常は音源に関して何らか(2)の仮定を置き，これにより立てられる規準をもとに両未知変数を最適推定する問題として定式化される。例えば， とする。また，A(ω) と s(ω, t) の間のスケールおよび符観測信号数が信号源数以上の場合には，独立成分分析が 号の任意性を除くため，有効な手法として知られ [1]，音源信号間の独立性を最∥ak (ω)∥22 = 1, a1,k (ω) ≥ 0, 1 ≤ k ≤ K(3)大化するように分離行列 (混合行列の逆行列) を推定することが目的となる。しかし，信号源の数が観測信号のを仮定する。ここで，k 番目の音源の成分を表す sk (ω, t)

数よりも多い劣決定の条件下では，たとえ混合過程が既に関して次の仮定を立てる。

知であったとしても式 (1) を満たす音源成分の解は無数に存在し，独立性の規準ではこの中から一意に解が決められない。このため，劣決定の条件では，音源に関してスパース性などの独立性よりさらに強い仮定が必要とな∗

仮定 1 (振幅成分の時間変化のコヒーレント性). 各音源信号を，振幅成分が周波数間でコヒーレントに時間変化する要素信号を I 個加算したものと見なす。

Underdetermined Blind Source Separation with Multichannel Complex NMF Assuming W-Disjoint Orthogonality of Sources. by Kazuma TAKEDA (University of Tsukuba), Hirokazu KAMEOKA, Hiroshi SAWADA,
Shoko ARAKI (NTT Corporation), Takeshi YAMADA, Shoji MAKINO (University of Tsukuba)

日本音響学会講演論文集- 801 -2011年3月

ある音源 k の i 番目の要素信号の周波数成分比をである。式 (2), (8), (10) より，未知パラメータ

Hki (ω) とし，これに対応する時刻 t におけるゲインおよび位相スペクトルをそれぞれ Uki (t), ϕik (ω, t) とする

θ := {ak (ω), πk , σ 2 (ω), νk2 }1≤k≤K,1≤ω≤Ωと，仮定 1 より，要素信号はixik (ω, t) = Hki (ω)Uki (t)ejϕk (ω,t)(4)

および k̂(ω, t) が既知の下での y(ω, t) の確率密度関数は，¡p(y(ω, t)|k̂(ω, t), θ) = NC A(ω)Zk̂(ω,t) µ(ω, t),¢A(ω)Zk̂(ω,t) ΛZk̂(ω,t) A(ω)H + σ 2 (ω)Iと表すことができる。ただし，X

Hki (ω) = 1(14)πk := {Hki (ω), Uki (t), ϕik (ω, t)}1≤i≤Ik ,1≤ω≤Ω,1≤t≤T(5)(15)ωとする。周波数成分比 Hki (ω) が時不変であることから， で与えられる。また，式 (8) より，s(ω, t) も既知の下でxi (ω, t) の振幅成分が周波数間でコヒーレントに時間変 は y(ω, t) の確率密度関数は，k化するような拘束を自動的に与えていることになっていp(y(ω, t)|s(ω, t), k̂(ω, t), θ)¡¢= NC A(ω)Zk̂(ω,t) s(ω, t), σ 2 (ω)Iる。これを用いて音源信号モデルを

sk (ω, t) =IkX

xik (ω, t) + ϵk (ω, t)(6)(16)で与えられる。

i=1と表すこととする。ここで，ϵk (ω, t) は音源ごとのモデル化誤差である。音源信号を以上のような要素信号の和で表す考え方は，複素 NMF [14] における観測信号のモデル化に基づくため，仮定 1 により立てられる音源モデルに基づく BSS を多チャンネル複素 NMF と呼ぶ。
次に，音源の混合過程に関して次の仮定を立てる。
仮定 2 (音源の時間周波数成分のスパース性). 各時間周EM アルゴリズムによる θ の最尤推定

3以下に，観測信号 y := {y(ω, t)}1≤ω≤Ω,1≤t≤T が与えられたもとで最尤の θ を推定するアルゴリズムを導く。y を不完全データと見なし，y および音源成分 s ={s(ω, t)}1≤ω≤Ω,1≤t≤T とアクティブな音源インデックスk̂ = {k̂(ω, t)}1≤ω≤Ω,1≤t≤T を合わせて完全データとして扱うと，log p(y|θ) を θ に関して最大化する問題は EM波数において単一音源の成分だけがアクティブとなる。 アルゴリズムによる不完全データ問題に帰着できる。このとき Q 関数は，完全データ c := {y, s, k̂} に対する θの対数尤度 log p(c|θ) に対し，観測データ y と θ = θ ′ が以上の仮定は，[2] では W-Disjoint Orthogonality (W与えられた下での s, k̂ に関する条件つき期待値DO) と呼ばれており，音声などのスパースな音源を対象­®
Q(θ, θ ′ ) = log p(c|θ) p(s,k̂|y,θ=θ′ )(17)とした劣決定 BSS に有効であることが知られ，広く用い­­®
®
られている。この仮定により，各時間周波数でアクティ= log p(c|θ) p(s|k̂,y,θ=θ′ ) p(k̂|y,θ=θ′ ) (18)
ブな音源インデックスを k̂(ω, t) とすると，式 (1) は®
®
θ ­­= log P (y|s, k̂, θ) p(s|k̂,y,θ=θ′ ) p(k̂|y,θ=θ′ )y(ω, t) = ak̂(ω,t) (ω)sk̂(ω,t) (ω, t) + n(ω, t)
(7)
­­®
®
+ log P (s|θ) p(s|k̂,y,θ=θ′ ) p(k̂|y,θ=θ′ ) (19)= A(ω)Zk̂(ω,t) s(ω, t) + n(ω, t)(8)と書き直すことができる。Zk̂(ω,t) は n 行 m 列の成分がδn,k̂(ω,t) · δm,k̂(ω,t) (ただし，δn,m は Kronecker デルタ)の K × K 行列である。
2.2

観測信号の統計モデル化と定義される。ただし，⟨f (ξ)⟩p(ξ) は分布 p(ξ) に従う ξξに関する f (ξ) の期待値を，= は ξ に関係する項のみの等号を表す。式 (19) の各項の具体形を以下に示す。まず，⟨⟨log p(y|s, k̂, θ)⟩p(s|k̂,y,θ=θ′ ) ⟩p(k̂|y,θ=θ′ ) は，­­以上の仮定に基づき観測信号の確率密度関数を導く。
音源 k のモデル化誤差を，分散が νk2 の白色 Gauss 性雑音ϵk (ω, t) ∼ NC (0, νk2 )ω(9)−と仮定すると，音源成分 s(ω, t) は¡¢s(ω, t) ∼ NC µ(ω, t), ΛKXmk̂(ω,t)¡ y¢tr Fk̂(ω,t)(ω, t)σ 2 (ω)(20)の形に書ける。ただし，¡¢Tµ(ω, t) := µ1 (ω, t), · · · , µK (ω, t)
IkXi
Hki (ω)Uki (t)ejϕk (ω,t)(11)

yFk̂(ω,t)(ω, t) :=y(ω, t)y(ω, t)H¡¢2
Λ := diag ν12 , · · · , νK

(21)

−A(ω)Zk̂(ω,t) s̄k̂(ω,t) (ω, t)y(ω, t)H(12)−Zk̂(ω,t) A(ω)H y(ω, t)s̄k̂(ω,t) (ω, t)H(13)+Zk̂(ω,t) A(ω)H A(ω)Zk̂(ω,t) Rk̂(ω,t) (ω, t)

i=1

日本音響学会講演論文集X

ω,t k̂(ω,t)=1(10)に従う。ただし，µk (ω, t) :=®
®
log p(y|s, k̂, θ) p(s|k̂,y,θ=θ′ ) p(k̂|y,θ=θ′ )Xθlog σ 2 (ω)= −T M- 802 -2011年3月

であり，M ステップ

3.2′s̄k̂(ω,t) (ω, t) :=µ (ω, t) + Gk̂(ω,t) (ω)y(ω, t)(22)

Rk̂(ω,t) (ω, t) :=s̄k̂(ω,t) (ω, t)s̄k̂(ω,t) (ω, t)H + Λ′
− Gk̂(ω,t) (ω)A′ (ω)Zk̂(ω,t) Λ′′′(23)M ステップでは，Q 関数が最大化，または減少しないことが保証されるように θ を更新すれば良い。まず，A(ω) 以外を固定した下で Q 関数を最大化する A(ω) は，A(ω) =Γ(ω)Σ(ω)−1XX
Γ(ω) :=y(ω, t)
mk̂ s̄k̂ (ω, t)H Zk̂

HGk̂(ω,t) (ω) :=Λ Zk̂(ω,t) A (ω)
(24)
¡ ′
¢
−1A (ω)Zk̂(ω,t) Λ′ Zk̂(ω,t) A′ (ω)H + σ ′2 (ω)ItΣ(ω) :=で あ る 。s̄k̂(ω,t) (ω, t) と Rk̂(ω,t) (ω, t) は そ れ ぞ れ

(30)

k̂

XXt(29)mk̂ Zk̂ Rk̂ (ω, t)Zk̂

(31)

k̂

⟨s(ω, t)⟩p(s|y,k̂,θ=θ′ ) と ⟨s(ω, t)s(ω, t)H ⟩p(s|y,k̂,θ=θ′ ) に対応した量であり，Gk̂(ω,t) (ω) は Wiener フィルタである。 で与えられる。なお，上記に従って A(ω) の更新を行う
場合，必ずしも式 (3) を満たさないため，更新後に適切なまた，mk̂(ω,t) は，正規化を行う必要がある。次に，音源パラメータ πk については，式 (26) より Q 関数の中で πk に関係する項が，mk̂(ω,t) := p(k̂(ω, t)|y(ω, t), θ = θ′ )(25)
X¯¯π 1¯[ŝ(ω, t)]k − µk (ω, t)¯2(32)Q(θ, θ′ ) =k 2
νk ω,tを簡略表記したものであり，(ω, t) において k̂ 番目の音K源がアクティブである確率を表す。なお，ξ ′ は θ = θ ′ のXŝ(ω,t):=mk̂(ω,t) s̄k̂(ω,t) (ω, t)(33)
ときの ξ の値を表すものとする。
k̂(ω,t)=1同様に，⟨⟨log P (s|θ)⟩p(s|k̂,y,θ=θ′ ) ⟩p(k̂|y,θ=θ′ ) は，と等しいため，この関数を大きくするように πk を更新す­­®
®
れば良い。ただし，[·]k はベクトルの k 番目の要素を表
log P (s|θ) p(s|k̂,y,θ=θ′ ) p(k̂|y,θ=θ′ )Xす。上述の式は，複素 NMF の目的関数 [14] と同形であθ= −Tlog |Λ|るため，[14] で提案した反復アルゴリズムと同様の方法
ωで πk を更新することで，Q(θ, θ ′ ) を πk に関して局所最KX X¢¡ −1 s(26) 大化することができる。導出の詳細は [14] に委ねるが，mk̂(ω,t) tr Λ Fk̂(ω,t) (ω, t)−X
ω,t k̂(ω,t)=1
Uki (t)|[ŝ(ω, t)]k |の形に書くことができ，F sk̂(ω,t)

X jHki (ω) ← Hki (ω) X tUki (t)Hk (ω)Ukj (t)(ω, t) は，tXs(ω, t) :=Rk̂(ω,t) (ω, t)Fk̂(ω,t)−µ(ω, t)s̄k̂(ω,t) (ω, t)H−s̄k̂(ω,t) (ω, t)µ(ω, t)+µ(ω, t)µ(ω, t)Hj

Hki (ω)|[ŝ(ω, t)]k |X jUki (t) ← Uki (t) X ωHki (ω)
Hk (ω)Ukj (t)Hω(27)で与えられる。以上より，Q 関数の中で θ に関係する項は式 (20) と式 (26) の和で表される。

(34)(35)

jϕik (ω, t) ← arg([ŝ(ω, t)]k )(36)のように繰り返し更新すれば式 (32) は局所最大化される。σ 2 (ω) および νk2 の更新則も，Q(θ, θ ′ ) の偏微分を 0と置くことで容易に得ることができるが，以下の実装で3.1

E ステップはこれらはいずれも固定値とした。

E ステップは，直前の M ステップで更新された θ をθ に代入し，θ′ に関係する変数に反映させるステップであるため，E ステップでは具体的には，s̄k̂(ω,t) (ω, t) とRk̂(ω,t) (ω, t) と mk̂(ω,t) が式 (22), (23) と式 (25) により′算出される。式 (25) はより具体的には，p(k̂(ω, t)|θ) =1/K とすると，3.3

上述の EM アルゴリズムにより得られた θ を用いて算出される ⟨s(ω, t)⟩p(s(ω,t)|y(ω,t),θ) を分離信号として出力することにする。実は，式 (33) がこれに対応している。
このことは，⟨s(ω, t)⟩p(s(ω,t)|y(ω,t),θ)′= ⟨⟨s(ω, t)⟩p(s(ω,t)|y(ω,t),k̂(ω,t),θ) ⟩p(k̂(ω,t)|y(ω,t),θ)mk̂(ω,t) :=p(k̂(ω, t)|y(ω, t), θ = θ )= PK

分離信号の獲得方法

p(y(ω, t)|k̂(ω, t), θ = θ′ )′′k′ =1 p(y(ω, t)|k̂(ω, t) = k , θ = θ )(28)=と書けるため，式 (15) の分布を用いて算出できる。

日本音響学会講演論文集= ⟨s̄k̂(ω,t) (ω, t)⟩p(k̂(ω,t)|y(ω,t),θ)KX

mk̂(ω,t) s̄k̂(ω,t) (ω, t) = ŝ(ω, t)(37)k̂(ω,t)=1- 803 -

2011年3月

Frequency (kHz)

2.0

1.5

1.0

0.50

図 1 音 源 1 (ピ ア ノ 音 に 対 応) の 基 底 ス ペ ク ト ル H1i (ω)(H11 (ω), · · · , H110 (ω) を横に並べて表示)
表1

分離前

分離後

-6.8
-6.1
2.9

17.7
15.4
14.0

音源 2: ボーカル音源 3: ドラム24

68

10Time (s)図 2 音源 1 (ピアノ音に対応) のアクティベーション U1i (t)(U11 (t), · · · , U110 (t) を縦に並べて表示)

各音源の分離前/後の SIR[dB]

音源 1: ピアノ0は提案手法によって大幅に上昇している。また，図 1, 2からは，倍音構造をもつ基底スペクトルが，時間変化しながら混ざり合うことでピアノの音が表現されているのが読み取れ，仮定 1 の通りに音源がモデル化されているにより確認できる。以上より，分離信号の時間周波数成ことが分かる。現段階では，混合行列 A(ω) が既知とい分 ⟨s(ω, t)⟩p(s(ω,t)|y(ω,t),θ) は E ステップで算出される。

う条件下でなら十分な音源分離性能が得られることが確

4動作実験
本章では提案法の音源分離性能の検証結果について述べる。提案法は，周波数領域 BSS における劣決定性とパーミュテーション問題に起因する解の非一意性を同時解消することを目指したものであるが，本実験ではトータルな性能評価に向けた予備検討として，劣決定性のみに対する提案法の効果を簡易に調べる目的で混合行列A(ω) を既知とした。
対象とする信号は，ピアノ，ボーカル，ドラムの 3 つの音源と，マイク間の距離を 1m，残響時間を 100ms として測定したインパルス応答を畳込み混合によって人工的に混合した 10 秒の楽音信号とした。また，観測信号の時間周波数成分は短時間 Fourier 変換 (標本化周波数 16kHz，フレーム長 128ms，オーバーラップ 64ms，Hanning 窓)により計算した。提案法の条件については，EM アルゴリズムの繰り返し数を 50 とし，各パラメータ θ の初期値を混合行列 A(ω) のみ正解を与え，その他の値は全て乱数で設定した。基底スペクトル数は I1 = · · · = IK = 10とした。また，分離信号は式 (36) により定義される推定値とした。獲得した分離信号は客観評価によって測定した。なお，客観評価値として，他音声の消し残りによる歪みを表す歪み具合を表す SIR(Source to InterferenceRatio) [15] を用いた。単位は dB であり，この歪み尺度

まとめ5

本稿では，劣決定ブラインド音源分離における，劣決定性およびパーミュテーション問題に起因する解の非一意性を同時に解消する手法として，音源の W-DO 性を仮定した多チャンネル複素 NMF を提案した。提案法による音源分離実験を通して，混合行列が既知であるなら高い精度で音源分離できることを確認した。今後の課題として，EM アルゴリズムにおけるパラメータの適切な初期値の選定方法が挙げられる。特に混合行列の初期値選定方法は最重要課題になると考えられる。また，提案法の優位性を示すために他の従来法との性能比較実験を行う予定である。

参考文献[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] Kurita et al., Proc. ICASSP’00, pp. 3140–3143, 2000.
[8] Murata et al., Neurocomputing, 41(1–4), pp. 1–24,
2001.の数値が高いほど性能が良い。
客観評価の値を表 1 に示す。ここでは，分離前の信号と各音源信号との SIR，提案手法による分離後の信号と各音源信号との SIR の値を比較した。また，図 1 に提案法で求まったピアノ音の基底スペクトル Hki (ω)，図 2 に各基底スペクトルのゲイン Uki (t) の例を示す。

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

[10] 小野, 音講論 (秋)’10, 2-10-7, pp. 581–582, 2010.
[11] 北野他, 音講論 (秋)’09, 2-4-15, pp. 645–646, 2009.
[12] Sawada et al., submitted to ICASSP’11 in 2010.
[13] Ozerov & Févotte, IEEE Trans. ASLP, 18(3), pp. 550–563, 2010.

提案手法によって混合信号を各々の音源信号に分離できたかを考察する。表 1 を見ると分かる様に，SIR の値

日本音響学会講演論文集

認できた。

[14] 亀岡他, 音講論 (秋)’08, 2-8-13, pp. 657–660, 2008.
[15] Vincent et al., Proc. ICA’07, pp. 552–559, 2007.

- 804 -2011年3月