语音增强原理之增益因子
上次关于语音增强的原理讲说了噪声估计问题，这次打算说下增益因子如何确定，也就是当噪声已知后，如何进行去噪的问题（把增益因子与带噪语音相乘即可）。这里主要说下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还是推荐大家最好都认真看一遍上面提到的资料。
