在大多数语音通信应用中，系统的输入端只能接收到被背景噪声干扰后的含噪语音，噪 声极大地干扰了语音通信的质量，降低了语音的清晰度和可懂性，对系统中编、解码等语音 处理模块产生不利的影响。
背景噪声抑制技术能从含噪语音中提取尽可能纯净的原始语音，此研究在语音处理领域 中属于“语音增强”范畴。噪声抑制有助于改善含噪语音信号的知觉质量，提高通信环境的 舒适感和服务质量；同时，噪声抑制作为语音编码过程的前处理模块能提高声码器在噪声环 境下压缩性能和稳定性；另外，此技术能有效提高语音识别系统在背景噪声环境下的稳健性。
大多数通信应用的语音输入模型都具有单通道语音输入和加性背景噪声的特点。在这种 输入模型中，观测到的含噪信号在时域或频域都可以表示成语音和噪声分量之和。在目前已 有的噪声抑制方法中，短时谱加权方法是最为主流的技术，此类方法将含噪语音做短时傅立 叶变换，根据各频率上语音和噪声分量设计增益系数，通过将此系数与含噪信号相乘而得到 抑制效果，在短时反傅立叶变换后得到处理后的语音。
谱相减方法(S.F.Boll，Suppression of acoustic noise in speech using spectral subtraction， IEEE Trans.ASSP，vol.ASSP-27，pp.113-120，Apr.，1979)是最典型的例子，其基本原理是在频 域中从含噪语音幅度中减去估计噪声幅度而相位不变从而实现抑制过程，由于此过程以谱加 权的方式实现，因而加权增益系数与信噪比关联。很多方法都遵循了这一思路，但在增益的 计算方法和噪声估计方法上各具特点。
如专利(Method and Apparatus for Suppression Noise in a Communication System，US patent 5,695,622)，该方法利用信噪比的修正指数和噪声能量来计算增益，并由各频带长时平均功率 谱来计算谱偏差从而进行噪声功率估计。
专利(Method and Device for Speech Enhancement in the Presence of Background Noise， WO2005064595)进一步对增益系数进行时域的平滑处理。
专利(Low frequency spectral enhancement system and method，US patent 6,233,549)，该方 法则在计算谱增益系数时强调对低频分量的增强。
一些谱加权方法将抑制噪声理解为估计原始语音的谱幅度并得到了更好的效果，常用的 估计准则包括最大似然ML(R.McAulay，Speech enhancement using a soft-decision noise suppression filter，IEEE Trans.A.S.S.P.，28，1980)，最小均方误差MMSE(Y.Ephraim，Speech enhancement using a minimum mean-square error short-time spectral amplitude estimator，IEEE Trans.A.S.S.P，32，1984)等。
其中MMSE估计方法最为常用，得到了不断改进，典型的如Y.Ephraim(Y.Ephraim，Speech enhancement using a minimum mean-square error log-spectral amplitude estimator，IEEE Trans. A.S.S.P，33，1985)使用MMSE准则估计所有频率上语音谱幅度的对数值，由此计算的谱加权 增益能得到更好效果。
专利(Core Estimator and Adaptive Gains from Signal To Noise Ratio in a Hybrid Speech Enhancement System，US patent 2002002455)采用软判决方式并考虑估计误差最小化和抑制带 来的语音失真之间的权衡，从而最终得到加权增益的计算方法。
MMSE方法的关键问题在于采用何种统计模型来反映频域语音信号的统计分布，现有方 法常用的模型是高斯模型，但它并不能很好地模拟真实情况；另一方面，噪声估计的准确性 和稳健度也是决定抑制方法性能的关键因素。

