2-P-22
複合自己回帰系による音声パワースペクトル密度モデルを用いたブラインド音源分離と残響除去∗☆濱村真理子 (NTT CS 研/九大院・芸工), 亀岡弘和, 吉岡拓也
ルルージョナトン, 柏野邦夫 (NTT CS 研)

1

1M T
源成分を並べたベクトルを Sk,n = (Sk,n, · · · , Sk,n) とはじめに
ブラインド音源分離 (BSS) やブラインド残響除去では，音源の成分と音源からマイクロホンまでの伝達特性する。本稿では [4] に倣い，多チャンネル有限インパルス応答によって表される分離システムがともに未知のもとで，マイクロホン入力信号から音源Sk,n =成分を復元することが目的であり，音源に関する何らかの仮定をもとに音源らしさの規準を定め，これにより立てられる目的関数を最小 (最大) 化する最適化問題として定式化されることが多い。音源らしさの規準は，どのような信号を生成しやすいか，という何らかの傾向を表した統計モデルを仮定することにより定義できる。どのような統計モデルを仮定すべきかは，対象とする音源信号の実際のふるまいや現象にいかに即しているか，数理的技法が馴染みやすい最適化問題に帰着するか，を考慮して判断することが重要である。

ズム [1] は，混合過程のインパルス応答が長くなると計算量が膨大になる問題があるのに対し，短時間 Fourier
変換 (STFT) 領域においてフレームが混合過程のインWk,τ Yk,n−τ ,(1)τ =0を用いる。ただし，Wk,τ , 0 ≤ τ ≤ nk は M ×M 行列の分m離フィルタである。さて，ここで，Sk,nを確率変数とし，′m
Sk,nと Skm′ ,n′ は (k, n, m) ̸= (k ′ , n′ , m′ ) のとき独立と仮mmmm (s定する。Sk,nの確率密度関数 fSk,nk,n ; θ ) は音源信号の統計モデルに該当し，θm はその分布を特徴づけるモデルパラメータである。この統計モデルの具体形が決まmmm (sれば Y := {Yk,n }k,n の確率密度関数は fSk,nk,n ; θ )を用いて自動的にBSS は音源の混合過程に関する仮定に応じてさまざまなアプローチに大別でき，それぞれに一長一短がある。
時間領域の畳み込み混合モデルに基づく BSS アルゴリnk∑

fY |Θ (Y |Θ) =∏k| det Wk,0 |∏( m m)
m |θ m S̄
fSk,nk,n|θ

t,mmと定義される。ただし，S̄k,nは分離信号∑τ Wk,τ

Yk,n−τ の m 番目のベクトル要素である。この同時密度関数は，未知パラメータ Θ := {{θm }m , {Wk,τ }k,τ }の尤度関数に対応するものであり，[4] のように最尤推

パルス応答より十分長い場合には時間領域の畳み込み混定の枠組でブラインド音源分離とブラインド残響除去を合は周波数領域の瞬時混合で近似することができ，効率同時に行うための目的関数となる。の良いアルゴリズムが実現可能となる [2]。しかし，混音源の統計モデルをどのように定めるかによっても合過程のインパルス応答が STFT のフレームの長さをBSS のアプローチは大別することができる。例えば音声超える残響環境下では，上記の近似の精度は必然的にの短時間信号のモデルとして自己回帰モデルが有名であ低下してしまう。一方近年の残響除去の研究では残響音るが，[5] や [4] では，自己回帰モデルに基づいて確率密声が STFT のサブバンドごとの音声信号と室内インパルス応答の畳み込みにより良く表せることが示されてmmm (s度関数 fSk,nk,n |θ ) を記述した，最尤法に基づく BSSの解法が提案されている。本稿では，[4] の BSS システ

おり [3]，高残響環境下では STFT 領域における畳み込ムの枠組に，自己回帰モデルに代わる新しい音声の統計み混合モデルをベースにした BSS [4] が性能および効率モデルを組み込むことで BSS の性能向上を図ることが性の面で有効であることが示されている。以下，このモ主目的である。

