webrtc中的噪声抑制之四：语音噪声概率计算

webrtc中的噪声抑制之四：语音噪声概率计算

回顾

webrtc的噪声抑制，通过quantile方法初始估计出噪声后，采用DD方法估算出先验和后验信噪比，然后根据此计算LRT均值特征，同时结合频谱平坦度和频谱差异，计算得出当前帧语音噪声概率，最后根据当前帧语音/噪声概率和估算出的先验/后验信噪比完成最终的噪声估计和维纳滤波。此前就维纳滤波和噪声估计做了学习研究，本文探讨对语音/噪声概率估计的研究方法。因为语音噪声概率计算其实属于语音激活检测的范畴。在深入webrtc算法之前，我们还是扫盲一下VAD-Voice Activity Detetion。

语音激活检测VAD-Voice Activity Detetion

语音激活检测VAD技术是语音通话的关键技术之一，传统的VAD方法如线谱频率、全带宽信号能量、低频带信号能量和过零率等，已经广泛用于电信领域，有兴趣可以看一下G.729B的VAD部分[1]。传统方法一般给出是非判决，高信噪比下表现很好，但低信噪比，以及话音起始或结束帧令传统方法无法逾越，所以在上个世纪70年代末期引入软判决方法，90年代引入了统计模型和最大似然比等概率的方法，极好的弥补了传统方法的缺陷。统计概率模型又往往采用高斯型概率分布，应该大学里都学过。所以下面我们先了解一下所谓的likelihood。

ML-Maximum Likelihood和LRT- Likelihood Ratio Test

首先看一下Likelihood or Likelihood function，一般翻译为似然函数，它是推理统计学（Statistical inference）一个重要的概念，很遗憾《推理统计学》没有在大学里学过，只能找一些网文来扫盲一下了。下面是baidu引用《概率统计》“王翠香编．北京大学出版社，2010”的介绍。随着机器学习和深度学习的蓬勃发展，这些都要变成必备知识了。

似然函数在推断统计学（Statistical inference）中扮演重要角色，如在最大似然估计和费雪信息之中的应用等等。
“似然性”与“或然性”或“概率”意思相近，都是指某种事件发生的可能性，但是在统计学中，“似然性”和“或然性”或
“概率”又有明确的区分。概率用于在已知一些参数的情况下，预测接下来的观测所得到的结果，而似然性则是用于在已
知某些观测所得到的结果时，对有关事物的性质的参数进行估计。

似然函数是一种关于统计模型中的参数的函数，表示模型参数中的似然性，给定输出x时，关于参数θ的似然函数L(θ|x)（在数值上）等于给定参数θ后变量X的概率：
$$L(\theta |x)=P(X=x|\theta )$$
在离散性概率分布的情况下，也写成如下的表达形式：
$$L(\theta |x)=p_{\theta }(x)=P(X=x)$$
$p(x)$表示X取x时的概率，但写成$P(X=x|\theta )$或$P(X=x;\theta )$不是条件概率，需要特别注意。似然函数的主要用法在于比较它相对取值，虽然这个数值本身不具备任何含义。例如，考虑一组样本，当其输出固定时，这组样本的某个未知参数往往会倾向于等于某个特定值，而不是随便的其他数，此时，似然函数是最大化的。
似然函数乘以一个正的常数之后仍然是似然函数，其取值并不需要满足归一化条件：
$$\sum \alpha L(\theta |x)̸=1;\alpha >0$$

最大似然估计（ML）

以离散型样本分布为例，假设分布概率函数$P=p(x;\theta )$, $x$是发生的样本，$\theta$是待估计参数，那么如果我们拥有了样本值：$x_1,x_2,.......,x_n$，定义
$$L(\theta )=L(x_1,x_2,.......,x_n;\theta )=\prod _{i=1}^{n}p(x_i;\theta )$$
为样本的似然函数。如果能得到
$$L(\theta )=L(x_1,x_2,.......,x_n;\hat{\theta })=MA{X}_{\theta \in \Theta }L(x_1,x_2,.......,x_n;\theta )$$
此时$\hat{\theta }$即被称为 参数$\theta$的极大似然估计。为了方便计算，对定义取对数，可以将成法变成加法。
$$lnL(x_1,x_2,.......,x_n;\theta )=ln\prod _{i=1}^{n}p(x_i;\theta )=\sum _{i=1}^{n}ln[p(x_i;\theta )]$$
能满足对数似然函数最大的$\hat{\theta }$等价于参数$\theta$的极大似然估计。这样就将问题进一步简化了。

似然比检验（LRT）

