Merge changes I7bbf776e,I1b827825

am: fe8b4a6579

* commit 'fe8b4a657979b49e1701bd92f6d5814a99e0b2be': (7237 commits)
  WIP: Changes after merge commit 'cb3f9bd'
  Make the nonlinear beamformer steerable
  Utilize bitrate above codec max to protect video.
  Enable VP9 internal resize by default.
  Filter overlapping RTP header extensions.
  Make VCMEncodedFrameCallback const.
  MediaCodecVideoEncoder: Add number of quality resolution downscales to Encoded callback.
  Remove redudant encoder rate calls.
  Create isolate files for nonparallel tests.
  Register header extensions in RtpRtcpObserver to avoid log spam.
  Make an enum class out of NetEqDecoder, and hide the neteq_decoders_ table
  ACM: Move NACK functionality inside NetEq
  Fix chromium-style warnings in webrtc/sound/.
  Create a 'webrtc_nonparallel_tests' target.
  Update scalability structure data according to updates in the RTP payload profile.
  audio_coding: rename interface -> include
  Rewrote perform_action_on_all_files to be parallell.
  Update reference indices according to updates in the RTP payload profile.
  Disable P2PTransport...TestFailoverControlledSide on Memcheck
  pass clangcl compile options to ignore warnings in gflags.cc
  ...

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.
+}