3-8-10
音声のスパース性と非負制約つき畳み込みモデルに基づくパワースペクトル領域残響除去∗○亀岡弘和, 中谷智広, 吉岡拓也 (NTT コミュニケーション科学基礎研究所)

1

序論

室内で発せられた音声信号を，話者から離れた位置にあるマイクロホンで収音する際，マイクロホン収音信号には不可避的に原音声信号の残響成分が混入し，これが原因となって音声の明瞭度が低下する。
収音信号から残響成分だけを除去し，音声の明瞭度を向上する技術は，実環境における音声認識システムの性能を向上する目的，テレビ会議システムなどの通信システムにおいて音声の明瞭度を向上する目的，録音音声を明瞭に再生する目的など幅広い用途に有用である。
音響信号から残響を除去する方法として，これまでさまざまなアプローチが提案されている。マイクロホンアレーを用いたアプローチでは，ターゲット音の到来方向を推定してそれ以外の到来方向からの音を抑制するもの [1, 2] や，事前に測定しておいた室内インパルス応答をもとに観測信号に逆フィルタ処理を施してターゲット音を復元するもの [3] 等がある。
一方，単一マイクロホン入力を対象としたアプロー
チでは，クリーン音声に関する仮定やモデル (調波性，自己回帰モデル，自己相関関数コードブック等) に基づいて，復元音声ができるだけクリーンな音声らしさを有するように室内インパルス応答の逆フィルタを推定するもの [4, 5, 6, 7] が検討されている。一般に室内インパルス応答は音源位置に応じて時々刻々と著しく変化することがあるため，これらのアプロー
チにおいては，短い観測信号からいかに頑健に逆フィルタを構成できるかが重要課題となっている。
これと並行した試みとして，サブバンドごとのパワーエンベロープに逆フィルタ処理を行うアプロー
チも検討されている [8, 9]。このアプローチは，室内インパルス応答のサブバンド信号の中でも音源位置に応じて著しく変化するのは特に位相の方であり，パワーエンベロープに関しては比較的変動しにくいという仮説を基礎としている。音源と室内インパルス応答のサブバンドパワーエンベロープ同士の畳み込みによるモデル化はパワースペクトルの加法性など様々な近似仮定の上に導かれるため，残響除去精度に関してはある程度の限界があることが予想されるが，上記の室内インパルス応答の逆フィルタ処理に基づくアプローチに比べて音源位置などの環境変化に対してある程度頑健に動作する可能性がある。
本稿では，この利点に着目してサブバンドごとのパワーエンベロープから残響成分を除去するアプロー
チに準拠し，非負値行列分解 (NMF)[10] と呼ぶ原理∗をヒントにした新しい解法を提案する。提案法は，音声および室内インパルス応答のサブバンドパワーエ
ンベロープを各離散時刻で自由度をつようにモデル化する ([8, 9] のように特定のパラメトリックな関数を仮定しない) 点，音声のスパース性を尺度に最適化規準を設計する点，両サブバンドパワーエンベロー
プの非負性を保証した乗法更新アルゴリズムと呼ぶ
効率的な反復計算法が導ける点，乗法更新アルゴリズムが FFT(Fast Fourier Transform) を効果的に利用してさらに高効率化ができる点，等が特徴である。

2

提案法の原理2.1

残響音声のサブバンドパワーエンベロープ

クリーン音声および室内インパルス応答の k 番目
(k = 1, · · · , K) のサブバンド信号をそれぞれ sk,t ,hk,t とすると，残響音声のサブバンド信号 xk,t は，各サブバンド信号が STFT 型のフィルタバンクからの出力である場合には，特定の仮定の下で
Xxk,t =sk,τ hk,t−τ(1)τのように近似できる [11]。ただし，t = 1, · · · , T は時刻に対応するインデックスである。式 (1) より，残響音声のサブバンドパワーエンベロープはXX¯¯¯xk,t ¯2 =(2)s∗k,t−τ h∗k,τ sk,t−τ 0 hk,τ 0=XXτ

ττ0¯

¯¯

¯s∗k,t−τ sk,t−τ 0 ¯hk,τ ¯¯hk,τ 0 ¯e−jφk,τ ejφk,τ 0τ0と表される。ただし，hk,τ = |hk,τ |ejφk,τ である。ここで，φk,τ を D = [−π, π) 上の一様分布に従う独立な確率変数とすると，|xk,t |2 の期待値は，£¤
E |xk,t |2ZX¯¯ ¯¯
1 jφk,τ −jφk,τ
¯sk,t−τ ¯2 ¯hk,τ ¯2dφk,τ
e=eD 2πτX¯¯¯¯s∗k,t−τ sk,t−τ 0 ¯hk,τ ¯¯hk,τ 0 ¯+
τ 6=τ 0Z

