社団法人 電子情報通信学会
THE INSTITUTE OF ELECTRONICS,
INFORMATION AND COMMUNICATION ENGINEERS

信学技報
TECHNICAL REPORT OF IEICE.

多チャンネル階乗隠れマルコフモデルによる音響情景分析のための統合的アプローチ
樋口 卓哉†

亀岡

弘和†,††† 東京大学大学院情報理工学系研究科〒 113–8654 東京都文京区本郷 7–3–1†† 日本電信電話株式会社 NTT コミュニケーション科学基礎研究所〒 243–0198 神奈川県厚木市森の里若宮 3–1
E-mail: †{higuchi,kameoka}@hil.t.u-tokyo.ac.jp
あらまし 本研究では，音源分離・音響イベント検出・残響除去・音源の到来方向推定という，音響情景分析に関する諸問題を取り扱う．これら音響情景分析に関する問題の根幹にあるのは，ブラインド音源分離の問題である．ブラインド音源分離の問題は不良設定問題であり，一般に音源に関して何らかの仮定を置かない限り，解を限定し解くことができない．本稿では，音響情景分析に関する諸問題が相互依存の関係にあることに着目し，音響情景に関する様々なパラメータによって観測信号を確率的にモデル化することで，ブラインド音源分離の問題における解を限定すると共に，パラメータ推論を通して統合的に音響情景分析を行う手法を提案する．
キーワード

非負値行列因子分解，ブラインド音源分離，残響除去，音響イベント検出，DOA

Unified approach for auditory scene analysis based on multichannel
factorial hidden Markov model
Takuya HIGUCHI† and Hirokazu KAMEOKA†,††† Graduate School of Information Science and Technology, The University of Tokyo
Hongo 7–3–1, Bunkyo-ku, Tokyo, 113–8654 Japan†† NTT Communication Science Laboratories, Nippon Telegraph and Telephone Corporation
Morinosatowakamiya 3–1, Atsugi-shi, Kanagawa, 243–0198 Japan
E-mail: †{higuchi,kameoka}@hil.t.u-tokyo.ac.jp
Abstract This paper deals with the problems of audio source separation, audio event detection, dereverberation
and DOA estimation. We focus on the fact that these problems are interdependent, and propose an unified approachfor these problem. We model the generative process of observed signals with parameters each of which corresponds
to a specific aspect of an audio source. Through the parameter inference, We can simultaneously salve the problems
of auditory scene analysis.
Key words Non-negative matrix factorization, blind source separation, dereverberation, audio event detection,
DOA

1. は じ め にが，それらの問題の根底にあるのは，ブラインド音源分離の問題である．ブラインド音源分離の問題とは，音源信号や音源か

本稿では，ブラインド音源分離や残響除去，音響イベント検らマイクまでの伝達特性が未知の場合に，複数の音源信号が混出，到来方向推定などの音響情景分析の諸問題を取り扱う．音合された観測信号から元の音源信号を推定する問題である．ブ

響情景分析技術の応用先として，聴覚障害者へのリアルタイム

ラインド音源分離の問題はいわゆる不良設定問題であり，一般

音響情景提示，音を利用した動画からの特定のコンテンツ抽出，的に，この問題を解くためには，音源信号に対して立てたなん

またロボットの音認識システムなどが考えられる．これらを実らかの仮定を基に最適化基準を立て，最適化問題を解く必要が現するためには音響情景分析に関する諸問題を解く必要がある

ある．当然のように，立てた仮定が成り立たない環境下ではそ—1—

の手法はうまく動作しないことが予想され，また立てた仮定が例えば，各音源の音響イベント (いつ鳴ったか) が分かってい弱すぎる場合でも，解を十分に制限することができず，音源分

れば，音源分離の問題は比較的容易になると考えられるが，一離を十分に行うことができない．従って，ブラインド音源分離

方で音源分離が十分うまく行えれば，各音源信号の音響イベン

を解くためには，対象とする環境下において成り立つ仮定を置

トを検出することは比較的容易になると考えられる．すなわち，

き，さらにその仮定が十分に解を制限することができる必要が

音響イベント検出の問題と音源分離の問題は相互依存の関係にある．

あるといえる．[6] の手法では，アクティベーションに対して音

例えば，観測信号数が音源数よりも多い優決定問題では，音源源の状態 (無音状態，音の立ち上がり，定常状態など) に依存し信号間の独立性を仮定して分離する独立成分分析 (Independent

た事前分布を仮定し，音量の大きな値のとりやすさを設定する

Component Analysis; ICA) が有用であることが知られており，ことで，音響イベント検出と音源分離を同時に実現していた．

音源信号間の独立性を最大化するように分離フィルタを推定す

また前述したように，従来の NMF では残響下では成り立たることが目的となる [1]．しかし，ICA では観測信号数が音源ない瞬時混合近似を用いていたため，残響下において精度良い

数よりも少ない劣決定問題を扱うことはできず，この場合は独音源分離を行うためには，事前に残響除去を行う必要がある．

立性よりもさらに強い仮定が必要である．

すなわち，音源分離性能は残響除去の精度に依存することにな

単チャンネルの観測信号に対するブラインド音源分離の有る．一方で，実環境においてどのように残響がかかるかは各音

