[人工智能] webrtc ns模块代码公式详细解读

开发: C++知识库 Java知识库 JavaScript Python PHP知识库人工智能区块链大数据移动开发嵌入式开发工具数据结构与算法开发测试游戏开发网络协议系统运维
教程: HTML教程 CSS教程 JavaScript教程 Go语言教程 JQuery教程 VUE教程 VUE3教程 Bootstrap教程 SQL数据库教程 C语言教程 C++教程 Java教程 Python教程 Python3教程 C#教程
数码: 电脑笔记本显卡显示器固态硬盘硬盘耳机手机 iphone vivo oppo 小米华为单反装机图拉丁

-> 人工智能 -> webrtc ns模块代码公式详细解读 -> 正文阅读

[人工智能]webrtc ns模块代码公式详细解读

总述

webrtc的降噪模块主要分为3个部分：模块初始化、噪声分析（analysis）、噪声抑制。模块初始化是在最开始为降噪模块设置参数以及初始化一些状态的值。噪声分析模块，主要进行噪声估计、speech/noise概率计算等。噪声抑制模块则根据前面计算的语音概率和噪声使用维纳滤波来抑制噪声。下面首先对一些相关公式进行介绍，再对整个算法流程进行分析。

模块初始化

webrtc的降噪代码主要写在noise_suppression.c中，开始降噪之前首先要进行一些模块初始化操作，实际执行初始化的函数为WebRtcNs_InitCore(Noise SuppressionC *self, uint32_t fs)。本算法仅支持采样率为8kHz或16kHz的音频。对于8k的音频，每次更新长度为blockLen=80的数据，与先前的数据拼接成anaLen=128长度的数据进行处理，我的理解即为帧移为80，帧长为128；若为16k的采样率，则上述长度变为原来的2倍。对anaLen长的数据求fft后取一半加一个点进行分析，也即magnLen = anaLen/2 + 1（因为实信号的fft是共轭对称的）。这部分代码本身有注释，这里就不进行过多分析

噪声分析

这一部分内容包括计算LRT均值特征、频谱平坦度、频谱模板相似度（差异度）等特征，然后利用这些特征计算语音概率，之后进行噪声频谱估计的更新。执行这一步的函数为WebRtcNs_AnalyzeCore(NoiseSuppressionC *self, const int16_t *speechFrame)，整个流程大致如下面框图所示：

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-exqEBhkS-1638020359765)(C:\Users\user\Desktop\Note\Note_pic\webrtc ns noise analysis.png)]

数据准备

每一次读入长为80个samples的数据，然后对分析帧进行更新，也就是将analyze的最新80个点换为读入的数据：

UpdateBuffer(speechFrame, self->blockLen, self->anaLen, self->analyzeBuf);

而后对analyzeBuf进行加窗并计算能量，如果能量为0则直接退出，不进行WebRtcNs_AnalyzeCore的后续计算：
```
energy = WindowingEnergy(self->window, self->analyzeBuf, self->anaLen, winData);
if (energy == 0.0) {
    return;
}
```
数据分析
- 傅里叶变换
  
  FFT(self, winData, self->anaLen, self->magnLen, real, imag, magn, lmagn, 1, &signalEnergy, &sumMagn);
  
  这一步还同时计算了信号能量，其值为频谱幅值平方的和
- 计算能量与幅值相关数据
```
//能量归一化，除以帧长（self->normMagnLen=1/MagnLen）
signalEnergy *= self->normMagnLen;   
self->signalEnergy = signalEnergy; 
//计算幅值之和
self->sumMagn = sumMagn;
```
分位数噪声估计（Quantile noise estimate）

此步骤的主要目的是对输入帧进行一个初始的噪声估计，使用“分位数噪声估计”。分位数噪声估计的基本思想是：对于某一个特定频点，将所有时间帧在该频点的幅度保存下来并排序。设置一个分位值，低于分位值的被视为噪声，高于分位值的被视为语音。但在实际计算中，为了保证实时性，不可能等收集几百帧数据后再排序取分位数，因而再webrtc算法中，采用了自适应的分数更新：设置一个对数分位值lquantile作为初始对数噪声，每当输入帧的幅值大于该分位值时，说明分位值偏小，使其增加，反之则使其减少：
```
if (lmagn[i] > self->lquantile[offset + i]) {
    self->lquantile[offset + i] += QUANTILE * delta * norm_counter_weight;
} else {
    self->lquantile[offset + i] -= (1.f - QUANTILE) * delta * norm_counter_weight;
}
```
可以看出delta和norm_counter_weight控制着更新步长。首先看delta：
```
if (self->density[offset + i] > 1.0) {
    delta = FACTOR / self->density[offset + i];
} else {
    delta = FACTOR;
}
```
density可以理解为概率密度，用来衡量每个频点分位数估计的准确程度，其更新如下：
```
if (fabsf(lmagn[i] - self->lquantile[offset + i]) < WIDTH) {
    self->density[offset + i] =
            self->counter[s] * self->density[offset + i] * norm_counter_weight + density_plus_weight;
}
```
即当lmgn和lquantile足够接近时，认为当前频点的噪声估计值lquantile比较精确，则增加此频点的概率密度。再回看前面的delta，当density大于1时，delta减小。也即：在当前频点的噪声估计值比较准确时（density>1），为了更精细化的估计，使更新步长减小（类似于神经网络训练的学习率逐渐下降）。

