多チャンネル Factorial hidden Markov model による
劣決定ブラインド音源分離と音響イベント検出の統合的アプローチ ∗◎樋口卓哉, 竹田裕史，中村友彦 (東大院情報理工)，亀岡弘和 (東大院情報理工，NTT CS 研)

1

はじめに

ブラインド音源分離の問題とは，音源信号や音源からマイクまでの伝達特性が未知の場合に，複数の音
源信号が混合された観測信号から音源信号を推定する問題である．一般的に，ブラインド音源分離の問題を解くためには，音源信号に対して立てたなんらかの仮定を基に最適化基準を立て，最適化問題を解く必要がある．例えば，観測信号数が音源数よりも多い優決定問題では，音源信号間の独立性を仮定して分離する独立成分分析 (Independent Component Analysis;
ICA) が有用であることが知られており，音源信号間の独立性を最大化するように分離フィルタを推定することが目的となる [1]。しかし，ICA では観測信号
数が音源数よりも少ない劣決定問題を扱うことはできず，この場合は独立性よりもさらに強い仮定が必要である。
単チャンネルの観測信号に対するブラインド音源分離の有効なアプローチとして，非負値行列因子分
解 (Non-negative Matrix Factorization; NMF) が知られている [2, 3]．この手法では観測信号のパワース
ペクトログラムを，2 つの非負値行列の積に分解する．
分解した各行列は，いくつかの基底パワースペクトルによって構成される基底行列と，それらの基底パワー
スペクトルの時変な音量を表すアクティベーション
によって構成されるアクティベーション行列となる．
ここで重要なのは，分解された各パワー基底スペク
トルが，観測信号の中で主となる要素，すなわち各音
源信号を表していると考えられることである．また，音源信号の空間的な情報も利用して音源分離を行う
ために，NMF を多チャンネルの音響信号へと拡張するアプローチがいくつか取られてきた [4, 5]．
しかし以上のアプローチでは，各音源信号のパワー
スペクトルが，1 つの基底パワースペクトルによって表現できることを仮定していた．しかし実際の音源
信号のパワースペクトルは時変であることが多く，1つの基底パワースペクトルで表現するのは不十分な場合が多い．例えば，無音状態，音の立ち上がり，定常状態などの音源信号の状態に応じて，音源は異なるパワースペクトルを持つ場合があり，それらを 1 つの基底スペクトルによって表現するのは不十分であると考えられる．また音量に着目してみても，無音状
態，有音状態では大きな値のとりやすさが異なると考えられる．このように音源は，背後にある音源の状態に応じて異なる信号を出力すると考えられる．
これらの観点から，無音状態，有音状態などの音源の状態を推定する問題 (音響イベント検出の問題) と，∗

ブラインド音源分離の問題は相互依存の関係にあると考えられ，本来同時に解かれるべきであるといえる．そこで本稿では，音源信号の時変な性質を，音源の状態を隠れ変数とする隠れマルコフモデル (hidden
Markov model; HMM) によってモデル化し，パラメー
タ推定を通してブラインド音源分離と音響イベント検出を統合的に行う手法を提案する．本稿ではこのモデルを多チャンネル factorial hidden Markov modelと呼ぶ．

2

多チャンネル factorial HMM

2.1

混合モデル

まず，観測信号の生成プロセスについて述べる．I個の音源信号が M 個のマイクロフォンで観測される場合を考える．ym (ωk , tl ) ∈ C を m 番目のマイクで観測された観測信号の周波数 ωk ，時刻 tl における時
間周波数成分，si (ωk , tl ) ∈ C を i 番目の音源信号の周波数 ωk ，時刻 tl における時間周波数成分とし，
1 ≤ k ≤ K と 1 ≤ l ≤ L をそれぞれ時間周波数領域
における周波数と時間のインデックスとする．ここで，室内インパルス応答長が時間周波数展開における時間窓長よりも十分に短い場合を仮定すると，瞬時混合近似を用いて観測信号は以下のように時間周波数領域において記述できる．

y(ωk , tl ) =I∑ai (ωk )si (ωk , tl ),

(1)

i=1

