diff options
Diffstat (limited to 'webrtc/modules/audio_processing/ns/ns_core.c')
| -rw-r--r-- | webrtc/modules/audio_processing/ns/ns_core.c | 1416 |
1 files changed, 1416 insertions, 0 deletions
diff --git a/webrtc/modules/audio_processing/ns/ns_core.c b/webrtc/modules/audio_processing/ns/ns_core.c new file mode 100644 index 0000000000..1d6091400e --- /dev/null +++ b/webrtc/modules/audio_processing/ns/ns_core.c @@ -0,0 +1,1416 @@ +/* + * Copyright (c) 2012 The WebRTC project authors. All Rights Reserved. + * + * Use of this source code is governed by a BSD-style license + * that can be found in the LICENSE file in the root of the source + * tree. An additional intellectual property rights grant can be found + * in the file PATENTS. All contributing project authors may + * be found in the AUTHORS file in the root of the source tree. + */ + +#include <assert.h> +#include <math.h> +#include <string.h> +#include <stdlib.h> + +#include "webrtc/common_audio/fft4g.h" +#include "webrtc/common_audio/signal_processing/include/signal_processing_library.h" +#include "webrtc/modules/audio_processing/ns/include/noise_suppression.h" +#include "webrtc/modules/audio_processing/ns/ns_core.h" +#include "webrtc/modules/audio_processing/ns/windows_private.h" + +// Set Feature Extraction Parameters. +static void set_feature_extraction_parameters(NoiseSuppressionC* self) { + // Bin size of histogram. + self->featureExtractionParams.binSizeLrt = 0.1f; + self->featureExtractionParams.binSizeSpecFlat = 0.05f; + self->featureExtractionParams.binSizeSpecDiff = 0.1f; + + // Range of histogram over which LRT threshold is computed. + self->featureExtractionParams.rangeAvgHistLrt = 1.f; + + // Scale parameters: multiply dominant peaks of the histograms by scale factor + // to obtain thresholds for prior model. + // For LRT and spectral difference. + self->featureExtractionParams.factor1ModelPars = 1.2f; + // For spectral_flatness: used when noise is flatter than speech. + self->featureExtractionParams.factor2ModelPars = 0.9f; + + // Peak limit for spectral flatness (varies between 0 and 1). + self->featureExtractionParams.thresPosSpecFlat = 0.6f; + + // Limit on spacing of two highest peaks in histogram: spacing determined by + // bin size. + self->featureExtractionParams.limitPeakSpacingSpecFlat = + 2 * self->featureExtractionParams.binSizeSpecFlat; + self->featureExtractionParams.limitPeakSpacingSpecDiff = + 2 * self->featureExtractionParams.binSizeSpecDiff; + + // Limit on relevance of second peak. + self->featureExtractionParams.limitPeakWeightsSpecFlat = 0.5f; + self->featureExtractionParams.limitPeakWeightsSpecDiff = 0.5f; + + // Fluctuation limit of LRT feature. + self->featureExtractionParams.thresFluctLrt = 0.05f; + + // Limit on the max and min values for the feature thresholds. + self->featureExtractionParams.maxLrt = 1.f; + self->featureExtractionParams.minLrt = 0.2f; + + self->featureExtractionParams.maxSpecFlat = 0.95f; + self->featureExtractionParams.minSpecFlat = 0.1f; + + self->featureExtractionParams.maxSpecDiff = 1.f; + self->featureExtractionParams.minSpecDiff = 0.16f; + + // Criteria of weight of histogram peak to accept/reject feature. + self->featureExtractionParams.thresWeightSpecFlat = + (int)(0.3 * (self->modelUpdatePars[1])); // For spectral flatness. + self->featureExtractionParams.thresWeightSpecDiff = + (int)(0.3 * (self->modelUpdatePars[1])); // For spectral difference. +} + +// Initialize state. +int WebRtcNs_InitCore(NoiseSuppressionC* self, uint32_t fs) { + int i; + // Check for valid pointer. + if (self == NULL) { + return -1; + } + + // Initialization of struct. + if (fs == 8000 || fs == 16000 || fs == 32000 || fs == 48000) { + self->fs = fs; + } else { + return -1; + } + self->windShift = 0; + // We only support 10ms frames. + if (fs == 8000) { + self->blockLen = 80; + self->anaLen = 128; + self->window = kBlocks80w128; + } else { + self->blockLen = 160; + self->anaLen = 256; + self->window = kBlocks160w256; + } + self->magnLen = self->anaLen / 2 + 1; // Number of frequency bins. + + // Initialize FFT work arrays. + self->ip[0] = 0; // Setting this triggers initialization. + memset(self->dataBuf, 0, sizeof(float) * ANAL_BLOCKL_MAX); + WebRtc_rdft(self->anaLen, 1, self->dataBuf, self->ip, self->wfft); + + memset(self->analyzeBuf, 0, sizeof(float) * ANAL_BLOCKL_MAX); + memset(self->dataBuf, 0, sizeof(float) * ANAL_BLOCKL_MAX); + memset(self->syntBuf, 0, sizeof(float) * ANAL_BLOCKL_MAX); + + // For HB processing. + memset(self->dataBufHB, + 0, + sizeof(float) * NUM_HIGH_BANDS_MAX * ANAL_BLOCKL_MAX); + + // For quantile noise estimation. + memset(self->quantile, 0, sizeof(float) * HALF_ANAL_BLOCKL); + for (i = 0; i < SIMULT * HALF_ANAL_BLOCKL; i++) { + self->lquantile[i] = 8.f; + self->density[i] = 0.3f; + } + + for (i = 0; i < SIMULT; i++) { + self->counter[i] = + (int)floor((float)(END_STARTUP_LONG * (i + 1)) / (float)SIMULT); + } + + self->updates = 0; + + // Wiener filter initialization. + for (i = 0; i < HALF_ANAL_BLOCKL; i++) { + self->smooth[i] = 1.f; + } + + // Set the aggressiveness: default. + self->aggrMode = 0; + + // Initialize variables for new method. + self->priorSpeechProb = 0.5f; // Prior prob for speech/noise. + // Previous analyze mag spectrum. + memset(self->magnPrevAnalyze, 0, sizeof(float) * HALF_ANAL_BLOCKL); + // Previous process mag spectrum. + memset(self->magnPrevProcess, 0, sizeof(float) * HALF_ANAL_BLOCKL); + // Current noise-spectrum. + memset(self->noise, 0, sizeof(float) * HALF_ANAL_BLOCKL); + // Previous noise-spectrum. + memset(self->noisePrev, 0, sizeof(float) * HALF_ANAL_BLOCKL); + // Conservative noise spectrum estimate. + memset(self->magnAvgPause, 0, sizeof(float) * HALF_ANAL_BLOCKL); + // For estimation of HB in second pass. + memset(self->speechProb, 0, sizeof(float) * HALF_ANAL_BLOCKL); + // Initial average magnitude spectrum. + memset(self->initMagnEst, 0, sizeof(float) * HALF_ANAL_BLOCKL); + for (i = 0; i < HALF_ANAL_BLOCKL; i++) { + // Smooth LR (same as threshold). + self->logLrtTimeAvg[i] = LRT_FEATURE_THR; + } + + // Feature quantities. + // Spectral flatness (start on threshold). + self->featureData[0] = SF_FEATURE_THR; + self->featureData[1] = 0.f; // Spectral entropy: not used in this version. + self->featureData[2] = 0.f; // Spectral variance: not used in this version. + // Average LRT factor (start on threshold). + self->featureData[3] = LRT_FEATURE_THR; + // Spectral template diff (start on threshold). + self->featureData[4] = SF_FEATURE_THR; + self->featureData[5] = 0.f; // Normalization for spectral difference. + // Window time-average of input magnitude spectrum. + self->featureData[6] = 0.f; + + // Histogram quantities: used to estimate/update thresholds for features. + memset(self->histLrt, 0, sizeof(int) * HIST_PAR_EST); + memset(self->histSpecFlat, 0, sizeof(int) * HIST_PAR_EST); + memset(self->histSpecDiff, 0, sizeof(int) * HIST_PAR_EST); + + + self->blockInd = -1; // Frame counter. + // Default threshold for LRT feature. + self->priorModelPars[0] = LRT_FEATURE_THR; + // Threshold for spectral flatness: determined on-line. + self->priorModelPars[1] = 0.5f; + // sgn_map par for spectral measure: 1 for flatness measure. + self->priorModelPars[2] = 1.f; + // Threshold for template-difference feature: determined on-line. + self->priorModelPars[3] = 0.5f; + // Default weighting parameter for LRT feature. + self->priorModelPars[4] = 1.f; + // Default weighting parameter for spectral flatness feature. + self->priorModelPars[5] = 0.f; + // Default weighting parameter for spectral difference feature. + self->priorModelPars[6] = 0.f; + + // Update flag for parameters: + // 0 no update, 1 = update once, 2 = update every window. + self->modelUpdatePars[0] = 2; + self->modelUpdatePars[1] = 500; // Window for update. + // Counter for update of conservative noise spectrum. + self->modelUpdatePars[2] = 0; + // Counter if the feature thresholds are updated during the sequence. + self->modelUpdatePars[3] = self->modelUpdatePars[1]; + + self->signalEnergy = 0.0; + self->sumMagn = 0.0; + self->whiteNoiseLevel = 0.0; + self->pinkNoiseNumerator = 0.0; + self->pinkNoiseExp = 0.0; + + set_feature_extraction_parameters(self); + + // Default mode. + WebRtcNs_set_policy_core(self, 0); + + self->initFlag = 1; + return 0; +} + +// Estimate noise. +static void NoiseEstimation(NoiseSuppressionC* self, + float* magn, + float* noise) { + size_t i, s, offset; + float lmagn[HALF_ANAL_BLOCKL], delta; + + if (self->updates < END_STARTUP_LONG) { + self->updates++; + } + + for (i = 0; i < self->magnLen; i++) { + lmagn[i] = (float)log(magn[i]); + } + + // Loop over simultaneous estimates. + for (s = 0; s < SIMULT; s++) { + offset = s * self->magnLen; + + // newquantest(...) + for (i = 0; i < self->magnLen; i++) { + // Compute delta. + if (self->density[offset + i] > 1.0) { + delta = FACTOR * 1.f / self->density[offset + i]; + } else { + delta = FACTOR; + } + + // Update log quantile estimate. + if (lmagn[i] > self->lquantile[offset + i]) { + self->lquantile[offset + i] += + QUANTILE * delta / (float)(self->counter[s] + 1); + } else { + self->lquantile[offset + i] -= + (1.f - QUANTILE) * delta / (float)(self->counter[s] + 1); + } + + // Update density estimate. + if (fabs(lmagn[i] - self->lquantile[offset + i]) < WIDTH) { + self->density[offset + i] = + ((float)self->counter[s] * self->density[offset + i] + + 1.f / (2.f * WIDTH)) / + (float)(self->counter[s] + 1); + } + } // End loop over magnitude spectrum. + + 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] = (float)exp(self->lquantile[offset + i]); + } + } + } + + self->counter[s]++; + } // End loop over simultaneous estimates. + + // Sequentially update the noise during startup. + if (self->updates < END_STARTUP_LONG) { + // Use the last "s" to get noise during startup that differ from zero. + for (i = 0; i < self->magnLen; i++) { + self->quantile[i] = (float)exp(self->lquantile[offset + i]); + } + } + + for (i = 0; i < self->magnLen; i++) { + noise[i] = self->quantile[i]; + } +} + +// Extract thresholds for feature parameters. +// Histograms are computed over some window size (given by +// self->modelUpdatePars[1]). +// Thresholds and weights are extracted every window. +// |flag| = 0 updates histogram only, |flag| = 1 computes the threshold/weights. +// Threshold and weights are returned in: self->priorModelPars. +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.0; + avgHistLrtCompl = 0.0; + avgSquareHistLrt = 0.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.0; + posPeak2SpecFlat = 0.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.0; + posPeak2SpecDiff = 0.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 ((fabs(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 ((fabs(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) / featureSum; + self->priorModelPars[6] = ((float)useFeatureSpecDiff) / featureSum; + + // 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. +} + +// Compute spectral flatness on input spectrum. +// |magnIn| is the magnitude spectrum. +// Spectral flatness is returned in self->featureData[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. +} + +// Compute prior and post SNR based on quantile noise estimation. +// Compute DD estimate of prior SNR. +// Inputs: +// * |magn| is the signal magnitude spectrum estimate. +// * |noise| is the magnitude noise spectrum estimate. +// Outputs: +// * |snrLocPrior| is the computed prior SNR. +// * |snrLocPost| is the computed post SNR. +static void ComputeSnr(const NoiseSuppressionC* self, + const float* magn, + const float* noise, + float* snrLocPrior, + 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->noisePrev[i] + 0.0001f) * self->smooth[i]; + // Post SNR. + snrLocPost[i] = 0.f; + if (magn[i] > noise[i]) { + snrLocPost[i] = magn[i] / (noise[i] + 0.0001f) - 1.f; + } + // DD estimate is sum of two terms: current estimate and previous estimate. + // Directed decision update of snrPrior. + snrLocPrior[i] = + DD_PR_SNR * previousEstimateStsa + (1.f - DD_PR_SNR) * snrLocPost[i]; + } // End of loop over frequencies. +} + +// Compute the difference measure between input spectrum and a template/learned +// noise spectrum. +// |magnIn| is the input spectrum. +// The reference/template spectrum is self->magnAvgPause[i]. +// Returns (normalized) spectral difference in 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.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->magnLen; + avgMagn /= self->magnLen; + + covMagnPause = 0.0; + varPause = 0.0; + varMagn = 0.0; + // Compute variance and covariance quantities. + for (i = 0; i < self->magnLen; i++) { + covMagnPause += (magnIn[i] - avgMagn) * (self->magnAvgPause[i] - avgPause); + varPause += + (self->magnAvgPause[i] - avgPause) * (self->magnAvgPause[i] - avgPause); + varMagn += (magnIn[i] - avgMagn) * (magnIn[i] - avgMagn); + } + covMagnPause /= self->magnLen; + varPause /= self->magnLen; + varMagn /= self->magnLen; + // Update of average magnitude spectrum. + self->featureData[6] += self->signalEnergy; + + avgDiffNormMagn = + varMagn - (covMagnPause * covMagnPause) / (varPause + 0.0001f); + // Normalize and compute time-avg update of difference feature. + avgDiffNormMagn = (float)(avgDiffNormMagn / (self->featureData[5] + 0.0001f)); + self->featureData[4] += + SPECT_DIFF_TAVG * (avgDiffNormMagn - self->featureData[4]); +} + +// Compute speech/noise probability. +// Speech/noise probability is returned in |probSpeechFinal|. +// |magn| is the input magnitude spectrum. +// |noise| is the noise spectrum. +// |snrLocPrior| is the prior SNR for each frequency. +// |snrLocPost| is the post SNR for each frequency. +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); + } +} + +// Update the noise features. +// Inputs: +// * |magn| is the signal magnitude spectrum estimate. +// * |updateParsFlag| is an update flag for parameters. +static void FeatureUpdate(NoiseSuppressionC* self, + const float* magn, + int updateParsFlag) { + // Compute spectral flatness on input spectrum. + ComputeSpectralFlatness(self, magn); + // Compute difference of input spectrum with learned/estimated noise spectrum. + ComputeSpectralDifference(self, magn); + // Compute histograms for parameter decisions (thresholds and weights for + // features). + // Parameters are extracted once every window time. + // (=self->modelUpdatePars[1]) + if (updateParsFlag >= 1) { + // Counter update. + self->modelUpdatePars[3]--; + // Update histogram. + if (self->modelUpdatePars[3] > 0) { + FeatureParameterExtraction(self, 0); + } + // Compute model parameters. + if (self->modelUpdatePars[3] == 0) { + FeatureParameterExtraction(self, 1); + self->modelUpdatePars[3] = self->modelUpdatePars[1]; + // If wish to update only once, set flag to zero. + if (updateParsFlag == 1) { + self->modelUpdatePars[0] = 0; + } else { + // Update every window: + // Get normalization for spectral difference for next window estimate. + self->featureData[6] = + self->featureData[6] / ((float)self->modelUpdatePars[1]); + self->featureData[5] = + 0.5f * (self->featureData[6] + self->featureData[5]); + self->featureData[6] = 0.f; + } + } + } +} + +// Update the noise estimate. +// Inputs: +// * |magn| is the signal magnitude spectrum estimate. +// * |snrLocPrior| is the prior SNR. +// * |snrLocPost| is the post SNR. +// Output: +// * |noise| is the updated noise magnitude spectrum estimate. +static void UpdateNoiseEstimate(NoiseSuppressionC* self, + const float* magn, + const float* snrLocPrior, + const float* snrLocPost, + 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. +} + +// Updates |buffer| with a new |frame|. +// Inputs: +// * |frame| is a new speech frame or NULL for setting to zero. +// * |frame_length| is the length of the new frame. +// * |buffer_length| is the length of the buffer. +// Output: +// * |buffer| is the updated buffer. +static void UpdateBuffer(const float* frame, + size_t frame_length, + size_t buffer_length, + float* buffer) { + assert(buffer_length < 2 * frame_length); + + memcpy(buffer, + buffer + frame_length, + sizeof(*buffer) * (buffer_length - frame_length)); + if (frame) { + memcpy(buffer + buffer_length - frame_length, + frame, + sizeof(*buffer) * frame_length); + } else { + memset(buffer + buffer_length - frame_length, + 0, + sizeof(*buffer) * frame_length); + } +} + +// Transforms the signal from time to frequency domain. +// Inputs: +// * |time_data| is the signal in the time domain. +// * |time_data_length| is the length of the analysis buffer. +// * |magnitude_length| is the length of the spectrum magnitude, which equals +// the length of both |real| and |imag| (time_data_length / 2 + 1). +// Outputs: +// * |time_data| is the signal in the frequency domain. +// * |real| is the real part of the frequency domain. +// * |imag| is the imaginary part of the frequency domain. +// * |magn| is the calculated signal magnitude in the frequency domain. +static void FFT(NoiseSuppressionC* self, + float* time_data, + size_t time_data_length, + size_t magnitude_length, + float* real, + float* imag, + float* magn) { + size_t i; + + assert(magnitude_length == time_data_length / 2 + 1); + + WebRtc_rdft(time_data_length, 1, time_data, self->ip, self->wfft); + + imag[0] = 0; + real[0] = time_data[0]; + magn[0] = fabsf(real[0]) + 1.f; + imag[magnitude_length - 1] = 0; + real[magnitude_length - 1] = time_data[1]; + magn[magnitude_length - 1] = fabsf(real[magnitude_length - 1]) + 1.f; + for (i = 1; i < magnitude_length - 1; ++i) { + real[i] = time_data[2 * i]; + imag[i] = time_data[2 * i + 1]; + // Magnitude spectrum. + magn[i] = sqrtf(real[i] * real[i] + imag[i] * imag[i]) + 1.f; + } +} + +// Transforms the signal from frequency to time domain. +// Inputs: +// * |real| is the real part of the frequency domain. +// * |imag| is the imaginary part of the frequency domain. +// * |magnitude_length| is the length of the spectrum magnitude, which equals +// the length of both |real| and |imag|. +// * |time_data_length| is the length of the analysis buffer +// (2 * (magnitude_length - 1)). +// Output: +// * |time_data| is the signal in the time domain. +static void IFFT(NoiseSuppressionC* self, + const float* real, + const float* imag, + size_t magnitude_length, + size_t time_data_length, + float* time_data) { + size_t i; + + assert(time_data_length == 2 * (magnitude_length - 1)); + + time_data[0] = real[0]; + time_data[1] = real[magnitude_length - 1]; + for (i = 1; i < magnitude_length - 1; ++i) { + time_data[2 * i] = real[i]; + time_data[2 * i + 1] = imag[i]; + } + WebRtc_rdft(time_data_length, -1, time_data, self->ip, self->wfft); + + for (i = 0; i < time_data_length; ++i) { + time_data[i] *= 2.f / time_data_length; // FFT scaling. + } +} + +// Calculates the energy of a buffer. +// Inputs: +// * |buffer| is the buffer over which the energy is calculated. +// * |length| is the length of the buffer. +// Returns the calculated energy. +static float Energy(const float* buffer, size_t length) { + size_t i; + float energy = 0.f; + + for (i = 0; i < length; ++i) { + energy += buffer[i] * buffer[i]; + } + + return energy; +} + +// Windows a buffer. +// Inputs: +// * |window| is the window by which to multiply. +// * |data| is the data without windowing. +// * |length| is the length of the window and data. +// Output: +// * |data_windowed| is the windowed data. +static void Windowing(const float* window, + const float* data, + size_t length, + float* data_windowed) { + size_t i; + + for (i = 0; i < length; ++i) { + data_windowed[i] = window[i] * data[i]; + } +} + +// 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. +} + +// Changes the aggressiveness of the noise suppression method. +// |mode| = 0 is mild (6dB), |mode| = 1 is medium (10dB) and |mode| = 2 is +// aggressive (15dB). +// Returns 0 on success and -1 otherwise. +int WebRtcNs_set_policy_core(NoiseSuppressionC* self, int mode) { + // Allow for modes: 0, 1, 2, 3. + if (mode < 0 || mode > 3) { + return (-1); + } + + self->aggrMode = mode; + if (mode == 0) { + self->overdrive = 1.f; + self->denoiseBound = 0.5f; + self->gainmap = 0; + } else if (mode == 1) { + // self->overdrive = 1.25f; + self->overdrive = 1.f; + self->denoiseBound = 0.25f; + self->gainmap = 1; + } else if (mode == 2) { + // self->overdrive = 1.25f; + self->overdrive = 1.1f; + self->denoiseBound = 0.125f; + self->gainmap = 1; + } else if (mode == 3) { + // self->overdrive = 1.3f; + self->overdrive = 1.25f; + self->denoiseBound = 0.09f; + self->gainmap = 1; + } + return 0; +} + +void WebRtcNs_AnalyzeCore(NoiseSuppressionC* self, const float* speechFrame) { + size_t i; + const size_t kStartBand = 5; // Skip first frequency bins during estimation. + int updateParsFlag; + float energy; + float signalEnergy = 0.f; + float sumMagn = 0.f; + float tmpFloat1, tmpFloat2, tmpFloat3; + float winData[ANAL_BLOCKL_MAX]; + float magn[HALF_ANAL_BLOCKL], noise[HALF_ANAL_BLOCKL]; + float snrLocPost[HALF_ANAL_BLOCKL], snrLocPrior[HALF_ANAL_BLOCKL]; + float real[ANAL_BLOCKL_MAX], imag[HALF_ANAL_BLOCKL]; + // Variables during startup. + float sum_log_i = 0.0; + float sum_log_i_square = 0.0; + float sum_log_magn = 0.0; + float sum_log_i_log_magn = 0.0; + float parametric_exp = 0.0; + float parametric_num = 0.0; + + // Check that initiation has been done. + assert(self->initFlag == 1); + updateParsFlag = self->modelUpdatePars[0]; + + // Update analysis buffer for L band. + UpdateBuffer(speechFrame, self->blockLen, self->anaLen, self->analyzeBuf); + + Windowing(self->window, self->analyzeBuf, self->anaLen, winData); + energy = Energy(winData, self->anaLen); + if (energy == 0.0) { + // We want to avoid updating statistics in this case: + // Updating feature statistics when we have zeros only will cause + // thresholds to move towards zero signal situations. This in turn has the + // effect that once the signal is "turned on" (non-zero values) everything + // will be treated as speech and there is no noise suppression effect. + // Depending on the duration of the inactive signal it takes a + // considerable amount of time for the system to learn what is noise and + // what is speech. + return; + } + + self->blockInd++; // Update the block index only when we process a block. + + FFT(self, winData, self->anaLen, self->magnLen, real, imag, magn); + + for (i = 0; i < self->magnLen; i++) { + signalEnergy += real[i] * real[i] + imag[i] * imag[i]; + sumMagn += magn[i]; + if (self->blockInd < END_STARTUP_SHORT) { + if (i >= kStartBand) { + tmpFloat2 = logf((float)i); + sum_log_i += tmpFloat2; + sum_log_i_square += tmpFloat2 * tmpFloat2; + tmpFloat1 = logf(magn[i]); + sum_log_magn += tmpFloat1; + sum_log_i_log_magn += tmpFloat2 * tmpFloat1; + } + } + } + signalEnergy /= self->magnLen; + self->signalEnergy = signalEnergy; + self->sumMagn = sumMagn; + + // Quantile noise estimate. + NoiseEstimation(self, magn, noise); + // Compute simplified noise model during startup. + if (self->blockInd < END_STARTUP_SHORT) { + // Estimate White noise. + self->whiteNoiseLevel += sumMagn / self->magnLen * self->overdrive; + // Estimate Pink noise parameters. + tmpFloat1 = sum_log_i_square * (self->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.f) { + tmpFloat3 = 0.f; + } + self->pinkNoiseNumerator += tmpFloat3; + tmpFloat2 = (sum_log_i * sum_log_magn); + tmpFloat2 -= (self->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.f) { + tmpFloat3 = 0.f; + } + if (tmpFloat3 > 1.f) { + tmpFloat3 = 1.f; + } + self->pinkNoiseExp += tmpFloat3; + + // Calculate frequency independent parts of parametric noise estimate. + if (self->pinkNoiseExp > 0.f) { + // Use pink noise estimate. + parametric_num = + expf(self->pinkNoiseNumerator / (float)(self->blockInd + 1)); + parametric_num *= (float)(self->blockInd + 1); + parametric_exp = self->pinkNoiseExp / (float)(self->blockInd + 1); + } + for (i = 0; i < self->magnLen; i++) { + // Estimate the background noise using the white and pink noise + // parameters. + if (self->pinkNoiseExp == 0.f) { + // Use white noise estimate. + self->parametricNoise[i] = self->whiteNoiseLevel; + } else { + // Use pink noise estimate. + 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 / (float)(self->blockInd + 1)); + noise[i] /= END_STARTUP_SHORT; + } + } + // Compute average signal during END_STARTUP_LONG time: + // used to normalize spectral difference measure. + if (self->blockInd < END_STARTUP_LONG) { + self->featureData[5] *= self->blockInd; + self->featureData[5] += signalEnergy; + self->featureData[5] /= (self->blockInd + 1); + } + + // Post and prior SNR needed for SpeechNoiseProb. + ComputeSnr(self, magn, noise, snrLocPrior, snrLocPost); + + FeatureUpdate(self, magn, updateParsFlag); + SpeechNoiseProb(self, self->speechProb, snrLocPrior, snrLocPost); + UpdateNoiseEstimate(self, magn, snrLocPrior, snrLocPost, noise); + + // Keep track of noise spectrum for next frame. + memcpy(self->noise, noise, sizeof(*noise) * self->magnLen); + memcpy(self->magnPrevAnalyze, magn, sizeof(*magn) * self->magnLen); +} + +void WebRtcNs_ProcessCore(NoiseSuppressionC* self, + const float* const* speechFrame, + size_t num_bands, + float* const* outFrame) { + // Main routine for noise reduction. + int flagHB = 0; + size_t i, j; + + float energy1, energy2, gain, factor, factor1, factor2; + float fout[BLOCKL_MAX]; + float winData[ANAL_BLOCKL_MAX]; + float magn[HALF_ANAL_BLOCKL]; + float theFilter[HALF_ANAL_BLOCKL], theFilterTmp[HALF_ANAL_BLOCKL]; + float real[ANAL_BLOCKL_MAX], imag[HALF_ANAL_BLOCKL]; + + // SWB variables. + int deltaBweHB = 1; + int deltaGainHB = 1; + float decayBweHB = 1.0; + float gainMapParHB = 1.0; + float gainTimeDomainHB = 1.0; + float avgProbSpeechHB, avgProbSpeechHBTmp, avgFilterGainHB, gainModHB; + float sumMagnAnalyze, sumMagnProcess; + + // Check that initiation has been done. + assert(self->initFlag == 1); + assert((num_bands - 1) <= NUM_HIGH_BANDS_MAX); + + const float* const* speechFrameHB = NULL; + float* const* outFrameHB = NULL; + size_t num_high_bands = 0; + if (num_bands > 1) { + speechFrameHB = &speechFrame[1]; + outFrameHB = &outFrame[1]; + num_high_bands = num_bands - 1; + flagHB = 1; + // Range for averaging low band quantities for H band gain. + deltaBweHB = (int)self->magnLen / 4; + deltaGainHB = deltaBweHB; + } + + // Update analysis buffer for L band. + UpdateBuffer(speechFrame[0], self->blockLen, self->anaLen, self->dataBuf); + + if (flagHB == 1) { + // Update analysis buffer for H bands. + for (i = 0; i < num_high_bands; ++i) { + UpdateBuffer(speechFrameHB[i], + self->blockLen, + self->anaLen, + self->dataBufHB[i]); + } + } + + Windowing(self->window, self->dataBuf, self->anaLen, winData); + energy1 = Energy(winData, self->anaLen); + if (energy1 == 0.0) { + // Synthesize the special case of zero input. + // Read out fully processed segment. + for (i = self->windShift; i < self->blockLen + self->windShift; i++) { + fout[i - self->windShift] = self->syntBuf[i]; + } + // Update synthesis buffer. + UpdateBuffer(NULL, self->blockLen, self->anaLen, self->syntBuf); + + for (i = 0; i < self->blockLen; ++i) + outFrame[0][i] = + WEBRTC_SPL_SAT(WEBRTC_SPL_WORD16_MAX, fout[i], WEBRTC_SPL_WORD16_MIN); + + // For time-domain gain of HB. + if (flagHB == 1) { + for (i = 0; i < num_high_bands; ++i) { + for (j = 0; j < self->blockLen; ++j) { + outFrameHB[i][j] = WEBRTC_SPL_SAT(WEBRTC_SPL_WORD16_MAX, + self->dataBufHB[i][j], + WEBRTC_SPL_WORD16_MIN); + } + } + } + + return; + } + + FFT(self, winData, self->anaLen, self->magnLen, real, imag, magn); + + if (self->blockInd < END_STARTUP_SHORT) { + for (i = 0; i < self->magnLen; i++) { + self->initMagnEst[i] += magn[i]; + } + } + + ComputeDdBasedWienerFilter(self, magn, theFilter); + + for (i = 0; i < self->magnLen; i++) { + // Flooring bottom. + if (theFilter[i] < self->denoiseBound) { + theFilter[i] = self->denoiseBound; + } + // Flooring top. + if (theFilter[i] > 1.f) { + theFilter[i] = 1.f; + } + if (self->blockInd < END_STARTUP_SHORT) { + theFilterTmp[i] = + (self->initMagnEst[i] - self->overdrive * self->parametricNoise[i]); + theFilterTmp[i] /= (self->initMagnEst[i] + 0.0001f); + // Flooring bottom. + if (theFilterTmp[i] < self->denoiseBound) { + theFilterTmp[i] = self->denoiseBound; + } + // Flooring top. + if (theFilterTmp[i] > 1.f) { + theFilterTmp[i] = 1.f; + } + // Weight the two suppression filters. + theFilter[i] *= (self->blockInd); + theFilterTmp[i] *= (END_STARTUP_SHORT - self->blockInd); + theFilter[i] += theFilterTmp[i]; + theFilter[i] /= (END_STARTUP_SHORT); + } + + self->smooth[i] = theFilter[i]; + real[i] *= self->smooth[i]; + imag[i] *= self->smooth[i]; + } + // Keep track of |magn| spectrum for next frame. + memcpy(self->magnPrevProcess, magn, sizeof(*magn) * self->magnLen); + memcpy(self->noisePrev, self->noise, sizeof(self->noise[0]) * self->magnLen); + // Back to time domain. + IFFT(self, real, imag, self->magnLen, self->anaLen, winData); + + // 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 = (float)sqrt(energy2 / (energy1 + 1.f)); + + // 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. + + Windowing(self->window, winData, self->anaLen, winData); + + // Synthesis. + for (i = 0; i < self->anaLen; i++) { + self->syntBuf[i] += factor * winData[i]; + } + // Read out fully processed segment. + for (i = self->windShift; i < self->blockLen + self->windShift; i++) { + fout[i - self->windShift] = self->syntBuf[i]; + } + // Update synthesis buffer. + UpdateBuffer(NULL, self->blockLen, self->anaLen, self->syntBuf); + + for (i = 0; i < self->blockLen; ++i) + outFrame[0][i] = + WEBRTC_SPL_SAT(WEBRTC_SPL_WORD16_MAX, fout[i], WEBRTC_SPL_WORD16_MIN); + + // For time-domain gain of HB. + if (flagHB == 1) { + // Average speech prob from low band. + // Average over second half (i.e., 4->8kHz) of frequencies spectrum. + avgProbSpeechHB = 0.0; + for (i = self->magnLen - deltaBweHB - 1; i < self->magnLen - 1; i++) { + avgProbSpeechHB += self->speechProb[i]; + } + avgProbSpeechHB = avgProbSpeechHB / ((float)deltaBweHB); + // If the speech was suppressed by a component between Analyze and + // Process, for example the AEC, then it should not be considered speech + // for high band suppression purposes. + sumMagnAnalyze = 0; + sumMagnProcess = 0; + for (i = 0; i < self->magnLen; ++i) { + sumMagnAnalyze += self->magnPrevAnalyze[i]; + sumMagnProcess += self->magnPrevProcess[i]; + } + avgProbSpeechHB *= sumMagnProcess / sumMagnAnalyze; + // Average filter gain from low band. + // Average over second half (i.e., 4->8kHz) of frequencies spectrum. + avgFilterGainHB = 0.0; + for (i = self->magnLen - deltaGainHB - 1; i < self->magnLen - 1; i++) { + avgFilterGainHB += self->smooth[i]; + } + avgFilterGainHB = avgFilterGainHB / ((float)(deltaGainHB)); + avgProbSpeechHBTmp = 2.f * avgProbSpeechHB - 1.f; + // Gain based on speech probability. + gainModHB = 0.5f * (1.f + (float)tanh(gainMapParHB * avgProbSpeechHBTmp)); + // Combine gain with low band gain. + gainTimeDomainHB = 0.5f * gainModHB + 0.5f * avgFilterGainHB; + if (avgProbSpeechHB >= 0.5f) { + gainTimeDomainHB = 0.25f * gainModHB + 0.75f * avgFilterGainHB; + } + gainTimeDomainHB = gainTimeDomainHB * decayBweHB; + // Make sure gain is within flooring range. + // Flooring bottom. + if (gainTimeDomainHB < self->denoiseBound) { + gainTimeDomainHB = self->denoiseBound; + } + // Flooring top. + if (gainTimeDomainHB > 1.f) { + gainTimeDomainHB = 1.f; + } + // Apply gain. + for (i = 0; i < num_high_bands; ++i) { + for (j = 0; j < self->blockLen; j++) { + outFrameHB[i][j] = + WEBRTC_SPL_SAT(WEBRTC_SPL_WORD16_MAX, + gainTimeDomainHB * self->dataBufHB[i][j], + WEBRTC_SPL_WORD16_MIN); + } + } + } // End of H band gain computation. +} |