再看norm_counter_weight，此变量实际上是counter的倒数。counter是一个长度为3的数组，初始值为[67,134,1]。程序中实际对每个频点设置了三个分位点（lquantile），三个counter[i]分别是三个分位点的计数器——每一次噪声估计都加一，达到200时重新置零。再回到前面的norm_counter_weght，则显然当累计的帧数越多时，此值越小，也即lquantile的更新步长越小，含义和之前类似：先进行大范围搜索，当数据逐渐增多时，缩小步长以达到收敛的目的。

噪声更新的循环问题：函数中有一个参数为self->updates，在前两百帧，每执行一次噪声估计，此参数加一：
```
if (self->updates < END_STARTUP_LONG) {
    self->updates++;
}
```
当self->updates小于200时，也即处于前两百帧时，每一帧都更新系统分位数（也即结构体self->quantile）：
```
if (self->updates < END_STARTUP_LONG) {
    for (i = 0; i < self->magnLen; i++) {
        self->quantile[i] = expf(self->lquantile[offset + i]);
    }
    memcpy(noise, self->quantile, self->magnLen * sizeof(*noise));
} else {
    memcpy(noise, self->quantile, self->magnLen * sizeof(*noise));
}
```
当self->updates大于200时，每66帧才更新一次系统分位数：
```
if (self->counter[s] >= END_STARTUP_LONG) {
    self->counter[s] = 0;
    if (self->updates >= END_STARTUP_LONG) {
        for (i = 0; i < self->magnLen; i++) {
            self->quantile[i] = expf(self->lquantile[offset + i]);
        }
    }
}
```
66的来历：counter是一个长度为3初始值为[67,134,1]的数组，每一帧处理后三个值都加一。则每过66帧，就会有一个counter[i]的值大于或等于END_STARTUP_LONG（=200），将此counter[i]的值置零，然后更新系统分位数。

输出噪声noise：从上上个代码块中的内容可以看出，无论如何，输出的噪声估计都是系统分位值。

几个值得注意的点：
- 前面所说的分位数都指：分位点所对应的噪声幅度水平，而不是一个表示比例的数。lquantile和quantile是指数与对数关系
- 变量offset只是为了使得分位数不用设置为三维数组。如counter[1]所对应的分位数便是lquantile[128: 256]，则调整offset即可取出不同counter对应的分位值
- 最后输出的噪声noise只是对噪声的一个初始估计，后面还会更新
一个简单的噪声模型（simplified noise model）的计算

当处理的帧数小于50的时候，为了更加准确地估计噪声，webrtc使用已经处理过的帧（index记为num_analyzed_frames）建立了一个简单的噪声模型。在处理index为 50-num_analyzed_frames 的帧数据时，便在分位数估计的基础之上使用此噪声模型进行噪声的估计，下面叙述流程：

在前五十帧时，计算一些需要的参数：
```
//这里END_STARTUP_SHORT= 50
if (self->blockInd < END_STARTUP_SHORT) {
    //这里舍弃前kStartBand个频点，也即低频处的若干个频点不参与计算
    for (i = kStartBand; i < self->magnLen; i++) {
        sum_log_i += self->log_lut[i];
        sum_log_i_square += self->log_lut_sqr[i];
        sum_log_magn += lmagn[i];
        sum_log_i_log_magn += self->log_lut[i] * lmagn[i];
    }
}
```
然后计算白噪声：
```
//参数overdrive是衡量去噪水平的量，和初始设定的去噪等级有关
self->whiteNoiseLevel += sumMagn * self->normMagnLen * self->overdrive;
```
计算粉红噪声相关参数，一个分子参数为， $i$ 表示当前帧的index：

$pink\_noise\_numerator=\displaystyle{\sum_0^{i}\max {(\frac{sum\_log\_i\_square*sum\_log\_magn-sum\_log\_i*sum\_log\_o\_log\_mgan }{sum\_log\_i\_square*(129-5)-sum\_log\_i*sum\_log\_i},0)}}$

其中 $(129 ? 5)$ 指忽略前5个频点（下同），一个指数参数为：

$pink\_noise\_exp=\displaystyle{\sum_0^i\max{(\min(\frac{sum\_log\_i*sum\_log\_magn-(129-5)*sum\_log\_i\_log\_magn}{sum\_log\_i\_square*(129-5)-sum\_log\_i*sum\_log\_i},1),0)}}$

粉红噪声参数利用上面两个参数继续计算如下：

$parametric\_num=(num\_analyzed\_frames+1.0)\times\exp{^{\frac{pink\_noise\_numerator}{num\_analyzed\_frames+1}}}$

$parametric\_exp=\displaystyle{\frac{pink\_noise\_exp}{num\_analyzed\_frames+1}}$

然后利用粉红噪声参数和白噪声参数进行背景噪声估计：

$parametric\_noise\_spectrum[i]=\begin{cases}white\_noise\_level,&pink\_noise\_exp=0\\\displaystyle{\frac{parametric\_num}{2^{parametric\_exp\times \log_2(use\_band)}}},&else\end{cases}$

其中：

$use\_band=\begin{cases}5,&i<5\\i&5\le i<129\end{cases}$

之后，根据噪声模型结合分位数噪声估计的方法及逆行噪声估计，注意，这些操作只在前50帧进行，超过50帧时，只进行上文提到的分位数噪声估计。

