1.基本原理：
采用维纳滤波器抑制估计出来的噪声
时域与频域状态下：Y = N + S (Y:原始信号,S:纯净声音,N:噪音)
根据中心极限定义，一般认为噪声和语音分布服从均值为0，方差为ui的正态分布
中心思想：从Y中估计噪声N，然后抑制N以得到语音
估计噪声方法： 对似然比函数进行改进，将多个语音/噪声分类特征合并到一个模型中形成一个多特征综合概率密度函数，对输入的每帧频谱进行分析。其可以有效抑制风扇/办公设备等噪声。
2.处理过程
1）加窗
for (i = 0; i < inst->anaLen; i++) {
winData[i] = inst->window[i] * inst->dataBuf[i];
energy1 += winData[i] * winData[i];
}采用hanning窗
汉宁窗公式如下：
特征：升余弦函数，消去高频干扰和漏能，适用于非周期性的连续信号
图特征如下：
2） 傅里叶变换
采用实傅里叶变换
WebRtc_rdft(inst->anaLen, 1, winData, inst->ip, inst->wfft);将原始数据进行时域到频域的变换，此处采用蝶形算法进行快速傅里叶变换
得到实部虚部与能量值
magn = sqrt(real*real+imag*imag)
imag[inst->magnLen - 1] = 0; //虚部
real[inst->magnLen - 1] = winData[1]; //实部
magn[inst->magnLen - 1] = (float)(fabs(real[inst->magnLen - 1]) + 1.0f); //增强幅度fTmp = real[i] * real[i];
fTmp += imag[i] * imag[i];
signalEnergy += fTmp;
magn[i] = ((float)sqrt(fTmp)) + 1.0f3）计算频谱平坦度
频谱平坦度介绍：
属于三个特征集中的一个特征，假设语音比噪音有更多的谐波，语音频谱往往在基频和谐波中出现峰值，噪音则相对比较平坦，即平坦度越高，噪音概率越高。
公式
几何平均值/算数平均值
代码为加速计算实现过程如下，采用取Log操作将几何平均值的计算转为加法
代码如下：
为了使特征差异不明显，对结果数据进行了取历史数据的平滑操作
void WebRtcNs_ComputeSpectralFlatness(NSinst_t* inst, float* magnIn) {
int i;
int shiftLP = 1; //option to remove first bin(s) from spectral measures
float avgSpectralFlatnessNum, avgSpectralFlatnessDen, spectralTmp;
// comute spectral measures
// for flatness 跳过第一个频点，即直流频点Den是denominator（分母）的缩写，avgSpectralFlatnessDen是上述公式分母计算用到的
avgSpectralFlatnessNum = 0.0;
avgSpectralFlatnessDen = inst->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
//TVAG是time-average的缩写，对于能量出现异常的处理。利用前一次平坦度直接取平均返回。 ？？
//修正：将本次结果处理为0后返回
for (i = shiftLP; i < inst->magnLen; i++) {
if (magnIn[i] > 0.0) {
avgSpectralFlatnessNum += (float)log(magnIn[i]);
} else {
inst->featureData[0] -= SPECT_FL_TAVG * inst->featureData[0];
return;
}
}
//normalize
avgSpectralFlatnessDen = avgSpectralFlatnessDen / inst->magnLen;
avgSpectralFlatnessNum = avgSpectralFlatnessNum / inst->magnLen;
//ratio and inverse log: check for case of log(0)
spectralTmp = (float)exp(avgSpectralFlatnessNum) / avgSpectralFlatnessDen;
//time-avg update of spectral flatness feature
inst->featureData[0] += SPECT_FL_TAVG * (spectralTmp - inst->featureData[0]);
// done with flatness feature
}
4）噪声初步估计
中心思想：采用分位数估计法
分位数估计的基本思想：通过选择数列中的特定值来进行类别划分
假设噪声在整个频谱中属于能量值较低的部分，对能量值进行排序后取较小的部分
为了方便计算采用对数分位数估计法。
噪声估计过程如下：
调用一次处理函数算一帧，代码中用blockInd标记
分位数估计过程如下：
代码实现过程如下：
对原始能量值取对数，进行对数分位数估计,lquantile存储计算的对数结果。
density为概率密度， 用来更新系数delta
基本思想：设定log（噪声）初始化为8.0 （log(mag)频谱最大值为10），设定0.3的概率密度，对lquantile进行加操作更新，
当大于原始频谱的时候则进行减操作更新，当估计与原始频谱值差距减小的时候则增加概率密度，减少每次对lq的删改操作。
使lq无限逼近原始频谱。
void WebRtcNs_NoiseEstimation(NSinst_t* inst, float* magn, float* noise) {
int i, s, offset;
float lmagn[HALF_ANAL_BLOCKL], delta;
//updates帧计数
if (inst->updates < END_STARTUP_LONG) {
inst->updates++;
}
for (i = 0; i < inst->magnLen; i++) {
lmagn[i] = (float)log(magn[i]);
}
// loop over simultaneous estimates 三阶段估计 计算三次
for (s = 0; s < SIMULT; s++) {
offset = s * inst->magnLen;
// newquantest(...)
for (i = 0; i < inst->magnLen; i++) {
// compute delta
if (inst->density[offset + i] > 1.0) {
delta = FACTOR * (float)1.0 / inst->density[offset + i];
} else {
delta = FACTOR;
}
// update log quantile estimate 自适应中位数估计
if (lmagn[i] > inst->lquantile[offset + i]) {
inst->lquantile[offset + i] += QUANTILE * delta
/ (float)(inst->counter[s] + 1);
} else {
inst->lquantile[offset + i] -= ((float)1.0 - QUANTILE) * delta
/ (float)(inst->counter[s] + 1);
}
// update density estimate
if (fabs(lmagn[i] - inst->lquantile[offset + i]) < WIDTH) {
inst->density[offset + i] = ((float)inst->counter[s] * inst->density[offset
+ i] + (float)1.0 / ((float)2.0 * WIDTH)) / (float)(inst->counter[s] + 1);
}
} // end loop over magnitude spectrum
//200帧以上的处理方法
if (inst->counter[s] >= END_STARTUP_LONG) {
inst->counter[s] = 0;
if (inst->updates >= END_STARTUP_LONG) {
for (i = 0; i < inst->magnLen; i++) {
inst->quantile[i] = (float)exp(inst->lquantile[offset + i]);
}
}
}
inst->counter[s]++;
} // end loop over simultaneous estimates
// Sequentially update the noise during startup
//200帧以内的处理方法 每次更新
if (inst->updates < END_STARTUP_LONG) {
// Use the last "s" to get noise during startup that differ from zero.
for (i = 0; i < inst->magnLen; i++) {
inst->quantile[i] = (float)exp(inst->lquantile[offset + i]);
}
}
for (i = 0; i < inst->magnLen; i++) {
noise[i] = inst->quantile[i];
}
}代码实现过程
前200帧时每次更新lq，之后当差距为67的counter中有一个满足200则更新(counter到达200后置零)。
前50帧的模型叠加
计算方法为将分位数估计的噪声与模型估计的噪声进行加权平均
代码如下:
if (inst->blockInd < END_STARTUP_SHORT) {
// Estimate White noise 总能量平均值
inst->whiteNoiseLevel += sumMagn / ((float)inst->magnLen) * inst->overdrive;
// Estimate Pink noise parameters
tmpFloat1 = sum_log_i_square * ((float)(inst->magnLen - kStartBand));
tmpFloat1 -= (sum_log_i * sum_log_i);
tmpFloat2 = (sum_log_i_square * sum_log_magn - sum_log_i * sum_log_i_log_magn);
tmpFloat3 = tmpFloat2 / tmpFloat1;
// Constrain the estimated spectrum to be positive
if (tmpFloat3 < 0.0f) {
tmpFloat3 = 0.0f;
}
inst->pinkNoiseNumerator += tmpFloat3;
tmpFloat2 = (sum_log_i * sum_log_magn);
tmpFloat2 -= ((float)(inst->magnLen - kStartBand)) * sum_log_i_log_magn;
tmpFloat3 = tmpFloat2 / tmpFloat1;
// Constrain the pink noise power to be in the interval [0, 1];
if (tmpFloat3 < 0.0f) {
tmpFloat3 = 0.0f;
}
if (tmpFloat3 > 1.0f) {
tmpFloat3 = 1.0f;
}
inst->pinkNoiseExp += tmpFloat3;
// Calculate frequency independent parts of parametric noise estimate.
if (inst->pinkNoiseExp == 0.0f) {
// Use white noise estimate
parametric_noise = inst->whiteNoiseLevel;
} else {
// Use pink noise estimate
parametric_num = exp(inst->pinkNoiseNumerator / (float)(inst->blockInd + 1));
parametric_num *= (float)(inst->blockInd + 1);
parametric_exp = inst->pinkNoiseExp / (float)(inst->blockInd + 1);
parametric_noise = parametric_num / pow((float)kStartBand, parametric_exp);
}
for (i = 0; i < inst->magnLen; i++) {
// Estimate the background noise using the white and pink noise parameters
if ((inst->pinkNoiseExp > 0.0f) && (i >= kStartBand)) {
// Use pink noise estimate
parametric_noise = parametric_num / pow((float)i, parametric_exp);
}
theFilterTmp[i] = (inst->initMagnEst[i] - inst->overdrive * parametric_noise);
theFilterTmp[i] /= (inst->initMagnEst[i] + (float)0.0001);
// Weight quantile noise with modeled noise
noise[i] *= (inst->blockInd);
tmpFloat2 = parametric_noise * (END_STARTUP_SHORT - inst->blockInd);
noise[i] += (tmpFloat2 / (float)(inst->blockInd + 1));
noise[i] /= END_STARTUP_SHORT;
}
}
代码推导过程如下
上图与下图对应关系
num为pinkNoiseNumerator
exp为pinkNoiseExp
上图的noise是分位数估计值,p_noise是模型估计值生成逐步衰减的粉红噪声模型，前50帧估计中随着帧数的增大 分位数估计占比逐渐明显
5）snr初步估计
通过上面计算的noise进行初步snr(sound noise ratio)估计
生成先验snr和后验snr
由于声音信号为线性关系，存在Y = N+S(Y:原始信号 N:噪声信号 S:语音信号)
先验snr在应用过程中采用上一帧的先验与本次的后验平滑值
关系和表示如下
代码如下
for (i = 0; i < inst->magnLen; i++) {
// post snr
snrLocPost[i] = (float)0.0;
if (magn[i] > noise[i]) {
snrLocPost[i] = magn[i] / (noise[i] + (float)0.0001) - (float)1.0;
}
// previous post snr
// previous estimate: based on previous frame with gain filter
previousEstimateStsa[i] = inst->magnPrev[i] / (inst->noisePrev[i] + (float)0.0001)
* (inst->smooth[i]);
// DD estimate is sum of two terms: current estimate and previous estimate
// directed decision update of snrPrior
snrLocPrior[i] = DD_PR_SNR * previousEstimateStsa[i] + ((float)1.0 - DD_PR_SNR)* snrLocPost[i];
// post and prior snr needed for step 2
} // end of loop over freqs
6）差异度计算
差异度为特征集中的其中一个
假设噪声频谱比较平坦，对应的差异度较低，用来衡量噪声的可能性
代码如下：
void WebRtcNs_ComputeSpectralDifference(NSinst_t* inst, float* magnIn) {
// avgDiffNormMagn = var(magnIn) - cov(magnIn, magnAvgPause)^2 / var(magnAvgPause)
int i;
float avgPause, avgMagn, covMagnPause, varPause, varMagn, avgDiffNormMagn;
avgPause = 0.0;
avgMagn = inst->sumMagn;
// compute average quantities
for (i = 0; i < inst->magnLen; i++) {
//conservative smooth noise spectrum from pause frames
avgPause += inst->magnAvgPause[i];
}
avgPause = avgPause / ((float)inst->magnLen);
avgMagn = avgMagn / ((float)inst->magnLen);
covMagnPause = 0.0;
varPause = 0.0;
varMagn = 0.0;
// compute variance and covariance quantities
//covMagnPause 协方差 varPause varMagn 各自的方差
for (i = 0; i < inst->magnLen; i++) {
covMagnPause += (magnIn[i] - avgMagn) * (inst->magnAvgPause[i] - avgPause);
varPause += (inst->magnAvgPause[i] - avgPause) * (inst->magnAvgPause[i] - avgPause);
varMagn += (magnIn[i] - avgMagn) * (magnIn[i] - avgMagn);
}
covMagnPause = covMagnPause / ((float)inst->magnLen);
varPause = varPause / ((float)inst->magnLen);
varMagn = varMagn / ((float)inst->magnLen);
// update of average magnitude spectrum
inst->featureData[6] += inst->signalEnergy;
//通过决定系数计算差异性
avgDiffNormMagn = varMagn - (covMagnPause * covMagnPause) / (varPause + (float)0.0001);
// normalize and compute time-avg update of difference feature
avgDiffNormMagn = (float)(avgDiffNormMagn / (inst->featureData[5] + (float)0.0001));
inst->featureData[4] += SPECT_DIFF_TAVG * (avgDiffNormMagn - inst->featureData[4]);
}代码对 magnAvgPause和MagIn进行了联合计算，MagIn为原始频谱 magnAvgPause为在声音概率较低的时候进行更新的参数
//语音概率较小 更新magnAvgPause
if (probSpeech < PROB_RANGE) {
inst->magnAvgPause[i] += GAMMA_PAUSE * (magn[i] - inst->magnAvgPause[i]);
}代码公式推导如下
代码中分别计算了两个变量的方差和协方差，从公式可以看出是决定系数的计算方法
决定系数介绍：在上面的线性回归过程中用来衡量数据的拟合程度，(类似的相关系数用来衡量非回归场景下的拟合程度)。
7）特征值筛选与权重计算
过程如下：
代码如下：
(flag=0)500帧内对特征值(lrt均值，平坦度，差异度)进行直方图分布计算(即将特征值(feature)除以某个特定的参数(bin_size)，就能得到声音频谱的整个特征值在某个固定范围内的分布情况(x轴:i))
i = feature/bin_size ==> feature = bin_size*i
(flag=1)达到500帧后对特征值进行收敛和筛选
处理过程就是首先对lrt特征进行起伏系数的计算以及值的限制
对平坦度和差异度两个系数提取最大出现次数对应的feature和feature数量，通过这个数据来判断这两个feature是否使用
lrt为必须使用的特征
得到三个特征的权重
void WebRtcNs_FeatureParameterExtraction(NSinst_t* inst, 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 = inst->featureData[3];
//flat_feature = inst->featureData[0];
//diff_feature = inst->featureData[4];
//update histograms
if (flag == 0) {
// LRT
if ((inst->featureData[3] < HIST_PAR_EST * inst->featureExtractionParams.binSizeLrt)
&& (inst->featureData[3] >= 0.0)) {
i = (int)(inst->featureData[3] / inst->featureExtractionParams.binSizeLrt);
inst->histLrt[i]++;
}
// Spectral flatness
if ((inst->featureData[0] < HIST_PAR_EST
* inst->featureExtractionParams.binSizeSpecFlat)
&& (inst->featureData[0] >= 0.0)) {
i = (int)(inst->featureData[0] / inst->featureExtractionParams.binSizeSpecFlat);
inst->histSpecFlat[i]++;
}
// Spectral difference
if ((inst->featureData[4] < HIST_PAR_EST
* inst->featureExtractionParams.binSizeSpecDiff)
&& (inst->featureData[4] >= 0.0)) {
i = (int)(inst->featureData[4] / inst->featureExtractionParams.binSizeSpecDiff);
inst->histSpecDiff[i]++;
}
}
// extract parameters for speech/noise probability
if (flag == 1) {
//lrt feature: compute the average over inst->featureExtractionParams.rangeAvgHistLrt
avgHistLrt = 0.0;
avgHistLrtCompl = 0.0;
avgSquareHistLrt = 0.0;
numHistLrt = 0;
for (i = 0; i < HIST_PAR_EST; i++) {
//由feature = i*bin_size 得出binMid=feature+0.5*0.1 即binMid约等于feature原值 hisLrt存储对应feature数量
binMid = ((float)i + (float)0.5) * inst->featureExtractionParams.binSizeLrt;
//即feature(lrt)<1的情况 binmid < 10(rangeAvfHistLrt:10)
if (binMid <= inst->featureExtractionParams.rangeAvgHistLrt) { avgHistLrt += inst->histLrt[i] * binMid;numHistLrt += inst->histLrt[i];
}
avgSquareHistLrt += inst->histLrt[i] * binMid * binMid;
avgHistLrtCompl += inst->histLrt[i] * binMid;
}
//类似lrt<1的情况下所有特征平均值 num存储特征总数量
if (numHistLrt > 0) {
avgHistLrt = avgHistLrt / ((float)numHistLrt);
}
//所有lrt下特征的平均值 此处分母用500
avgHistLrtCompl = avgHistLrtCompl / ((float)inst->modelUpdatePars[1]);
//特征平方的平均值
avgSquareHistLrt = avgSquareHistLrt / ((float)inst->modelUpdatePars[1]);
//起伏特征：所有特征平方平均值-特征小于1的情况平均值*所有特征平均值
fluctLrt = avgSquareHistLrt - avgHistLrt * avgHistLrtCompl;
// get threshold for lrt feature: thresFluctLrt 0.05
//接下来就是系数的值域限制
if (fluctLrt < inst->featureExtractionParams.thresFluctLrt) {
//very low fluct, so likely noise 起伏系数？
inst->priorModelPars[0] = inst->featureExtractionParams.maxLrt;
} else {
inst->priorModelPars[0] = inst->featureExtractionParams.factor1ModelPars
* avgHistLrt;
// check if value is within min/max range
if (inst->priorModelPars[0] < inst->featureExtractionParams.minLrt) {
inst->priorModelPars[0] = inst->featureExtractionParams.minLrt;
}
if (inst->priorModelPars[0] > inst->featureExtractionParams.maxLrt) {
inst->priorModelPars[0] = inst->featureExtractionParams.maxLrt;
}
}
// done with lrt feature
//
// for spectral flatness and spectral difference: compute the main peaks of histogram
maxPeak1 = 0;
maxPeak2 = 0;
posPeak1SpecFlat = 0.0;
posPeak2SpecFlat = 0.0;
weightPeak1SpecFlat = 0;
weightPeak2SpecFlat = 0;
// peaks for flatness 平整度最大值和次大值
for (i = 0; i < HIST_PAR_EST; i++) {
binMid = ((float)i + (float)0.5) * inst->featureExtractionParams.binSizeSpecFlat;
if (inst->histSpecFlat[i] > maxPeak1) {
// Found new "first" peak
maxPeak2 = maxPeak1;
weightPeak2SpecFlat = weightPeak1SpecFlat;
posPeak2SpecFlat = posPeak1SpecFlat;
maxPeak1 = inst->histSpecFlat[i];
weightPeak1SpecFlat = inst->histSpecFlat[i];
posPeak1SpecFlat = binMid;
} else if (inst->histSpecFlat[i] > maxPeak2) {
// Found new "second" peak
maxPeak2 = inst->histSpecFlat[i];
weightPeak2SpecFlat = inst->histSpecFlat[i];
posPeak2SpecFlat = binMid;
}
}
//compute two peaks for spectral difference 差异度最大值和次大值
maxPeak1 = 0;
maxPeak2 = 0;
posPeak1SpecDiff = 0.0;
posPeak2SpecDiff = 0.0;
weightPeak1SpecDiff = 0;
weightPeak2SpecDiff = 0;
// peaks for spectral difference
for (i = 0; i < HIST_PAR_EST; i++) {
binMid = ((float)i + (float)0.5) * inst->featureExtractionParams.binSizeSpecDiff;
if (inst->histSpecDiff[i] > maxPeak1) {
// Found new "first" peak
maxPeak2 = maxPeak1;
weightPeak2SpecDiff = weightPeak1SpecDiff;
posPeak2SpecDiff = posPeak1SpecDiff;
maxPeak1 = inst->histSpecDiff[i];
weightPeak1SpecDiff = inst->histSpecDiff[i];
posPeak1SpecDiff = binMid;
} else if (inst->histSpecDiff[i] > maxPeak2) {
// Found new "second" peak
maxPeak2 = inst->histSpecDiff[i];
weightPeak2SpecDiff = inst->histSpecDiff[i];
posPeak2SpecDiff = binMid;
}
}
// for spectrum flatness feature
useFeatureSpecFlat = 1;
// merge the two peaks if they are close
//inst->featureExtractionParams.limitPeakSpacingSpecFlat = 0.1
// inst->featureExtractionParams.limitPeakWeightsSpecFlat = (float)0.5
//平坦度 最大次大满足需求后叠加取均值 以及判断是否使用
if ((fabs(posPeak2SpecFlat - posPeak1SpecFlat)
< inst->featureExtractionParams.limitPeakSpacingSpecFlat)
&& (weightPeak2SpecFlat
> inst->featureExtractionParams.limitPeakWeightsSpecFlat
* weightPeak1SpecFlat)) {
weightPeak1SpecFlat += weightPeak2SpecFlat;
posPeak1SpecFlat = (float)0.5 * (posPeak1SpecFlat + posPeak2SpecFlat);
}
//reject if weight of peaks is not large enough, or peak value too small
if (weightPeak1SpecFlat < inst->featureExtractionParams.thresWeightSpecFlat
|| posPeak1SpecFlat < inst->featureExtractionParams.thresPosSpecFlat) {
useFeatureSpecFlat = 0;
}
// if selected, get the threshold
if (useFeatureSpecFlat == 1) {
// compute the threshold
inst->priorModelPars[1] = inst->featureExtractionParams.factor2ModelPars
* posPeak1SpecFlat;
//check if value is within min/max range 限制大小在：0.1-0.95
if (inst->priorModelPars[1] < inst->featureExtractionParams.minSpecFlat) {
inst->priorModelPars[1] = inst->featureExtractionParams.minSpecFlat;
}
if (inst->priorModelPars[1] > inst->featureExtractionParams.maxSpecFlat) {
inst->priorModelPars[1] = inst->featureExtractionParams.maxSpecFlat;
}
}
// done with flatness feature
// for template feature
useFeatureSpecDiff = 1;
// merge the two peaks if they are close
//差异度 取极大次大值合并 判断是否使用
if ((fabs(posPeak2SpecDiff - posPeak1SpecDiff)
< inst->featureExtractionParams.limitPeakSpacingSpecDiff)
&& (weightPeak2SpecDiff
> inst->featureExtractionParams.limitPeakWeightsSpecDiff
* weightPeak1SpecDiff)) {
weightPeak1SpecDiff += weightPeak2SpecDiff;
posPeak1SpecDiff = (float)0.5 * (posPeak1SpecDiff + posPeak2SpecDiff);
}
// get the threshold value
inst->priorModelPars[3] = inst->featureExtractionParams.factor1ModelPars
* posPeak1SpecDiff;
//reject if weight of peaks is not large enough
if (weightPeak1SpecDiff < inst->featureExtractionParams.thresWeightSpecDiff) {
useFeatureSpecDiff = 0;
}
//check if value is within min/max range
if (inst->priorModelPars[3] < inst->featureExtractionParams.minSpecDiff) {
inst->priorModelPars[3] = inst->featureExtractionParams.minSpecDiff;
}
if (inst->priorModelPars[3] > inst->featureExtractionParams.maxSpecDiff) {
inst->priorModelPars[3] = inst->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 < inst->featureExtractionParams.thresFluctLrt) {
useFeatureSpecDiff = 0;
}
// select the weights between the features
// inst->priorModelPars[4] is weight for lrt: always selected
// inst->priorModelPars[5] is weight for spectral flatness
// inst->priorModelPars[6] is weight for spectral difference
featureSum = (float)(1 + useFeatureSpecFlat + useFeatureSpecDiff);
inst->priorModelPars[4] = (float)1.0 / featureSum;
inst->priorModelPars[5] = ((float)useFeatureSpecFlat) / featureSum;
inst->priorModelPars[6] = ((float)useFeatureSpecDiff) / featureSum;
// set hists to zero for next update
if (inst->modelUpdatePars[0] >= 1) {
for (i = 0; i < HIST_PAR_EST; i++) {
inst->histLrt[i] = 0;
inst->histSpecFlat[i] = 0;
inst->histSpecDiff[i] = 0;
}
}
} // end of flag == 1
}得到以下参数
8）语音/噪音概率估计
主要是计算在带噪语音和特征条件下是语音的概率为
P(H1 | YF) (H1：语音状态 Y：原始信号 F：特征集)
利用的原理如下：
1)贝叶斯公式：假设A和B为独立分布的时间，存在P(AB)=p(A|B)*P(B) （同时发生AB的概率等于B的发生概率*B发生条件下A发生的概率）在这里假设噪音跟语音是独立分布的两种事件
2)语音的线性关系：设定H0为噪音 H1为语音 Y为原始信号 存在 P(H0|YF) + P(H1|YF) = 1 , P(H0|Y) + P(H1|Y) =1
推导过程如下
P(H0|F) = 1-P(H1|F)
根据推导出的公式得出求出以下两个值 就可以计算出P(H1|YF)
简化为下面的公式
求解似然比系数采用联合高斯概率密度分布的操作
高斯概率密度介绍如下：
这里假设噪音和语音符合均值为0的概率密度分布
对应似然比计算如下
代码如下：
void WebRtcNs_SpeechNoiseProb(NSinst_t* inst, float* probSpeechFinal, float* snrLocPrior,
float* snrLocPost) {
int i, 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;
widthPrior1 = (float)2.0 * WIDTH_PR_MAP; //width for pause region:
// lower range, so increase width in tanh map
widthPrior2 = (float)2.0 * WIDTH_PR_MAP; //for spectral-difference measure
//threshold parameters for features
threshPrior0 = inst->priorModelPars[0];
threshPrior1 = inst->priorModelPars[1];
threshPrior2 = inst->priorModelPars[3];
//sign for flatness feature
sgnMap = (int)(inst->priorModelPars[2]);
//weight parameters for features
weightIndPrior0 = inst->priorModelPars[4];
weightIndPrior1 = inst->priorModelPars[5];
weightIndPrior2 = inst->priorModelPars[6];
// compute feature based on average LR factor
// this is the average over all frequencies of the smooth log lrt
//通过似然比与高斯密度概率函数求解lrt均值
logLrtTimeAvgKsum = 0.0;
for (i = 0; i < inst->magnLen; i++) {
tmpFloat1 = (float)1.0 + (float)2.0 * snrLocPrior[i];
tmpFloat2 = (float)2.0 * snrLocPrior[i] / (tmpFloat1 + (float)0.0001);
besselTmp = (snrLocPost[i] + (float)1.0) * tmpFloat2;
inst->logLrtTimeAvg[i] += LRT_TAVG * (besselTmp - (float)log(tmpFloat1)- inst->logLrtTimeAvg[i]);
logLrtTimeAvgKsum += inst->logLrtTimeAvg[i];
}
//以上计算得到极大似然比
//取极大似然比的平均值
logLrtTimeAvgKsum = (float)logLrtTimeAvgKsum / (inst->magnLen);
inst->featureData[3] = logLrtTimeAvgKsum;接下来求特征条件下的语音概率
用到了上面特征筛选中计算出的权重与特征值
代码如下
这里采用了取反正切sigmod函数用来限制feature的值域(反正切的值域为[-1,1])
方法如下：
代码如下
//以下操作对三个特征值进行sigmod函数变换 将值控制在同一范围内
// sigmod函数采用tanh（单调增 值域【-1，1】 代码中控制在[-0.5,0.5]中间）
// 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 = (float)0.5 * ((float)tanh(widthPrior *
(logLrtTimeAvgKsum - threshPrior0)) + (float)1.0);
//spectral flatness feature
tmpFloat1 = inst->featureData[0];
widthPrior = widthPrior0;
//use larger width in tanh map for pause regions 使用flatnesss
if (sgnMap == 1 && (tmpFloat1 > threshPrior1)) {
widthPrior = widthPrior1;
}
//feature0>feature1 不使用flatness
if (sgnMap == -1 && (tmpFloat1 < threshPrior1)) {
widthPrior = widthPrior1;
}
// compute indicator function: sigmoid map
indicator1 = (float)0.5 * ((float)tanh((float)sgnMap *
widthPrior * (threshPrior1 - tmpFloat1)) + (float)1.0);
//for template spectrum-difference
tmpFloat1 = inst->featureData[4];
widthPrior = widthPrior0;
//use larger width in tanh map for pause regions
if (tmpFloat1 < threshPrior2) {
widthPrior = widthPrior2;
}
// compute indicator function: sigmoid map
indicator2 = (float)0.5 * ((float)tanh(widthPrior * (tmpFloat1 - threshPrior2))
+ (float)1.0);
//combine the indicator function with the feature weights
indPrior = weightIndPrior0 * indicator0 + weightIndPrior1 * indicator1 + weightIndPrior2
* indicator2;
// done with computing indicator function
//compute the prior probability
inst->priorSpeechProb += PRIOR_UPDATE * (indPrior - inst->priorSpeechProb);接下来取计算好的两个概率计算最终的语音概率 得到特征条件以及原始信号条件下的语音概率
//final speech probability: combine prior model with LR factor:
gainPrior = ((float)1.0 - inst->priorSpeechProb) / (inst->priorSpeechProb + (float)0.0001);
for (i = 0; i < inst->magnLen; i++) {
invLrt = (float)exp(-inst->logLrtTimeAvg[i]);
invLrt = (float)gainPrior * invLrt;
probSpeechFinal[i] = (float)1.0 / ((float)1.0 + invLrt);
}9）根据语音概率的结果更新之前估计的噪音(反馈)
由P(H0|YF) = 1-P(H1|YF)得出
代码如下：计算类似上面噪声初始化估计的过程
//更新完语音概率后 更新noise
gammaNoiseTmp = NOISE_UPDATE;
for (i = 0; i < inst->magnLen; i++) {
probSpeech = probSpeechFinal[i];
probNonSpeech = (float)1.0 - probSpeech;
// temporary noise update:
// use it for speech frames if update value is less than previous
noiseUpdateTmp = gammaNoiseTmp * inst->noisePrev[i] + ((float)1.0 - gammaNoiseTmp)
* (probNonSpeech * magn[i] + probSpeech * inst->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
//语音概率较小 更新magnAvgPause
if (probSpeech < PROB_RANGE) {
inst->magnAvgPause[i] += GAMMA_PAUSE * (magn[i] - inst->magnAvgPause[i]);
}
// noise update
if (gammaNoiseTmp == gammaNoiseOld) {
noise[i] = noiseUpdateTmp;
} else {
noise[i] = gammaNoiseTmp * inst->noisePrev[i] + ((float)1.0 - gammaNoiseTmp)
* (probNonSpeech * magn[i] + probSpeech * inst->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
// done with step 2: noise update
//
// STEP 3: compute dd update of prior snr and post snr based on new noise estimate
//filter:比重
for (i = 0; i < inst->magnLen; i++) {
// post and prior snr
currentEstimateStsa = (float)0.0;
if (magn[i] > noise[i]) {
currentEstimateStsa = magn[i] / (noise[i] + (float)0.0001) - (float)1.0;
}
// DD estimate is sume of two terms: current estimate and previous estimate
// directed decision update of snrPrior
snrPrior = DD_PR_SNR * previousEstimateStsa[i] + ((float)1.0 - DD_PR_SNR)
* currentEstimateStsa;
// gain filter
tmpFloat1 = inst->overdrive + snrPrior;
tmpFloat2 = (float)snrPrior / tmpFloat1;
theFilter[i] = (float)tmpFloat2;
} // end of loop over freqs
// done with step310）维纳滤波
维纳滤波介绍如下：
由于上面已经计算出Pri即先验概率的值 所以可以直接用来求最终的结果
期望信号 = 原始信号*维纳系数
维纳系数计算的代码如下： override设置为1
for (i = 0; i < inst->magnLen; i++) {
// post and prior snr
currentEstimateStsa = (float)0.0;
if (magn[i] > noise[i]) {
currentEstimateStsa = magn[i] / (noise[i] + (float)0.0001) - (float)1.0;
}
// DD estimate is sume of two terms: current estimate and previous estimate
// directed decision update of snrPrior
snrPrior = DD_PR_SNR * previousEstimateStsa[i] + ((float)1.0 - DD_PR_SNR)
* currentEstimateStsa;
// gain filter
tmpFloat1 = inst->overdrive + snrPrior;
tmpFloat2 = (float)snrPrior / tmpFloat1;
theFilter[i] = (float)tmpFloat2;
} // end of loop over freqs
// done with step3通过维纳滤波系数以及原始信号计算语音信号 代码如下：
for (i = 0; i < inst->magnLen; i++) {
// flooring bottom
if (theFilter[i] < inst->denoiseBound) {
theFilter[i] = inst->denoiseBound;
}
// flooring top
if (theFilter[i] > (float)1.0) {
theFilter[i] = 1.0;
}
if (inst->blockInd < END_STARTUP_SHORT) {
// flooring bottom
if (theFilterTmp[i] < inst->denoiseBound) {
theFilterTmp[i] = inst->denoiseBound;
}
// flooring top
if (theFilterTmp[i] > (float)1.0) {
theFilterTmp[i] = 1.0;
}
// Weight the two suppression filters
theFilter[i] *= (inst->blockInd);
theFilterTmp[i] *= (END_STARTUP_SHORT - inst->blockInd);
theFilter[i] += theFilterTmp[i];
theFilter[i] /= (END_STARTUP_SHORT);
}
// smoothing
#ifdef PROCESS_FLOW_0
inst->smooth[i] *= SMOOTH; // value set to 0.7 in define.h file
inst->smooth[i] += ((float)1.0 - SMOOTH) * theFilter[i];
#else
inst->smooth[i] = theFilter[i];
#endif
real[i] *= inst->smooth[i];
imag[i] *= inst->smooth[i];
}将计算好的数据存储到prev中用于下次的计算迭代平滑
将数据进行一些收敛判断，进行反傅里叶变换，就得出最终的语音信号
for (i = 0; i < inst->magnLen; i++) {
inst->noisePrev[i] = noise[i];
inst->magnPrev[i] = magn[i];
}
// back to time domain
winData[0] = real[0];
winData[1] = real[inst->magnLen - 1];
for (i = 1; i < inst->magnLen - 1; i++) {
winData[2 * i] = real[i];
winData[2 * i + 1] = imag[i];
}
//反傅里叶变换 幅度谱生成原始数据
WebRtc_rdft(inst->anaLen, -1, winData, inst->ip, inst->wfft);
参考:书籍《实时语音处理实践指南》
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