1 −jφk,τ
dφk,τ
e2πDX¯¯2 ¯¯2¯sk,t−τ ¯ ¯hk,τ ¯=Z

1 jφk,τ 0edφk,τ 02πD(3)τと書ける。以上のように，残響音声のサブバンドパ
ワーエンベロープは，期待値の意味で，クリーン音声と室内インパルス応答のサブバンドパワーエンベロープの畳み込みで表される。

Power spectrum domain dereverberation based on speech sparseness and non-negative convolution model.
by KAMEOKA Hirokazu, NAKATANI Tomohiro and YOSHIOKA Takuya (NTT Communication ScienceLaboratories)

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

2008年9月

2.2

音声スパース性に基づく残響除去の問題設定

表記の簡単のため，以後 Sk,t ≡ |sk,t | , Hk,t ≡|hk,t |2 と置く。2.1 節の議論より，残響音声の k 番目のサブバンドのパワーエンベロープ Xk,t を
XXk,t ≡Sk,τ Hk,t−τ(4)2τのようにモデル化する。ここで，S, H のスケールの任P意性を除くため，便宜的に t Hk,t = 1 を仮定する。
音声には時間周波数領域においてエネルギーがまばらにしか存在しない (スパース性) という性質があるが，式 (4) のモデルにおいて，スパース性の性質をSk,t を推定するための手がかりにすることが提案法の狙いである。具体的には，観測信号のサブバンドパワーエンベロープが Yk,t のとき，Xk,t ' Yk,t となるような非負の Sk,t , Hk,t の中で，できるだけ Sk,t がスパースになるような解を求めるのがここでの目的である。今，Y と X との間に，Yk,t = Xk,t + ²k,t(5)なる関係があるとする。このモデル化誤差 ²k,t は，2.1
節において立てられた近似仮定に起因するあらゆる誤差を含む。この ²k,t を N (0, σ 2 ) に従う Gauss 性
白色雑音と仮定すると，Y ≡ (Yk,t )K×T に対するS ≡ (Sk,t )K×T , H ≡ (Hk,t )K×T の尤度は，
Y 122√
P (Y |S, H, σ 2 ) =e−(Yk,t −Xk,t ) /2σ(6)2πσk,tとなる。ここで S の事前確率 P (S) を一般化正規分布
pYS1− k,tbp¡¢
P (S) =(7)e2Γ 1 + p1 bk,tとし，P (Hk,1 , · · · , Hk,T ) を k ごとに独立な一様分布
P( t Hk,t = 1 の制約があるため，正確にはパラメー
タがすべて 1 の Dirichlet 分布) とすると，事後確率は

P (S, H, σ 2 |Y ) ∝ P (Y |S, H, σ 2 )P (S)(8)と表される。ただし，p, b は一般化正規分布の形状を規定する定数であり，0 < p < 2 のとき P (S) は優ガウス的となり，スパース性を測るための適切な尺度になる [12]。式 (8) より，log P (S, H, σ 2 |Y ) を S, H に関して最大化する問題は，p, σ 2 , b をあらかじめ定めておく定数とすると，
X¡X¢2f (S, H) ≡Yk,t − Xk,t + 2λ|Sk,t |p (9)k,t2.3

非負性を保証する乗法更新アルゴリズム

本節では，非負値行列分解の効率的な解法として知られる乗法更新アルゴリズム [10] をヒントにし，S,H の非負性を保証しながら f (S, H) を反復的に小さくできるアルゴリズムを導出する。
まず，S の乗法更新式を導く。1 ステップ前の S とH の更新値をそれぞれ S 0 , H 0 として，0S 0 k,t Hk,t−τmk,t,τ = X
0
S 0 k,t Hk,t−τ(11)τとすると，f (S, H 0 ) ≤+X³0X (mk,t,τ Yk,t − Sk,t Hk,t−τ)2mk,t,τ

k,t,τ
0pSk,tp−2

20 p0 pSk,t+ 2|Sk,t| − p|Sk,t|´(12)k,tが成り立つ。上不等式右辺を f¯(S) と置く。証明は省略するが，f¯(S) を最小化するように Sk,t を更新すれ∂ f¯(S)ば f (S, H) の非増加性が保証される。 ∂Sk,τ = 0 を解
くと，f (S, H) の非増加性が保証される更新式
X0Hk,t−τYk,t0X
Sk,τ = Sk,τt00Hk,t−τXk,t+ λp|S 0 k,τ |p−1(13)tを得る。ただし，0=Xk,tX

