2-9-8
複合自己回帰系に基づく音響信号の生成モデル∗○亀岡弘和, 柏野邦夫
(NTT コミュニケーション科学基礎研究所)

1計モデルとしよう。ただし，θ(m) は統計モデルのパラ序論

(m0 )

(m)

観測データがある確率分布に従って生成されたものと見なせるとき，その確率分布の族を一般に統計モデルという。本稿では，音声などを対象とした様々なブラインメータである。まず，Sk,n と Sk0 ,n0 (m 6= m0 ) の独立性
より，Sk,n の確率密度関数は，θ = {θ(m) }m とすると，
Y(m)fS (m) (sk,n ; θ(m) )(3)
fSk,n (sk,n ; θ) =k,nm

ド信号処理の問題への応用を念頭に置いた，音響信号のである。次に，式 (2) より Ŝk,n の確率密度関数は新しい統計モデルについて論じる。
ブラインド音源分離やブラインド残響除去では，音源

fSbk,n (ŝk,n ) = | det Wk |fSk,n (Wk ŝk,n ; θ)の成分と音源からマイクロホンまでの伝達特性がともに未知のもとで，マイクロホン入力信号から音源成分を復元することが目的となる。明らかにこれは不良設定問題であるため，音源に関する何らかの仮定をもとに音源らしさの規準を定め，これにより立てられる目的関数を最(4)(m0 )(m)と表すことができ，さらに Sk,n と Sk0 ,n0 (n 6= n0 ) の独T
T
立性より，Ŝk = (ŜTk,1 , · · · , Ŝk,N ) と置くと，Y
fSbk (ŝk ) =fSbk,n (ŝk,n )(5)n

小 (最大) 化する最適化問題として定式化されることがTTを得る。次に，Yk = (Yk,1, · · · , Yk,N)T と置けば，式多い。例えば，独立成分分析に基づくブラインド音源分(1) は Gk Yk = Ŝk と書けるため Yk の確率密度関数は，離では，音源間の信号の統計的独立性を測るための一つ

fYk (yk ) = | det Gk |fSbk (Gk yk )の指標である信号の非 Gauss 性がしばしば音源らしさを表す規準として用いられる。このように，音源らしさの規準は，どのような信号を生成しやすいか，という何らかの傾向を表した統計モデルを仮定することにより定義できる。どのような統計モデルを仮定すべきかは，それが対象とする音源信号の実際のふるまいや現象にいかに即しているか，数理的技法が馴染みやすい最適化問題で与えられる。ただし，I −Gk,1
 .
 .
 .
Gk = 

−Gk,nlに帰着するか，を考慮して判断することが重要である。
今，M 個の音源の信号を M 個のマイクロホンで収音する状況を考え，これらの収音信号を短時間 Fourier 変換 (STFT) により展開した時間周波数成分を観測データと見なすことにする。ここで，l 番目のマイクロホンで観..

.

..

.

···
..
.

OO..

−Gk,1..
.
−Gk,nl

..

(1)

(M )

応するインデックス) とし，Yk,n = (Yk,n , · · · , Yk,n )Tとする。M 個の音源成分を並べたベクトルを Sk,n =(1)(M )(Sk,n , · · · , Sk,n )T と表すと，Yk,n と Sk,n は近似的に
Yk,n =nlX

Gk,t Yk,n−t + Ŝk,n..

.
···

.
−Gk,1(7)

I(m)

(m0 )

であり，故に det Gk = 1 である。最後に，Sk,n と Sk0 ,n0(k 6= k 0 ) の独立性より，Y1 , · · · , YK の同時密度関数はY
fY1 ,··· ,YK (y1 , · · · , yK ) =fYk (yk )(8)kと書かれる。以上から分かるように，式 (3)〜(8) より，(m)

統計モデル fS (m) (sk,n ; θ(m) ) の具体形を決めれば，自k,n

動的に観測信号 Y1 , · · · , YK の同時密度関数が定義される。この同時密度関数は，未知パラメータ θ, {Wk },(1){Gk,t } の尤度関数に対応するものであり，[1] のように最尤推定の枠組でブラインド音源分離とブラインド残響(2)除去を同時に行うための統一的な目的関数となる。

