上次關於語音增強的原理講說了噪聲估計問題，這次打算說下增益因子如何確定，也就是當噪聲已知后，如何進行去噪的問題（把增益因子與帶噪語音相乘即可）。這里主要說下MMSE濾波，順帶說下譜減法、維納濾波。當然也有其它方式來實現語音增強的，比如基於矩陣分解原理的子空間法、基於自適應濾波器的降噪，有的方法icoolmedia比較清楚，有的也還在學習之中，同時也歡迎各位朋友就不足之處批評指正。
先說下經典的譜減法。我們還是假設帶噪語音y(n)由純凈語音x(n)和加性噪聲d(n)組成，其時域表示與頻域表示為
\[\begin{array}{l}
y(n) = x(n) + d(n) \\
Y(\omega ) = X(\omega ) + D(\omega ) \\
\end{array}\]
在語音增強領域中，最常使用的是頻域功率譜形式，那么，帶噪語音的功率譜可以表示為Y與其共軛相乘，展開可得
\[\begin{array}{l}
|Y(\omega ){|^2} = [X(\omega ) + D(\omega )][{X^*}(\omega ) + {D^*}(\omega )] \\
= |X(\omega ){|^2} + |D(\omega ){|^2} + X(\omega ){D^*}(\omega ) + {X^*}(\omega )D(\omega ) \\
= |X(\omega ){|^2} + |D(\omega ){|^2} + 2{\mathop{\rm Re}\nolimits} \{ X(\omega ){D^*}(\omega )\} \\
\end{array}\]
展開式的第三項被稱為交叉項，當純凈語音與加性噪聲不相關時，交叉項為0，那么，當我們已經估計出噪聲信號的功率譜時，純凈語音信號的估值就可以表示為
\[|\hat X(\omega ){|^2} = |Y(\omega ){|^2} - |\hat D(\omega ){|^2}\]
根據線性濾波理論，可以將這個濾波過程建模為
\[|\hat X(\omega ){|^2} = {H^2}(\omega )|Y(\omega ){|^2}\]
綜合以上兩式，H可以表示為
\[H(\omega ) = \sqrt {\frac{{|\hat X(\omega ){|^2}}}{{|Y(\omega ){|^2}}}} = \sqrt {\frac{{|Y(\omega ){|^2} - |\hat D(\omega ){|^2}}}{{|Y(\omega ){|^2}}}} = \sqrt {1 - \frac{{|\hat D(\omega ){|^2}}}{{|Y(\omega ){|^2}}}} \]
這里的H，就是線性濾波系統的傳遞函數，在語音增強領域，通常也稱為增益函數/抑制函數、或者增益因子/抑制因子，都是是同一個意思。另外，要注意，上面的這個增益因子表示只是一個理想的過程。因為當從帶噪語音中減去估計噪聲后，總會遺留一些或長或短的小譜峰，這些譜峰比較影響聽感。這種現象就是出現了音樂噪聲。因此，如果譜減法要實際使用時，必須做如下改變。
- 當對噪聲估計過高時，就出現了過估計現象，很可能估計出來的噪聲功率大於帶噪語音的功率，這時，不能簡單的把純凈語音的功率置0，而應該設置與噪聲相關的一個譜下限值。設置譜下限的目的在於縮小遺留的小譜峰的差異，控制殘留噪聲的多少和音樂噪聲的大小
- 人為的對噪聲多減去一部分，這樣做的目的是為了盡可能的減少小譜峰的遺留，避免音樂噪聲。
經過這兩方面的改變，譜減法具有如下形式：
\[|\hat X(\omega ){|^2} = \left\{ {\begin{array}{*{20}{c}}
{|Y(\omega ){|^2} - \alpha |\hat D(\omega ){|^2}} \\
{\beta |\hat D(\omega ){|^2}} \\
\end{array}} \right.\]
這里當存在語音時，用第一個式子，當沒有語音存在時，就用下面的式子，其中，alpha就是過減因子，取為一個大於1的值，beta為一個遠小於1的值，具體取值范圍請參考Berouti等人的論文：Enhancement of speech corrupted by acoustic noise，這里不再詳加敘述。
這時，當存在語音時譜減法的增益因子就變為
\[H(\omega ) = \sqrt {\frac{{|Y(\omega ){|^2} - \alpha |\hat D(\omega ){|^2}}}{{|Y(\omega ){|^2}}}} = \sqrt {\frac{{\gamma (\omega ) - \alpha }}{{\gamma (\omega )}}} \]
這里的gamma是后驗信噪比，為帶噪語音與噪聲的功率之比。譜減法增益因子就說完了。頻域維納濾波的增益因子可以參考我以前寫的博客一個頻域語音降噪算法實現及改進方法中的內容，里面有詳細的推導過程，這里就不詳加敘述了。下面重點說下MMSE降噪算法是如何確定增益因子的。
MMSE估計器用在語音增強之中，就是在貝葉斯准則下估計出來的純凈語音頻幅度與實際幅度的均方誤差最小，因此，也可以稱為貝葉斯MSE。而要做到這一點，我們可以充分利用帶噪語音的先驗信息來提高估計的准確性。即，假定我們已知信號的噪聲DFT系數的概率密度，在此基本上，充分利用這種已知的先驗信息，提高估計的准確性。因此，貝葉斯MSE用公式表示如下：
\[Bmse({\hat X_k}) = E[{({X_k} - {\hat X_k})^2}] = \int {\int {{{({X_k} - {{\hat X}_k})}^2}} } p(Y,{X_k})dYd{X_k}\]
我們來推導下使貝葉斯MSE最小的估計量，首先應用貝葉斯原理，聯合概率密度可以寫為：\[p(Y,{X_k}) = p({X_k}|Y)p(Y)\]
所以
\[Bmse({{\hat X}_k}) = \int {\left[ {\int {{{({X_k} - {{\hat X}_k})}^2}} p({X_k}|Y)d{X_k}} \right]} p(Y)dY\]
對中括號中的積分求導
\[\begin{array}{l}
\frac{d}{{d{{\hat X}_k}}}\int {{{({X_k} - {{\hat X}_k})}^2}} p({X_k}|Y)d{X_k} = \int {\frac{d}{{d{{\hat X}_k}}}} {({X_k} - {{\hat X}_k})^2}p({X_k}|Y)d{X_k} \\
= \int { - 2(} {X_k} - {{\hat X}_k})p({X_k}|Y)d{X_k} \\
= - 2\int {{X_k}p({X_k}|Y)d{X_k}} + 2{{\hat X}_k}\int {p({X_k}|Y)d{X_k}} \\
\end{array}\]
令等式等於0，得
\[{{\hat X}_k} = \int {{X_k}p({X_k}|Y)d{X_k}} = E[{X_k}|Y] = E[{X_k}|Y({\omega _0})Y({\omega _1})...Y({\omega _{N - 1}})]\]
在此我們假設傅里葉變換系數之音是統計獨立的。因此上式可以表示為
\[{{\hat X}_k} = E[{X_k}|Y({\omega _0})Y({\omega _1})...Y({\omega _{N - 1}})] = E[{X_k}|Y({\omega _k})] = \int {{X_k}p({X_k}|Y({\omega _k}))d{X_k}} \]
可以看到，要想得到MMSE估計器，我們首先需要計算純凈語音第k個分量的后驗概率密度函數，它可以通過貝葉斯准則得到：
\[p({X_k}|Y) = \frac{{p(Y({\omega _k})|{X_k})p({X_k})}}{{p(Y({\omega _k}))}} = \frac{{p(Y({\omega _k})|{X_k})p({X_k})}}{{\int {p(Y({\omega _k})|{x_k})p({x_x})d{x_k}} }}\]
這里xk是隨機變量Xk的實際值。把上面這個后驗概率密度函數表達式代入我們推導出來的MMSE估計器中
\[{{\hat X}_k} = E[{X_k}|Y({\omega _k})] = \int\limits_0^\infty {{x_k}p({x_k}|Y({\omega _k}))d{x_k}} = \frac{{\int\limits_0^\infty {{x_k}p(Y({\omega _k})|{x_k})p({x_k})} d{x_k}}}{{\int\limits_0^\infty {p(Y({\omega _k})|{x_k})p({x_x})d{x_k}} }} = \frac{{\int\limits_0^\infty {\int\limits_0^{2\pi } {{x_k}p(Y({\omega _k})|{x_k},{\theta _k})p({x_k},{\theta _k})d{\theta _k}d{x_k}} } }}{{\int\limits_0^\infty {\int\limits_0^{2\pi } {p(Y({\omega _k})|{x_k},{\theta _k})p({x_x},{\theta _k})d{\theta _k}d{x_k}} } }}\]
其中
\[\begin{array}{l}
p(Y({\omega _k})|{x_k},{\theta _k}) = \frac{1}{{\pi {\lambda _d}(k)}}\exp \left\{ { - \frac{1}{{{\lambda _d}(k)}}|Y({\omega _k}) - X({\omega _k}){|^2}} \right\} \\
p({x_x},{\theta _k}) = \frac{{{x_k}}}{{\pi {\lambda _k}(k)}}\exp \left\{ { - \frac{{x_k^2}}{{{\lambda _k}(k)}}} \right\} \\
\end{array}\]
代入MMSE估計器中，我們最終得到MMSE幅度譜估計器（推導過程請參考：語音增強-理論與實踐中的附錄B）
\[{{\hat X}_k} = \frac{{\sqrt {{v_k}} }}{{{\gamma _k}}}\Gamma (1.5)\Phi ( - 0.5,1; - {v_k}){Y_k}\]
其中，Γ(.)為伽馬函數，Φ（a,b;c)為合流超幾何函數，ξ為先驗信噪比、最后一個式子為后驗信噪比。
\[\begin{array}{l}
{v_k} = \frac{{{\xi _k}}}{{1 + {\xi _k}}}{\gamma _k} \\
{\xi _k} = \frac{{{\lambda _x}(k)}}{{{\lambda _d}(k)}} \\
{\gamma _k} = \frac{{Y_k^2}}{{{\lambda _d}(k)}} \\
\end{array}\]
最后，把合流超幾何函數寫成貝塞爾函數的形式，我們就得到了最終的MMSE估計器的表達式：
\[{{\hat X}_k} = \frac{{\sqrt \pi }}{2}\frac{{\sqrt {{v_k}} }}{{{\gamma _k}}}\exp \left( { - \frac{{{v_k}}}{2}} \right)\left[ {(1 + {v_k}){I_0}\left( {\frac{{{v_k}}}{2}} \right) + {v_k}{I_1}\frac{{{v_k}}}{2}} \right]{Y_k}\]
如果我們定義：
\[G({\xi _k},{\gamma _k}) = \frac{{{{\hat X}_k}}}{{{Y_k}}} = \frac{{\sqrt \pi }}{2}\frac{{\sqrt {{v_k}} }}{{{\gamma _k}}}\exp \left( { - \frac{{{v_k}}}{2}} \right)\left[ {(1 + {v_k}){I_0}\left( {\frac{{{v_k}}}{2}} \right) + {v_k}{I_1}\frac{{{v_k}}}{2}} \right]\]
的話，這里G就是我們要求的MMSE幅度估計器的增益。
另外想說一下，MMSE估計的推導思路我弄明白了，主要是通過參考《語音增強-理論與實踐》、《統計信號處理基礎-估計與檢測理論》這兩本書做到的，但關於合流超幾何函數與貝塞爾函數的推導內容還沒完全搞明白，如果不是對理論推導過程非常感興趣的話，這里也沒有必要深究，只要會使用這個結果就行了。
使用MMSE做語音增強，經典的出處應該是Speech enhancement using minimum mean-square error這篇論文，但里面講的並不詳細，這里盡可能的給出能讓大家理解流程的推導。當然，如果感興趣的話，icoolmedia還是推薦大家最好都認真看一遍上面提到的資料。