多チャンネル 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つの基底パワースペクトルで表現するのは不十分な場合が多い．また音量に着目してみても，無音状態，
有音状態では大きな値のとりやすさが異なると考えられる．そこで我々は以前，音源信号の時変な性質
を，音源の状態を隠れ変数とする隠れマルコフモデ
ル (hidden Markov model; HMM) によってモデル化
し，パラメータ推定を通してブラインド音源分離と音響イベント検出を統合的に行う手法を提案した [6]．
しかし従来のモデルでは，室内インパルス応答長が時間周波数展開における時間窓長よりも十分に短い場合を仮定しており，残響がある場合を考慮していない．従って従来の手法をそのまま利用して，残響を考慮した精度良い音源分離を行うためには，従来法とは別に，残響除去を行う必要がある．すなわち，音源分離性能は残響除去の精度に依存することになる．
一方で，実環境においてどのように残響がかかるか
は各音源信号によって異なるため，残響除去の精度は
音源分離性能に依存する．
このような観点から，残響除去の問題と音源分離の問題は相互依存の関係にあると考えられ，本来同時に解かれるべきであるといえる．そこで本稿では，
従来のモデル [6] を拡張し，残響下での観測信号を時
間周波数領域における畳み込み混合の形で表現し，パラメータ推定を通して音源分離・音響イベント検出・残響除去を統合的に行う手法を提案する．本稿ではこのモデルを多チャンネル factorial hidden Markov
model と呼ぶ．

∗

2

多チャンネル factorial HMM

2.1 混合モデル
まず，観測信号の生成プロセスについて述べる．I個の音源信号が M 個のマイクロフォンで観測される場合を考える．ym (t) ∈ R を m 番目のマイクで観測される観測信号，si (t) ∈ R を i 番目の音源信号だとすると，観測信号は時間領域において以下のように記述できる．
y(t) =I ∑∞∑ai (τ )si (t − τ ).

(1)

i=1 τ =0

ここで y(t) = (y1 (t), . . . , yM (t))T ∈ RM ， ai (t) =(ai,1 (t), . . . , ai,M (t))T ∈ RM である． ai,m (t) は i 番
目の音源と m 番目のマイクの間の室内インパルス応答を表す．ここで，室内インパルス応答長が時間周波数展開における時間窓長よりも十分に短い場合，瞬時混合近似を用いて観測信号は以下のように時間周波数領域において記述できる．

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

(2)

i=1

ただし， y(ωk , tl ) = (y1 (ωk , tl ), . . . , yM (ωk , tl ))T ∈CM ， ai (ωk ) = (ai,1 (ωk ), . . . , ai,M (ωk ))T ∈ CM である．ここで ym (ωk , tl ) ∈ C は m 番目のマイクで観測された観測信号の周波数 ωk ，時刻 tl における時間
周波数成分であり，si (ωk , tl ) ∈ C は i 番目の音源信号の周波数 ωk ，時刻 tl における時間周波数成分である．
ai (ωk ) は i 番目の音源信号に対する周波数 ωk における伝達周波数特性を表す．1 ≤ k ≤ K と 1 ≤ l ≤ Lはそれぞれ時間周波数領域における周波数と時間のインデックスである．しかし残響がある場合には一般に，室内インパルス応答長は時間窓長に対して十分に短いとはいえず，式 (2) の瞬時混合近似は成り立たない．そこで，時間周波数領域における畳み込み混合の形で観測信号を近似する．

y(ωk , tl ) ≈

I ∑
N
∑

ai (ωk , tn )si (ωk , tl − tn ).

(3)

i=1 n=0

ここで ai (ωk , n) は i 番目の音源信号に対する周波数
ωk における伝達周波数特性の時刻 tn の成分であり，
0 ≤ n ≤ N は伝達周波数特性の時間周波数領域における時間インデックスである．ここで ai (ωk , tn ) は i
番目の音源信号が時間周波数領域において tn だけ先の時刻にどれだけ影響を与えるか，すなわち時間周
波数領域においてどれだけ残響がかかるかを表している．表記の都合上，以下では ωk ，tl ，tn を k ，l，nの添え字でそれぞれ表す．

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

Unified approach for source separation, dereverberation and source activity detection based on multichannel factorial hidden Markov model. by HIGUCHI, Takuya (The University of Tokyo), KAMEOKA
Hirokazu (The University of Tokyo, NTT)(4)

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