代码为：
```
if (self->pinkNoiseExp == 0.f) {
    for (i = 0; i < self->magnLen; i++) {
        //如果粉红噪声为0，使用白噪声参数估计背景噪声参数
        self->parametricNoise[i] = self->whiteNoiseLevel;
        // Weight quantile noise with modeled noise.
        noise[i] *= (self->blockInd);
        tmpFloat2 = self->parametricNoise[i] * (END_STARTUP_SHORT - self->blockInd);
        noise[i] += tmpFloat2 * norm;
        noise[i] *= norm_end;
    }
} else {
    // 粉红噪声相关参数计算
    parametric_num = expf(self->pinkNoiseNumerator * norm);
    parametric_num *= (float) (self->blockInd + 1);
    parametric_exp = self->pinkNoiseExp * norm;
    for (i = 0; i < self->magnLen; i++) {
        //否则使用使用粉红噪声参数估计背景噪声参数
        float use_band = (float) (i < kStartBand ? kStartBand : i);
        self->parametricNoise[i] = parametric_num / powf(use_band, parametric_exp);
        // Weight quantile noise with modeled noise.
        noise[i] *= (self->blockInd);
        tmpFloat2 = self->parametricNoise[i] * (END_STARTUP_SHORT - self->blockInd);
        noise[i] += tmpFloat2 * norm;
        noise[i] *= norm_end;
    }
}
```

计算信噪比（SNR）

使用ComputeSnr()函数计算先验信噪比和后验信噪比，它们可以用来计算前文所述的LRT参数，进而计算语音存在概率，计算如下：

static void ComputeSnr(const NoiseSuppressionC *self,
                       const float *magn,
                       const float *noise,
                       float *snrLocPrior, float *logSnrLocPrior,
                       float *snrLocPost) {
    size_t i;

    for (i = 0; i < self->magnLen; i++) {
        // Previous post SNR.
        // Previous estimate: based on previous frame with gain filter.
        float previousEstimateStsa = (self->magnPrevAnalyze[i] * self->smooth[i]) / (self->noisePrev[i] + epsilon);
        // Post SNR.
        snrLocPost[i] = 0.f;
        if (magn[i] > noise[i]) {
            snrLocPost[i] = (magn[i] - noise[i]) / (noise[i] + epsilon);
        }
        // DD estimate is sum of two terms: current estimate and previous estimate.
        // Directed decision update of snrPrior.
        snrLocPrior[i] = 2.f * (
                DD_PR_SNR * previousEstimateStsa + (1.f - DD_PR_SNR) * snrLocPost[i]);
        logSnrLocPrior[i] = log1pf(snrLocPrior[i]);
    }  // End of loop over frequencies.
}

数据更新

首先使用FeatureUpdate()函数对频谱坦度、LRT和频谱模板相似度三个特征进行更新，分为两种情况：

a. 帧数小于200时：即还处于噪声阶段，关键之处在于初始噪声的估计，只有模板相似度会进行更新

b. 帧数大于200时：可能有语音存在，初始的噪声估计也已经完成，三个特征都进行更新
- 频谱平坦度计算
  
  计算原理见式 $(1 ? 9)$ ，在FeatureUpdate()函数中调用ComputeSpectralFlatness()函数计算频谱平坦度：
  
  $spectral\_tmp=\displaystyle{\frac{e^{avg\_spec\_flatness\_num}}{avg\_spect\_flatness\_denom}}\\ =\frac{e^{\frac{\sum_{i=0}^{128}\ln(signal\_spectrum[i])}{129}}}{\frac{signal\_spectral\_sum-signal\_spectrum[0]}{129}}=\frac{\prod_k|Y_k|^{\frac{1}{129}}}{\frac{1}{129}\sum_{k=0}^{129}Y_k}$
  
  然后进行时间平滑更新并储存至self->featureData[0]中，具体代码如下:
```
static void ComputeSpectralFlatness(NoiseSuppressionC *self,
                                    const float *magnIn, const float *logmagnIn) {
    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;
    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 += logmagnIn[i];
        } else {
            self->featureData[0] -= SPECT_FL_TAVG * self->featureData[0];
            return;
        }
    }
    // Normalize.
    avgSpectralFlatnessDen = avgSpectralFlatnessDen * self->normMagnLen;
    avgSpectralFlatnessNum = avgSpectralFlatnessNum * self->normMagnLen;

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

    // Time-avg update of spectral flatness feature.
    self->featureData[0] += SPECT_FL_TAVG * (spectralTmp - self->featureData[0]);
    // Done with flatness feature.
}
```
- 频谱模板相似度计算
  
  这个特征可以衡量输入频谱和算法学习到的噪声的差异，原理见 $(1 ? 10)$ 。代码中首先计算噪声以及输入信号的平均值和方差，接着计算输入信号和噪声的协方差，然后：
  
  $spectral\_diff=signal\_variance-\displaystyle{\frac{covariance^2}{noise\_variance+0.0001}=\sigma_Y-\frac{Cov(Y,N)^2}{\sigma_N}}$
  
  之后对其进行归一化，然后类似地，做时间平滑，最后储存在self->featureData[4]中，具体代码为：