似然比检验是基于上述理念，进一步引入某个假设，来证明新引入假设是否成立的方法。首先，假设原来的条件已经令似然函数取得最大值，新引入的假设不会超过这个最大似然函数值，但如果约束条件有效，有约束的最大值应当接近无约束的最大值，两者的比值接近于1。否则远远小于1。继续上一节的定义，我们要检验一个无效假设 $H_0:\theta ={\theta }_0$，反之则$H_1:\theta ̸={\theta }_0$。我们可以先列出似然比函数：
$$\lambda =\frac{L_0(\theta )=L(x_1,x_2,.......,x_n;{\theta }_0)}{L(\theta )=L(x_1,x_2,.......,x_n;\theta )}$$
这里的 $\lambda$范围在（0，1）之间，LRT的定义同样是取自然对数：
$$LRT=-2ln\lambda =2[ln(L(x_1,x_2,.......,x_n;\theta ))-ln(L_0(x_1,x_2,.......,x_n;{\theta }_0))]$$
关于似然函数和似然比的概念简单介绍这么多，因为在噪声概率统计的时候Webrtc用到了LRT的方法。结合具体场景我们可以加深对LRT的理解。

LRT均值特征

假设语音主要受加性噪声干扰，每一帧可能有两种情况
$$\begin{array}{rl}& H_0:SpeechAbsent----X=N\\ & H_1:SpeechPresent---X=N+S\end{array}$$
其中$S,N,X$分别表示L维语音、噪声和受噪语音的傅立叶变换系数向量，为了在频域上分析，引入下标k，$S_k,N_k,X_k$则表示对应第k个频域元素。我们采用高斯统计模型，回忆一下pdf函数的形式：
$$f(x)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{(x-\mu {)}^2}{2{\sigma }^2}}$$
而语音信号每一帧DFT系数采用了渐进正态分布（ asymptotically independent Gaussian random）模型，对于之前假设的条件概率分布定义如下：
p ( X ∣ H 0 ) = ∏ k = 0 L − 1 1 π λ N ( k ) e x p { − ∣ X k ∣ 2 λ N ( k ) } p ( X ∣ H 1 ) = ∏ k = 0 L − 1 1 π [ λ N ( k ) + λ S ( k ) ] e x p { − ∣ X k ∣ 2 λ N ( k ) + λ S ( k ) } \begin{aligned} &p(X|H_0)= \prod_{k=0}^{L-1} \frac{1}{\pi\lambda_N(k)} exp\lbrace-\frac{|X_k|^2}{\lambda_N(k)} \rbrace\\ &p(X|H_1)= \prod_{k=0}^{L-1} \frac{1}{\pi[\lambda_N(k)+\lambda_S(k)]} exp\lbrace-\frac{|X_k|^2}{\lambda_N(k)+\lambda_S(k)} \rbrace \end{aligned} ​p(X∣H0​)=k=0∏L−1​πλN​(k)1​exp{−λN​(k)∣Xk​∣2​}p(X∣H1​)=k=0∏L−1​π[λN​(k)+λS​(k)]1​exp{−λN​(k)+λS​(k)∣Xk​∣2​}​
定义${\lambda }_N(k)和{\lambda }_S(k)$是$N_k和S_k$的方差。而对第k个频率分量的LRT定义如下：
Λ k = p ( X ∣ H 0 ) p ( X ∣ H 1 ) = 1 1 + ξ k e x p { γ k ξ k 1 + ξ k } ξ k = λ S ( k ) λ N ( k ) − − − − − − p r i o r i S N R γ k = ∣ X k ∣ 2 λ N ( k ) − − − − − − p o s t S N R \begin{aligned} &\Lambda_k=\frac{p(X|H_0)}{p(X|H_1)}=\frac{1}{1+\xi_k}exp\lbrace\frac{\gamma_k\xi_k}{1+\xi_k} \rbrace\\ &\xi_k=\frac{\lambda_S(k)}{\lambda_N(k)} ------ priori SNR\\ &\gamma_k=\frac{|X_k|^2}{\lambda_N(k)} ------ postSNR\\ \end{aligned} ​Λk​=p(X∣H1​)p(X∣H0​)​=1+ξk​1​exp{1+ξk​γk​ξk​​}ξk​=λN​(k)λS​(k)​−−−−−−prioriSNRγk​=λN​(k)∣Xk​∣2​−−−−−−postSNR​
对于0均值的统计信号，方差就是功率估计，由此再次用信噪比的概念简化分析流程。先验信噪比是未知的，其中理想情况下，$\widehat{{\xi }_k^{(ML)}}={\gamma }_k-1$。ML也表示最大似然意义下的先验信噪比。所以根据后验信噪比的ISD（Itakura–Saito distortion）准则。
$$log{\hat{\Lambda }}_k^{(ML)}={\gamma }_k-log{\gamma }_k-1$$
上面的算法据说更倾向于$H_1$，或者为了减少偏差，利用直接判决DD法估算先验信噪比。
$${\hat{\xi }}_k(n{)}^{(DD)}=\alpha \frac{{\hat{\Lambda }}_k^2(n-1)}{{\lambda }_n(k,n-1)}+(1-\alpha )P[{\gamma }_k(n)-1]$$
这个和之前的直接判决法应该是等价的，具体性能可以参考论文[2]。
不过WebRtc中的信噪比是频谱幅值之比，具体计算过程还没有完全和上述公式对应起来，博文[3]的推导我也没有完全看明白，此处略有遗憾。