y k,l |a1:I,k,0:N , σ1:I,k,l−N :l∼ NC (y k,l ; 0,∑2
C i,k,n σi,k,l−n), (5)

i,n

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

2.3 多チャンネル NMF[4, 5] における生成モデル
通常の NMF のモデルにおいては ([7] 参照)，音源
信号のパワースペクトルはスケールを除いて時不変であることが仮定されていた． 上記のモデルにこの仮定を組み込むと，
2
σi,k,l= wi,k hi,l ,(6)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 ).

(8)

こ こ で Categorical(x; y) = yx で あ り，ρd =(ρd,1 , . . . , ρd,D ) は状態 d から各状態 1, . . . , D への遷移確率を表し， ρ = (ρd,d′ )D×D は遷移行列である．
ここで，hi,l が，zi,l によって異なるハイパーパラメー
タを持つガンマ分布に従うと仮定すると，(9)となる．ここで α1:D と β1:D はそれぞれガンマ分
布の形状パラメータとスケールパラメータであり，α−1 −x/βe
Gamma(x; α, β) = x Γ(α)β
である．zi,l が無音状
α態に対応するときは hi,l は小さな値をとってほしい
ので，小さな値をとる確率が高くなるようにガンマ分布のハイパーパラメータを設定し，zi,l が有音状態に対応するときは一様分布に近くなるようにガンマ分布のハイパーパラメータを設定すればよい．状態 d
である i 番目の音源の基底パワースペクトルを wi,k,dと表すとすると，時刻 l における i 番目の音源信号の(10)

y k,l の生成モデルは各音源信号による HMM の足し合わせとなるので，最終的な生成モデルは factorial
HMM とみなせる．観測信号の最終的な生成モデル
は a1:I,k,0:N ，w1:I,k,1:D ， h1:I,l−N :l ，z1:I,l−N :l が既知の条件下で，式 (8)，式 (9) と合わせて以下のように書ける．
y k,l |a1:I,k,0:N , w1:I,k,1:D , h1:I,l−N :l , z1:I,l−N :l∑∼ NC (y k,l ; 0,C i,k,n wi,k,zi,l−n hi,l−n ). (11)i,n3パラメータ推定アルゴリズム

3.1 目的関数
ここでは，上記の生成モデルに対する，補助関数法に基づくパラメータ推定アルゴリズムについて述べる．
モデルにおける推定したい変数は W = w1:I,1:K,1:D ,H = h1:I,1:L , C = C 1:I,1:K,0:N ，Z = z1:I,1:L である．上記の変数の集合を Θ で表す．以下では ρ は実
験的に定められた定数とする．パラメータ推定のためには事後確率p(Θ|Y ) =2.4 HMM を用いた音源信号の生成モデル
上記のように，多チャンネル NMF はそれぞれの音
源信号のパワースペクトルはスケールを除いて時不変であることが仮定されていた．しかし，多くの音源のパワースペクトルはその “状態” に応じて異なると考えられる．そこで音源信号の時変な性質を表現するため，ここで各音源信号のパワースペクトルの時
系列と音量を，HMM を用いて表現する．
時刻 l における i 番目の音源の状態を表す隠れ変数 zi,l ∈ {1, . . . , D} を導入し，状態の時系列
zi,1 , . . . , zi,L がマルコフ連鎖に従うと仮定する．

hi,l |zi,l ∼ Gamma(hi,l ; αzi,l βzi,l ),si,k,l |wi,k,1:D , hi,l , zi,l∼ NC (si,k,l ; 0, wi,k,zi,l hi,l ).

(7)

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

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

パワースペクトルは zi,l に依存し，si,k,l の生成モデ
ルは以下のように書ける．

p(Y , Θ)
,
p(Y )

(12)

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

p(Y , Θ) ∝p(Y |Θ)p(H|Z)p(Z).

(13)

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

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

(14)

Θ

式 (12)，(13)，(14) から，この最適化問題は以下のよ
うに書き直せる．
(
Θ̂ = argmax log p(Y |Θ)Θ)+ log p(H|Z) + log p(Z) .
(15)
また，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

(16)
ここで X̂ k,l =∑i,n C i,k,n wi,k,zi,l−n hi,l−n である．