00Hk,t−τSk,τ(14)τである。式 (13) のとおり，S の更新値が 1 ステップ
前の S 0 と更新係数との積となるため，このような形の更新式を乗法更新式という [10]。また，式 (13) より，S 0 および H 0 の要素がすべて非負値であれば更新係数は非負となるため，S の要素はすべて非負値0に更新される (Sk,t= 0 であれば Sk,t = 0 となる)。
次に H の乗法更新式を導く。同様に，1 ステップ
前の S と H の更新値をそれぞれ S 0 , H 0 として，式(11) を用いると，f (S 0 , H) ≤0X (mk,t,τ Yk,t − Sk,tHk,t−τ )2mk,t,τk,t,τ+ 2λk,tX

0 p|Sk,t|(15)k,tを最小化することと等しい。S, H はいずれも非負な量であるので，Sk,τ ≥ 0, Hk,τ ≥ 0 となる制約のもとで f (S, H) を最小化する解を求めたい。ただし，λは p, σ 2 , b に依存して決まる適当な定数であり，モデル化誤差に対するコストと第 2 項のスパース性コストの比重に相当する。以上より，最適化問題は以下のようにまとめられる。

が成り立つ。上不等式右辺を fe(H) と置く。同様に

minimize f (S, H)
X
Hk,t = 1, Hk,t ≥ 0, Sk,t ≥ 0
subject toを 得 る 。た だ し ，上 記 更 新 式 の 導 出 に お い て は ，Pt Hk,t = 1 の拘束は考慮していないので，式 (16)の更新後に規格化する必要がある。

t

日本音響学会講演論文集(10)∂ fe(H)∂Hk,τ = 0 を解くと，乗法更新式

X

t0XHk,τ = Hk,τ

0Sk,t−τ
Yk,t00Sk,t−τ
Xk,t

(16)

t- 706 -2008年9月

2.4

畳み込み NMF との関係

250

200

JX

TX

t→

Ht S

20(21)(22)

tt→
S0 Y T0
Ht = Ht ¯ I ¯ t→
00T

80

100120140160

180200

Frame #

250

Frequency bin #200150100

5020406080100120140160

180200

Frame #(c) 従来法 [11] による残響除去音声
250

Frequency bin #

200

150

100

50

20406080100120140160180200

Frame #合成残響に対する提案法の動作結果例

FFT による乗法更新式の高速計算

式 (13), (16) の更新式には，S 0 と H 0 の畳み込み，H と Y の相関関数，H 0 と X 0 の相関関数，S 0 と Yの相関関数，S 0 と X 0 の相関関数の計算を要するが，これらはすべて FFT を使って非常に効率良く計算することができる。提案法は畳み込み NMF の特殊型であるがゆえに 1 回の更新を式 (22), (23) のとおりに計算するよりはるかに効率的に計算できるのである。

3

動作実験

提案法をさまざまな条件の残響音声に対して適用し動作確認を行った。ここにその動作結果を従来法によるそれとともにいくつかを例示する。以下，従来法は [11] の手法を単一チャネル入力で行ったものをさす。実験データの標本化および周波数分析の条件はTable 1 に示すとおりとした。また，提案法の条件に関しては，各パラメータは

p = 1.2, λ =(23)S Xと書かれる。ただし，(·)p は行列要素ごとに p 乗する演算子とする。また，¯ は Hadamard 積 (行列要素ごとの積をとる演算)，行列の除法は要素ごとの商をとる演算とする。

日本音響学会講演論文集

600は式 (4) の行列表現となっており，式 (17) と対比すると同型のモデルであることが分かる。ただし，式 (17)に対し，Ht が対角行列である点で式 (4) の残響音声モデルは畳み込み NMF モデルの特殊型と言える。
式 (13), (16) の乗法更新式はそれぞれ同様に行列表現にすると，
Xt←
H0t Ytt←
p−1
H0t X0 +λpS040(b) 合成残響音声 (RT=0.5s)

2.5(20)t=1

S = S0 ¯ X

180(d) 提案法による残響除去音声と置くと，

X=

160

100

Fig. 1···140150HK,t
S1,T.. . SK,T12050となる。Wj = (w1,j , w2,j , · · · , wI,j ) とすると，wi,1 ,wi,2 , · · · , wi,J を並べたものが i 番目の 2 次元配列の基底ということになる。
ここで，式 (4) を以上の表記の流儀に従って表す。
Ht を，H1,t , H2,t , · · · , HK,t を対角要素にした行列
H1,t0
..
Ht ≡ (19)

.

···

100200

込み NMF である。ただし， (·) は行列の成分をすべて j − 1 個分右にシフトする演算子とする。例えば，!
!
ÃÃ
2→1→0 1 2 3
1 2 3 4,, A =A =0 5 6 75 6 7 8!
!
ÃÃ
4→
3→
0 0 0 10 0 1 2(18)
, A =
A =
0 0 5 60 0 0 5