効なアプローチとして，非負値行列因子分解 (Non-negative

源の位置によって異なるため，残響除去の精度は音源分離性能

Matrix Factorization; NMF) が知られている [2] [3]．この手法
では観測信号のパワースペクトログラムを，2 つの非負値行列の積に分解する．分解した各行列は，いくつかの基底パワース
ペクトルによって構成される基底行列と，それらの基底パワー
スペクトルの時変な音量を表すアクティベーションによって構成されるアクティベーション行列となる．ここで重要なのは，分解された各基底パワースペクトルが，観測信号の中で主となる要素，すなわち各音源信号を表していると考えられることである．また，音源信号の空間的な情報も利用して音源分離を行うために，NMF を多チャンネルの音響信号へと拡張するアプ
ローチがいくつか取られてきた [4] [5]．
しかし NMF(あるいはその多チャンネル観測信号への拡張手法) で立てられていた仮定が，実環境において必ずしも成り立
つとは限らない．例えば NMF では，観測信号を構成する限られた数の基底パワースペクトルが，それぞれ各音源信号を表していることを仮定していた．しかし実際の音源信号のパワース
ペクトルは時変であることが多く，1 つの基底パワースペクト
ルで表現するのは不十分な場合がある．そこで我々は以前，音源信号の時変なパワースペクトルを隠れマルコフモデルで表現し，多チャンネル NMF のモデルと統合することで，音響イ
ベント検出と音源分離を統合的に行う手法 [6] を提案した．また，[4] [5] では室内インパルス応答長が時間周波数展開における時間窓長に対して十分に短いことを仮定し，瞬時混合近似を用いていたが，残響下ではこの近似は成り立たないことが多い．

に依存する．すなわち，残響除去の問題と音源分離の問題は相互依存の関係にあるといえる．[7] の手法では，各音源ごとに残響成分を推定しながら音源分離を行うことによって，残響除去と音源分離を統合的に行う手法となっていた．
さらに，音源の到来方向推定の問題も音源分離の問題と相互依存の関係にあるといえる．点音源を仮定すると，一般的に，音源信号の直接波に対応する空間相関行列は，音源の到来方向に応じてある特定の構造を持つことが知られている．この事実に基づき，空間相関行列をモデル化する一つの手法として，音源の取りうるあらゆる到来方向を隠れ変数として，到来方向に基づいて空間相関行列をモデル化する手法が提案されている [8] [9]．また別のモデル化として，音源の取りうるあらゆる到来方向に基づいて陽に記述された空間相関行列の重みつき和として，空間相関行列をモデル化する方法 [10] も提案されている．音源の到来方向に応じて空間相関行列がある特定の構造を持つということは，音源信号の到来方向が分かっていれば，音源分離の手がかりに成りうることを意味しており，また当然ながら，音源分離が十分うまくいけば音源の到来方向推定は比較的容易であると考えられる．すなわち，音源の到来方向推定の問題と音源分離の問題もまた，相互依存の関係にあるのである．
そこで本稿では，多チャンネル NMF [5] における観測信号の生成モデルを拡張し，音響イベント・残響・音源の到来方向を表すパラメータによって観測信号の生成プロセスを確率的に記述することによって，パラメータ推論を通して音響イベント検この問題に対して我々は，観測信号を時間周波数領域における出・残響除去・音源の到来方向推定・音源分離を統合的に行う

畳み込み混合で表現することで，時間周波数展開におけるフ手法を提案する．本稿ではこのモデルを多チャンネル階乗隠れ

レーム外の残響をモデル化し，[6] の手法と組み合わせることにマルコフモデルと呼ぶ．

よって，残響除去・音響イベント検出・音源分離を統合的に行う手法を提案した [7]．
以上の説明は，従来の NMF において立てられていた仮定のうち，実環境では成り立たなかった仮定 (強すぎた仮定) を，より実環境下の音響信号に則した形に拡張した，という側面での説明であるが，一方でこれらの手法は，音響イベント検出や残響除去の問題という，音源分離の問題と相互依存にある問題
を，音源分離の問題と統合的に解くことによって，従来の NMF
ベースの手法では利用されなかった，音響イベントや残響成分をブラインド音源分離の問題を解く手がかりとしているという側面も存在する．

2. 音響イベント検出と音源分離の統合的アプローチ
2. 1 多チャンネル階乗隠れマルコフモデル
2. 1. 1 瞬時混合近似を用いた混合モデル
まず，観測信号の生成プロセスについて述べる．I 個の音
源信号が M 個のマイクロフォンで観測される場合を考える．
ym (ωk , tl ) ∈ C を m 番目のマイクで観測された観測信号の周
波数 ωk ，時刻 tl における時間周波数成分，si (ωk , tl ) ∈ C を i
番目の音源信号の周波数 ωk ，時刻 tl における時間周波数成分とし，1 <= L をそれぞれ時間周波数領域にお
=l<
=Kと1<
=k<

—2—

ける周波数と時間のインデックスとする．ここで，室内インパ

離の手がかりとすることのできるモデル化を行う．

ルス応答長が時間周波数展開における時間窓長よりも十分に短まず，多くの音源のパワースペクトルは，無音状態，音の立い場合を仮定すると，瞬時混合近似を用いて観測信号は以下のち上がり，定常状態などその音源の状態 (音響イベント) に応じように時間周波数領域において記述できる．