デルをベースに BSS の問題を具体的に示す。まず，M個の音源の信号を M 個のマイクロホンで収音する状況を考え，これらの収音信号を STFT により展開した時間周波数成分を観測データと見なすことにしよう。ここl(たで，l 番目のマイクロホンで観測される信号を Yk,n

だし，k, n は周波数，時刻に対応するインデックス) とM T1) とする。また，M 個の音, · · · , Yk,nし，Yk,n = (Yk,n∗

複合自己回帰系による音声の統計モデル

2自己回帰モデルは，白色の Gauss 性雑音を入力とした全極システムからの出力と捉えることができ，入力とフィルタをそれぞれ声帯振動波形と声道伝達特性に対応させることで音声生成過程の一つの簡潔なモデルと見なすことができる。しかし，駆動信号を白色雑音で仮定すBlind source separation and dereverberation using speech power spectral density model derived from composite autoregressive system. by HAMAMURA Mariko (NTT Communication Science Laboratories/Kyushu
University), KAMEOKA Hirokazu, YOSHIOKA Takuya, LE ROUX Jonathan and KASHINO Kunio (NTTCommunication Science Laboratories)

日本音響学会講演論文集- 699 -2010年9月

ることは，駆動信号が一様なパワースペクトル密度 (以の混ざり具合により各フレームの音声スペクトルが特徴後，PSD) をもつことを仮定していることに相当し，高づけられる。n 番目のフレームにおいて，i 番目の駆動信い基本周波数になるほど適合しなくなる。従って，駆動号 PSD と j 番目の全極型フィルタの組み合せから生成さ信号の PSD 自体も推定すべきパラメータとして扱いたれる音声要素のゲインを Uni,j とすると，当該音声要素のいが，白色性に代わる何らかの新たな制約が必要であi,ji,j T
DFTXni,j = (X1,n, · · · , XK,n) は多次元複素正規分布
i,ji,ji,ji,jNC (0, Λi,j)に従う。
Λは
Λnnn = diag(λ1,n , · · · , λK,n )の対角共分散行列で，対角要素はる。そこで，発話全体を通して声の高さや音素の種類が限られているという音声の特徴に着目し，高々数種類の声道伝達特性と高々数種類の駆動信号 PSD の組み合せで音声を表現できるという仮定の下考案したのが以前提

λi,jk,n =Hki Uni,jj|A (ej2πk/K )|2(5)

案した複合自己回帰系 [6] に基づく音声信号のモデルでAj (z) = 1 − aj1 z −1 · · · − ajP z −P(6)ある。以下，このモデルについて具体的に説明する。

′′

ここからは特定の音声信号について議論することとで与えられる。ここで，Xni,j と Xni ′,j は (i, j) ̸= (i′ , j ′ )し，音源成分のインデックス m は省略する。 ここで nまたは n ̸= n′ のとき独立と仮定すると，Xn1,1 , · · · , XnI,J番目のフレームにおける信号 {xn [t]}Tt=1 を，P 次の自の和，すなわち音声信号の STFT Sn ∈ CK はやはり正己回帰過程

規分布に従い，xn [t] =P∑

ap xn [t − p] + ϵn [t]Sn =(2)∑∑ip=1Φn =からの標本値系列と仮定する。ただし，ϵ[t] は声帯駆動()
Xni,j ∼ NC 0, Φn(7)Λi,jn(8)j∑∑ij

信号に対応し，平均 0 の Gauss 性雑音としているが，必が言える。式 (8) は，正規分布する独立な確率変数の和ずしも白色雑音とは限らない点を強調しておく。xn =は各々の正規分布の分散の和を分散とする正規分布に従(xn [1], · · · , xn [K])T ∈ RK と定義し，その離散 Fourier
変換 (以後，DFT) を Xn = (X1,n , · · · , XK,n )T ∈ CKとする。式 (2) より，DFT の線形性，および，ϵn [t] の定常性，Gauss 性から，Xn は，平均が 0，分散共分散()行列が Λ = diag λ1 , · · · , λK の多次元複素正規分布

NC (0, Λ) に従う。ただし，で与えられる。ただし，θ = {Hki , ajp , Uni,j }k,n,i,j は以上の統計モデルにおける未知パラメータである。Φn の対角要素の ϕk,n は複合系により生成される不規則信号のHk Un,|A(ej2πk/K )|2(3)A(z) = 1 − a1 z −1 · · · − aP z −P(4)λk,n =うことを意味する。以上より，Sk,n の統計モデルは()
1|sk,n |2fSk,n (sk,n |θ) =exp −(9)πϕk,nϕk,nである。H1 , · · · , HK は駆動信号 ϵ[t] の PSD であり，スペクトル微細構造に対応する。一方で，1/|A(ej2π1/K )|2 ,1/|A(ej2π2/K )|2 , · · · , 1/|A(ej2πK/K )|2 は全極型伝達関数の PSD であり，スペクトル包絡に対応する。
以上をもとに，複合自己回帰系を構築することができる。複合自己回帰系は，I 種類の駆動信号 PSD と J 種類の全極型フィルタにより構成される複合系であり，i 番目の駆動信号 PSD と j 番目の全極型伝達関数をそれぞれHki , 1/Aj (ej2πk/K ) とする。この複合系は，すべての駆動信号 PSD と全極型フィルタを組み合わせた計 I × J 種類

PSD であり，ϕk,n =∑∑ijHki Uni,j|Aj (ej2πk/K )|2(10)である。ところで，J = 1, P = 0 の特殊ケースにおい∑て，ϕk,n は ϕk,n = i Hki Uni,1 のとおり行列積の形となり，この下での式 (9) は，非負値行列因子分解 (NMF)のある統計的解釈に基づいて導かれた尤度関数 [8] と同じ形となる。よって，複合自己回帰系による音源モデルを用いた全体の BSS システムは多チャンネル NMF[9]と構造的に深い関連があることが示唆される。

最適化アルゴリズム

3の音声要素を生成することができる。ここで例えば，各観測信号 Y が与えられたとき，未知パラメータ Θ =フレームで単一の音声要素だけがアクティブになると仮{{θm }m , {Wk,τ }} についての事後確率 fΘ|Y (Θ; Y ) ∝fY |Θ (Y ; Θ) fΘ (Θ) の対数定した場合，この複合系の出力信号の統計モデルは [7] と同様の形態の混合正規分布により記述することができる。

L(θ) := log fY |Θ (Y |Θ) + log fΘ (Θ)

これに対し，複合自己回帰系では，これらの計 I ×J 種類(11)の音声要素の和を音声信号と考える。各音声要素のゲイを最大化する θ を探索するアルゴリズムについて述べンがフレーム n ごとに唯一異なる値を取り得るパラメーる。この最適化問題は解析的に求めることはできないが，タであり，この値を通して決まる I × J 種類の音声要素

次の 3 つの最大化ステップ

日本音響学会講演論文集- 700 -2010年9月

(S1)θm← argmax L(θ) for all m,θm(S2)

Wk,0(S3)

k{Wk,τ }nτ =1← argmax L(θ) for all k
Wk,0← argmax L(θ) for all kを得る。ここで，不完全データ S̄k,n と完全データ Xk,nに S̄k,n = 1T Xk,n (ただし，1 = (1, · · · , 1)T ) なる多対[]H一関係があるため，E Xk,n Xk,n|Sk,n = S̄k,n , θ = θ′は具体的にnk{Wk,τ }τ =1を繰り返すことにより局所最適解を求めることができる。もし，各音源の PSD の推定値 ϕk,n が分かれば (S2)と (S3) の更新は [4] で述べられているプロセスで行うことができるため，(S1) の部分最適化解の探索アルゴリ∑ズムについてのみ以下に示す。(S1) は τ Wk,τ Yk,n−τmの m 番目のベクトル要素 S̄k,nを θm に関して∑mm
mm |θ m (S̄log fSk,nk,n |θ ) + log fθ m (θ )(12)k,nを m ごとに最大化することと等しい。よってここより再び特定の m についてのみ考えるので表記の簡単化のため m は省略する。
式 (9) より，log fSk,n (S̄k,n |θ) は |S̄k,n |2 と ϕk,n との[]H
E Xk,n Xk,n|Sk,n = S̄k,n , θ = θ′ =()−1 T ′Λ′k,n − Λ′k,n 1 1T Λ′k,n 1
1 Λk,n +( T ′)−1 ( T ′)−1 T ′2 ′|S̄k,n | Λk,n 1 1 Λk,n 1
1 Λk,n 1
1 Λk,n (16)と書ける。ただし Λ′k,n は Λk,n に θ = θ′ を代入したものを表す。θ に依らない項をまとめて c とすると上式は[]
j j2πk/K 2∑∑Ψi,j)|k,n |A (e′i i,jQ(θ, θ ) = −log Hk Un +Hki Uni,jk,n i,j∑∑−αUni,j ＋ｃ(17)ni,jと書き直せる。ただし，i,ji,j)λ′ k,n ( ′λ′ k,ni,j′ i,jΨk,n = ′ϕ k,n − λ k,n + ′ |S̄k,n |2ϕ k,nϕ k,n板倉斎藤距離と定数項を除いて等しい。第二項の事前確率は，音声要素のアクティベーションを抑制するための(18)

スパース正則化項として利用する。以下では，指数分布

は各要素信号の PSD 推定値であり，E ステップでは θ′(片側 Laplace 分布)の更新に伴いこの値が更新されることになる。

fθ (θ) =∏α exp(−αUni,j )(13)以上より，各パラメータの M ステップ更新則が導ける。Q 関数を Hki と Uni,j に関して偏微分して 0 と置くと，i,j,nを仮定する。この事前確率のスパース化効果は α が大きいほど高くなる。このスパース正則化によって，各フレームの音声スペクトル |S̄k,n |2 の中に高い頻度で生起するスペクトル微細構造ないしスペクトル包絡構造ほど駆動信号 PSD ないし全極型フィルタの推定結果となっHki =1 ∑ ∑ i,j j j2πk/K 2 i,jΨk,n |A (e)| /UnNJ n j

αUni,j 2 + KUni,j −∑ i,j/Ψk,n |Aj (ej2πk/K )|2 Hki = 0,(20)kて表れやすい仕組となる。
式 (12) を最大化する θ は解析的に求めることはできな
1,1I,J T
いが，以下に示すように Xk,n = (Xk,n, · · · , Xk,n) をが得られる。同様に，aj1 , · · · , ajP に関して偏微分して 0と置き，連立させると，Yule-Walker 方程式完全データと見なして Expectation-Maximization (EM)rpj =アルゴリズムにより効率良く局所解を探索することがで加さえ保証されていれば全体のアルゴリズムの収束性は′i,j′′ ′ ′保証される。Xk,nと Xki ′,j,n′ は (k, n, i, j) ̸= (k , n , i , j )

P∑

(21)

を得る。ただし，rpj は，インデックス j の全ての音声要素のスペクトル包絡の平均を取ったものの逆 DFTのとき独立であるため，完全データ Xk,n の対数尤度は

rpj =

具体的に

∑ ∑ Ψi,jk,nk

log fX|θ (X|θ)∑[()]
H=−log det πΛk,n + tr Λ−1k,n Xk,n Xk,nj
ajq rp−q(p = 1, · · · , P )q=1

きる。なお，(S1)〜(S3) の各ステップで目的関数の単調増′

Uni,j ≥ 0(19)i i,jn,i Hk Un

epj2πk/K(22)である。以上より，aj1 , · · · , ajP の更新値は，式 (21) を(14)Levinson-Durbin アルゴリズムで解くことで得られる。
なお，駆動信号 PSD と全極型フィルタとの間にスケーと書ける。上式に対し，S̄k,n および θ = θ′ の条件付期ルの任意性があるため，各反復の H の更新後に Hki ←∑Hki / k′ Hki ′ のような正規化ステップを設けることにすk,n

待値を取り，log fθ (θ) を加えると，Q 関数∑[Q(θ, θ′ ) = log fθ (θ) −log det πΛk,nk,n([H′+ tr Λ−1k,n E Xk,n Xk,n |Sk,n = S̄k,n , θ = θ

日本音響学会講演論文集る。このステップを含めることで厳密にはアルゴリズムの収束性は保証されなくなるが，実験を通してこのこと])]が収束性能に大きく影響を及ぼすような場面は特に見受(15)けられなかった。

- 701 -

2010年9月

表1Source#1#24

入力 SIR(dB)#1-0.6
0.6

Channel#2#3-0.3-0.1
0.30.1

表2#40.6-0.6

Source#1#2

I と J の各設定条件における出力 SIR(dB)

I=2
J = 12J = 15
18.0
16.8
11.7
10.9

実験結果

I=5J =8J = 1019.9
18.5
14.1
13.9

I = 10J =8J = 1217.7
17.6
13.6
14.5

Baseline1
11.6
11.0

Baseline2
17.2
13.9

参考文献

ATR 音声データベースの音声データ (女性話者の音声信号で標本化周波数は 16kHz) を用いて提案法による残響環境下でのブラインド音源分離の実験を行った。混合過程のインパルス応答は，残響時間が 0.5s の可変残響室にて計測したた 2 入力 4 出力のものを用いた。各音源信号の入力信号対妨害比 (Signal-to-Interface ratio, SIR)
は表 1 に示すとおりである。また，音源 1 と音源 2 の直接音対残響音比 (Direct-to-reverberation ratio, DDR)
はそれぞれ 5.7dB，4.7dB であった。マイクロホン入力
信号の時間周波数成分はポリフェーズフィルタバンクにより算出し，STFT で言うところのフレーム長 32ms 相当，フレームシフト 8ms 相当の分析条件とした。分離システムにおけるフィルタ長 nk は，Fk < 0.8 において

nk = 25，0.8 ≤ Fk < 1.5 において nk = 20，1.5 ≤ Fk < 3において nk = 15，Fk ≥ 3 において nk = 10 とした。ただし，Fk は k 番目の周波数ビンの kHz 単位での周波数[1] S. Douglas, H. Sawada and S. Makino, “Natural gradient multichannel blind deconvolution andspeech separation using causal FIR ﬁlters,” IEEE
Trans. Speech, Audio Process., vol. 13, no. 1, pp.
92–104, 2005.
[2] P. Smaragdis, “Blind separation of convolved mixtures in the frequency domain,” Neur. Comp., vol.
22, pp. 21–34, 1998.
[3] T. Nakatani, T. Yoshioka, K. Kinoshita,
M. Miyoshi and B.-H. Juang,“Blind speech
dereverberation with multi-channel linear prediction based on short time Fourier transform
representation,”in Proc. Int’l Conf. Acoust.,
Speech, Signal Process., 2008, pp. 85–88.
[4] T. Yoshioka, T. Nakatani, M. Miyoshi and H. G.
Okuno, “Blind separation and dereverberation ofspeech mixtures by joint optimization,” acceptedfor publication in IEEE Trans. Audio, Speech, Language Process., 2010.

[5] S. Dégerine and A. Zaı̈di, “Separation of an instantaneous mixture of Gaussian autoregressive sources
P は 12，スパース正則化パラメータ α は K/100 とし，
by the exact maximum likelihood approach,” IEEE
(S1) から (S3) の反復回数は 3 回，(S1) の内部で実行すTrans. Signal Processing, vol. 52, no. 6, pp. 1499–
る EM アルゴリズムの反復回数は 100 回とした。
1512, 2004.
提案法を様々な I と J の組み合わせの設定条件で行っ [6] H. Kameoka and K. Kashino, “Composite Autoregressive System for Sparse Source-Filter Represenた際の各々の分離性能の結果を，ベースライン手法によtation of Speech,” In Proc. 2009 IEEE Internaる分離性能の結果と併せて表 2 に示す。Baseline1 およtional Symposium on Circuits and Systems (ISび Baseline2 はそれぞれ澤田らの手法 [10] と吉岡らのCAS2009), pp. 2477-2480, 2009.
手法 [4] をさす。Baseline1 については時間周波数解析 [7] L. Benaroya, F. Bimbot and R. Gribonval, “Auをフレーム長 256ms，フレームシフト 64ms の条件で，
dio source separation with a single sensor,” IEEE
Trans. Audio Speech Language Processing, vol. 14,
Baseline1 についてはフレーム長 16ms，フレームシフト
no. 1, pp. 191–199, 2006.
8ms の条件でそれぞれ行った。実験結果より，提案法に[8] C. Févotte, N. Bertin and J.-L. Durrieu, “Nonnegよる分離精度は I と J が 5 と 8 のとき最も高く，両ベー
ative matrix factorization,with the Itakura-Saito
スライン手法を凌駕する結果を得た。
divergence. With application to music,analysis,”
Neural Comput., vol. 21, no. 3, pp. 793–830, Mar.
2009.
5 まとめ
[9] A. Ozerov and C. Févotte, “Multichannel nonneg本稿では，我々が以前提案した複合自己回帰系と呼ぶ
ative matrix factorization in convolutive mixtures
for audio source separation,” IEEE Trans. Audio,
音声信号の統計モデルを既存の BSS システムに組み込Speech,Language Process., vol. 18, no. 3, pp. 550–
むことで BSS の性能向上を図ることを目的とし，実験
563, March 2010.
によりその一有効性を示すことができた。また，[4] の[10] H. Sawada, S. Araki and S. Makino, “Measuring
手法と同様，残響環境下においても良い分離精度が得ら
dependence of binwise separated signals for perれることが分かった。
mutation alignment in frequency-domain BSS,” inProc. Int’l Symp. Circ., Syst., 2007, pp. 3247–3250.
値を表す。複合自己回帰系における各全極モデルの次数

日本音響学会講演論文集- 702 -2010年9月