```
static void ComputeSpectralDifference(NoiseSuppressionC *self,
                                      const float *magnIn) {
    // avgDiffNormMagn = var(magnIn) - cov(magnIn, magnAvgPause)^2 /
    // var(magnAvgPause)
    size_t i;
    float avgPause, avgMagn, covMagnPause, varPause, varMagn, avgDiffNormMagn;

    avgPause = 0;
    avgMagn = self->sumMagn;
    // Compute average quantities.
    for (i = 0; i < self->magnLen; i++) {
        // Conservative smooth noise spectrum from pause frames.
        avgPause += self->magnAvgPause[i];
    }
    avgPause *= self->normMagnLen;
    avgMagn *= self->normMagnLen;

    covMagnPause = 0;
    varPause = 0;
    varMagn = 0;
    // Compute variance and covariance quantities.
    for (i = 0; i < self->magnLen; i++) {
        const float avgPauseDiff = self->magnAvgPause[i] - avgPause;
        const float avgMagnDiff = magnIn[i] - avgMagn;
        covMagnPause += avgMagnDiff * avgPauseDiff;
        varPause += avgPauseDiff * avgPauseDiff;
        varMagn += avgMagnDiff * avgMagnDiff;
    }
    covMagnPause *= self->normMagnLen;
    varPause *= self->normMagnLen;
    varMagn *= self->normMagnLen;
    // Update of average magnitude spectrum.
    self->featureData[6] += self->signalEnergy;

    avgDiffNormMagn =
            varMagn - (covMagnPause * covMagnPause) / (varPause + epsilon);
    // Normalize and compute time-avg update of difference feature.
    avgDiffNormMagn = avgDiffNormMagn / (self->featureData[5] + epsilon);
    self->featureData[4] +=
            SPECT_DIFF_TAVG * (avgDiffNormMagn - self->featureData[4]);
}
```
- 更新直方图（histogram）
  
  直方图（histogram）中有参数的阈值和权重，每500帧提取一次，用于更新模型参数。计算分位两部分，当帧数不是500的整数倍时（对应代码中self->modelUpdatePars[3]>0），只对直方图内的参数进行更新；当帧数为500的整数倍时（对应代码中self->modelUpdatePars[3]==0），更新整个模型的参数，并将直方图中的参数重置。直方图更新时，会更新频谱坦度、LRT、模板相似度三个参数。
  
  参数更新调用的是FeatureParameterExtraction()函数，其有两个参数，第一个参数为主结构体，第二个参数为flag，表示是否到达500帧。其运算过程如下：
  
  ? a. 当flag为0时，只更新直方图本身，即更新LRT, flatness, difference各三个参数，分别储存在： self->featureData[3]，self->featureData[0]，self->featureData[4]中
  
  ? b. 当flag为1时，更新整个模型的参数：
  
  ? i) LRT特征：利用了平均的思想，计算如下
  
  ? $bin\_mid=(i+0.5)*kBinSizeLrt$
  
  ? $average=\displaystyle{\frac{\sum_{i=0}^9[Lrt\_histogram[i]\times(i+0.5)\times kBinSizeLrt]}{\sum_{i=0}^9 lrt_histogram[i]}}$
  
  ? $average\_squared=\displaystyle{\frac{\sum_{i=0}^{999}(Lrt\_histogram[i]\times[(i+0.5)\times kBinSizeLrt]^2)}{kFeatureUpdateWindowSize}}$
  
  ? $average\_compl=\displaystyle{\frac{\sum_{i=0}^{999}(Lrt\_histogram[i]\times(i+0.5)\times kBinSizeLrt)}{kFeatureUpdateWindowSize}}$
  
  ? 然后计算LRT的波动程度：
  
  ? $low\_lrt\_fluctuations=average\_squarerd-average*average\_compl$
  
  ? 然后根据此波动程度，更新模型的LRT阈值：
  
  ? $prior\_model\_lrt = \begin{cases}kMaxLrt=1.0,&low\_lrt\_fluctuation=true(波动较小类似噪声)\\1.2\times average,&else(限定1.2*average\in[0.2,1])\end{cases}$
  
  ? ii)频谱坦度特征：对于频谱平坦度特征，先计算其两个主要峰值以及峰值的位置，如果最高峰和第二高峰相邻，并且第二峰峰值超过最高峰峰值的一半，则认为最高峰的位置在这两个峰之间，峰值为这两个峰的峰值之和。之后，计算权值（flatness在三个特征中所占的权重），然后更新整个模型的 flatness参数
  
  ? iii)频谱相似度特征：与平坦度特征类似，都是求峰值
  
  ?
  
  ? 之后计算每个特征的权重，然后将直方图重置。
  
  FeatureParameterExtraction()代码如下：