3.2 補助関数法に基づく最適化アルゴリズム
目的関数 L(Θ) を Θ について解析的に最大化することは難しいが，補助関数法に基づく反復計算に
よって局所最適となる Θ を求めることができる [5]．
今回の最適化問題に適用するため，まず L(Θ) =maxΛ L+ (Θ, Λ) を満たす補助関数 L+ (Θ, Λ) を設計する．ここで Λ は補助変数である．紙面の都合上詳細は省略するが，Θ ← argmaxΘ L+ (Θ, Λ) と Λ ←argmaxΛ L+ (Θ, Λ) を繰り返すことにより，L(Θ) を局所最大化できる．このとき，Θ と Λ について解析的に最適化可能な L+ (Θ, Λ) を設計することが重要である．今回の場合は，以下のように補助関数 L+ (Θ, Λ)を設計できる．
L(Θ)≥L+ (Θ, Λ)( (Fig. 1 作成した混合音のスペクトログラム (残響時
−1H
1 ∑ ∑ tr(y k,l y k,l Ri,k,l,n C i,k,n Ri,k,l,n )間 600 ms).
=−
2
w
h
i,k,z
i,l−n
i,l−n
i,n
k,l
)
) これらを 0 と置くことで，以下の更新則が導ける．
wi,k,zi,l+ tr(U −1
k,l C i,k,n )wi,k,zi,l−n hi,l−n + log |U k,l | − Mvu∑)tr(y k,l+n y k,l+n H Ri,k,l+n,n C −1u∑(i,k,n Ri,k,l+n,n )
u l,n
hi,l
+
(αzi,l − 1) log hi,l − hi,l /βzi,l − αzi,l log βzi,l
t
,
←
∑
−1
i,l
l,n tr(U k,l+n C i,k,n )hi,l+ log p(Z).
(17)
(22)
∑
hi,l
ここで Ri,k,l,n と U k,l は i,n Ri,k,l,n = I を満たす√
(αzi,l − 1) + (αzi,l − 1)2 + µi,l νi,l
エルミート正定値行列であり， R と U の集合を Λ←,(23)で表す. tr(·) は行列のトレースを表す． 式 (17) の等νi,l号成立条件は
ただし−1

Ri,k,l,n = C i,k,n wi,k,zi,l−n hi,l−n X̂ k,l ,U k,l = X̂ k,l ,

(18)
(19)

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

1. R と U について L+ を最大化.

k,n

(24)
νi,l
=

∑

tr(U −1k,l+n C i,k,n )wi,k,zi,l + 2/βzi,l ,

(25)

k,n

2. W ，H ，C ，Z について L+ を最大化.
ステップ 1 における R と U の更新では，式 (18) と式 (19) を使えばよい．ステップ 2 では，L+ をそれぞれの変数で変微分して 0 となるものを求めることで，更新則が導ける．W と H に関する L+ の変微分はそれぞれ以下のようになる．

∂L+∂wi,k,zi,l(−1
H1 ∑ tr(y k,l+n y k,l+n Ri,k,l+n,n C i,k,n Ri,k,l+n,n )=
22
wi,k,z
hi,l
i,l
l,n
)
− tr(U −1(20)k,l+n C i,k,n )hi,l ,である．C についての L+ の変微分は以下のようになる．
−1
H∑( C −1∂L+i,k,n Ri,k,l,n y k,l y k,l Ri,k,l,n C i,k,n=∂C i,k,nwi,k,zi,l−n hi,l−nl)− U −1
w
h
.
(26)
i,k,z
i,l−n
i,l−n
k,l
これを 0 と置くと，以下の Riccati 方程式が得られる．

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

Ai,k,n =∑B i,k,n = C i,k,n ((−1
H1 ∑ tr(y k,l+n y k,l+n Ri,k,l+n,n C i,k,n Ri,k,l+n,n )=2
wi,k,zi,l h2i,lk,n)−1− tr(U k,l+n C i,k,n )wi,k,zi,l
(21)

(27)

−1

wi,k,zi,l−n hi,l−n X̂ k,l ,

l

∂L+
∂hi,l

+ (αzi,l − 1)/hi,l − 1/βzi,l .

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

l
−1

−1

X̂ k,l y k,l y k,l H X̂ k,l )C i,k,n ,(28)である．紙面の都合上詳細は省略するが，この Riccati
方程式を解くことで C の更新則が得られる．
L+ を Z の関数と見ると，各音源に対応する HMMの対数事後確率の和となっている．従って最適な状態の時系列 zi,1 , . . . , zi,L をビタビアルゴリズムによって効率的に，各 i ごとに個別に求めることができる．
以上の更新則から W ，H ，C ，Z を反復計算により求めることは，音源分離・音響イベント・残響除去の問題を画一的最適化規準に基づき協調的に解いていることに相当している．

Table 1 提案法と従来法によって得られた SIR.
SIR [dB]
提案法
従来法

ベル
43.35
32.37

携帯電話
32.05
30.15

ホッチキス
7.19
-8.77

象には [6] の手法を用いた．この手法は提案法において N = 0 と置いたものに相当する．分離音 ŷ i,k,l は
ウィーナーフィルタ−1ŷ i,k,l = wi,k,zi,l hi,l C i,k,0 X̂ k,l y k,l .

(29)に よって 得 た ．客 観 評 価 基 準 と し て ，signal-tointerference ratio (SIR)[9] を用いた．高い SIR は高い音源分離性能を表す．分離前の SIR はベル, 携帯
電話の着信音，ホッチキスの音でそれぞれ -16.61,
16.58，-39.17 [dB] であった．
Table 1 に提案法と従来法で得られた SIR の値を示す．提案法によって得られた分離音の SIR の平均は従来法の SIR の平均を 9.61 [dB] 上回った．Fig. 2 に無響下でのホッチキスの音のスペクトログラム (上)，従来法によって得られた分離音のスペクトログラム (中上)，提案法によって得られた分離音のスペクトログ
ラム (中下)，音響イベント検出結果 (下) を示す．黒がその時刻で推定された状態を表す．提案法によって残響が除去されていると共に，音響イベント検出がおおむね正しく行われていることがわかる．

5

Fig. 2 無響下でのホッチキスの音のスペクトログラ
ム (上)，従来法によって得られた分離音のスペクトロ
グラム (中上)，提案法によって得られた分離音のス
ペクトログラム (中下)，音響イベント検出結果 (下)を示す．黒がその時刻で推定された状態を表す．

4

評価実験

提案法の音源分離・残響除去性能と音響イベント検
出性能の評価のために，実験を行った．RWCP データ
ベース非音声ドライソース [8] の中のホッチキスの音，
ベルの音，携帯電話の着信音に対して同じく RWCP
データベース [8] のインパルス応答 (残響時間 600 ms，
マイク間距離 2.83 cm，マイクの数 M = 2) を畳み込
み，人工的に残響下での多チャンネルの混合信号を作成した．Fig. 1 に作成した混合音のスペクトログラ
ムを示す．サンプリング周波数は 32 kHz とした．フ
レーム長 16 ms，フレームシフト長 8 ms で STFT を行い，時間周波数展開を行った．HMM の状態数 D
は 2 とした．α1 と β1 を 1，10−2 とそれぞれ設定し，α2:3 と β2:3 を 1 と 1020 と設定することで，d = 1 を無音状態とみなした．
あらかじめ無響下で録音されたホッチキスの音，ベルの音，携帯電話の着信音に対して N = 0 としてそれぞれ提案法を適用し基底スペクトルと遷移確率を学習した後，混合音に対して N = 20 として提案法を適用した．
C 1:I,1:K,0 の初期値については，対角成分√を 1/ M ，それ以外の成分を 0 とした．
√C 1:I,1:K,1:Nの初期値については，対角成分を 10−3 / M ，それ以外の成分を 0 とした．パラメータ推定アルゴリズム
はまず N = 0 として 100 回反復し，その後徐々に Nを増やしながら全体として 150 回反復した．比較対

おわりに

本稿ではブラインド音源分離・残響除去・音響イベ
ント検出を統合的に行う手法を提案した．残響下での観測信号を，時間周波数領域における畳み込み混合の形でモデル化し，また音源信号の状態を隠れ変数とする 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] 樋口 他，音講論 (秋)，2014．
[7] C. Févotte, N. Bertin, and J.L. Durrieu, Neuralcomputation, vol. 21, no. 3, pp. 793–830, 2009.
[8] S. Nakamura et al., in Proc. LREC ’00, pp. 965–
968, 2000.
[9] E. Vincentet al., IEEE Trans. ASLP, pp. 1462–
1469, 2006.