本发明的目的在于提供一种基于多统计模型和最小均方误差的背景噪声抑制方法，以很 好的模拟真实情况，提高噪声估计的准确性和稳健性。
本发明具体是这样实现的：
一种基于多统计模型和最小均方误差的背景噪声抑制方法，包括以下步骤：
步骤1、对当前输入帧的语音信号进行短时傅立叶变换；
步骤2、利用上帧保留的各频率上纯语音幅度方差估计和噪声幅度方差估计，以及当前 输入帧中语音信号各频率分量，计算当前输入帧中语音信号各频率分量的实部和虚部的估计；
步骤3、计算当前输入帧各频率分量的先验语音不存在概率，据此进一步修正步骤2计 算的实部和虚部的估计；
步骤4、根据当前输入帧的修正后的实部和虚部的估计，计算纯语音幅度方差估计并保 留给下一帧使用；
步骤5、计算当前输入帧的似然比，判决当前输入帧是否为纯噪声帧，如是，则更新噪 声幅度方差估计；
步骤6、采用短时傅立叶反变换和叠接相加得到噪声抑制后的语音。
所述步骤2进一步包括：
分别采用拉普拉斯分布和伽马分布来模拟噪声和语音频谱分量的实部和虚部概率密度分 布；
利用实部和虚部概率密度分布根据最小均方误差准则分别计算实部和虚部的条件期望 值，作为实部和虚部的估计。
所述噪声和语音频谱分量的实部和虚部概率密度分布如下式表示：
$\left\{\begin{array}{c}p\left(N_R\right)=\frac{1}{{\lambda }_n}\exp (-\frac{2\left|N_R\right|}{{\lambda }_n})\\ p\left(N_I\right)=\frac{1}{{\lambda }_n}\exp (-\frac{2\left|N_I\right|}{{\lambda }_n})\end{array}\right.$
$\left\{\begin{array}{c}p\left(S_R\right)=\frac{\sqrt[4]{3}}{2\sqrt{{\mathrm{\pi \lambda }}_s}\sqrt[4]{2}}{\left|S_R\right|}^{-\frac{1}{2}}\exp (-\frac{\sqrt{3}\left|S_R\right|}{\sqrt{2}{\lambda }_s})\\ p\left(S_I\right)=\frac{\sqrt[4]{3}}{2\sqrt{{\mathrm{\pi \lambda }}_s}\sqrt[4]{2}}{\left|S_I\right|}^{-\frac{1}{2}}\exp (-\frac{\sqrt{3}\left|S_I\right|}{\sqrt{2}{\lambda }_s})\end{array}\right.$
上式中NR、NI、SR和SI分别代表噪声、语音频谱分量的实部和虚部，λn和λs分别代表 噪声和语音频谱分量的方差。
所述实部条件的期望值为：
$E\left[S_R(k,l)\right|Y_R(k,l)]=\frac{\sqrt[4]{1.5}}{2{\lambda }_n(k,l)\sqrt{{\mathrm{\pi \lambda }}_s(k,l)}p\left(Y_R(k,l)\right)}{\int }_{-\infty }^{\infty }{S}_R(k,l)·{\left|S_R(k,l)\right|}^{-0.5}.$
$\exp [-\frac{2|Y_R(k,l)-S_R(k,l)|}{{\lambda }_n(k,l)}]·\exp [-\frac{\sqrt{1.5}\left|S_R(k,l)\right|}{{\lambda }_s(k,l)}]{\mathrm{dS}}_R(k,l)$
$=\frac{\sqrt[4]{1.5}}{2{\lambda }_n(k,l)\sqrt{{\mathrm{\pi \lambda }}_s(k,l)}p\left[Y_R(k,l)\right]}\left\{\frac{2}{3}\exp \right[\frac{-2Y_R(k,l)}{{\lambda }_n(k,l)}\left]Y_R^{3/2}(k,l)\Phi \right[\mathrm{1.5,2.5},-2Y_R(k,l)G_1]$
$+\exp \left[\frac{-\sqrt{1.5}{Y}_R(k,l)}{{\lambda }_s(k,l)}\right]{\left({2G}_2\right)}^{-1.5}\psi [-0.5,-\mathrm{0.5,2}{G}_2{Y}_R(k,l)]$
$-0.856·\exp \left[\frac{-2Y_R(k,l)}{{\lambda }_n(k,l)}\right]{\left({2G}_2\right)}^{-1.5}\}$
上式中YR(k，l)表示含噪信号Y在l时刻的第k个频率分量的实部，Ф(a，b，z)＝M(a，b，z)表 示第一类合流超几何函数，Ψ(a，b，z)同样表示超几何函数，可由Ф(a，b，z)计算出来，
$G_1=\frac{\sqrt{1.5}{\lambda }_n+2{\lambda }_s}{2{\lambda }_s{\lambda }_n}$ 和 $G_2=\frac{\sqrt{1.5}{\lambda }_n-2{\lambda }_s}{2{\lambda }_s{\lambda }_n};$
所述虚部条件的期望值：
$E\left[S_I(k,l)\right|Y_I(k,l)]=\frac{\sqrt[4]{1.5}}{2{\lambda }_n(k,l)\sqrt{{\mathrm{\pi \lambda }}_s(k,l)}p\left(Y_I(k,l)\right)}{\int }_{-\infty }^{\infty }{S}_I(k,l)·{\left|S_I(k,l)\right|}^{-0.5}.$
$\exp [-\frac{2|Y_I(k,l)-S_I(k,l)|}{{\lambda }_n(k,l)}]·\exp [-\frac{\sqrt{1.5}\left|S_I(k,l)\right|}{{\lambda }_s(k,l)}]{\mathrm{dS}}_I(k,l)$
$=\frac{\sqrt[4]{1.5}}{2{\lambda }_n(k,l)\sqrt{{\mathrm{\pi \lambda }}_s(k,l)}p\left[Y_I(k,l)\right]}\left\{\frac{2}{3}\exp \right[\frac{-2Y_I(k,l)}{{\lambda }_n(k,l)}\left]Y_I^{3/2}(k,l)\Phi \right[\mathrm{1.5,2.5},-2Y_I(k,l)G_1]$
$+\exp \left[\frac{-\sqrt{1.5}{Y}_I(k,l)}{{\lambda }_s(k,l)}\right]{\left({2G}_2\right)}^{-1.5}\psi [-0.5,-\mathrm{0.5,2}{G}_2{Y}_I(k,l)]$
$-0.856·\exp \left[\frac{-2Y_I(k,l)}{{\lambda }_n(k,l)}\right]{\left({2G}_2\right)}^{-1.5}\}$
上式中YI(k，l)表示含噪信号Y在l时刻的第k个频率分量的虚部。
所述步骤3中进一步包括：
计算当前输入帧平方幅度的和，计算信噪比后，进行时域递归平滑，然后计算全局概率， 为每个频率计算信噪比后，进行时域递归平滑，然后计算局部概率，所述先验语音不存在概 率等于1与全局概率与局部概率乘积的差；
利用先验语音不存在概率和语音存在不确定假设修正实部和虚部的估计。
所述修正的实部估计为：
$E\left[S_R\right|Y_R]=\frac{\Gamma (k,l)}{1+\Gamma (k,l)}E[S_R|Y_R,H_1]=$
$\frac{\sqrt[4]{1.5}}{2{\lambda }_n(k,l)\sqrt{{\mathrm{\pi \lambda }}_s(k,l)}}\{2·\exp [-\frac{2Y_R(k,l)}{{\lambda }_n(k,l)}\left]\sqrt{Y_R(k,l)}\Phi \right[\mathrm{0.5,1.5},-2G_1{Y}_R(k,l)]$
$+\exp [-\frac{1.5}{{\lambda }_s(k,l)}{Y}_R(k,l)]\frac{1}{\sqrt{2G_2}}\Psi \left[\mathrm{0.5,0.5,2}{G}_2{Y}_R(k,l)\right]$
$+\exp [-\frac{{2Y}_R(k,l)}{{\lambda }_n(k,l)}]\sqrt{\frac{\pi }{2G_2}}\}·\frac{\Gamma (k,l)}{1+\Gamma (k,l)}$
上式中 $\Gamma (k,l)=\frac{p\left(Y_R(k,l)\right|H_1)}{p\left(Y_R(k,l)\right|H_0)}·\frac{[1-P\left(H_0\right)]}{P\left(H_0\right)},$ p(YR(k，l)|H1)表示语音存在时的YR(k，l)概 率密度分布，而p(YR(k，l)|H0)则表示只有噪声时的概率密度；
所述修正的虚部估计为：
$E\left[S_R\right|Y_R]=\frac{\Gamma (k,l)}{1+\Gamma (k,l)}E[S_R|Y_R,H_1].$
所述步骤1的短时傅立叶变换后还包括：
检测步骤，对输入的一个或者多个单频音混合的信号音，计算所有频率分量的能量和， 而后求出2个平方幅度的最大值，并从能量和中减取最大值的和，若最大值和大于能量和且 彼此接近，则判决当前的输入为信号音调，不进行任何抑制处理。
本发明采用多个统计模型来分别拟合语音和噪声频域分量的统计分布，因而能更准确地 逼近实际应用中语音和噪声的真实分布；由于考虑了含噪语音存在不确定性对抑制过程的影 响，能获得更高的抑制效果；采用了最大似然比方法进行VAD检测，据此进行噪声功率谱估 计，估计过程更为准确、稳健；采用全局加局部方法来进行先验语音不存在概率的估计；采 用特殊流程来避免对单频和多频音信号的不良影响，不会影响DTMF、传真音调的检测，具 有高的噪声抑制效率和较低的计算复杂度，适合用于各种语音通信系统。
附图说明
图1是本发明所述方法的框架图；
图2是本发明所述方法中运动检测步骤的流程图。

大多数通信应用的语音输入模型都具有单通道语音输入和加性背景噪声的特点，本发明 涉及此模型下的噪声抑制问题。针对噪声抑制问题，本发明提出一个基于多统计模型的自适 应滤波方法。图1所示的是整个方法的框架原理。本发明使用短时傅立叶变换将输入信号变 换到频域，而后利用上一帧获取的参数计算当前输入帧中语音信号各频率分量的实部和虚部 的估计，然后计算语音存在概率并修正语音信号估计，在更新当前参数估计后，利用短时傅 立叶反变换得到抑制后的语音。
具体实施方式的步骤依次为以下七个小节：
1.语音信号的时-频分析
语音和背景噪声都具有高度非平稳特点，单一傅立叶变换不能反映信号随时间变化的频 谱信息，如语音的时变共振峰以及噪声的功率谱等，因而所有的语音噪声抑制必须采用时频 分析方法。短时傅立叶变换(Short-time Fourier，STFT)是最重要的时频分析方法。
STFT过程首先采用分析窗(analysis window)对当前语音数据加权，分析窗函数只在其 支撑内不为0。本发明中，分析窗函数支撑为L时，语音帧长度l为L的25％。STFT对窗加 权后的语音数据进行离散傅立叶变换。比较L和l可见，本发明中相邻分析窗有3/4重叠，此 过程如图2所示。STFT过程如式(1)，其中N即为分析窗长度，而w(n)为分析窗函数。
$Y(k,l)=\underset{n=0}{\overset{N-1}{\Sigma }}y(n+\mathrm{lN})w\left(n\right)\exp [-j\left(\frac{2\pi }{N}\right)\mathrm{nk}]---\left(1\right)$
对Y(k，l)进行降噪处理后得到 进行短时傅立叶反变换(STIFT)和叠接相加(OLA) 方法可得到处理后的语音信号 由于语音降噪中谱加权系数是时变的，因而STIFT必须 采用与h(n)双正交的合成窗(synthesis window)(J.Wexler，Discrete Gabor expansions，Signal Processing，Nov，1990)。
2.语音谱幅度系数的统计模型
目前基于MMSE估计的谱幅度加权抑制方法基本上都使用高斯分布(Gaussian Distribution)来建立噪声和语音的各频谱分量的概率分布模型，此模型主要优势在于数学处 理的方便，实际上并不能准确描述语音和噪声频谱分量的分布。
若假定语音频谱分量的实、虚部服从高斯分布，采用MMSE作为估计准则会得到一个线 性估计子，且此估计子是一个实值的滤波器，也就是说，它得到的MMSE相位估计实际上依 然等于含噪语音的相应频谱分量的相位(Y.Ephraim，Speech enhancement using a minimum mean-square error short-time spectral amplitude estimator，IEEE Trans.A.S.S.P，32，1984)。
研究和试验都表明，语音频谱分量的实部和虚部都更符合伽马分布(Gamma Distribution)，同时采用伽马分布得到的MMSE估计子是高度非线性、复数值的滤波器，这 将会得到更好的噪声抑制性能。因此，本发明分别采用拉普拉斯分布和伽马分布来模拟噪声 和语音频谱分量的实部和虚部概率密度分布，具体分布如下式(2)、(3)所示。
$\left\{\begin{array}{c}p\left(N_R\right)=\frac{1}{{\lambda }_n}\exp (-\frac{2\left|N_R\right|}{{\lambda }_n})\\ p\left(N_I\right)=\frac{1}{{\lambda }_n}\exp (-\frac{2\left|N_I\right|}{{\lambda }_n})\end{array}\right.---\left(2\right)$
$\left\{\begin{array}{c}p\left(S_R\right)=\frac{\sqrt[4]{3}}{2\sqrt{{\mathrm{\pi \lambda }}_s}\sqrt[4]{2}}{\left|S_R\right|}^{-\frac{1}{2}}\exp (-\frac{\sqrt{3}\left|S_R\right|}{\sqrt{2}{\lambda }_s})\\ p\left(S_I\right)=\frac{\sqrt[4]{3}}{2\sqrt{{\mathrm{\pi \lambda }}_s}\sqrt[4]{2}}{\left|S_I\right|}^{-\frac{1}{2}}\exp (-\frac{\sqrt{3}\left|S_I\right|}{\sqrt{2}{\lambda }_s})\end{array}\right.---\left(3\right)$
上式中NR、NI、SR和SI分别代表噪声、语音频谱分量的实部和虚部，λn和λs分别代表 噪声和语音频谱分量的方差，t时刻分析窗的语音、噪声的第k个频谱分量的实部和虚部分别 表示为SR(t，k)、SI(t，k)和NR(t，k)、NI(t，k)，其概率密度分布分别为相应方差λs和λn的拉 普拉斯和伽马分布随机数。
由于语音的高度非平稳性，对每个k，不同时刻t得到的分布参数 {λs(1，k)，λs(2，k)，λs(3，k)…}和{λn(1，k)，λn(2，k)，λn(3，k)…}更应该理解为一个随机序列，本发明 需要从含噪语音中在线估计这些随机序列。
3.谱幅度系数的最小均方估计
人的听觉系统对频域变化更为敏感，因而从含噪语音的频域分量中估计纯语音的频域分 量能得到更好的结果。随机信号的估计要求分布模型以及误差测度已知，均方误差(MSE) 是最常用的估计准则，它要求计算出的估计信号和纯语音信号的平方误差的期望值最小。本 发明不使用高斯分布模型，所有的估计子都应包括实部和虚部的估计(R.Martin，Speech Enhancement using MMSE Short Time Spectral Estimation with Gamma Distributed Speech Priors， Proceeding of IEEE ICASSP，May，2002)。
含噪信号Y在l时刻的第k个频率分量表示为Y(k，l)＝YR(k，l)+jYI(k，l)，它包括噪声和语 音分量，即有Y(k，l)＝[SR(k，l)+DR(k，l)]+j[SI(k，l)+DI(k，l)]。MMSE估计问题可以归结为 在已知观测值Y(k，l)的条件下估计 $\hat{S}(k,l)={\hat{S}}_R(k,l)+j{\hat{S}}_I(k,l)$ 使得误差最小 由信号估计理论可知信号的条件最小均方误差值是信号的条件期望，即 $\hat{S}(k,l)=E\left[S(k,l)\right|Y(0,l),Y(1,l),...].$ 考虑到FFT频谱系数间的无关性假定和实部、虚部无关 性假定，最终可得MMSE估计子如式(4)所示：
$\hat{S}(k,l)=E\left[S_R(k,l)\right|Y_R(k,l)]+\mathrm{jE}[S_I(k,l)\left|Y_I(k,l)\right]---\left(4\right)$
依据上一节给出的拉普拉斯分布和伽马分布，利用实部和虚部概率密度分布进一步分别 计算实部和虚部的条件期望值，可得如(5)式所示的结果(R.Martin，Speech Enhancement using MMSE Short Time Spectral Estimation with Gamma Distributed Speech Priors，Proceeding of IEEE ICASSP，May，2002)。
$E\left[S_R(k,l)\right|Y_R(k,l)]=\frac{\sqrt[4]{1.5}}{2{\lambda }_n(k,l)\sqrt{{\mathrm{\pi \lambda }}_s(k,l)}p\left(Y_R(k,l)\right)}{\int }_{-\infty }^{\infty }{S}_R(k,l)·{{|S}_R(k,l)|}^{-0.5}.$
$\exp [-\frac{2|Y_R(k,l)-S_R(k,l)|}{{\lambda }_n(k,l)}]·\exp [-\frac{\sqrt{1.5}\left|S_R(k,l)\right|}{{\lambda }_s(k,l)}]{\mathrm{dS}}_R(k,l)$
$=\frac{\sqrt[4]{1.5}}{2{\lambda }_n(k,l)\sqrt{{\mathrm{\pi \lambda }}_s(k,l)}p\left[Y_R(k,l)\right]}\left\{\frac{2}{3}\exp \right[\frac{-2Y_R(k,l)}{{\lambda }_n(k,l)}\left]Y_R^{3/2}(k,l)\Phi \right[\mathrm{1.5,2.5},-2Y_R(k,l)G_1]$
$+\exp \left[\frac{-\sqrt{1.5}{Y}_R(k,l)}{{\lambda }_s(k,l)}\right]{\left(2G_2\right)}^{-1.5}\psi [-0.5,-\mathrm{0.5,2}{G}_2{Y}_R(k,l)]$
$-0.856·\exp \left[\frac{-2Y_R(k,l)}{{\lambda }_n(k,l)}\right]{\left(2G_2\right)}^{-1.5}\}---\left(5\right)$
上式中YR(k，l)表示含噪信号Y在l时刻的第k个频率分量的实部，Φ(a，b，z)＝M(a，b，z)表 示第一类合流超几何函数，Ψ(a，b，z)同样表示超几何函数，可由Φ(a，b，z)计算出来，
$G_1=\frac{\sqrt{1.5}{\lambda }_n+2{\lambda }_s}{2{\lambda }_s{\lambda }_n}$ 和 $G_2=\frac{\sqrt{1.5}{\lambda }_n-2{\lambda }_s}{2{\lambda }_s{\lambda }_n};$
其中p(YR(k，l))表示观察含噪语音频率分量实部的概率密度，有
$p\left[Y_R(k,l)\right]=\frac{\sqrt[4]{1.5}}{2{\lambda }_n(k,l)\sqrt{{\mathrm{\pi \lambda }}_s(k,l)}}{\int }_{-\infty }^{\infty }{\left|S_R(k,l)\right|}^{-0.5}\exp (-\frac{2|Y_R(k,l)-S_R(k,l)|}{{\lambda }_n(k,l)})$
$\exp (-\frac{\sqrt{1.5}\left|S_R(k,l)\right|}{{\lambda }_s(k,l)})d{S}_R$
$=\frac{\sqrt[4]{1.5}}{2{\lambda }_n(k,l)\sqrt{{\mathrm{\pi \lambda }}_s(k,l)}}\{2·\exp [-\frac{2Y_R(k,l)}{{\lambda }_n(k,l)}\left]\sqrt{Y_R(k,l)}\Phi \right[\mathrm{0.5,1.5},-2G_1{Y}_R(k,l)]$
$+\exp [-\frac{\sqrt{1.5}}{{\lambda }_s(k,l)}{Y}_R(k,l)]\frac{1}{\sqrt{2G_2}}\Psi [0.5,\mathrm{0.5,2}{G}_2{Y}_R(k,l)]$
$+\exp [-\frac{2Y_R(k,l)}{{\lambda }_n(k,l)}]\sqrt{\frac{\pi }{2G_2}}\}---\left(6\right)$
式中Φ[a，b，z]＝M·(a，b，z)表示第一类合流超几何函数，可利用级数求和计算，而
$\Psi [a,b,z]=\frac{\Gamma (1-b)}{\Gamma (a-b+1)}M(a,b,z)+\frac{\Gamma (b-1)}{\Gamma \left(a\right)}{z}^{1-b}M(a-b+\mathrm{1,2}-b,z)---\left(7\right)$
(I.S.Gradshteyn，Table of Inergrals，Series，and Products，1994)。
同理可得虚部的条件期望值，作为虚部的估计
$E\left[S_I(k,l)\right|Y_I(k,l)]=\frac{\sqrt[4]{1.5}}{2{\lambda }_n(k,l)\sqrt{{\mathrm{\pi \lambda }}_s(k,l)}p\left(Y_I(k,l)\right)}{\int }_{-\infty }^{\infty }{S}_I(k,l)·{\left|S_I(k,l)\right|}^{-0.5}.$
$\exp [-\frac{2|Y_I(k,l)-S_I(k,l)|}{{\lambda }_n(k,l)}]·\exp [-\frac{\sqrt{1.5}\left|S_I(k,l)\right|}{{\lambda }_s(k,l)}]{\mathrm{dS}}_I(k,l)$
$=\frac{\sqrt[4]{1.5}}{2{\lambda }_n(k,l)\sqrt{{\mathrm{\pi \lambda }}_s(k,l)}p\left[Y_I(k,l)\right]}\left\{\frac{2}{3}\exp \right[\frac{-2Y_I(k,l)}{{\lambda }_n(k,l)}\left]Y_I^{3/2}(k,l)\Phi \right[\mathrm{1.5,2.5},-2Y_I(k,l)G_1]$
$+\exp \left[\frac{-\sqrt{1.5}{Y}_I(k,l)}{{\lambda }_s(k,l)}\right]{\left({2G}_2\right)}^{-1.5}\psi [-0.5,-\mathrm{0.5,2}{G}_2{Y}_I(k,l)]$
$-0.856·\exp \left[\frac{-2Y_I(k,l)}{{\lambda }_n(k,l)}\right]{\left(2G_2\right)}^{-1.5}\}---\left(8\right)$
式中YI(k，l)表示含噪信号Y在l时刻的第k个频率分量的虚部，其他参见前式说明；
其中p[YI(k，l)]计算方法与如式(6)同理。
在实际语音通信中，各分布的参数往往不可能先验的知道，必须从含噪数据中估计，估 计方法见下文。
可以发现，式(5)和(8)利用了两个独立、服从拉普拉斯和伽马分布的随机变量的联 合分布，也就是说估计子实际上包括语音和噪声分量，也即假定了语音的无条件存在。实际 的语音通信环境中的含噪语音信号包括大量停顿带来的静音，这不仅包括句子之间的过渡甚 至还包括音节之间的暂停，因而在含噪信号中语音是不确定的，它只是依概率存在。因而上 述假定在实际应用中是不正确的：在输入数据中语音存在着大量暂停而背景噪声则始终存在。
4.语音存在的不确定性
由于(5)式和(8)式只是当前输入帧处于语音存在状态下的估计，本发明根据语音存 在不确定性对此做进一步扩展。含噪语音模型Y(k，l)＝S(k，l)+D(k，l)假定语音始终存在于输 入数据，若用H0和H1分别表示语音存在与否后，更准确的含噪模型为：
$Y\left(w\right)=\left\{\begin{array}{cc}S\left(w\right)+D\left(w\right),& H_1\\ D\left(w\right),& H_0\end{array}\right.---\left(9\right)$
为表达方便，本小节各表达式省去下标。考虑语音存在不确定性后的MMSE估计 E[SR|YR]应该改写成式(10)，其中
E[SR|YR]＝E[SR|YR，H1]P(H1|YR)
                                                 (10)
               +E[SR|YR，H0]P(H0|YR)
E[SR|YR，H0]表示语音不存在时由YR得到的MMSE估计，显然语音不存在时不可能得到语音 的估计，该项应该为0，因此考虑语音存在不确定性时的语音信号的实部估计为：
E[SR|YR]＝E[SR|YR，H1]P(H1|YR)                 (11)
计算(11)式要求后验概率P(H1(k)|Yk)已知，这可通过贝叶斯法则来计算，即：
$P\left(H_1\right|Y_R)=\frac{p\left(Y_R(k,l)\right|H_1)P\left(H_1\right)}{p\left(Y_R(k,l)\right|H_1)P\left(H_1\right)+p\left(Y_R\right|H_0)P\left(H_0\right)}$
$=\frac{\Gamma (k,l)}{1+\Gamma (k,l)}---\left(12\right)$
上式中 $\Gamma (k,l)=\frac{p\left(Y_R(k,l)\right|H_1)}{p\left(Y_R(k,l)\right|H_0)}·\frac{[1-P\left(H_0\right)]}{P\left(H_0\right)},$ p(YR(k，l)|H1)的计算如式(6)所示，而 p(YR(k，l)|H0)则表示只有噪声时的条件概率密度，其计算如式(2)所示。令P(H0)＝q表示 先验语音不存在概率。显然，q是先验未知的，如何得到q的估计 见第6小节。
由此得到了本发明的MMSE估计：
$E\left[S_R\right|Y_R]=\frac{\Gamma (k,l)}{1+\Gamma (k,l)}E[S_R|Y_R,H_1]=$
$\frac{\sqrt[4]{1.5}}{2{\lambda }_n(k,l)\sqrt{{\mathrm{\pi \lambda }}_s(k,l)}}\{2·\exp [-\frac{2Y_R(k,l)}{{\lambda }_n(k,l)}\left]\sqrt{Y_R(k.l)}\Phi \right[\mathrm{0.5,1.5},-2G_1{Y}_R(k,l)]$
$+\exp [-\frac{1.5}{{\lambda }_s(k,l)}{Y}_R(k,l)]\frac{1}{\sqrt{2G_2}}\Psi \left[\mathrm{0.5,0.5,2}{G}_2{Y}_R(k,l)\right]$
$+\exp [-\frac{2Y_R(k,l)}{{\lambda }_n(k,l)}]\sqrt{\frac{\pi }{2G_2}}\}·\frac{\Gamma (k,l)}{1+\Gamma (k,l)}---\left(13\right)$
上式中 $\Gamma (k,l)=\frac{p\left(Y_R(k,l)\right|H_1)}{p\left(Y_R(k,l)\right|H_0)}·\frac{[1-P\left(H_0\right)]}{P\left(H_0\right)},$ p(YR(k，l)|H1)表示语音存在时的YR(k，l)概 率密度分布，而p(YR(k，l)|H0)则表示只有噪声时的概率密度。
语音信号修正的虚部的估计与实部同理。
$E\left[S_R\right|Y_R]=\frac{\Gamma (k,l)}{1+\Gamma (k,l)}E[S_R|Y_R,H_1]---\left(14\right)$
5.谱分量方差的估计方法
本发明在估计纯语音谱系数时，要求语音和噪声谱分量的方差已知，但实际环境中，这 两个参数是不知道，而且没有任何先验知识，只能从含噪语音中估计。考虑到实际环境中语 音和噪声的非平稳性，方法应跟踪这些参数的变化。本发明使用上一帧抑制处理后语音信号 各频率分量的幅度方差作为λs(k，l)的估计 λn(k，l)的估计则比较复杂，本发明使用 VAD模块来判断当前输入帧是否纯噪声帧，如果是纯噪声帧则更新噪声参数，。这种硬判决 估计方法认为输入语音信号在语音-噪声和纯噪声两种状态中切换，噪声方差的估计只应在纯 噪声态中进行。在实际通信环境中，语音和噪声往往表现出高度非平稳特点，时变的统计特 性使得这种硬判决方法表现得更为稳健，因而此类方法获得了广泛使用。
本发明在上节基础上提出了一种基于似然比(Likelihood Ratio)的VAD方法，在式(6) 计算的概率密度基础上进一步通过似然函数的比较来判决当前输入帧是否为纯噪声帧。由于 每个频域分量包括近似独立的实部和虚部，其似然比同时包括实部和虚部，其形如(15)式 定义，且各频谱分量彼此不相关。
$\Lambda (k,l)=\frac{p\left[Y_R(k,l)\right|H_1\left]p\right[Y_I(k,l)\left|H_1\right]}{p\left[Y_R(k,l)\right|H_0\left]p\right[Y_I(k,l)\left|H_0\right]}---\left(15\right)$
上式中p(YR|H1)和p(YR|H0)的计算分别如式(6)和(2)所示，p(YI|H1)和p(YI|H0) 的计算同理。
由于VAD必须整帧进行，因而整帧的似然比为：
$\log \left[\Lambda \left(l\right)\right]=\frac{1}{K}\underset{k=0}{\overset{K-1}{\Sigma }}\log \left[\Lambda (k,l)\right]---\left(16\right)$
据此，本发明的VAD判决过程如式(17)所示，
$\left\{\begin{array}{c}{H}_0,\mathrm{if}\left[\log \left(\Lambda \right)\right]<{\theta }_{\Lambda }\\ H_1,\mathrm{if}\left[\log \left(\Lambda \right)\right]>{\theta }_{\Lambda }\end{array}---\left(17\right)\right.$
当[log(Λ)]＜θΛ时，判决作H1判决，即语音-噪声帧，否则作H0判决，即纯噪声帧。若当 前输入帧判决为纯噪声帧，则如式(18)所示计算噪声的各个谱分量。
${\hat{\sigma }}_n^2(k,l)=\left\{\begin{array}{cc}{\hat{\sigma }}_n^2(k,l-1),& H_0\\ {\alpha }_{\sigma }{\hat{\sigma }}_n^2(k,l-1)+(1-{\alpha }_{\sigma }){\sigma }_n^2(k,l),& H_1\end{array}\right.---\left(18\right)$
从式(13)可以观察到计算抑制结果还需要估计纯语音的各频率分量的幅度方差，本发 明直接采用上一帧的滤波语音作为纯语音的估计，从而得到这一参数的估计。
6.先验语音不存在概率的估计
对噪声抑制器式(13)和(14)来说，先验语音不存在概率是一个重要参数。在实际应 用中，此参数不但先验未知而且随时间和频率不同而变化，因而必须逐频率在线估计。本发 明提出如下的估计方法。
首先如式(19)所示计算当前输入帧平方幅度的和：
$A_{\mathrm{Sum}}^2\left(l\right)=\underset{k=0}{\overset{K-1}{\Sigma }}{A}^2(k,l)---\left(19\right)$
在计算信噪比 $\eta \left(l\right)=\frac{A_{\mathrm{Sum}}^2}{{\sigma }_n^2}$ 后，如式(20)所示进行时域递归平滑，
η(l)＝βη η(l-1)+(1-βη)η(l)    (20)
其中βη＝0.9。然后如式(21)所示计算全局概率Pglob(l)：
$P_{\mathrm{glob}}\left(l\right)=\left\{\begin{array}{cc}0,& {\eta }_{\min }\ge \bar{\eta }\left(l\right)\\ \frac{\log \left[\bar{\eta }\left(l\right)\right]-\log {\eta }_{\min }}{\log \frac{{\eta }_{\max }}{{\eta }_{\mathrm{mn}}}},& {\eta }_{\min }\le \bar{\eta }\left(l\right)\le {\eta }_{\max }\\ 1,& {\eta }_{\max }\le \bar{\eta }\left(l\right)\end{array}\right.---\left(21\right)$
其中ηmax和ηmin是经验常数，分别为-3dB和-11dB。
在为每个频率计算信噪比 $\gamma (k,l)=\frac{A^2(k,l)}{{\sigma }_n^2}$ 后，如式(22)所示进行时域递归平滑，
γ(k，l)＝βγ γ(k，l-1)+(1-βγ)γ(k，l)    (22)
其中βγ＝0.9。然后如式(23)所示计算局部概率Ploc(k，l)，
$P_{\mathrm{loc}}(k,l)=\left\{\begin{array}{cc}0,& {\gamma }_{\min }\ge \bar{\gamma }\left(l\right)\\ \frac{\log \left[\bar{\gamma }(k,l)\right]-\log {\eta }_{\min }}{\log \frac{{\gamma }_{\max }}{{\gamma }_{\min }}},& {\gamma }_{\min }\le \bar{\gamma }\left(l\right)\le {\gamma }_{\max }\\ 1,& {\gamma }_{\max }\le \bar{\gamma }\left(l\right)\end{array}\right.---\left(23\right)$
其中ηmax和ηmin是经验常数，分别为-1dB和-9dB。最终得到的先验语音不存在概率为
$\hat{q}(k,l)=1-P_{\mathrm{loc}}(k,l)P_{\mathrm{glob}}(k,l)---\left(24\right)$
分析式(24)可知，本发明充分利用了语音信号相邻帧间的时域相关性，结合考虑了当前输 入帧的全局语音不存在和各频率局部语音分量不存在的可能性，估计过程具有更好的稳健性。
7.ITU-G.160协议的要求
本发明主要用于各种语音增强设备(Voice Enhancement Device，VED)，以提高语音通信 的质量，但是很多时候网络中还传输DTMF音、传真音等信号音调。显然，任何噪声抑制算 法在处理过程中不能对这些信号音产生不利影响，对此ITU-G.160协议提出了明确的要求。 这些信号音都有一个或2个单频音混合而成，在频域上表现为在一个或多个频率上具有尖峰。 本发明的步骤中包括短时傅立叶变换，在此基础上检测这些尖峰是很方便的，可以很容易的 把信号音调区分出来。为满足G.160协议的要求，本发明在短时傅立叶变换之后加入一个检 测环节，对输入进行判断，如果发现只有一个或多个明显的尖峰，则判为信号音调，不进行 任何抑制处理。在判决过程中，本发明计算所有频率分量的能量和，而后求出2个平方幅度 的最大值，并从能量和中减取最大值的和，若最大值和大于能量和且彼此接近，则判当前输 入为信号音调。