2-1-1
バイアス付き MMSE 短時間振幅スペクトル推定法の理論解析およびミュージカルノイズフリー雑音抑圧への拡張∗☆中井駿介, 宮崎亮一，猿渡洋，中村哲 (奈良先端大)，近藤多伸 (ヤマハ)

1

はじめに

ここで，式 (2) 中の変数 ν(f, τ ) は

近年，音声認識システムやテレビ会議システムといった音声通信に関するアプリケーションが増加している．しかし，これらのシステムでは周囲の雑音の影響により収録音声の品質が劣化するという問題がある．この問題を解決するため，周囲の雑音を抑圧し目的音声を強調する音声強調法の開発が重要な課題の一つとして挙げられている．これまで，代表的な非線形音声強調法として，スペクトル減算法 (spectralsubtraction: SS) [1]，ウィーナフィルタ (Wiener filtering: WF) [2] および最小平均二乗誤差規範短時間振幅スペクトル（minimum mean-square error shorttime spectral amplitude: MMSE-STSA）推定器 [3]などが提案されている．しかし，これらの手法ではミュージカルノイズと呼ばれる非線形処理特有の人工的な音の歪みが生じ，致命的な品質劣化を引き起こすという問題がある．
我々は先行研究において，ミュージカルノイズ発生量を定量的に評価するため，高次統計量を用いた評価尺度を提案している [4]．また，この評価尺度に基づいて SS や WF における内部パラメータを設定し，弱い雑音抑圧処理をくり返すことでミュージカルノイズが全く発生しない，ミュージカルノイズフリー雑音抑圧が可能であることを理論解析により明らかにしている [5, 6]．しかし，これまで MMSE-STSA 推定器におけるミュージカルノイズフリー雑音抑圧については検討がなされていない．
そこで，本稿ではミュージカルノイズフリー雑音抑圧を目的として，MMSE-STSA 推定器における事前SNR 推定に改良を加えたバイアス付き MMSE-STSA推定器を提案し，提案手法に関する理論解析を行う．
そして理論解析の結果より，提案手法によりミュージカルノイズフリー雑音抑圧が可能であることを明らかにする．

2

関連研究

2.1

MMSE-STSA 推定器 [3]MMSE-STSA 推定器は，真の音声と推定音声の振幅スペクトルの誤差を最小とする音声強調法であり，以下の式で表される．

ν(f, τ ) =ξ(f, τ )γ(f, τ )1 + ξ(f, τ )(3)である．また，f は周波数ビン，τ は時間フレームイン
デックス，Y (f, τ ) は推定音声，X(f, τ ) = S(f, τ ) +N (f, τ ) は観測信号，S(f, τ ) および N (f, τ ) はそれぞれ真の音声成分および雑音成分である．さらに，
ξ(f, τ ) および γ(f, τ ) はそれぞれ事前 SNR および事後 SNR であり，以下の式で定義される．

ξ(f, τ ) =E[|S(f, τ )|2 ]/E[|N (f, τ )|2 ]

(4)

2(5)

2γ(f, τ ) =|X(f, τ )| /E[|N (f, τ )| ]2

式 (4) および (5) における E[|N (f, τ )| ] は，一般に非音声区間の雑音パワースペクトルの期待値 PN̂ (f )により推定することができる．しかし，真の音声成分である E[|S(f, τ )|2 ] は直接推定することができない．
そのため，一般に，事前 SNR ξ(f, τ ) は以下で定義される判定帰還型推定 (decision directed: DD) 法を用いて逐次推定される．