频谱平坦度

频谱平坦度也叫做维纳熵，即各个频率分量的几何平均数与算术平均数的比值， 网上有定义，下面直接出公式了。
$$\begin{gathered} Flatnes{s}_{magn}=\frac{\sqrt[N]{{\prod }_{n=0}^{N-1}{Y}_k(m)}}{\frac{{\sum }_{n=0}^{N-1}{Y}_k(m)}{N}} \\ =\frac{exp(\frac{1}{N}{\sum }_{n=0}^{N-1}ln[Y_k(m)])}{\frac{1}{N}{\sum }_{n=0}^{N-1}{Y}_k(m)} \end{gathered}$$
直观的说如果每个分量都相等，比值为1；分量值相当，比值接近于1，这时认为是白噪声；分量差异大的话，比值接近0，一般认为有语音产生。下面是函数。

static void ComputeSpectralFlatness(NoiseSuppressionC* self,
                                    const float* magnIn) {
  size_t i;
  size_t shiftLP = 1;  // Option to remove first bin(s) from spectral measures.
  float avgSpectralFlatnessNum, avgSpectralFlatnessDen, spectralTmp;

  // Compute spectral measures.
  // For flatness.
  avgSpectralFlatnessNum = 0.0;
  avgSpectralFlatnessDen = self->sumMagn;
  for (i = 0; i < shiftLP; i++) {
    avgSpectralFlatnessDen -= magnIn[i];
  }
  // Compute log of ratio of the geometric to arithmetic mean: check for log(0)
  // case.
  for (i = shiftLP; i < self->magnLen; i++) {
    if (magnIn[i] > 0.0) {
      avgSpectralFlatnessNum += (float)log(magnIn[i]);
    } else {
      self->featureData[0] -= SPECT_FL_TAVG * self->featureData[0];
      return;
    }
  }
  // Normalize.
  avgSpectralFlatnessDen = avgSpectralFlatnessDen / self->magnLen;
  avgSpectralFlatnessNum = avgSpectralFlatnessNum / self->magnLen;

  // Ratio and inverse log: check for case of log(0).
  spectralTmp = (float)exp(avgSpectralFlatnessNum) / avgSpectralFlatnessDen;

  // Time-avg update of spectral flatness feature.
  self->featureData[0] += SPECT_FL_TAVG * (spectralTmp - self->featureData[0]);
  // Done with flatness feature.
}

频谱偏差

这个是将输入的频谱和估算的频谱模板进行对比来评价语音和噪声概率的方法。是一种谱估计方法，暂时没追本溯源到出处，分析欠奉。

噪声概率计算

根据上述三个指标，经过时间平滑和tanh影射，最后算出当前语音帧的噪声概率。