```
static void FeatureParameterExtraction(NoiseSuppressionC *self, int flag) {
    int i, useFeatureSpecFlat, useFeatureSpecDiff, numHistLrt;
    int maxPeak1, maxPeak2;
    int weightPeak1SpecFlat, weightPeak2SpecFlat, weightPeak1SpecDiff,
            weightPeak2SpecDiff;

    float binMid, featureSum;
    float posPeak1SpecFlat, posPeak2SpecFlat, posPeak1SpecDiff, posPeak2SpecDiff;
    float fluctLrt, avgHistLrt, avgSquareHistLrt, avgHistLrtCompl;

    // 3 features: LRT, flatness, difference.
    // lrt_feature = self->featureData[3];
    // flat_feature = self->featureData[0];
    // diff_feature = self->featureData[4];

    // Update histograms.
    if (flag == 0) {
        // LRT
        if ((self->featureData[3] <
             HIST_PAR_EST * self->featureExtractionParams.binSizeLrt) &&
            (self->featureData[3] >= 0.0)) {
            i = (int) (self->featureData[3] /
                       self->featureExtractionParams.binSizeLrt);
            self->histLrt[i]++;
        }
        // Spectral flatness.
        if ((self->featureData[0] <
             HIST_PAR_EST * self->featureExtractionParams.binSizeSpecFlat) &&
            (self->featureData[0] >= 0.0)) {
            i = (int) (self->featureData[0] /
                       self->featureExtractionParams.binSizeSpecFlat);
            self->histSpecFlat[i]++;
        }
        // Spectral difference.
        if ((self->featureData[4] <
             HIST_PAR_EST * self->featureExtractionParams.binSizeSpecDiff) &&
            (self->featureData[4] >= 0.0)) {
            i = (int) (self->featureData[4] /
                       self->featureExtractionParams.binSizeSpecDiff);
            self->histSpecDiff[i]++;
        }
    }

    // Extract parameters for speech/noise probability.
    if (flag == 1) {
        // LRT feature: compute the average over
        // self->featureExtractionParams.rangeAvgHistLrt.
        avgHistLrt = 0;
        avgHistLrtCompl = 0;
        avgSquareHistLrt = 0;
        numHistLrt = 0;
        for (i = 0; i < HIST_PAR_EST; i++) {
            binMid = ((float) i + 0.5f) * self->featureExtractionParams.binSizeLrt;
            if (binMid <= self->featureExtractionParams.rangeAvgHistLrt) {
                avgHistLrt += self->histLrt[i] * binMid;
                numHistLrt += self->histLrt[i];
            }
            avgSquareHistLrt += self->histLrt[i] * binMid * binMid;
            avgHistLrtCompl += self->histLrt[i] * binMid;
        }
        if (numHistLrt > 0) {
            avgHistLrt = avgHistLrt / ((float) numHistLrt);
        }
        avgHistLrtCompl = avgHistLrtCompl / ((float) self->modelUpdatePars[1]);
        avgSquareHistLrt = avgSquareHistLrt / ((float) self->modelUpdatePars[1]);
        fluctLrt = avgSquareHistLrt - avgHistLrt * avgHistLrtCompl;
        // Get threshold for LRT feature.
        if (fluctLrt < self->featureExtractionParams.thresFluctLrt) {
            // Very low fluctuation, so likely noise.
            self->priorModelPars[0] = self->featureExtractionParams.maxLrt;
        } else {
            self->priorModelPars[0] =
                    self->featureExtractionParams.factor1ModelPars * avgHistLrt;
            // Check if value is within min/max range.
            if (self->priorModelPars[0] < self->featureExtractionParams.minLrt) {
                self->priorModelPars[0] = self->featureExtractionParams.minLrt;
            }
            if (self->priorModelPars[0] > self->featureExtractionParams.maxLrt) {
                self->priorModelPars[0] = self->featureExtractionParams.maxLrt;
            }
        }
        // Done with LRT feature.

        // For spectral flatness and spectral difference: compute the main peaks of
        // histogram.
        maxPeak1 = 0;
        maxPeak2 = 0;
        posPeak1SpecFlat = 0;
        posPeak2SpecFlat = 0;
        weightPeak1SpecFlat = 0;
        weightPeak2SpecFlat = 0;

        // Peaks for flatness.
        for (i = 0; i < HIST_PAR_EST; i++) {
            binMid = (i + 0.5f) * self->featureExtractionParams.binSizeSpecFlat;
            if (self->histSpecFlat[i] > maxPeak1) {
                // Found new "first" peak.
                maxPeak2 = maxPeak1;
                weightPeak2SpecFlat = weightPeak1SpecFlat;
                posPeak2SpecFlat = posPeak1SpecFlat;

                maxPeak1 = self->histSpecFlat[i];
                weightPeak1SpecFlat = self->histSpecFlat[i];
                posPeak1SpecFlat = binMid;
            } else if (self->histSpecFlat[i] > maxPeak2) {
                // Found new "second" peak.
                maxPeak2 = self->histSpecFlat[i];
                weightPeak2SpecFlat = self->histSpecFlat[i];
                posPeak2SpecFlat = binMid;
            }
        }

        // Compute two peaks for spectral difference.
        maxPeak1 = 0;
        maxPeak2 = 0;
        posPeak1SpecDiff = 0;
        posPeak2SpecDiff = 0;
        weightPeak1SpecDiff = 0;
        weightPeak2SpecDiff = 0;
        // Peaks for spectral difference.
        for (i = 0; i < HIST_PAR_EST; i++) {
            binMid =
                    ((float) i + 0.5f) * self->featureExtractionParams.binSizeSpecDiff;
            if (self->histSpecDiff[i] > maxPeak1) {
                // Found new "first" peak.
                maxPeak2 = maxPeak1;
                weightPeak2SpecDiff = weightPeak1SpecDiff;
                posPeak2SpecDiff = posPeak1SpecDiff;

                maxPeak1 = self->histSpecDiff[i];
                weightPeak1SpecDiff = self->histSpecDiff[i];
                posPeak1SpecDiff = binMid;
            } else if (self->histSpecDiff[i] > maxPeak2) {
                // Found new "second" peak.
                maxPeak2 = self->histSpecDiff[i];
                weightPeak2SpecDiff = self->histSpecDiff[i];
                posPeak2SpecDiff = binMid;
            }
        }

        // For spectrum flatness feature.
        useFeatureSpecFlat = 1;
        // Merge the two peaks if they are close.
        if ((fabsf(posPeak2SpecFlat - posPeak1SpecFlat) <
             self->featureExtractionParams.limitPeakSpacingSpecFlat) &&
            (weightPeak2SpecFlat >
             self->featureExtractionParams.limitPeakWeightsSpecFlat *
             weightPeak1SpecFlat)) {
            weightPeak1SpecFlat += weightPeak2SpecFlat;
            posPeak1SpecFlat = 0.5f * (posPeak1SpecFlat + posPeak2SpecFlat);
        }
        // Reject if weight of peaks is not large enough, or peak value too small.
        if (weightPeak1SpecFlat <
            self->featureExtractionParams.thresWeightSpecFlat ||
            posPeak1SpecFlat < self->featureExtractionParams.thresPosSpecFlat) {
            useFeatureSpecFlat = 0;
        }
        // If selected, get the threshold.
        if (useFeatureSpecFlat == 1) {
            // Compute the threshold.
            self->priorModelPars[1] =
                    self->featureExtractionParams.factor2ModelPars * posPeak1SpecFlat;
            // Check if value is within min/max range.
            if (self->priorModelPars[1] < self->featureExtractionParams.minSpecFlat) {
                self->priorModelPars[1] = self->featureExtractionParams.minSpecFlat;
            }
            if (self->priorModelPars[1] > self->featureExtractionParams.maxSpecFlat) {
                self->priorModelPars[1] = self->featureExtractionParams.maxSpecFlat;
            }
        }
        // Done with flatness feature.

        // For template feature.
        useFeatureSpecDiff = 1;
        // Merge the two peaks if they are close.
        if ((fabsf(posPeak2SpecDiff - posPeak1SpecDiff) <
             self->featureExtractionParams.limitPeakSpacingSpecDiff) &&
            (weightPeak2SpecDiff >
             self->featureExtractionParams.limitPeakWeightsSpecDiff *
             weightPeak1SpecDiff)) {
            weightPeak1SpecDiff += weightPeak2SpecDiff;
            posPeak1SpecDiff = 0.5f * (posPeak1SpecDiff + posPeak2SpecDiff);
        }
        // Get the threshold value.
        self->priorModelPars[3] =
                self->featureExtractionParams.factor1ModelPars * posPeak1SpecDiff;
        // Reject if weight of peaks is not large enough.
        if (weightPeak1SpecDiff <
            self->featureExtractionParams.thresWeightSpecDiff) {
            useFeatureSpecDiff = 0;
        }
        // Check if value is within min/max range.
        if (self->priorModelPars[3] < self->featureExtractionParams.minSpecDiff) {
            self->priorModelPars[3] = self->featureExtractionParams.minSpecDiff;
        }
        if (self->priorModelPars[3] > self->featureExtractionParams.maxSpecDiff) {
            self->priorModelPars[3] = self->featureExtractionParams.maxSpecDiff;
        }
        // Done with spectral difference feature.

        // Don't use template feature if fluctuation of LRT feature is very low:
        // most likely just noise state.
        if (fluctLrt < self->featureExtractionParams.thresFluctLrt) {
            useFeatureSpecDiff = 0;
        }

        // Select the weights between the features.
        // self->priorModelPars[4] is weight for LRT: always selected.
        // self->priorModelPars[5] is weight for spectral flatness.
        // self->priorModelPars[6] is weight for spectral difference.
        featureSum = (float) (1 + useFeatureSpecFlat + useFeatureSpecDiff);
        self->priorModelPars[4] = 1.f / featureSum;
        self->priorModelPars[5] = ((float) useFeatureSpecFlat) * self->priorModelPars[4];
        self->priorModelPars[6] = ((float) useFeatureSpecDiff) * self->priorModelPars[4];

        // Set hists to zero for next update.
        if (self->modelUpdatePars[0] >= 1) {
            for (i = 0; i < HIST_PAR_EST; i++) {
                self->histLrt[i] = 0;
                self->histSpecFlat[i] = 0;
                self->histSpecDiff[i] = 0;
            }
        }
    }  // End of flag == 1.
}
```
计算映射函数和语音存在概率