て異なると考えられるので，時間変化する音源の状態に応じて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. 1. 2 観測信号の生成プロセス
次に，式 (1) に基づいて観測信号の生成プロセスを確率的に記述する．まず，音源信号が区分的に定常であることを仮定し，2の複素正規分布に
各時間周波数点で si,k,l が平均 0，分散 σi,k,l従うとすると，音源信号の生成プロセスは，
2
si,k,l |σi,k,l ∼ NC (si,k,l ; 0, σi,k,l)，l における i 番目の音源の状態を表す隠れ変数 zi,l ∈ {1, . . . , Q}を導入し，状態の時系列 zi,1 , . . . , zi,L がマルコフ連鎖に従うと仮定すると，
zi,l |zi,l−1 ∼ Categorical(zi,l ; ρzi,l−1 )，(6)と 書 け る ．こ こ で Categorical(x; y) = yx で あ り，ρq =(ρq,1 , . . . , ρq,Q ) は状態 q から各状態 1, . . . , Q への遷移確率を表し， ρ = (ρq,q′ )Q×Q は遷移行列である．状態 q である i
番目の音源の基底パワースペクトルを wi,k,q と表すとすると，
時刻 l における i 番目の音源信号のパワースペクトルは zi,l に依存し，si,k,l の生成モデルは以下のように書き直せる．
si,k,l |wi,k,1:Q , hi,l , zi,l(2)

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

2と書き表せる．ここで σi,k,lは周波数 k ，時刻 l における i 番(7)

目の音源のパワースペクトル密度を表す．式 (1) と式 (2) から,さらに，音源の音量に着目すると，無音状態と有音状態では

a1:I,k と σ1:I,k,l が既知の条件下で観測信号 y k,l は同じく複素正規分布に従う．

当然音量の大きな値の取りやすさが異なると考えられるので，

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

(3)

i

音量もまた音源の状態に依存して異なる振る舞いをするといえる．そこで，hi,l が，zi,l によって異なるハイパーパラメータを持つガンマ分布に従うと仮定すると，

hi,l |zi,l ∼ Gamma(hi,l ; αzi,l βzi,l ),ただし C i,k = ai,k aH
i,k は i 番目の音源に対する周波数 k における空間相関行列と呼ばれる行列であり，また，NC (x; µ, Σ) ∝exp(−(x − µ)H Σ−1 (x − µ)) である．
2. 1. 3 多チャンネル NMF [4] [5] における生成モデル
通常の NMF のモデルにおいては ( [11] 参照)，音源信号のパワースペクトルはスケールを除いて時不変であることが仮定されていた． 上記のモデルにこの仮定を組み込むと，2
σi,k,l= wi,k hi,l ,(4)2となる．ここで σi,k,lは時不変の基底パワースペクトル wi,k と(8)となる．ここで α1:Q と β1:Q はそれぞれガンマ分布の形状パラメータとスケールパラメータであり，Gamma(x; α, β) =xα−1 e−x/β
Γ(α)β αである．zi,l が無音状態に対応するときは hi,l は小さな値をとってほしいので，小さな値をとる確率が高くなるよ
うにガンマ分布のハイパーパラメータを設定し，zi,l が有音状態に対応するときは一様分布に近くなるようにガンマ分布のハ
イパーパラメータを設定すればよい．具体的には，無音状態に対応する状態では α と β をそれぞれ 1，10−2 などと設定し，有音状態のときはそれぞれ 1，1020 などと設定すればよい．こ時変な音量を表す hi,l の積の形で表現されている．これによりのような設定により，音響イベント検出と音源分離を協調的に

si,k,l の生成モデルは以下のように書き直せる．

解くことが可能となる．

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

(5)

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

2. 1. 4 音響イベントに基づく音源信号の生成モデル [6]上記のように，多チャンネル NMF ではそれぞれの音源信号のパワースペクトルはスケールを除いて時不変であることが仮定されていた．この仮定により，解を限定しブラインド音源分
離の問題を解く手がかりとしていたのである．しかし多くの場
合，音源信号は時変なパワースペクトルを持つため，NMF で仮定されている時不変なパワースペクトルでは，音源信号を十分に表現することができない場合がある．そこで，音響イベン
トに基づいて音源信号の生成プロセスを記述することで，時変なパワースペクトルを表現すると共に，音響イベントを音源分

観測信号の最終的な生成モデルは a1:I,k ，w1:I,k,1:Q ， h1:I,l ，
z1:I,l が既知の条件下で，式 (6)，式 (8) と合わせて以下のように書ける．

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

(9)i

観測信号の生成モデルは，各音源信号に関するパラメータ

w1:I,1:K,1:Q ， h1:I,1:L や，各音源信号の音響イベントに関するパラメータ z1:I,1:L などに基づいて記述されているため，この生成モデルに基づいて最適なパラメータを求めることは，音源分離の問題と音響イベント検出の問題を統合的に解くことに相当している．
2. 2 パラメータ推論アルゴリズム
2. 2. 1 目 的 関 数
ここでは，上記の生成モデルに対する，補助関数法に基づくパ—3—

ラメータ推論アルゴリズムについて述べる．モデルにおける推定
したい変数は W = w1:I,1:K,1:Q , H = h1:I,1:L , C = C 1:I,1:K ，
Z = z1:I,1:L である．上記の変数の集合を Θ で表す．以下では
ρ は実験的に定められた定数とする．我々の目的は以下の式を満たす Θ̂ を求めることである．