t=1Ŝk,n = Wk−1 Sk,n.

..



.

(l)

測される信号を Yk,n (ただし，k, n は周波数，時刻に対(6)(m)以上を踏まえ，本稿では，fS (m) (sk,n ; θ(m) ) に該当す
k,nの関係によって結ばれる [1]。ただし，Gk,t は周波数ごる音響信号の新しい統計モデルを提案する。また，提との残響除去フィルタのフィルタ係数行列 (nl はフィル

案する統計モデルに基づき，本報告では，最もシンプタ長)，Wk は周波数領域の瞬時混合に対する分離行列

ルなケースである Gk,1 = · · · = Gk,nl = O, Wk = I,(m)(m)に対応する。さて，ここで，Sk,n を確率変数とし，Sk,n0(m )と Sk0 ,n0 は (k, n, m) 6= (k 0 , n0 , m0 ) のとき独立と仮定す(m)(m)る。そして確率密度関数 fS (m) (sk,n ; θ(m) ) を Sk,n の統k,n∗

M = 1 の場合の θ に関する最適化問題を解く基本アルゴリズムを導く。従って，以後は (m) の表記は無意味なので省略する。

Statistical model for acoustic signals based on composite autoregressive system. by KAMEOKA Hirokazu(NTT Corporation), KASHINO Kunio (NTT Corporation)

日本音響学会講演論文集- 671 -

2009年3月

2

音響信号の統計モデル
音響信号のスペクトルに見られる微細構造と包絡構造は相補的な量であり，前者は，周期性の有無の情報や周{i, j} 番目の要素信号の n 番目のフレームにおけるゲインを Uni,j で表すと，2.1 節での議論のとおり，その STFT
i,ji,j T
KXi,jは，n = (X1,n , · · · , XK,n ) ∈ C

期性を有する場合にはその周期ないし基本周波数の情i,jXi,jn ∼ NC (0, Λn )報を含み，後者は音素や音色に相当する情報を含んでい¡ i,ji,j ¢のように分散共分散行列が Λi,jn = diag λ1,n , · · · , λK,nる。音声や楽音などの意思伝達や表現の媒体として使われる音響信号の多くは，微細構造と包絡構造を個別に時の多次元複素正規分布に従う。ただし，間変化させて多様なスペクトルを形成しながらこれら相らに，音声 (楽音) においては，声の高さ (音階) や音素

2.1

Hki Uni,j|Aj (ej2πk/K )|2(13)Aj (z) = 1 − aj1 z −1 · · · − ajP z −P(14)λi,jk,n =補的な情報を時間に載せて人間に伝達している。またさ(音色) の範囲や種類は限られるため，限られた種類の微細構造ないし包絡構造が個別に明滅しながら時変スペクトルが構成されていると仮定できそうである。
本章では，以上の考えに基づき，音響信号を I 種類のパワースペクトル密度 (以後，PSD) をもつ Gauss 性雑音の駆動信号と J 種類の全極型フィルタからなる複合系から生成されたものと捉えてモデル化する。

i0 ,j 0である。ここで，Xi,jn と Xn0i

Φn =jXX準備

x[t] =PX

Λi,jn(16)jが言える。式 (16) は，正規分布する独立な確率変数の和は各々の正規分布の分散の和を分散とする正規分布にap x[t − p] + ²[t](9)従うことを意味する。複合系のパラメータをまとめてp=1[ ©

θ=からの標本値系列と仮定する。ここで，²[t] は，自己相関関数が h[t] の定常な Gauss 性雑音とし，必ずしも白色雑音とは限らない点を強調しておく。さて，

x[1], · · · , x[K] の離散 Fourier 変換 (以後，DFT) を X =(X1 , · · · , XK )T ∈ CK とすると，式 (9)，DFT の線形性，および，²[t] の定常性と Gauss 性により，X は，平¡¢均が 0，分散共分散行列が Λ = diag λ1 , · · · , λK の多次元複素正規分布 NC (0, Λ) に従う。ただし，Hk,|A(ej2πk/K )|2(10)A(z) = 1 − a1 z −1 · · · − aP z −P

(11)

λk =

は (i, j) 6= (i0 , j 0 ) またはI,Jn 6= n0 のとき独立と仮定すると，X1,1n , · · · , Xn の和
信号の STFT Sn ∈ CK はやはり正規分布に従い，XX¡¢Sn =Xi,j(15)n ∼ NC 0, Φni

まず，ある短時間信号 {x[t]} を，P 次の自己回帰過程(12)

Hki , ajp , Uni,jª(17)k,n,i,jとすると，以上より，Sk,n の確率密度関数が具体的にµ¶
1|sk,n |2fSk,n (sk,n ; θ) =exp −(18)πφk,nφk,nと定義されたことになる。ここで，Φn の対角要素の φk,nは，複合系から生成される不規則信号の PSD であり，

φk,n =XXij

Hki Uni,j|Aj (ej2πk/K )|2(19)で与えられる。以上より，式 (18) が複合自己回帰系に基づく音響信号の統計モデルを与える。

である。H1 , · · · , HK は h[1], · · · , h[K] の DFT，すなわ
ち駆動信号 ²[t] の PSD であり，スペクトル微細構造に対応する。一方で，1/|A(ej2πk/K )|2 は全極型伝達関数の PSD であり，スペクトル包絡に対応する。

2.23最適化アルゴリズム

3.1

問題設定

ここでは，1 章で示した問題において，Gk,1 = · · · =

複合自己回帰系

ここで，上述のとおり，I 種類の駆動信号 PSD と J

Gk,nl = O, Wk = I, M = 1 の場合，すなわち，2 章で仮定した統計モデル fSk,n (sk,n ; θ) の実現値が直接観測種類の全極型フィルタにより構成される複合系によりされるケースを考える。従って，観測信号 y1,1 , · · · , yK,N

音響信号の STFT(フレーム内の信号の DFT) をモデルが与えられたとき，化する。この複合系は，異なる PSD をもつ合計 I × J

log

種類の要素信号をフレームごとに個別のゲインによっfSk,n (yk,n ; θ) + log fθ (θ)(20)

k,nてアクティベートし，それらを重畳したものを生成する機構をもつ。さて，i 番目の駆動信号 PSD と j 番目

Yを最大化する θ を求める最大事後確率推定 (MAP) 問題の全極型伝達関数をそれぞれ Hki , 1/Aj (ej2πk/K ) とし， を考える。式 (18) より，log fSk,n (yk,n ; θ) は観測パワー

日本音響学会講演論文集- 672 -2009年3月

スペクトル |yk,n |2 と PSD モデル φk,n との板倉斎藤距離と定数項を除いて等しい。一方，第二項の事前確率は，同時にアクティベートできる要素信号の個数に関し，観測データに依らない何らかの傾向を統計モデルに持たせる用途に用いる。ここでは，逆ガンマ分布µ¶Y 1βfθ (θ) ∝exp−i,j αUni,ji,j,n Unと書き直せる。ただし，"¢#i,j ¡λ0 k,n |yk,n |2 − φ0k,ni,j0 i,jΨk,n = λ k,n 1 +φ02k,nは各要素信号の PSD 推定値であり，E ステップでは θ0の更新に伴いこの値が更新されることになる。

(21)

以上より，各パラメータの M ステップ更新則が導ける。Q 関数を Hki と Uni,j に関して偏微分して 0 と置くと，を仮定する。この事前確率はアクティベーションをスi,j)|
1 X X Ψk,n |A (ei,jNJ n jUn"#j j2πk/K 2X Ψi,j)|
1k,n |A (ei,jβ+Un =K +αHki

パース化する効果をもち，その効果は α が大きいほど高として扱いたい場面では α を大きく取れば良く，逆に，音楽のように複数の音が同時に発音しうる音響信号の統計モデルとして扱いたい場面では小さく取れば良い。

EM アルゴリズム立であるため，完全データ Ŷ = (Ŷ1,1 T , · · · , ŶK,N T )Tの対数尤度は具体的に
iXh−1HΛŷlog fY(ŷ;θ)
=−log
detπΛ
+ŷk,nbk,n k,n k,nk,n¡¢i
Hlog det πΛk,n + tr Λ−1k,n ŷk,n ŷk,n(26)(27)が得られる。同様に，aj1 , · · · , ajP に関して偏微分して 0

rpj =jajq rp−q(p = 1, · · · , P )(28)を得る。ただし，rpj は，rpj =#"X X X Ψi,jk,nepj2πk/Ki U i,jHnknik(29)である。以上より，aj1 , · · · , ajP の更新値は

 j  ja1r0
 .   .
 . = .
 .   .
ajPrPj −1と書ける。上式に対し，Yk,n = yk,n および θ = θ0 のときの条件付期待値を取り，log fθ (θ) を加えると，Q 関数

XhPX
q=1k,nQ(θ, θ0 ) = log fθ (θ) −2と置き，連立させると，Yule-Walker 方程式

1,1I,J T(Xk,n, · · · , Xk,n) を完全データと見なして EM アルゴ
リズム [2] により解くことができる。
0 0i,j00 0 0
Xk,nと Xki 0,j,n0 は (k, n, i, j) 6= (k , n , i , j ) のとき独=−j2πk/Kk

以上の MAP 推定の局所最適解は，要素信号 Ŷk,n =XhjHki =くなる。従って，fSk,n (sk,n ; θ) を単一音声の統計モデル

3.2(25)···
..
.
···

−1  j jr1−Pr1 ..   .. .   . r0jrPj(30)を解くことで得られる。これは，Levinson-Durbin アルlog det πΛk,nk,n¡£¤¢iH0+ tr Λ−1Eŷŷ|Y=y;θ(22)k,nk,nk,nk,nk,nを得る。ここで，不完全データ (観測データ)Yk,n と完全£¤
データ Ŷk,n に Yk,n = CŶk,n (ただし，C = 1, · · · , 1 )ゴリズムにより高速に計算できる。

4提案法の主要な特徴
1. P = 0, J = 1, fθ (θ) = const. において，P =0 ⇒|Aj (ej2πk/K )|2 = 1 より式 (26), (27) はなる関係があるため，

i,1

Hki =

¤
£
H|Yk,n = yk,n ; θ0 =(23)E ŷk,n ŷk,n¡¢−1Λ0k,n − Λ0k,n CT CΛ0k,n CT
CΛ0k,n +¡¢−1 ¡¢−1CΛ0k,n|yk,n |2 Λ0k,n CT CΛ0k,n CT
CΛ0k,n CT(31)kと書ける。これは，板倉斎藤距離最小化規範の非負値行列因子分解 (NMF) の乗法更新則 [3] と等価である。これよりアルゴリズム面で提案法は NMF とである。ただし，Λ0k,n は Λk,n に θ = θ0 を代入したものを表す。θ に依らない項をまとめて c とすると上式は"#j j2πk/K 2XXΨi,j)|k,n |A (e0i i,jQ(θ, θ ) = −log Hk Un +Hki Uni,jk,n i,j¸XX·β−
α log Uni,j + i,j + c(24)Unn i,ji,1

1 X Ψk,n
1 X Ψk,n, Uni,1 =i,1N n UnKHki

関係が深いことが示唆される。

2. 白色雑音を駆動とする通常の自己回帰モデルは，基本周波数が高い音声には適合しないことが知られる
が，これは駆動信号の白色性の仮定からの乖離が顕著になるためである。一方，提案モデルでは，駆動信号の白色性は仮定せず，その自己相関関数 (または PSD) 自体が推定すべきパラメータとなる。

日本音響学会講演論文集- 673 -2009年3月

Frequency bin #2502001501005050100150200250300350400450

Frame #

Frequency bin #250200150100 (TGSWGPE[KPFGZ

50

図250100150200250300350400





(TGSWGPE[KPFGZ

推定した駆動信号 PSD(左) と全極型フィルタ (右)450

Frame #

図1

表1

観測信号 (上段) と推定した PSD モデル (下段)

3. [1] で用いられる時変全極モデルは自己回帰パラメー
タがフレームごとに自由度をもつが，提案モデルでは全フレームにわたって高々J 組の自己回帰パラメータしか自由度をもたない。また，駆動信号のPSD パラメータの自由度も全フレームにわたって高々I である。これらの拘束が，系と入力の両特性が未知のもとで系を同定する手がかりを与える。

絡をフラットにしたものを，スケールを等化してから平均化する操作となっており，このことから
i,JΨi,1k,1 , · · · , Ψk,N に含まれる共通の微細構造を抽出しようとする働きが理解される。

I \J

1355

3.07

3.54

4.91

103.78

4.63

5.81

15

5.18

6.32

6.65

まざまな I と J の条件のもとで観測信号のパワースペクトルと推定した PSD モデルとの間の一致度を信号対雑音比 (SNR)

4. 式 (26) は，E ステップで更新された要素信号の PSD
i,J
推定値 Ψi,1k,n , · · · , Ψk,n を前ステップで更新された
全極型伝達関数の PSD で除算し，スペクトル包I, J の各条件でのモデル化誤差

PSNR(dB) = 10 log10 Pk,n2k,n |yk,n |¯¯p¯|yk,n | − φk,n ¯2により測定した。表 1 はこの結果を示したものである。

6まとめ
本稿では，音声などを対象とした様々なブラインド信号処理の問題への応用を念頭に置き，I 種類の駆動信号

5. 式 (29) は，E ステップで更新された要素信号の PSD の PSD と J 種類の全極型フィルタによって構成されるI,j1I
推定値 Ψ1,j複合系 (複合自己回帰系) に基づく音響信号の新しい統k,n , · · · , Ψk,n を Hk , · · · , Hk で除算し，微細構造をフラットにしたものを平均化した PSD を 計モデルを提案した。この統計モデルの実現値が直接観逆 Fourier 変換する操作となっている。式 (30) は， 測できる状況において，観測データが与えられたもとでこの結果求まる自己相関関数をもつ仮想的な信号に のモデルパラメータの最大事後確率推定行う基本アルゴ
対し，従来の全極型の最尤スペクトル法を適用して リズムを EM アルゴリズムにより導いた。このアルゴいることになる。
リズムは，板倉齋藤距離を規準とした非負値行列因子分5解と同形の駆動信号 PSD の更新則と，Yule-Walker 方動作実験

程式に基づく自己回帰パラメータの更新則からなる。

ATR 音声データベースの音声データ (女性話者) を用いて提案法の基本動作の確認を行った。yk,n は STFT(標今後は，提案モデルをブラインド音源分離，ブラインド残響除去の問題に本格適用していく予定である。

本化周波数 16kHz， フレーム長 64ms，フレーム周期

参考文献32ms，Hanning 窓) により計算した。
図 2 は，図 1 上段の観測信号に対して I = 10, J = 5の条件のもとで推定した駆動信号 PSD と全極型フィルタを示したものである。また，図 1 下段は推定した PSD
モデル φk,n を図示したものである。図 2 を見ると，Hkiの推定値には調波構造らしい特徴的な構造が表れていることが確認できる。提案する統計モデルがどの程度音声信号に合致したものとなっているかを確認するため，さ日本音響学会講演論文集[1] 吉岡, 中谷, 三好, ブラインド音源分離と残響除去の統合のための一手法, 音講論 (秋), 703–704, 2008.
[2] M. Feder, E. Weinstein, Parameter estimationof superimposed signals using the EM algorithm,
IEEE Trans. ASSP, 36(4), pp. 477–489, 1988.
[3] C. Févotte, N. Bertin, J.-L. Durrieu, Nonnegative
matrix factorization with the Itakura-Saito divergence with application to music analysis, Tech.
Rep. TELECOM ParisTech 2008D006, 2008.

- 674 -2009年3月