ただし y(ωk , tl ) = (y1 (ωk , tl ), . . . , yM (ωk , tl ))T ∈CM である. ai (ωk ) は i 番目の音源信号に対する周波数 ωk における伝達周波数特性を表す．表記の都合上，以下では ωk ，tl を k ，l の添え字でそれぞれ表す．

2.2 観測信号の生成プロセス
次に，式 (1) に基づいて観測信号の生成プロセスを確率的に記述する．まず，音源信号 si,k,l が平均 0，分2
散 σi,k,lの複素正規分布に従うと仮定する．
2
si,k,l |σi,k,l ∼ NC (si,k,l ; 0, σi,k,l
).

(2)

2
ここで σi,k,lは周波数 k ，時刻 l における i 番目の音源のパワースペクトル密度を表す．式 (1) と式 (2) から,
a1:I,k と σ1:I,k,l が既知の条件下で観測信号 y k,l は同
じく複素正規分布に従う．
∑
2
y k,l |a1:I,k , σ1:I,k,l ∼ NC (y k,l ; 0,C i,k σi,k,l), (3)
i

Unified approach for underdetermined blind source separation and source activity detection based on
multichannel factorial hidden Markov model. by HIGUCHI, Takuya, TAKEDA, Hirofumi, NAKAMURA,
Tomohiko (The University of Tokyo), KAMEOKA Hirokazu (The University of Tokyo, NTT)

ただし C i,k = ai,k aH
i,k は i 番目の音源に対する周波数 k における空間相関行列と呼ばれる行列であり，また，NC (x; µ, Σ) ∝ exp(−(x − µ)H Σ−1 (x − µ)) で
ある．

2.3

多チャンネル NMF[4, 5] における生成モデル

通常の NMF のモデルにおいては ([6] 参照)，音源
信号のパワースペクトルはスケールを除いて時不変であることが仮定されていた． 上記のモデルにこの仮定を組み込むと，
2
σi,k,l= wi,k hi,l ,

si,k,l |wi,k,1:D , hi,l , zi,l(4)

2となる．ここで σi,k,lは時不変の基底パワースペクト
ル wi,k と時変な音量を表す hi,l の積の形で表現されている．これにより si,k,l の生成モデルは以下のように書き直せる．

si,k,l |wi,k , hi,l ∼ NC (si,k,l ; 0, wi,k hi,l ).

対応するときは一様分布に近くなるようにガンマ分
布のハイパーパラメータを設定すればよい．具体的には，無音状態に対応する状態では α と β をそれぞれ 1，10−2 などと設定し，有音状態のときはそれぞれ 1，1020 などと設定すればよい．このような設定により，音響イベント検出と音源分離を協調的に解くことが可能となる．状態 d である i 番目の音源の基底パ
ワースペクトルを wi,k,d と表すとすると，時刻 l における i 番目の音源信号のパワースペクトルは zi,l に依存し，si,k,l の生成モデルは以下のように書ける．

(5)∼ NC (si,k,l ; 0, wi,k,zi,l hi,l ).

y k,l の生成モデルは各音源信号による HMM の足し合わせとなるので，最終的な生成モデルは factorial
HMM とみなせる．観測信号の最終的な生成モデルは
a1:I,k ，w1:I,k,1:D ， h1:I,l ，z1:I,l が既知の条件下で，式 (6)，式 (7) と合わせて以下のように書ける．

式 (5) による Y = {y k,l }k,l の生成モデルは，NMFの多チャンネル観測信号への自然な拡張とみなすことができ，このモデルに基づく BSS のアプローチは多チャンネル NMF[4, 5] と呼ばれている．

2.4

HMM を用いた音源信号の生成モデル

zi,l |zi,l−1 ∼ Categorical(zi,l ; ρzi,l−1 ).

(9)iパラメータ推定アルゴリズム

3.1 目的関数
ここでは，上記の生成モデルに対する，補助関数法に基づくパラメータ推定アルゴリズムについて述べる．
モデルにおける推定したい変数は W = w1:I,1:K,1:D ,H = h1:I,1:L , C = C 1:I,1:K ，Z = z1:I,1:L である．
上記の変数の集合を Θ で表す．以下では ρ は実験的に定められた定数とする．パラメータ推定のためには事後確率
p(Θ|Y ) =p(Y , Θ)
,
p(Y )

(10)