+)∑((αzi,l − 1) log hi,l − hi,l /βzi,l − αzi,l log βzi,li,l+ log p(Z).

(15)

ここで Ri,k,l と U k,l は∑i Ri,k,l = I を満たすエルミート正定値行列であり， R と U の集合を Λ で表す. tr(·) は行列のトΘ̂ = argmax log p(Θ|Y ).

(10)

Θ

レースを表す． 式 (15) の等号成立条件は−1

ここで Y = y 1:K,1:L は多チャンネル観測信号の時間周波数成Ri,k,l = C i,k wi,k,zi,l hi,l X̂ k,l ,(16)分の集合である．ベイズの定理から，

U k,l = X̂ k,l ,(17)p(Θ|Y ) =p(Y , Θ)
,
p(Y )

(11)

となる．
以上から，L は次の 2 つのステップを繰り返すことによってであり，また 2.1 節で定義された条件つき確率を用いて，同時局所最大化できる．

確率 p(Y , Θ) は以下のように記述できる．

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

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

（ 2 ） W ，H ，C ，Z について L+ を最大化.

(12)

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

)+ log p(H|Z) + log p(Z) .

数，各音源のパラメータごとの和の形で書かれているために，ことができる．
ステップ 1 における R と U の更新では，式 (16) と式 (17)(13)

ここで，2.1 節で記述した生成モデルから，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 ). (14)
2
2
2
k,l

ただし X̂ k,l =

複雑な形をしていたが，補助関数 L+ (Θ, Λ) は各時間，各周波ステップ 2 において並列計算で各パラメータごとに更新を行う(Θ̂ = argmax log p(Y |Θ)Θ

元の目的関数 L(Θ) では，各パラメータがお互いに関係しあうを使えばよい．ステップ 2 では，L+ をそれぞれの変数で変微分して 0 となるものを求めることで，更新則が導ける．W とH に関する 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,zi,l hi,ll)− tr(U −1(18)k,l C i,k )hi,l ,(−1
H
1 ∑ tr(y k,l y k,l Ri,k,l C i,k Ri,k,l )
∂L+
=
∂hi,l
2
wi,k,zi,l h2i,lk)− tr(U −1C)wi,ki,k,zi,lk,l∑i C i,k wi,k,zi,l hi,l である．

2. 2. 2 補助関数法に基づく最適化アルゴリズム
式 (13)，(14) をみると，今回の目的関数は各パラメータ同士がお互いに関係しあう複雑な形をしており，一般的に最適化が困難であるが，補助関数法に基づく反復計算によって局所最適となる Θ を求めることができる．補助関数法による，目的関数 F (Θ) の最大化問題の最適化アルゴリズムでは，まず補助変数 Λ を導入し，F (Θ) = maxΛ F + (Θ, Λ) を満たす補助関数 F + (Θ, Λ) を設計する．そして，Θ ← argmaxΘ F + (Θ, Λ)と Λ ← argmaxΛ F + (Θ, Λ) を 交 互 に 繰 り 返 し ，補 助 関 数
F + (Θ, Λ) の最大化を行うことで，間接的に元の目的関数 F (Θ)
を最大化する．ここで重要なのは，Θ について最大化しやすいF + (Θ, Λ) を設計することである．
今 回 の 最 適 化 問 題 に 補 助 関 数 法 を 適 用 す る た め ，ま ず
L(Θ) = maxΛ L+ (Θ, Λ) を満たす補助関数 L+ (Θ, Λ) を設計す
る．今回の場合は，以下のように補助関数 L+ (Θ, Λ) を設計できる．

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

(19)

これらを 0 と置くことで，以下の更新式が導ける．
vu ∑ tr(y y H R−1ui,k,l C i,k Ri,k,l )k,l k,l
u l
hi,l
t
,
wi,k,zi,l ←∑−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

(20)

ただし

µi,l =+>=L (Θ, Λ)( (−1
H1 ∑ ∑ tr(y k,l y k,l Ri,k,l C i,k Ri,k,l )=−2
wi,k,zi,l hi,l
i
k,l
)
)
+
log
|U
|
−
M
+ tr(U −1
h
C
)w
k,l
i,l
i,k
i,k,z
k,l
i,l

∑ tr(y k,l y k,l H Ri,k,l C −1
i,k Ri,k,l )kνi,l =
L(Θ)

(21)

∑

,

(22)

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

wi,k,zi,lkである．C についての L+ の変微分は以下のようになる．
−1
H∑( C −1∂L+i,k Ri,k,l y k,l y k,l Ri,k,l C i,k=∂C i,kwi,k,zi,l hi,ll)h.
− U −1
w
i,l
i,k,z
k,l
i,l

(24)

—4—

ける伝達周波数特性の時刻 tn の成分であり，0 <
=n<
=N は

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

伝達周波数特性の時間周波数領域における時間インデックスでC i,k Ai,k C i,k = B i,k ,(25)において tn だけ先の時刻にどれだけ影響を与えるか，すなわ

ただし

Ai,k =

ある．ここで ai (ωk , tn ) は i 番目の音源信号が時間周波数領域∑

ち時間周波数領域においてどれだけ残響がかかるかを表して−1

wi,k,zi,l hi,l X̂ k,l ,(26)

