webrtc中的噪声抑制之一：频域维纳滤波

webrtc中的噪声抑制之一：频域维纳滤波

前言

在开源的噪声抑制算法中,webrtc ns是很有名的,社区里也有很多分享的文章,但要么深要么浅,还有一些误导读者的,所以趁着移植项目的机会,从盲人摸象到庖丁解牛的学习一番这里面的算法原理和工程实现。
WebRtc Ns模块采用的是频域维纳滤波的方法，结合VAD检测得到前验信噪比和后验信噪比，算出频域维纳滤波器的系数，在频域实现了噪声的滤除。该模块有定点和浮点两部分代码，其中定点代码相对比较复杂，而浮点代码主要分WebRtcNs_AnalyzeCore和WebRtcNs_ProcessCore两部分模块，先分析当前帧，再对当前帧进行滤波，结构清晰，方便阅读和学习，所以本文研究浮点代码的实现。

频域维纳滤波器

使用维纳滤波进行语音降噪的过程，其实是把降噪过程视为一个线性时不变系统，当带噪语音通过这个系统时，在均方误差最小化准则下，使得系统的输出与期望的纯净语音信号最接近的过程。
实际当中我们能观察到的信号$y(n)$是含有噪声的。 假设带噪语音信号可以表述为
$$y(n)=s(n)+n(n)$$
$s(n)$为语音信号，$d(n)$为加性噪声。维纳滤波方法就是涉及一个数字滤波器$h(n)$，当输入$y(n)$经过滤波器后，滤波器的输出可以最大程度的逼近$s(n)$。
$$\hat{s}(n)=y(n)\ast h(n)=\sum _{k=0}^{K-1}y(n-k)h(k)$$
上式进行DFT得到
$$\hat{S}(\omega )=H(\omega )Y(\omega )$$
在任意频率${\omega }_k$，我们可以计算出误差估计
$$E({\omega }_k)=S({\omega }_k)-\hat{S}({\omega }_k)=S({\omega }_k)-H({\omega }_k)Y({\omega }_k)$$
由此定义频域均方误差的代价函数
$$\begin{array}{rl}{J}_2& =E[|E({\omega }_k){|}^2]\\ & =E\{[S({\omega }_k)-H({\omega }_k)Y({\omega }_k){]}^{\ast }[S({\omega }_k)-H({\omega }_k)Y({\omega }_k)]\}\\ & =E[|S({\omega }_k){|}^2]-H({\omega }_k)E[S^{\ast }({\omega }_k)Y({\omega }_k)]-H^{\ast }({\omega }_k)E[S({\omega }_k)Y^{\ast }({\omega }_k)]+|H({\omega }_k){|}^{2}E[|Y({\omega }_k){|}^2]\\ & =E[|S({\omega }_k){|}^2]-H({\omega }_k)P_{ys}({\omega }_k)-H^{\ast }({\omega }_k)P_{sy}({\omega }_k)+|H({\omega }_k){|}^2{P}_{yy}({\omega }_k)\end{array}$$
公式的最后一行定义了$P_{yy}({\omega }_k)=E[|Y({\omega }_k){|}^2]$，叫做信号的功率谱。$P_{ys}({\omega }_k)=E[S^{\ast }({\omega }_k)Y({\omega }_k)]$ 定义为输入信号$y(n)$和期望信号$s(n)$的互功率谱，$P_{sy}({\omega }_k)=E[S({\omega }_k)Y({\omega }_k{)}^{\ast }]$是互功率谱的另外一种形式，$P_{ys}({\omega }_k)=P_{sy}^{\ast }({\omega }_k)$。对$J_2$求复导数，并令其等于0，得到如下结果：（共轭求偏导这块有点没理解）
$$\frac{\partial J_2}{\partial H({\omega }_k)}=H^{\ast }({\omega }_k)P_{yy}({\omega }_k)-P_{ys}({\omega }_k)=[H({\omega }_k)P_{yy}({\omega }_k)-P_{sy}({\omega }_k){]}^{\ast }=0$$
由此得到$H({\omega }_k)$的维纳解：
$$H({\omega }_k)=\frac{P_{sy}({\omega }_k)}{P_{yy}({\omega }_k)}$$
假设噪声为加性噪声，且与语音信号不相关
$$\begin{gathered} P_{sy}({\omega }_k)=P_{ss}({\omega }_k) \\ {P}_{yy}({\omega }_k)=P_{xx}({\omega }_k)+P_{nn}({\omega }_k) \end{gathered}$$
上式代入维纳解，我们得到：
$$H({\omega }_k)=\frac{P_{ss}({\omega }_k)}{P_{ss}({\omega }_k)+P_{nn}({\omega }_k)}$$
我们定义这种形式为“频域维纳滤波器”。如果得到$H({\omega }_k)$，通过频域相乘，很容易得到$\hat{S}(\omega )=H(\omega )Y(\omega )$。我们观察频域维纳滤波器的计算公式，只涉及到功率谱的计算，看上去也比时域要简单许多，但是语音信号的短时平稳特性，令我们求真实的输入信号功率谱很麻烦，噪声功率谱也不容易得到。所以需要继续推导寻找近似逼近的方法。我们定义${\xi }_{{\omega }_k}={\xi }_k=\frac{P_{ss}({\omega }_k)}{P_{nn}({\omega }_k)}$为频率${\omega }_k$处的先验信噪比$SN{R}_{pre}$，即信号没有被噪声干扰的信噪比。${\gamma }_{{\omega }_k}={\gamma }_k=\frac{P_{yy}({\omega }_k)}{P_{nn}({\omega }_k)}$为频率${\omega }_k$处的后验信噪比$SN{R}_{post}$，即信号引入加性噪声后的信噪比。则$H({\omega }_k)$可以写成着两种信噪比的表达方法。
$$H({\omega }_k)=\frac{{\xi }_k}{1+{\xi }_k}=1-\frac{1}{{\gamma }_k}$$
至此，对频域维纳滤波器的求解变成了求解信号的先验或者后验信噪比问题。看上去难度似乎大大降低了，然而实际上对两种信噪比的求解也非常困难。所以有的算法希望利用两种信噪比来平滑计算结果，所以引入一个平滑因子 $\alpha$，导出 ：
$$\begin{array}{rl}{\xi }_i(k)& =\alpha {\xi }_i(k)+(1-\alpha ){\xi }_i(k)\\ & =\alpha {\xi }_i(k)+(1-\alpha )({\gamma }_i(k)-1)\\ & \approx \alpha {\xi }_{i-1}(k)+(1-\alpha )({\gamma }_i(k)-1)\end{array}$$
因为计算当前帧的先验信噪比是一个非因果事件，所以利用上次滤波后的先验信噪比很显然是一种流行的做法，那么最后我们得出先验信噪比的估计值为：
$$\widehat{{\xi }_i}(k)=\alpha {\xi }_{i-1}(k)+(1-\alpha )({\gamma }_i(k)-1)$$
进而我们得出当前帧的维纳滤波器
$$\hat{H}(k)=\frac{\widehat{{\xi }_i}(k)}{1+\widehat{{\xi }_i}(k)}$$
最后利用此公式实现频域的维纳滤波：
$$\hat{S}(\omega )=H(\omega )Y(\omega )$$
可以看出，这和时域从维纳解导出自适应滤波的过程还是有些差异的。从公式我们观察两点结论：