を計算することが必要となる．ここで Y = y 1:K,1:L
は多チャンネル観測信号の時間周波数成分の集合である．ここで 2 章で定義された条件つき確率を用いて，同時確率 p(Y , Θ) は以下のように記述できる．

(6)

こ こ で Categorical(x; y) = yx で あ り，ρd =(ρd,1 , . . . , ρd,D ) は状態 d から各状態 1, . . . , D への遷移確率を表し， ρ = (ρd,d′ )D×D は遷移行列である．
ここで，hi,l が，zi,l によって異なるハイパーパラメー
タを持つガンマ分布に従うと仮定すると，

hi,l |zi,l ∼ Gamma(hi,l ; αzi,l βzi,l ),y k,l |a1:I,k , w1:I,k,1:D , h1:I,l , z1:I,l∑∼ NC (y k,l ; 0,C i,k wi,k,zi,l hi,l ).

3

上記のように，多チャンネル NMF はそれぞれの音
源信号のパワースペクトルはスケールを除いて時不変であることが仮定されていた．しかし，多くの音源のスペクトルはその “状態” に応じて異なると考えられる．ここでいう “状態” とは，例えば，ピアノの音であれば “アタック”，“ディケイ”，“サステイン”，“
リリース” などを指す．また無音状態や有音状態についても音源の状態とみなすことができ，音源信号の音
量も，音源の状態に依存すると捉えることができる．
このような音源信号の時変な性質を表現するため，ここで各音源信号のパワースペクトルの時系列と音量を，HMM を用いて表現するモデルを提案する．
時刻 l における i 番目の音源の状態を表す隠れ変数 zi,l ∈ {1, . . . , D} を導入し，状態の時系列
zi,1 , . . . , zi,L がマルコフ連鎖に従うと仮定する．

(8)

(7)

となる．ここで α1:D と β1:D はそれぞれガンマ分
布の形状パラメータとスケールパラメータであり，α−1 −x/βe
Gamma(x; α, β) = x Γ(α)β
である．zi,l が無音状
α態に対応するときは hi,l は小さな値をとってほしいので，小さな値をとる確率が高くなるようにガンマ分
布のハイパーパラメータを設定し，zi,l が有音状態にp(Y , Θ) ∝p(Y |Θ)p(H|Z)p(Z).

(11)

目的関数は L(Θ) = log p(Θ|Y ) と定義でき，我々の目的は以下の式を満たす Θ̂ を求めることである．

Θ̂ = argmax log p(Θ|Y ).

(12)

Θ

式 (10)，(11)，(12) から，この最適化問題は以下のよ
うに書き直せる．

(
Θ̂ = argmax log p(Y |Θ)Θ)+ log p(H|Z) + log p(Z) .

(13)

また，2 章で記述した生成モデルから，log p(Y |Θ) は
以下のように書ける．

log p(Y |Θ)
∑ M
−1
1
1
(− log 2π − log |X̂ k,l | − y k,l H X̂ k,l y k,l ).
=
2
2
2
k,l

(14)
ここで X̂ k,l =3.2∑

i C i,k wi,k,zi,l hi,l である．

補助関数法に基づく最適化アルゴリズム

目的関数 L(Θ) を Θ について解析的に最大化することは難しいが，補助関数法に基づく反復計算に
よって局所最適となる Θ を求めることができる [5]．
今回の最適化問題に適用するため，まず L(Θ) =maxΛ L+ (Θ, Λ) を満たす補助関数 L+ (Θ, Λ) を設計する．ここで Λ は補助変数である．紙面の都合上詳細は省略するが，Θ ← argmaxΘ L+ (Θ, Λ) と Λ ←argmaxΛ L+ (Θ, Λ) を繰り返すことにより，L(Θ) を局所最大化できる．このとき，Θ と Λ について解析的に最適化可能な L+ (Θ, Λ) を設計することが重要である．今回の場合は，以下のように補助関数 L+ (Θ, Λ)を設計できる．

それぞれ以下のようになる．
(
−1
H
∂L+
1 ∑ tr(y k,l y k,l Ri,k,l C i,k Ri,k,l )=2∂wi,k,zi,l2
wi,k,z
hi,l
i,l
l
)
−1
− tr(U k,l C i,k )hi,l ,
(18)