いる．観測信号の最終的な生成モデルは a1:I,k,0:N ，w1:I,k,1:Q ，∑−1−1
B i,k = C i,k (
wi,k,zi,l hi,l X̂ k,l y k,l y k,l H X̂ k,l )C i,k , (27)h1:I,l−N :l ，z1:I,l−N :l が既知の条件下で，式 (6)，式 (8) と合わ
せて以下のように書き直せる．

llである．以下の方法でこの Riccati 方程式を解くことで C の更新則が得られる [5]．まず，以下の 2M × 2M の行列に対して

y k,l |a1:I,k,0:N , w1:I,k,1:Q , 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 ).

(32)

i,n

固有値分解を行う．

[

0−B i,k

]
−Ai,k
.
0

この生成モデルに基づいて最適なパラメータを求めることは，(28)

ことに相当している [7]．

ここで e1,i,k . . . eM,i,k を負の固有値に対応する固有ベクトルだとし，2M 次元の固有ベクトルを m = 1 . . . M において以下のように分解する．

[

f m,i,kem,i,k =g m,i,k .

](29)

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

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

残響除去・音響イベント検出・音源分離の問題を統合的に解く(30)

ただし F i,k = [f 1,i,k , . . . , f M,i,k ]，Gi,k = [g 1,i,k , . . . , g M,i,k ]である．

3. 2 パラメータ推論アルゴリズム
3. 2. 1 目 的 関 数
式 (32) の生成モデルにより，2.2.1 節の目的関数における式∑(14) の X̂ k,l を，X̂ k,l = i,n C i,k,n wi,k,zi,l−n hi,l−n と置き換えたものが，今回のモデルにおける目的関数に相当する．これ
により，L(Θ) は，各時刻のパラメータ同士が互いに関係しあう形となってしまうが，この場合も補助関数法を適用し，各時刻ごとのパラメータの和の形で表された補助関数を設計することで，効率的な最適化アルゴリズムを導出できる．
3. 2. 2 補助関数法に基づく最適化アルゴリズム
この場合の補助関数は，以下のように設計できる．
L(Θ)

L+ は各音源のパラメータごとの和の形で書かれているために，L+ を Z の関数と見ると，各音源に対応する隠れマルコフ
モデルの対数事後確率の和とみなすことができる．従って最適な状態の時系列 zi,1 , . . . , zi,L を動的計画法によって効率的に，各音源ごとに個別に求めることができる．
以上の更新則から W ，H ，C ，Z を反復計算により求める
ことは，ブラインド音源分離の問題と音響イベント検出の問題を画一的最適化規準に基づき協調的に解いていることに相当している．
2. 3 動作実験例
以上の提案法により劣決定音源分離実験を行った結果，多チャンネル NMF [5] と比較して，高い音源分離性能を示した．
また，音響イベント検出もある程度の精度で行うことができた．
実験の詳細は [6] を参照してほしい．

3. 残響除去との統合的アプローチ+>=L (Θ, Λ)( (−1
H1 ∑ ∑ tr(y k,l y k,l Ri,k,l,n C i,k,n Ri,k,l,n )=−2
wi,k,zi,l−n hi,l−n
i,n
k,l
)
)
+ log |U k,l | − M+ tr(U −1
k,l C i,k,n )wi,k,zi,l−n hi,l−n)∑(+(αzi,l − 1) log hi,l − hi,l /βzi,l − αzi,l log βzi,li,l+ log p(Z).
ここで Ri,k,l,n と U k,l は(33)∑

i,n Ri,k,l,n = I を満たすエルミー

ト正定値行列であり， R と U の集合を Λ で表す. 式 (33) の等号成立条件は−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 ,

(34)
(35)

3. 1 畳み込み混合による観測信号の混合モデル
式 (1) による混合モデルでは，瞬時混合近似を仮定していた．
しかし残響がある場合には一般に，室内インパルス応答長は時間窓長に対して十分に短いとはいえず，瞬時混合近似は成り立
たない．そこで，時間周波数領域における畳み込み混合の形でとなる．観測信号の生成プロセスを畳み込み混合の形で表現し観測信号を近似する．

計算によって求めることができると共に，音響イベントを表す

y(ωk , tl ) ≈

I ∑
N
∑

ていたがゆえに，元の目的関数は各時刻のパラメータが他の時
刻のパラメータと関係しあう複雑な形をしていたが，設計した補助関数は，各時刻ごとのパラメータの和の形で表現されていることが分かる．これにより，各時刻ごとのパラメータを並列
隠れ変数 zi,1:L の時系列を，動的計画法により効率的に求める

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

(31)

i=1 n=0

ここで ai (ωk , n) は i 番目の音源信号に対する周波数 ωk におことが可能となる．
紙面の都合上詳細は省略するが，2.2.2 節の最適化アルゴリ
ズムと同様，各パラメータの更新則は以下の形で求まる．

—5—

wi,k,zi,lvu ∑ tr(y−1Huk,l+n y k,l+n Ri,k,l+n,n C i,k,n Ri,k,l+n,n )
u l,n
h
i,l
←t
, (36)∑−1l,n tr(U k,l+n C i,k,n )hi,l√(αzi,l − 1) + (αzi,l − 1)2 + µi,l νi,l
hi,l ←
,
(37)
νi,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 )k,nνi,l =