ˆ τ ) =αγ̂(f, τ − 1)G2 (f, τ − 1)
ξ(f,+(1 − α)max[0, γ̂(f, τ )−1] (0 ≤ α < 1) (6)ここで，α は DD 法のパラメータである忘却係数，γ̂(f, τ ) = |X(f, τ )|2 /PN̂ (f ) は推定事後 SNR である．

2.2

高次統計量に基づくミュージカルノイズ発生量の評価尺度 [4]ミュージカルノイズ発生量を定量的に評価するために，以下で定義されるカートシス比 (kurtosis ratio)
を対象信号のパワースペクトルの非音声区間に適用する [4]．

kurtosis ratio =kurtproc
kurtorg(7)

ここで kurtproc は信号処理後のカートシス，kurtorg
は信号処理前のカートシスである．また，カートシスは以下で定義される．
µ4
kurt = 2(8)
µ∫ 2∞
µm =
xm p(x)dx(9)−∞

Y (f, τ ) =G(f, τ )X(f, τ )
(1)
()
√ √
ν(f,τ))G(f, τ ) = 2π γ(f,τ ) exp − ν(f,τ2[()()]))(1+ν(f, τ )) I0 ν(f,τ+ν(f, τ )I1 ν(f,τ22(2)∗ここで，µm は m 次モーメント，p(x) はパワースペ
クトルドメインの信号 x の確率密度関数（p.d.f.）である．なお，カートシス比の値が 1 の場合は，ミュージカルノイズが全く発生していないことを表し，値が大きいほど信号処理によってミュージカルノイズが発生したことを表す．

Theoretical analysis on biased MMSE short-time spectral amplitude and its extension to musical-noisefree speech enhancement, by Shunsuke Nakai, Ryoichi Miyazaki, Hiroshi Saruwatari, Satoshi Nakamura(Nara Institute of Science and Technology) and Kazunobu Kondo (YAMAHA).

日本音響学会講演論文集- 615 -2013年9月

2.3

ここで，式 (13) 中の変数 Z(m) はミュージカルノイズフリー雑音抑圧 [5]近年，Miyazaki らによって，ミュージカルノイズフリー雑音抑圧が提案されている [5]．ミュージカル
ノイズフリー雑音抑圧は，信号処理前後におけるカートシスの不動条件および雑音抑圧量 (noise reductionrate: NRR) の増加条件を満たすことで達成される．
ここで，NRR は信号処理前後の SNR の差で定義される．

NRR = SNRout − SNRin(10)ここで，SNRin および SNRout はそれぞれ信号処理前後の SNR である．また，非線形音声強調法における雑音の減算量が音声の減算量に比べはるかに大きいと仮定すると，NRR は以下のように表すことができる．
µ1NRR ≃ 10log10 ′(11)µ1
ここで，µ1 および µ′1 はそれぞれ信号処理前後の雑音パワースペクトルの 1 次モーメントである．Miyazakiらは，SS や WF においてカートシス不動条件およびNRR 増加条件を満たすパラメータを数理解析により明らかにしている．

3

バイアス付き MMSE-STSA 推定器

3.1

バイアス付き事前 SNR 推定

MMSE-STSA 推定器では，式 (4) に示す事前 SNR
推定にバイアス値を与えることでミュージカルノイズ発生量を減らすことができる [2]．本稿ではミュージカルノイズ発生量を減らすため，以下に示すバイ
アス付き事前 SNR 推定を導入する．
ˆ τ ) =αγ̂(f, τ − 1)G2 (f, τ − 1)
ξ(f,+(1 − α)max[ε, γ̂(f, τ )−1](12)

ここで，ε は最小許容値とする．以降，式 (12) に示すバイアス付き事前 SNR 推定を導入した MMSE-STSA推定器をバイアス付き MMSE-STSA 推定器と呼ぶ．

3.2

従来の MMSE-STSA 推定器の理論解析 [8]およびその改善

2.2 節および 2.3 節で示したようにカートシス比および NRR の理論値を求めるためには，信号処理前後の雑音パワースペクトルの 1，2 および 4 次モーメン
トが必要となる．そこで，以降では非音声区間における m 次モーメントを導出する．
MMSE-STSA 推定器は DD 法による無限の再帰処理を含んでおり，信号処理後のモーメントを解析的な形で表す事ができない．そこで，Breithaupt らの事
前 SNR の近似 [7] を用いた解析手法が提案されている [8]．文献 [8] では，まず MMSE-STSA 推定器による信号処理後のキュムラントを求め，キュムラントからモーメントへの変換を行う．MMSE-STSA 推定器の出力パワースペクトルのキュムラントは以下の式で与えられる [8]．
(π )m(13)κm[|Y |2 ] = (1 − α) ηθ Z(m)κm[ξml ]4

日本音響学会講演論文集

exp(−mλξ )1 − exp(−mλξ )()1 4λξ =ln·α πZ(m) =(14)(15)である．また，κm[|Y |2 ] は MMSE-STSA 推定器の出力パワースペクトルの m 次キュムラント，κm[ξml ] は
ξml (f, τ ) = max[0, γ̂(f, τ )−1] の m 次キュムラントである．さらに，η および θ はそれぞれ雑音パワースペ
クトルをガンマ分布でモデリングした場合の形状母数および尺度母数である．ガンマ分布は以下の式で与えられる．

p(x) =xη−1 exp(−x/θ)θη Γ(η)(16)

ここで，ξml (f, τ ) は減算係数が 1 の正規化された SSとみなせるため，ξml (f, τ ) のモーメントは容易に求めることができる [9]．文献 [8] では，ξml (f, τ ) のモーメントをキュムラントへ変換し式 (13) に代入することで，MMSE-STSA 推定器の出力パワースペクトルのキュムラントを求める．
この解析手法では，MMSE-STSA 推定器におけるカートシス比および NRR の大まかな傾向を得ることができる．しかし，カートシス比の理論値が実測値と大きく異なるという問題がある．そこで，本稿では式(15) に若干の改善を加える．
()
ρ 4
λξ =ln·(17)α π
ここで，ρ はスムージングパラメータである．ρ は事
前 SNR の近似 [7] おける平滑化の程度を調整するパラメータで，ρ を大きくすることで平滑化を抑えることができる．以降では，実データと整合をとるため
ρ = 2 とし，式 (17) を用いるものとする．

3.3

バイアス付き MMSE-STSA 推定器の理論解析

3.2 節に示した解析手法をバイアス付き MMSESTSA 推定器に適用すると，出力パワースペクトルのキュムラントは以下の式で与えられる．
(π )m
κm[|Ybias |2 ] = (1 − α) ηθ Z(m)κm[ξbias ] (18)4
ここで，κm[|Ybias |2 ] はバイアス付き MMSE-STSA 推定器の出力パワースペクトルの m 次キュムラント，
κm[ξbias ] は ξbias (f, τ ) = max[ε, γ̂(f, τ )−1] の m 次キュムラントである．式 (18) における ξbias (f, τ ) には最小許容値 ε が存在するため，ξml (f, τ ) のように正規化された SS とみなすことはできず，モーメントを求めることができない．そこで，本節では ξbias (f, τ ) のモーメントの導出を行う．
まず ξbias (f, τ ) の p.d.f. について考える．ここでは非音声区間を仮定しているため，γ̂(f, τ ) は正規化された雑音パワースペクトル系列となる．よって，γ̂(f, τ )の p.d.f. は以下の式で表される (Fig. 1 (a) 参照)．

- 616 -p[γ̂] (x) =xη−1 exp(−ηx)(1/η)η Γ(η)(19)

2013年9月

次に γ̂(f, τ ) − 1 の p.d.f. を考えると，ここでの減算はパワーがゼロ方向への平行移動にあたる．よって，
γ̂(f, τ ) − 1 の p.d.f. は以下の式で表される (Fig. 1 (b)参照)．

(x + 1)η−1 exp(−η(x + 1))(1/η)η Γ(η)

p[γ̂−1] (x) =0∫ ε+11(ηx)η−1 exp(−ηx)dx(1/η)Γ(η) 01
Γ1 (η, η(ε + 1))
(22)
=
Γ(η)
∫bとなる．ただし，Γ1 (a, b) = 0 z a−1 exp(−z)dz は第一種不完全ガンマ関数である．よって式 (21) は以下のように書き直すことができる．
(x < ε) 01(x = ε)p[ξbias ] (x) =Γ(η) Γ1 (η, η(ε + 1))
 (x+1)η−1 exp(−η(x+1)) (x > ε)=(1/η)η Γ(η)(23)
式 (9) および (23) より，ξbias (f, τ ) のモーメントは以下の式で与えられる．
∫ ∞µm[ξbias ] =xm p[ξbias ] (x)−∞∫ ∞m (x + 1)=

xε1+εΓ1 (η, η(ε + 1))
Γ(η)
m(24)

ここで，式 (24) の右辺第一項において t = η(x + 1)と変数変換を行うと∫ ∞(x + 1)η−1 exp(−η(x + 1))
xmdx(1/η)η Γ(η)ε∫ ∞
1=(t − η)m tη−1 exp(−t)dtm Γ(η)
ηη(ε+1)1 ∑Γ(m+1)(−η)lmη Γ(η)
Γ(l+1)Γ(m−l+1)
l=0∫ ∞·tη+m−1 exp(−t)dtm=

η(ε+1)1 ∑Γ(m+1)(−η)lmη Γ(η)
Γ(l+1)Γ(m−l+1)m=l=0· Γ2 (η+m−l, η(ε+1))

日本音響学会講演論文集

Fig. 1(25)

P.d.f. after
flooring

P.d.f. after
subtraction

-1-1

ε

P.d.f. deformation.

∫∞となる．ただし，Γ2 (a, b) = b z a−1 exp(−z)dz は第二種不完全ガンマ関数である．また式 (25) の導出には以下に示す二項定理を用いた．
(t + a)m =m∑l=0al

Γ(m + 1)tm−l
Γ(l + 1)Γ(m − l + 1)(26)

ここで式 (24) および (25) で与えられる ξbias (f, τ )のモーメント µm[ξbias ] をキュムラントへ変換すると

κ1[ξbias ] = µ1[ξbias ]κ2[ξbias ] = µ2[ξbias ] − µ1[ξbias ]κ3[ξbias ] = µ3[ξbias ] − 3µ2[ξbias ] µ1[ξbias ] + 2µ32[ξbias ]
κ4[ξbias ] = µ4[ξbias ] − 4µ3[ξbias ] µ1[ξbias ] − 3µ22[ξbias ]+12µ2[ξbias ] µ21[ξbias ] − 6µ41[ξbias ](27)となる．さらに，式 (27) を式 (18) に代入し，キュムラントからモーメントへ変換すると，バイアス付き
MMSE-STSA 推定器の出力パワースペクトルのモーメントが以下の式で得られる．

µ1 = κ1[|Ybias |2 ]µ = κ
2
22[|Ybias |2 ] + κ1[|Ybias |2 ]
µ4 = κ4[|Ybias |2 ] + 4κ3[|Ybias |2 ] κ1[|Ybias |2 ] + 3κ22[|Ybias |2 ]
+6κ2[|Ybias |2 ] κ21[|Ybias |2 ] + κ41[|Ybias |2 ](28)η−1exp(−η(x + 1))dx(1/η)η Γ(η)

(c)

(b)
Normalized
gamma distribution(20)

最後に max [ε, γ̂(f, τ ) − 1] による p.d.f. 変化を考えると，ε より小さい成分が全て ε の位置に積み上げられる (Fig. 1 (c) 参照)．よって，ξbias (f, τ ) の p.d.f. は以下の式で表すことができる．
(x < ε) 0
∫εp(x)dx(x = ε)p[ξbias ] (x) =(21) −1 [γ̂−1]p[γ̂−1] (x)(x > ε)∫ε∫ ε+1
ここで， −1 p[γ̂−1] (x)dx は 0 p[γ̂] (x)dx に等しいため，
∫ ε∫ ε+1p[γ̂−1] (x)dx =p[γ̂] (x)dx−1(a)

4

評価実験

4.1

実験条件

ここでは，バイアス付き MMSE-STSA 推定器においてミュージカルノイズフリー状態が存在することを検証するため，式 (28) を用いてカートシス比および NRR の理論解析を行う．DD 法の忘却係数は 0.98，雑音パワースペクトルの p.d.f. の形状母数 η は白色ガウス雑音を示す 1.0 とした．これらの条件の下，最小許容値 ε を 0 から 1.0 まで変化させた．
また，実際にミュージカルノイズフリー状態が存在することを確認するため，実音源を用いた客観評価実験を行った．音声信号には新聞記事読み上げコーパスを用いた．また，音声信号および白色ガウス雑音信号を等 SNR で混合し，バイアス付き MMSE-STSA推定器を用いて処理を行い，各評価値を算出した．混合信号の先頭 5 秒には雑音のみが存在する非音声区間を設け，非音声区間の雑音パワースペクトルの平均値を推定雑音とした．

- 617 -2013年9月

参考文献

2.0

ε=0Kurtosis ratio

Suppression of acoustic noise in[1] S. F. Boll,speech using spectral subtraction, IEEE Transactions on Acoustics, Speech and Signal Processing, vol.27, no.2, pp.113–120, 1979.

(ConventionalMMSE-STSAestimator)

Musical-noise-free condition

1.5

1.0

[2] P. C. Loizou, Speech Enhancement Theory andPractice, CRC Press, Taylor & Francis Group,FL, 2007.

0.5

ε=1018222026

24

Noise reduction rate [dB]Fig. 2 Theoretical behavior for biased MMSESTSA estimator.
2.5(ConventionalMMSE-STSAestimator)

Kurtosis ratio

2.0

ε=0

[4] Y. Uemura, Y. Takahashi, H. Saruwatari,
K. Shikano, K. Kondo, Automatic optimization scheme of spectral subtraction based onmusical noise assessment via higher-order statistics, Proceedings of International Workshop on
Acoustic Echo and Noise Control, 2008.

Musical-noise-free condition

1.5[3] Y. Ephraim, and D. Malah, Speech enhancement using a minimum mean-square error shorttime spectral amplitude estimator, IEEE Transactions on Acoustics, Speech and Signal Processing, vol.32, no.6, pp.1109–1121, 1984.

1.0

ε=10.510111213141516

Noise reduction rate [dB]Fig. 3 Experimental result for biased MMSESTSA estimator.
4.2

実験結果

理論解析の結果を Fig. 2 に示す．Figure 2 より，εに関してミュージカルノイズ発生量と雑音抑圧性能はトレードオフの関係にあることが分かる．また，正の NRR においてカートシス比が 1 となる点が存在することより，カートシス不動条件および NRR 増加条件を満たす ε が存在することが分かる．すなわち，バイアス付き MMSE-STSA 推定器においてミュージカルノイズフリー雑音抑圧が可能であることを理論的に明らかにした．
次に実音源を用いた客観評価実験の結果を Fig. 3に示す．Figure 3 より，実測値においても ε を適切に設定することで，ミュージカルノイズフリー雑音抑圧を達成できることが確認できた．

5まとめ

本稿ではミュージカルノイズフリー雑音抑圧を目的として，バイアス付き MMSE-STSA 推定器を提案し，その理論解析を行った．理論解析の結果より，バイアス付き MMSE-STSA 推定器にはミュージカルノイズフリー状態が存在することを示した．また実音源を用いた評価実験により，実際にミュージカルノイズフリー雑音抑圧を達成していることを確認した．
謝辞 本研究の一部は，科学技術振興機構・戦略的創造研究推進事業（CREST）の支援を受けた．

日本音響学会講演論文集[5] R. Miyazaki, H. Saruwatari, T. Inoue, Y. Takahashi, K. Shikano, K. Kondo, Musical-noisefree speech enhancement based on optimized
iterative spectral subtraction, IEEE Transactions on Audio, Speech, and Language Processing, vol.20, no.7, pp.2080–2094, 2012.
[6] R. Miyazaki, H. Saruwatari, K. Shikano,K. Kondo, Musical-noise-free speech enhancement based on iterative Wiener filtering, Proceedings of IEEE International Symposium onSignal Processing and Information Technology,
2012.
[7] C. Breithaupt and R. Martin, Analysis of thedecision-directed SNR estimator for speech enhancement with respect to low-SNR and transient conditions, IEEE Transactions on Audio,Speech, and Language Processing, vol.19, no.2,pp.277–289, 2011.
[8] H. Saruwatari, S. Kanehara, R. Miyazaki,
K. Shikano, K. Kondo,
Musical noise analysis for Bayesian minimum mean-square errorspeech amplitude estimators based on higherorder statistics, Proceedings of Interspeech2013,
in press 2013[9] T. Inoue, H. Saruwatari, Y. Takahashi,
K. Shikano, K. Kondo, Theoretical analysis ofmusical noise in generalized spectral subtraction
based on higher-order statistics, IEEE Transactions on Audio, Speech and Language Processing,vol.19, no.6, pp.1770–1779, 2011.

- 618 -2013年9月