static void SpeechNoiseProb(NoiseSuppressionC* self,
                            float* probSpeechFinal,
                            const float* snrLocPrior,
                            const float* snrLocPost) {
  size_t i;
  int sgnMap;
  float invLrt, gainPrior, indPrior;
  float logLrtTimeAvgKsum, besselTmp;
  float indicator0, indicator1, indicator2;
  float tmpFloat1, tmpFloat2;
  float weightIndPrior0, weightIndPrior1, weightIndPrior2;
  float threshPrior0, threshPrior1, threshPrior2;
  float widthPrior, widthPrior0, widthPrior1, widthPrior2;

  widthPrior0 = WIDTH_PR_MAP;
  // Width for pause region: lower range, so increase width in tanh map.
  widthPrior1 = 2.f * WIDTH_PR_MAP;
  widthPrior2 = 2.f * WIDTH_PR_MAP;  // For spectral-difference measure.

  // Threshold parameters for features.
  threshPrior0 = self->priorModelPars[0];
  threshPrior1 = self->priorModelPars[1];
  threshPrior2 = self->priorModelPars[3];

  // Sign for flatness feature.
  sgnMap = (int)(self->priorModelPars[2]);

  // Weight parameters for features.
  weightIndPrior0 = self->priorModelPars[4];
  weightIndPrior1 = self->priorModelPars[5];
  weightIndPrior2 = self->priorModelPars[6];

  // Compute feature based on average LR factor.
  // This is the average over all frequencies of the smooth log LRT.
  logLrtTimeAvgKsum = 0.0;
  for (i = 0; i < self->magnLen; i++) {
    tmpFloat1 = 1.f + 2.f * snrLocPrior[i];
    tmpFloat2 = 2.f * snrLocPrior[i] / (tmpFloat1 + 0.0001f);
    besselTmp = (snrLocPost[i] + 1.f) * tmpFloat2;
    self->logLrtTimeAvg[i] +=
        LRT_TAVG * (besselTmp - (float)log(tmpFloat1) - self->logLrtTimeAvg[i]);
    logLrtTimeAvgKsum += self->logLrtTimeAvg[i];
  }
  logLrtTimeAvgKsum = (float)logLrtTimeAvgKsum / (self->magnLen);
  self->featureData[3] = logLrtTimeAvgKsum;
  // Done with computation of LR factor.

  // Compute the indicator functions.
  // Average LRT feature.
  widthPrior = widthPrior0;
  // Use larger width in tanh map for pause regions.
  if (logLrtTimeAvgKsum < threshPrior0) {
    widthPrior = widthPrior1;
  }
  // Compute indicator function: sigmoid map.
  indicator0 =
      0.5f *
      ((float)tanh(widthPrior * (logLrtTimeAvgKsum - threshPrior0)) + 1.f);

  // Spectral flatness feature.
  tmpFloat1 = self->featureData[0];
  widthPrior = widthPrior0;
  // Use larger width in tanh map for pause regions.
  if (sgnMap == 1 && (tmpFloat1 > threshPrior1)) {
    widthPrior = widthPrior1;
  }
  if (sgnMap == -1 && (tmpFloat1 < threshPrior1)) {
    widthPrior = widthPrior1;
  }
  // Compute indicator function: sigmoid map.
  indicator1 =
      0.5f *
      ((float)tanh((float)sgnMap * widthPrior * (threshPrior1 - tmpFloat1)) +
       1.f);

  // For template spectrum-difference.
  tmpFloat1 = self->featureData[4];
  widthPrior = widthPrior0;
  // Use larger width in tanh map for pause regions.
  if (tmpFloat1 < threshPrior2) {
    widthPrior = widthPrior2;
  }
  // Compute indicator function: sigmoid map.
  indicator2 =
      0.5f * ((float)tanh(widthPrior * (tmpFloat1 - threshPrior2)) + 1.f);

  // Combine the indicator function with the feature weights.
  indPrior = weightIndPrior0 * indicator0 + weightIndPrior1 * indicator1 +
             weightIndPrior2 * indicator2;
  // Done with computing indicator function.

  // Compute the prior probability.
  self->priorSpeechProb += PRIOR_UPDATE * (indPrior - self->priorSpeechProb);
  // Make sure probabilities are within range: keep floor to 0.01.
  if (self->priorSpeechProb > 1.f) {
    self->priorSpeechProb = 1.f;
  }
  if (self->priorSpeechProb < 0.01f) {
    self->priorSpeechProb = 0.01f;
  }

  // Final speech probability: combine prior model with LR factor:.
  gainPrior = (1.f - self->priorSpeechProb) / (self->priorSpeechProb + 0.0001f);
  for (i = 0; i < self->magnLen; i++) {
    invLrt = (float)exp(-self->logLrtTimeAvg[i]);
    invLrt = (float)gainPrior * invLrt;
    probSpeechFinal[i] = 1.f / (1.f + invLrt);
  }
}

参考资料：
[1]：《G.729B中话音激活检测算法的改进及仿真》 胡斌 等 桂林电子科技大学
[2]：《A Statistical Model-Based Voice Activity Detection》Jongseo Sohn, Nam Soo Kim, Wonyong Sung, IEEE.
[3]： 《WebRTC之noise suppression算法》https://blog.csdn.net/shichaog/article/details/52514816
[4]：https://blog.csdn.net/s09094031/article/details/84673428 单通道噪声抑制算法总结
[5]：https://wenku.baidu.com/view/f0e81ab5b8f67c1cfad6b865.html《基于改进噪声谱估计的单通道语音增强研究》
[6]：《语音信号处理》第三版 赵力等 机械工业出版社