wi,k,zi,l∑, (38)

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

(39)

k,n

である．C i,k,n の更新については，以下の Riccati 方程式を章と同様の手順で解けばよい．

C i,k,n Ai,k,n C i,k,n = B i,k,n ,

(40)

ただし

Ai,k,n =∑−1

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

(41)

l

∑
B i,k,n = C i,k,n (wi,k,zi,l−n hi,l−nl−1−1X̂ k,l y k,l y k,l H X̂ k,l )C i,k,n ,(42)である．Z の更新については，前述したように，L+ を各音源，

図 1

無響下で録音された音源信号のスペクトログラム (上)，2 章の手法によって得られた分離音のスペクトログラム (中上)，本章の提案法によって得られた残響除去済み分離音のスペクトログ

各時刻ごとのパラメータの和の形になるように設計したことに

ラム (中下)，提案法による音響イベント検出結果 (下)．音響イ

よって，最適な音源の状態 zi,1 , . . . , zi,L を動的計画法によってベント検出結果は，黒がその時刻に推定された状態を表す．

効率的に，各音源ごとに個別に求めることができる．

3. 3 動作実験例
以上の提案法により残響下で録音された観測信号に対して，
教師あり音源分離・残響除去・音響イベント検出を試みた結果，
2 章で説明した [6] の手法と比較して，高い音源分離性能を示した．図 1 に無響下で録音された音源信号のスペクトログラム(上)，2 章の手法 [6] によって得られた分離音のスペクトログラ
ム (中上)，本章の提案法によって得られた残響除去済み分離音のスペクトログラム (中下)，提案法による音響イベント検出結
果 (下) を示す．実験の詳細は [7] を参照してほしい．

ここで ȷ は虚数単位，B [m] はマイクロフォン間の距離，C[m/s] は音速である．i 番目の音源の到来方向 θi が既知の場合
では，直接波に対応する空間相関行列は J(θi , ωk ) と等しくなることが期待される．しかし実際には，時間周波数展開におけるフレーム内に含まれる残響成分によって，C i,k,0 は J (θi , ωk )と等しくならないことが予想されるばかりでなく，音響信号か
ら直接音源の到来方向を観測することはできない．そこで，[10]の手法と同様に，まず到来方向を O 個の到来方向 ϑ1 , . . . , ϑOに離散化し，さらに重み定数 di,1 . . . di,O を導入することで，空4. 音源の到来方向推定との統合的アプローチ

間相関行列を以下のようにモデル化する．

4. 1 音源の到来方向に基づく空間相関行列のモデル化
上記のモデルでは，空間相関行列に対して何の仮定も置かれておらず，自由度の高いモデルとなっていたがゆえに，推定した空間相関行列が望ましくない局所解に陥ってしまうことがあった．しかし一般的に，音源信号の直接波に対応する実際の空間

C i,k,0 =相関行列は，点音源を仮定すると，音源の到来方向に応じてある特定の構造を持つことが知られている．従って，音源の到来
方向に基づいて空間相関行列をモデル化することが可能である．
マイクロフォンの数 M = 2 の場合では，方向 θ(0 <=θ<= π) にある音源の空間相関行列は，以下のように陽に記述可能である．
[
]
[
]∗
1
J (θ, ω) = ȷωB cos θ/C 1 eȷωB cos θ/C
(43)
e

∑

di,o J(ϑo , ωk ).

(44)

o

∑
ここで di,1 . . . di,O は非負値の値であり， o di,o = 1 を満たす．このモデル化により，実際の伝達特性にある程度即した形
で，空間相関行列の自由度を制限することができる．さらに，

di,1 . . . di,O を推定することによって，音源の到来方向推定を同時に行うことが可能となる．なぜなら，比較的大きな値として推定された重み定数に対応した到来方向が，音源の到来方向であると期待されるからである．
観 測 信 号 の 最 終 的 な 生 成 モ デ ル は a1:I,k,0:N ，w1:I,k,1:Q ，
h1:I,l−N :l ，z1:I,l−N :l ，d1:I,1:O が既知の条件下で，式 (6)，式(8) と合わせて以下のように書き直せる．

—6—

y k,l |a1:I,k,0:N , w1:I,k,1:Q , h1:I,l−N :l , z1:I,l−N :l , d1:I,1:O∑∼ NC (y k,l ; 0,di,o J(ϑo , ωk )wi,k,zi,l hi,l
+

∑

+

∑

λq,i,l log p(Z),

ここで Ri,k,l,n,q ，U k,l はi,oC i,k,n wi,k,zi,l−n hi,l−n ).

(45)

∑

i,n Ri,k,l,n = I を満たすエルミー

ト正定値行列であり，λq,i,l は∑q λq,i,l = 1 を満たす非負値のスカラー値である．式 (47) の等号成立条件は，Ri,k,l,n,q ，U k,l

i,n=0
|

この生成モデルに基づいて最適なパラメータを求めることは，については音源の到来方向推定・残響除去・音響イベント検出・音源分離−1

Ri,k,l,n,q = C i,k,n wi,k,zi,l−n =q hi,l−n X̂ k,l ,の問題を統合的に解くことに相当している．

4. 2 パラメータ推論アルゴリズム
4. 2. 1 目 的 関 数∑
3.2.1 節の目的関数において，C i,k,0 を o di,o J(ϑo , ωk )と置き換え，また音源の状態数 Q を増やした場合に，i 番目の音源の音響イベントを表す zi,1:L が望ましくない局所解に落ち