由SpeechNoiseProb()进行计算，先进行映射的计算然后再计算语音存在概率
- 映射
  
  语音存在概率时利用映射函数对特征映射后得到的，典型的映射函数有：
  
  $\sinh x=\displaystyle{\frac{e^x-e^{-x}}{2}},\ \ \ \cosh x=\displaystyle{\frac{e^x+e^{-x}}{2}},\ \ \ \tanh x=\displaystyle{\frac{\sinh x}{cosh x}=\frac{e^x-e^{-x}}{e^x+e^{-x}}}$
  
  LRT的映射为：
  
  $indicator0=0.5\times[tanh(width\_prior\times(model.lrt-prior\_model.lrt))+1]$
  
  flattness的映射为：
  
  $indicator1=0.5\times[tanh(1\times width\_prior\times(prior\_model.flatness\_threshold-model.spectral\_flatness))+1]$
  
  频谱相似度的映射为：
  
  $indicator2=0.5\times[tanh(width\_prior\times(model.spectral\_diff-prior\_model.template\_diff\_threshold))+1]$
  
  最后将不同的映射加权求和：
  
  $ind\_prior=prior\_model.lrt\_weighting\times indicator0\\+prior\_model.flatness\_weighing\times indicator1\\ +prior\_model.difference\_weighting\times indicator2\\ =1\times indicator0+0\times indicator1+0\times indicator2$
  
  flatness和difference的权重实际为0，这也就是说，语音存在概率的判断完全依赖于LRT