1. 通过该公式我们将降噪问题转化为求解信噪比问题，降低了问题复杂度。
2. 频域维纳滤波需要准确的估算出先验信噪比和后验信噪比，这两个值得准确程度和收敛速度决定了滤波器以及整个降噪算法的性能。

WebRtc中的WienerFilter

以下是代码

// Estimate prior SNR decision-directed and compute DD based Wiener Filter.
// Input:
//   * |magn| is the signal magnitude spectrum estimate.
// Output:
//   * |theFilter| is the frequency response of the computed Wiener filter.
static void ComputeDdBasedWienerFilter(const NoiseSuppressionC* self,
                                       const float* magn,
                                       float* theFilter) {
  size_t i;
  float snrPrior, previousEstimateStsa, currentEstimateStsa;

  for (i = 0; i < self->magnLen; i++) {
    // Previous estimate: based on previous frame with gain filter.
    previousEstimateStsa = self->magnPrevProcess[i] /
                           (self->noisePrev[i] + 0.0001f) * self->smooth[i];
    // Post and prior SNR.
    currentEstimateStsa = 0.f;
    if (magn[i] > self->noise[i]) {
      currentEstimateStsa = magn[i] / (self->noise[i] + 0.0001f) - 1.f;
    }
    // DD estimate is sum of two terms: current estimate and previous estimate.
    // Directed decision update of |snrPrior|.
    snrPrior = DD_PR_SNR * previousEstimateStsa +
               (1.f - DD_PR_SNR) * currentEstimateStsa;
    // Gain filter.
    theFilter[i] = snrPrior / (self->overdrive + snrPrior);
  }  // End of loop over frequencies.
}

从代码的核心行可以看出这个和上一篇算法推算的原理一致，是利用两个信噪比综合得出当前系数的办法。具体细节处理还有些差异，暂时没有对比，当我们完成对整个vad统计下信噪比的计算之后，再回头看看能否理解每一步公式的细节。

小结

本篇学习和分析了WebRtc中的频域维纳滤波器，整个Ns都是围绕着这个来工作的，本文参考了引用里提到的书籍和几篇博文，特此感谢。文中有纰漏的地方非常欢迎纠错。

引用

1.《MATLAB在语音信号分析与合成中的应用》宋之用 北京航空航天大学出版社
2.《一个频域语音降噪算法实现及改进方法》 https://www.cnblogs.com/icoolmedia/p/weiner_audio_ns.html
3.《WebRTC之noise suppression算法》https://blog.csdn.net/shichaog/article/details/52514816
4.《Webrtc NS模块算法》https://blog.csdn.net/qq_28882043/article/details/80885240
5.《webrtc 单通道降噪算法（ANS）简析》https://blog.csdn.net/audio_algorithm/article/details/80812408