U k,l = X̂ k,l .

ことを考える．今回の目的関数は，∑
L(Θ) = log
p(Y |W , H, C, D, Z)
Z∑
Z

p(H|Z) + log∑ 1∑=(− log|X̂ k,l |2Z∑λq,i,l = Fq,i,l Bq,i,l /∑Fq,i,l Bq,i,l ,∑∑Fq,i,l−1 ρzi,l−1 =q′ ,zi,l =q ,Z

Bq′ ,i,l+1 p(Θ|zi,l+1 = q ′ )ρzi,l =q,zi,l+1 =q′ , (52)

ただし

(46)

Z

わりに d1:I,1:O を推定する必要がある．目的関数をみると，ここでもまた，各到来方向に対応する重み定数 d1:I,1:O が互いにp(Θ|zi,l = q)
(
( (
−1
H
1 ∑ ∑ tr(y k,l y k,l Ri,k,l,0,q J i,k,o Ri,k,l,0,q )∝ exp −
2di,o wi,k,zi,l =q hi,l
i,o
k,l
)
−1
+ tr(U k,l J i,k,o )di,o wi,k,zi,l =q hi,l+により，各音源ごと，各到来方向ごとの和の形になるように補助関数を設計し，効率的な最適化アルゴリズムを導出すること∑ ( tr(y k,l+n y k,l+n H Ri,k,l+n,n,q C −1
i,k,n Ri,k,l+n,n,q )

wi,k,zi,l =q hi,l))+ tr(U −1C)w
hi,k,zi,l =q i,lk,l+n i,k,ni,n=0|

関係しあう形となっているが，これも補助関数法を用いること

+

∑(

))
(αzi,l =q − 1) log hi,l − hi,l /βzi,l =q − αzi,l =q log βzi,l =q

i,l

(53)

ができる．

4. 2. 2 補助関数法に基づく最適化アルゴリズム
2, 3 章の最適化アルゴリズムでは，i 番目の音源の音響イベントを表す zi,1:L の最適な時系列をビタビアルゴリズムにより求めていたが，本章では各時刻でどの状態にあるらしいかを求めるアルゴリズムを導出する．この場合の補助関数は，以下のように設計できる．

L(Θ)+>=L (Θ, Λ)((tr(y k,l y k,l H Ri,k,l,0,q J −11∑ ∑i,k,o Ri,k,l,0,q )
λq,i,l=−h2dwi,o i,k,qi,l i,l
i,q,o
k,l
)
+ tr(U −1
k,l J i,o )di,o wi,k,qi,l hi,l(tr(y k,l y k,l H Ri,k,l,n,q C −1
i,k,n Ri,k,l,n,q )
wi,k,qi,l−n hi,l−n
i,q,n=0
|
)
)
h
+
log
|U
|
−
M
+ tr(U −1
C
)w
k,l
i,k,n
i,k,qi,l−n i,l−n
k,l
+

+

∑
i,l,q

(51)

q′

q′

p(Z)

∑ (∑
となる．ここで X̂ k,l =o di,o J(ϑo , ωk )wi,k,zi,l hi,l +
i
)
∑
C
w
h
であり，今回は
C 1:I,1:K,0 の代
i,k,n i,k,zi,l−n i,l−n
n=0
|

∑

(50)

q

Fq,i,l = p(Θ|zi,l = q)

k,l

Z

(49)

ズムを用いて，

Bq,i,l =∑−1
1
exp y k,l H X̂ k,l y k,l )
− log
2
Z
∑
∑
+ log
p(H|Z) + log
p(Z),

(48)

である．さらに λq,i,l に関しては Forward-Backward アルゴリてしまうことを防ぐため，目的関数を Z について周辺化する

+ log

(47)

q

λq,i,l()λq,i,l (αqi,l − 1) log hi,l − hi,l /βqi,l − αqi,l log βqi,lである．di,o の更新則は，2.2.2 節の最適化アルゴリズムと同様の手順により，以下の形で求まる．

vu∑
tr(y k,l y k,l H Ri,k,l,0 J −1R)ui,k,o i,k,l,0
u k,l,q λq,i,l
wi,k,z =q hi,li,ldi,o ← t ∑. (54)−1
k,l,q λq,i,l tr(U k,l J i,k,o )wi,k,zi,l =q hi,l4. 3 評 価 実 験
提案法の音源分離・到来方向推定性能の評価のために，実験を行った．ATR 音声データベース [12] の中の音声 (男性 1 人，
女性 2 人) に RWCP データベース [13] のインパルス応答 (残響時間 380 ms，マイク間距離 11.48 cm，マイクの数 M = 2)を畳み込み，人工的に残響下での多チャンネルの混合信号を作成した．音源の到来方向はそれぞれ π/6，π/2，13π/18 [rad]である．音声を代えて 10 個の混合信号を作成し，実験に用いた．サンプリング周波数は 16 kHz とした．フレーム長 64 ms，
フレームシフト長 16 ms で STFT を行い，時間周波数展開を行った．HMM の状態数 Q は 5 とした．α1 と β1 を 1，10−1とそれぞれ設定し，α2:10 と β2:10 を 1 と 1010 と設定すること
で，d = 1 を無音状態とみなした．D の初期値は図 2(上) のよ—7—.