- 语音存在概率计算
  
  $prior\_speech\_prob\_=prior\_speech\_prob\_+0.1\times(ind\_prior-prior\_speech\_prob\_)$
  
  $prior\_speech\_prob\_=\max(\min(prior\_speech\_prob\_,1),0.01)$
  
  $gain\_prior=\displaystyle{\frac{1.0-prior\_speech\_prob\_}{prior\_speech\_prob\_}}$
  
  $inv\_lrt[i]=e^{model.avg\_log\_lrt[i]}$
  
  $speech\_probability\_[i]=\displaystyle{\frac{1}{1.0+gain\_prior\times inv\_lrt[i]}}$
噪声估计：根据语音存在概率对噪声进行平滑更新，具体计算由UpdateNoiseEstimate()函数完成，原理在前文“相关算法公式分析”有介绍，即：

$|\hat{N}_k(m)|=\gamma_n |\hat{N}_k(m-1)|+(1-\gamma_n)[P(H_1|Y_k(m))\hat{N}_k(m-1)+P(H_0|Y_k(m))Y_k(m)]$

具体计算为：

$noise\_update\_tmp=\gamma \times prev\_noise\_spectrum\_[i]+(1-\gamma)\{(1-speech\_probability[i])\\ \times signal\_spectrum[i] +speech\_probability[i]\times prev\_noise\_spectrum[i]\}$

WebRTC中的 $\gamma$ 取为0.9，上市表明，当当前帧为语音时，主要使用过去帧的噪声来估计当前帧噪声

之后更新 $\gamma$ ：

$\gamma_{old}=kNoiseUpdate=0.9$

$\gamma=\begin{cases}0.99,&prob\_speech>k\_ProbRange\\kNoiseUpdate,&else\end{cases}$

当语音概率小于0.2时：

$conservative\_noise\_spectrum\_[i]+=0.05\times(signal\_spectrum[i]-conservative\_noise\_spectrum\_[i])$

然后更新噪声频谱，当 $\gamma=\gamma_{dd}=0.9$ 时：

$noise\_spectrum\_[i]=noise\_update\_tmp=\gamma_{old}\times prev\_noise\_spectrum\_[i]+(1-\gamma_{old})\times \\ \{(1-speech\_probability[i])\times signal\_spectrum[i]+speech\_probability[i]\times prev\_noise\_spectrum\_[i]\}$

当 $\gamma=0.99\neq\gamma_{old}$ 时：

$noise\_spectrum\_[i]=noise\_update\_tmp=\gamma\times prev\_noise\_spectrum\_[i]+(1-\gamma)\times \\ \{(1-speech\_probability[i])\times signal\_spectrum[i]+speech\_probability[i]\times prev\_noise\_spectrum\_[i]\}$

且 $noise\_spectrum\_[i]=\min(noise\_spectrum\_[i],noise\_update\_tmp)$

代码如下：
```
static void UpdateNoiseEstimate(NoiseSuppressionC *self,
                                const float *magn,
                                float *noise) {
    size_t i;
    float probSpeech, probNonSpeech;
    // Time-avg parameter for noise update.
    float gammaNoiseTmp = NOISE_UPDATE;
    float gammaNoiseOld;
    float noiseUpdateTmp;

    for (i = 0; i < self->magnLen; i++) {
        probSpeech = self->speechProb[i];
        probNonSpeech = 1.f - probSpeech;
        // Temporary noise update:
        // Use it for speech frames if update value is less than previous.
        noiseUpdateTmp = gammaNoiseTmp * self->noisePrev[i] +
                         (1.f - gammaNoiseTmp) * (probNonSpeech * magn[i] +
                                                  probSpeech * self->noisePrev[i]);
        // Time-constant based on speech/noise state.
        gammaNoiseOld = gammaNoiseTmp;
        gammaNoiseTmp = NOISE_UPDATE;
        // Increase gamma (i.e., less noise update) for frame likely to be speech.
        if (probSpeech > PROB_RANGE) {
            gammaNoiseTmp = SPEECH_UPDATE;
        }
        // Conservative noise update.
        if (probSpeech < PROB_RANGE) {
            self->magnAvgPause[i] += GAMMA_PAUSE * (magn[i] - self->magnAvgPause[i]);
        }
        // Noise update.
        if (gammaNoiseTmp == gammaNoiseOld) {
            noise[i] = noiseUpdateTmp;
        } else {
            noise[i] = gammaNoiseTmp * self->noisePrev[i] +
                       (1.f - gammaNoiseTmp) * (probNonSpeech * magn[i] +
                                                probSpeech * self->noisePrev[i]);
            // Allow for noise update downwards:
            // If noise update decreases the noise, it is safe, so allow it to
            // happen.
            if (noiseUpdateTmp < noise[i]) {
                noise[i] = noiseUpdateTmp;
            }
        }
    }  // End of freq loop.
}
```