S1,1
 ..
S≡ .
SK,180250(17)j→60

Frame #となるように W1 , · · · , WJ と U を求めるのが畳みとし，S を40(a) クリーン音声 (女性話者)

j=10100

20

j→

Wj U

150

50

Frequency bin #Y'Frequency bin #以上のアルゴリズムは畳み込み NMF[13] の特殊型と解釈できる。通常の NMF では 1 次元配列の基底を想定するのに対し，2 次元配列の基底を想定するのが畳み込み NMF の基本的な考え方である。具体的には，観測行列を Y とすると，

XE2, E=Yk,t × 10−8pEk,tに設定し，反復計算回数は 20 とした。
図 1〜4 に，各種実験データに対して提案法および従来法を適用して得られた残響除去音声 (提案法の場
合，Sk,t ) のスペクトログラムを，観測信号のスペクト- 707 -2008年9月

250

200200Frequency bin #Frequency bin #

250

150

100

50150

1005050100150200

250100200300

Frame #

500

600

400

500600200

Frequency bin #

250

200

Frequency bin #400

Frame #
250(a) クリーン音声 (女性話者)

150

150

10050100

100

50200300

Frame #
50100150200

Fig. 3 移動音源 (男性話者) の実環境収録音声 (上段) に250

Frame #

対する提案法の動作結果 (下段)(b) 雑音つき合成残響音声
250250200Frequency bin #Frequency bin #20015010050150

1005050100150200250

Frame #

50150200250300

350400450500

300350400450500250250

Frequency bin #200200Frequency bin #100

Frame #(c) 従来法による残響除去音声

150

100150

10050

5050

100

150200250

Frame #
50

100

150

200

250250

Frame #(c) 提案法による残響除去音声
Frequency bin #Fig. 2200

雑音重畳合成残響に対する提案法の動作結果

150

10050

Table 1 標本化・周波数分析条件

50100

150200

250

300

350

400

450500

Frame #

標本化周波数STFT 条件

フレーム長
フレーム周期
窓関数

8kHz
64ms
32ms
Hanning 窓

Fig. 4 移動音源 (女性話者) の実環境収録音声 (上段) に対する従来法の動作結果 (中段) と提案法の動作結果 (下段)

ログラムと併せて図示する。図 1 は，ATR 単語デー
タベースの女性話者音声に，可変残響室で測定したインパルス応答 (残響時間 (RT) は 0.5s) を畳み込んで作成した合成残響音声データを，図 2 は，同じくATR 単語データベースの女性話者音声に，同様の残響時間のインパルス応答を畳み込んで合成した信号にさらに信号対雑音比が 30dB となるように定常の白色 Gauss 雑音を重畳させた雑音つき残響音声データを，図 3, 4 は，それぞれ，男性話者と女性話者が可変残響室内を移動しながら発した音声を実際に収録した実環境音声データを対象とした結果を示している。

4

参考文献[1] Flanagan et al., J. Acoust. Soc. Am., 78, 1508–1518, 1985.

[2] Schmidt, IEEE Trans. AP, 34(3), 276–280, 1986.
[3] Miyoshi and Kaneda, IEEE Trans. ASSP, 36(2),145–152, 1988.

[4] 中谷他, 信学論, J88-D-II(3), 509–520, 2005.
[5] Kinoshita et al., Proc. ICASSP’06, 1, 817–820,2006.

おわりに

本稿では，音源位置などの環境変化に対して頑健に動作する残響除去法の実現を目的とし，サブバンドごとのパワーエンベロープから残響成分を除去する問題に対し，非負値行列分解 (NMF) の原理をヒントにした新しい解法を提案した。残響音声のサブバ

日本音響学会講演論文集

ンドパワーエンベロープのモデルについての妥当性，音声のスパース性に基づく最適化規準，そのもとで導かれるモデルの制約つき最適化アルゴリズム，等を中心に論じ，さまざまな実験データに対する提案法の動作結果例を示した。

[6] Yoshioka et al., Proc. IWAENC’06, 2006.
[7] Nakatani et al., Proc. ICASSP’07, 193–197, 2007.
[8] 広林他, 信学論, J81-A(10), 1323–1330, 1998.
[9] Unoki et al., ICASSP’03, 1, 840–843, 2003.
[10] Lee and Seung, Nature, 401, 788–791, 1999.
[11] Nakatani et al., Proc. ICASSP’08, 85–88, 2008.
[12] Karvanen et al., Proc. ICA’03, 125–130, 2003.
[13] Smaragdis et al., Proc. ICA’04, 494–499, 2004.

- 708 -2008年9月