表 1 分離処理前の混合音の SDR/SIR の平均値と，提案法と [7] の手法によって得られた分離音の SDR/SIR の平均値．

SDR [dB] SIR [dB]
提案法

-4.17

[7] の手法

-6.49

1.94

分離処理前

-40.33

-4.22

5.90

図 2 音源ごとに色分けされた，D の初期値 (上) と提案法による Dの推定結果例 (下)．

うに設定した．C 1:I,1:K,1:N の初期値については，対角成分を√10−1 / M ，それ以外の成分を 0 とした．W の初期値は乱数で与え，H の初期値は一様とし，λq,i,l の初期値は 1/Q とした．
パラメータ推定アルゴリズムはまず N = 0 としてその後徐々に N を N = 3 となるまで増やしながら全体として 30 回反復した．比較対象には [7] の手法を用いた．分離音 ŷ i,k,l はウィー
ナーフィルタ−1ŷ i,k,l = wi,k,zi,l hi,l C i,k,0 X̂ k,l y k,l .

(55)

によって得た．客観評価基準として，signal-to-distortion/interference ratio (SDR/SIR) [14] を用いた．高い SDR/SIR
は高い音源分離性能を表す．
表 1 に提案法と従来法で得られた SDR/SIR の平均値を示す．提案法によって得られた分離音の SDR と SIR は従来法を上回った．図 2(下) に D の推定結果例を示す．音源の到来方向
がおおむね推定できていることが分かる．

5. お わ り に
本稿では，音響情景分析に関する様々な諸問題を取り扱った．
それらの諸問題が相互依存の関係にあることに着目し，音響情

ponent Analysis, John Wiley & Sons, 2001.
[2] D. D. Lee, and H. S. Seung, “Learning the parts of objects
with nonnegative matrix factorization,” Nature, vol. 401,
pp.788–791, 1999.
[3] P. Smaragdis, and J. C. Brown, “Non-negative matrix factorization for polyphonic music transcription,” WASPAA
2003, pp. 177–180, Oct. 2003.
[4] A. Ozerov, and C. Févotte, “Multichannel nonnegative matrix factorization in convolutive mixtures for audio source
separation,” IEEE Trans. Audio, Speech and Language
Processing, vol. 18, no. 3, pp. 550-563, Mar.2010.
[5] H. Sawada, H. Kameoka, S. Araki and N. Ueda, “Eﬃcientalgorithms for multichannel extensions of Itakura-Saito nonnegative matrix factorization,” ICASSP 2012, pp. 261–264,
2012.
[6] T. Higuchi, H. Takeda, T. Nakamura and H. Kameoka, “A
unified approach for underdetermined blind signal separation and source activity detection by multichannel factorial
hidden Markov models,” Interspeech 2014, pp. 850–854,
2014.
[7] T. Higuchi and H. Kameoka, “Joint audio source separationand dereverberation based on multichannel factorial hidden
Markov model,” MLSP 2014.
[8] H. Kameoka, M. Sato, T. Ono, N. Ono, S. Sagayama, “Blindseparation of infinitely many sparse sources,” IWAENC
2012, H-09, Sep. 2012.
[9] T. Higuchi, N. Takamune, T. Nakamura, H. Kameoka,”Underdetermined blind separation and tracking of moving
sources based on DOA-HMM,” ICASSP 2014, pp. 32153219, May 2014.
[10] J. Nikunen and T. Virtanen, “Multichannel audio separation by direction of arrival based spatial covariance modeland non-negative matrix factorization,” ICASSP 2014, pp.
6727–6731, May 2014.
[11] C. Févotte, N. Bertin, and J. -L. Durrieu, “Nonnegative
matrix factorization with the itakura-saito divergence. withapplication to music analysis,” Neural Computation, vol.
21, no. 3, 2009.
[12] A. Kurematsu, K. Takeda, Y. Sagisaka, S. Katagiri,
H. Kuwabara, and K. Shikano,, K. Takeda, Y. Sagisaka,
S. Katagiri, H. Kuwabara, and K. Shikano, “ATR Japanese
speech database as a tool of speech recognition and synthesis,” Speech Communication, pp. 357–363, 1990.
[13] S. Nakamura, K. Hiyane, F. Asano, T. Nishiura, and T. Yamada, “Acoustical sound database in real environments forsound scene understanding and hands-free speech recognition,” LREC 2000, pp. 965–968, 2000.
[14] E. Vincent, R. Gribonval, and C. Févotte, “Performancemeasurement in blind audio source separation,” IEEE
Transactions on Audio, Speech, and Language Processing,pp. 1462–1469, 2006.

景を表す様々なパラメータによって観測信号の生成プロセスを確率的に記述することで，パラメータ推論を通して統合的な音響情景分析を行うアプローチを提案した．さらに，本提案モデルにおける対数事後確率は，様々なパラメータが互いに関係しあう複雑な形をしているが，補助関数法を用いることによって，目的関数を最適化する効率的な反復アルゴリズムを導出した．
本稿は [6] と [7] の内容をまとめた上で，到来方向推定を同時に行うために，[7] の手法を拡張したものである．

6. 謝

辞

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

献[1] A. Hyvärinen, J. Karhunen, and E. Oja, Independent Com-

—8—