噪声抑制

当完成噪声估计后，便可以使用维纳滤波对噪声进行抑制，进而实现降噪，这一部分主要是由WebRtcNs_Process()函数调用WebRtcNs_ProcessCore()函数实现的，主要步骤如下：

数据形成

每次将当前10ms对应的160个数据点和前一帧最后的96个数据点拼起来，形成长度为256的当前帧（为了方便做FFT）
加窗并在滤波前计算能量
FFT
计算信号幅度谱
计算维纳滤波器增益，下面前三步主要由ComputeDdBasedWienerFilter()实现
- 根据过去帧计算
  
  $prev\_tsa=\displaystyle{\frac{spectrum\_prev\_process\_[i]}{prev\_noise\_spectrum[i]\times filter\_[i]}}$
- 当前帧数据的计算：当信号幅度谱大于噪声的时。计算一个先验SNR
  
  $current\_tsa=\begin{cases}\displaystyle{\frac{signal\_spectrum[i]}{noise\_spectrum[i]}-1},&signal\_spectrum[i]>noise\_spectrum[i]\\ 0,&else\end{cases}$
- 根据过去和当前的来计算滤波器增益
  
  $filter\_[i]=\displaystyle{\frac{snr\_prior}{suppresion\_params\_.over\_subtraction\_factor+snr\_prior}}$
  
  $filter\_[i]=\max(\min(filter\_[i],1.0),suppression\_params\_.minimum\_attenuating\_gain)$
- 前50帧的情况：当帧数小于50时，认为还处于噪声阶段。此时，维纳滤波器是通过结合当前帧数据以及估计的噪声模型来计算的，具体计算如下：
  
  $filter\_initial=\displaystyle{\frac{\sum_0^{50}signal\_spectrum[i]-suppression\_params\_over.over\_subtraction\_factor\times parametric\_noise\_spectrum[i]}{\sum_0^{50}signal\_spectrum[i]}}$
  
  $filter\_[i]=\displaystyle{\frac{filter\_[i]\times num\_analyzed\_frames+filter\_initial\times(50-num\_analyzed\_frames)}{50}}$

进行维纳滤波

将频谱和维纳滤波器相乘

self->smooth[i] = theFilter[i];
real[i] *= self->smooth[i];
imag[i] *= self->smooth[i];

时域处理
- 通过IFFT将数据变换回时域
- 计算一个调整因子（scale factor），这一步只在200帧时候进行
  
  首先计算滤波后的能量：
  
  $energy\_after\_filtering=\sum_{i=0}^{256}extended\_frame[i]^2$
  
  然后加上要进行IFFT的窗（即FFT窗的逆窗），然后根据降噪的效果，对尺度因子进行计算：
  
  如果处于前200帧，则 $gain\_adjustment[ch]=1.0$ ，否则：
  
  $gain=\displaystyle{\sqrt{\frac{energy\_after\_filtering}{energy\_before\_filtering+1.0}}}$
  
  $scale\_factor1=\begin{cases}1.0,&gain\le kBLim=0.5\\ \displaystyle{\frac{1}{gain}},&gain\times[1+1.3\times(gain-kBLim)]>1,\\1+1.3\times(gain-kBLim),&else\end{cases}$
  
  $scale\_factor2=\begin{cases}1-0.3\times[kBLim-\max(gain,suppression\_params\_.minimum\_attenuating\_gain)],&gain<kBLim\\ 1.0,&else\end{cases}$
  
  $gain\_adjustments=prior\_speech\_probability\times scale\_factor1+\\(1-prior\_speech\_probability)\times scale\_factor2$
  
  $gain\_adjustments$ 即为factor,代码如下：
```
// Scale factor: only do it after END_STARTUP_LONG time.
    factor = 1.f;

    if (self->gainmap == 1 && self->blockInd > END_STARTUP_LONG) {
        factor1 = 1.f;
        factor2 = 1.f;
        energy2 = Energy(winData, self->anaLen);
        gain = sqrtf(energy2 / (energy1 + epsilon) + epsilon_squ);

        // Scaling for new version.
        if (gain > B_LIM) {
            factor1 = 1.f + 1.3f * (gain - B_LIM);
            if (gain * factor1 > 1.f) {
                factor1 = 1.f / gain;
            }
        }
        if (gain < B_LIM) {
            // Don't reduce scale too much for pause regions:
            // attenuation here should be controlled by flooring.
            if (gain <= self->denoiseBound) {
                gain = self->denoiseBound;
            }
            factor2 = 1.f - 0.3f * (B_LIM - gain);
        }
        // Combine both scales with speech/noise prob:
        // note prior (priorSpeechProb) is not frequency dependent.
        factor = self->priorSpeechProb * factor1 +
                 (1.f - self->priorSpeechProb) * factor2;
    }  // Out of self->gainmap == 1.
```
- synthesis：将factor应用于此步骤
```
for (i = 0; i < self->anaLen; i++) {
        self->syntBuf[i] += factor * winData[i] * self->window[i];
    }
```
- 高频带增益的处理
  
  通常来讲不必要，略。