(
−1
H
∂L+
1 ∑ tr(y k,l y k,l Ri,k,l C i,k Ri,k,l )
=
∂hi,l
2
wi,k,zi,l h2i,lk)C)w− tr(U −1
i,k,zi,lk,l i,k+ (αzi,l − 1)/hi,l − 1/βzi,l . (19)
これらを 0 と置くことで，以下の更新式が導ける．
vu∑
tr(y k,l y k,l H Ri,k,l C −1
ui,k Ri,k,l )
u l
hi,l
t
wi,k,zi,l ←, (20)∑−1l tr(U k,l C i,k )hi,l√(αzi,l − 1) + (αzi,l − 1)2 + µi,l νi,l
hi,l ←
,
νi,l
(21)
ただし

µi,l =∑ tr(y k,l y k,l H Ri,k,l C −1
i,k Ri,k,l )
wi,k,zi,l

k

,

(22)

∑
νi,l =
tr(U −1(23)
k,l C i,k )wi,k,zi,l + 2/βzi,l ,≥L+ (Θ, Λ)
k
( (
−1
H
1 ∑ ∑ tr(y k,l y k,l Ri,k,l C i,k Ri,k,l )である．C についての L+ の変微分は以下のように=−
2whi,k,zi,l i,l
i
なる．
k,l
)
)
−1
H
∑( C −1∂L+i,k Ri,k,l y k,l y k,l Ri,k,l C i,k+ tr(U −1
k,l C i,k )wi,k,zi,l hi,l + log |U k,l | − M=∂C i,kwi,k,zi,l hi,l
)
l
∑(
)
+
(αzi,l − 1) log hi,l − hi,l /βzi,l − αzi,l log βzi,l−1− U k,l wi,k,zi,l hi,l .
(24)
i,l
L(Θ)

+ log p(Z).

(15)∑

ここで Ri,k,l と U k,l は i Ri,k,l = I を満たすエルミート正定値行列であり， R と U の集合を Λ で表す. tr(·) は行列のトレースを表す． 式 (15) の等号成
立条件は
−1

これを 0 と置くと，以下の Riccati 方程式が得られる．

C i,k Ai,k C i,k = B i,k ,ただし

Ai,k =

∑

(25)

−1

wi,k,zi,l hi,l X̂ k,l ,lRi,k,l = C i,k wi,k,zi,l hi,l X̂ k,l ,(16)

U k,l = X̂ k,l ,

(17)

となる．
以上から，L は次の 2 つのステップを繰り返すこと
によって局所最大化できる．

1. R と U について L+ を最大化.
2. W ，H ，C ，Z について L+ を最大化.
ステップ 1 における R と U の更新では，式 (16) と式 (17) を使えばよい．ステップ 2 では，L+ をそれぞれの変数で変微分して 0 となるものを求めることで，更新則が導ける．W と H に関する L+ の変微分は∑B i,k = C i,k (
wi,k,zi,l hi,ll
−1
−1X̂ k,l y k,l y k,l H X̂ k,l )C i,k ,(26)である．以下の方法でこの Riccati 方程式を解くこと
で C の更新則が得られる [5]．まず，以下の 2M ×2Mの行列に対して固有値分解を行う．
[
]
0
−Ai,k
.
(27)
−B i,k0
ここで e1,i,k . . . eM,i,k を負の固有値に対応する固有ベクトルだとし，2M 次元の固有ベクトルを m =1 . . . M において以下のように分解する．
[]
f m,i,kem,i,k =(28)
g m,i,k .

ルゴリズムは 100 回反復した．比較対象には [5] の手法を用いた．分離音 ŷ i,k,l はウィーナーフィルタ−1ŷ i,k,l = wi,k,zi,l hi,l C i,k X̂ k,l y k,l ,Fig. 1 ベルの音のスペクトログラム (上), 提案法によって得られた分離音 (中), 音源の状態推定結果 (下).
黒がその時刻で推定された状態を表す.
Table 1 提案法と従来法によって得られた SDR の平均値と標準偏差.
SDR(±SD) [dB]

ベル

笛

ホッチキス

提案法
従来法

19.92(±11.50)
13.33(±8.22)

31.21(±6.69)
23.37(±4.57)

16.81(±10.09)
8.28(±5.98)

ここで f m,i,k と g m,i,k は M 次元のベクトルである．
C i,k の更新則は以下のように得られる．

C i,k ← Gi,k F −1
i,k .

(29)

た だ し F i,k= [f 1,i,k , . . . , f M,i,k ]，Gi,k=[g 1,i,k , . . . , g M,i,k ] である．
L+ を Z の関数と見ると，各音源に対応する HMMの対数事後確率の和となっている．従って最適な状態の時系列 zi,1 , . . . , zi,L をビタビアルゴリズムによって効率的に，各 i ごとに個別に求めることができる．
以上の更新則から W ，H ，C ，Z を反復計算により求めることは，ブラインド音源分離の問題と音響
イベント検出の問題を画一的最適化規準に基づき協調的に解いていることに相当している．

4

評価実験

提案法の音源分離性能と音響イベント検出性能の評
価のために，実験を行った．RWCP データベース非音
声ドライソース [7] の中のホッチキス，ベル，笛の音に対して同じく RWCP データベース [7] のインパルス応
答 (残響時間 0 ms，マイク間距離 5.85 cm，マイクの数
M = 2) を畳み込み，人工的に多チャンネルの混合信号を作成した．サンプリング周波数は 32 kHz とした．
フレーム長 16 ms，フレームシフト長 8 ms で STFTを行い，時間周波数展開を行った．HMM の状態数 D
は 3 とした．α1 と β1 を 1，10−2 とそれぞれ設定し，α2:3 と β2:3 を 1 と 1020 と設定することで，d = 1 を無音状態とみなした．遷移確率 ρ は ρ1 = (0.9, 0.1, 0)，
ρ2 = (0, 0.5, 0.5)，ρ3 = (0.5, 0, 0.5) とそれぞれ設定√した．C の初期値については，対角成分を 1/ M ，それ以外の成分を 0 とした．W と H については，まずランダムな初期値から単チャンネルの IS-ダイバージェンス規準の NMF を行い，結果として得られた
W と H を初期値とした．以上をランダムな W とH の初期値を変えて 10 回行った．パラメータ推論ア(30)に よって 得 た ．客 観 評 価 基 準 と し て ，signal-todistortion ratio (SDR)[8] を用いた．高い SDR は高い音源分離性能を表す．分離前の SDR はベル, 笛，
ホッチキスの音でそれぞれ-8.36, 6.64，-12.78 [dB] であった．
Table 1 に 10 回の施行で得られた SDR の平均値と標準偏差を示す．提案法によって得られた分離音の SDR の平均は従来法の SDR の平均を 7.65 [dB] 上
回った．Fig. 1 にベルの音のスペクトログラム (上)，提案法によって得られた分離音のスペクトログラム(中)，状態推定結果 (下) を示す．黒がその時刻で推定された状態を表す．音響イベント検出がおおむね正しく行われていることがわかる．

5

おわりに

本稿ではブラインド音源分離と音響イベント検出を統合的に行う手法を提案した．音源信号の状態を隠れ変数とする HMM を用いて各音源信号の生成モデルを記述し，設計した全体の生成モデルに基づくパラメータ推定を通して，ブラインド音源分離と音響イベント検出を協調的に行うことができるのが本
手法のポイントである．

6

謝辞

本研究は JSPS 科研費 26730100 の助成を受けたも
のです。

参考文献
[1] A. Hyvärinen et al., John Wiley & Sons, 2001.
[2] D. D. Lee, and H. S. Seung, Nature, vol. 401,
pp.788–791, 1999.
[3] P. Smaragdis, and J. C. Brown, in Proc. WASPAA 2003, Oct. 2003, pp. 177–180.
[4] A. Ozerov, and C. Févotte, IEEE Trans. Audio,
Speech and Language Processing, vol. 18, no. 3,
pp. 550-563, Mar.2010.
[5] H. Sawada et al., in Proc. ICASSP, pp. 261-264,
2012.
[6] C. Févotte, N. Bertin, and J.L. Durrieu, Neuralcomputation, vol. 21, no. 3, pp. 793–830, 2009.
[7] S. Nakamura et al., in Proc. LREC ’00, pp. 965–
968, 2000.
[8] E. Vincentet al., IEEE Trans. ASLP, pp. 1462–
1469, 2006.