/*
 *  Copyright (c) 2011 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 <string.h>
#include <math.h>
//#include <stdio.h>
#include <stdlib.h>
#include "noise_suppression.h"
#include "ns_core.h"
#include "windows_private.h"
#include "fft4g.h"
#include "signal_processing_library.h"

// Set Feature Extraction Parameters
void WebRtcNs_set_feature_extraction_parameters(NSinst_t* inst) {
  //bin size of histogram
  inst->featureExtractionParams.binSizeLrt      = (float)0.1;
  inst->featureExtractionParams.binSizeSpecFlat = (float)0.05;
  inst->featureExtractionParams.binSizeSpecDiff = (float)0.1;

  //range of histogram over which lrt threshold is computed
  inst->featureExtractionParams.rangeAvgHistLrt = (float)1.0;

  //scale parameters: multiply dominant peaks of the histograms by scale factor to obtain
  // thresholds for prior model
  inst->featureExtractionParams.factor1ModelPars = (float)1.20; //for lrt and spectral diff
  inst->featureExtractionParams.factor2ModelPars = (float)0.9;  //for spectral_flatness:
  // used when noise is flatter than speech

  //peak limit for spectral flatness (varies between 0 and 1)
  inst->featureExtractionParams.thresPosSpecFlat = (float)0.6;

  //limit on spacing of two highest peaks in histogram: spacing determined by bin size
  inst->featureExtractionParams.limitPeakSpacingSpecFlat = 
      2 * inst->featureExtractionParams.binSizeSpecFlat;
  inst->featureExtractionParams.limitPeakSpacingSpecDiff =
      2 * inst->featureExtractionParams.binSizeSpecDiff;

  //limit on relevance of second peak:
  inst->featureExtractionParams.limitPeakWeightsSpecFlat = (float)0.5;
  inst->featureExtractionParams.limitPeakWeightsSpecDiff = (float)0.5;

  // fluctuation limit of lrt feature
  inst->featureExtractionParams.thresFluctLrt = (float)0.05;

  //limit on the max and min values for the feature thresholds
  inst->featureExtractionParams.maxLrt = (float)1.0;
  inst->featureExtractionParams.minLrt = (float)0.20;

  inst->featureExtractionParams.maxSpecFlat = (float)0.95;
  inst->featureExtractionParams.minSpecFlat = (float)0.10;

  inst->featureExtractionParams.maxSpecDiff = (float)1.0;
  inst->featureExtractionParams.minSpecDiff = (float)0.16;

  //criteria of weight of histogram peak  to accept/reject feature
  inst->featureExtractionParams.thresWeightSpecFlat = (int)(0.3
      * (inst->modelUpdatePars[1])); //for spectral flatness
  inst->featureExtractionParams.thresWeightSpecDiff = (int)(0.3
      * (inst->modelUpdatePars[1])); //for spectral difference
}

// Initialize state
int WebRtcNs_InitCore(NSinst_t* inst, WebRtc_UWord32 fs) {
  int i;
  //We only support 10ms frames

  //check for valid pointer
  if (inst == NULL) {
    return -1;
  }

  // Initialization of struct
  if (fs == 8000 || fs == 16000 || fs == 32000) {
    inst->fs = fs;
  } else {
    return -1;
  }
  inst->windShift = 0;
  if (fs == 8000) {
    // We only support 10ms frames
    inst->blockLen = 80;
    inst->blockLen10ms = 80;
    inst->anaLen = 128;
    inst->window = kBlocks80w128;
    inst->outLen = 0;
  } else if (fs == 16000) {
    // We only support 10ms frames
    inst->blockLen = 160;
    inst->blockLen10ms = 160;
    inst->anaLen = 256;
    inst->window = kBlocks160w256;
    inst->outLen = 0;
  } else if (fs == 32000) {
    // We only support 10ms frames
    inst->blockLen = 160;
    inst->blockLen10ms = 160;
    inst->anaLen = 256;
    inst->window = kBlocks160w256;
    inst->outLen = 0;
  }
  inst->magnLen = inst->anaLen / 2 + 1; // Number of frequency bins

  // Initialize fft work arrays.
  inst->ip[0] = 0; // Setting this triggers initialization.
  memset(inst->dataBuf, 0, sizeof(float) * ANAL_BLOCKL_MAX);
  WebRtc_rdft(inst->anaLen, 1, inst->dataBuf, inst->ip, inst->wfft);

  memset(inst->dataBuf, 0, sizeof(float) * ANAL_BLOCKL_MAX);
  memset(inst->syntBuf, 0, sizeof(float) * ANAL_BLOCKL_MAX);

  //for HB processing
  memset(inst->dataBufHB, 0, sizeof(float) * ANAL_BLOCKL_MAX);

  //for quantile noise estimation
  memset(inst->quantile, 0, sizeof(float) * HALF_ANAL_BLOCKL);
  for (i = 0; i < SIMULT * HALF_ANAL_BLOCKL; i++) {
    inst->lquantile[i] = (float)8.0;
    inst->density[i] = (float)0.3;
  }

  for (i = 0; i < SIMULT; i++) {
    inst->counter[i] = (int)floor((float)(END_STARTUP_LONG * (i + 1)) / (float)SIMULT);
  }

  inst->updates = 0;

  // Wiener filter initialization
  for (i = 0; i < HALF_ANAL_BLOCKL; i++) {
    inst->smooth[i] = (float)1.0;
  }

  // Set the aggressiveness: default
  inst->aggrMode = 0;

  //initialize variables for new method
  inst->priorSpeechProb = (float)0.5; //prior prob for speech/noise
  for (i = 0; i < HALF_ANAL_BLOCKL; i++) {
    inst->magnPrev[i]      = (float)0.0; //previous mag spectrum
    inst->noisePrev[i]     = (float)0.0; //previous noise-spectrum
    inst->logLrtTimeAvg[i] = LRT_FEATURE_THR; //smooth LR ratio (same as threshold)
    inst->magnAvgPause[i]  = (float)0.0; //conservative noise spectrum estimate
    inst->speechProbHB[i]  = (float)0.0; //for estimation of HB in second pass
    inst->initMagnEst[i]   = (float)0.0; //initial average mag spectrum
  }

  //feature quantities
  inst->featureData[0] = SF_FEATURE_THR;  //spectral flatness (start on threshold)
  inst->featureData[1] = (float)0.0;      //spectral entropy: not used in this version
  inst->featureData[2] = (float)0.0;      //spectral variance: not used in this version
  inst->featureData[3] = LRT_FEATURE_THR; //average lrt factor (start on threshold)
  inst->featureData[4] = SF_FEATURE_THR;  //spectral template diff (start on threshold)
  inst->featureData[5] = (float)0.0;      //normalization for spectral-diff
  inst->featureData[6] = (float)0.0;      //window time-average of input magnitude spectrum

  //histogram quantities: used to estimate/update thresholds for features
  for (i = 0; i < HIST_PAR_EST; i++) {
    inst->histLrt[i] = 0;
    inst->histSpecFlat[i] = 0;
    inst->histSpecDiff[i] = 0;
  }

  inst->blockInd = -1; //frame counter
  inst->priorModelPars[0] = LRT_FEATURE_THR; //default threshold for lrt feature
  inst->priorModelPars[1] = (float)0.5;      //threshold for spectral flatness:
  // determined on-line
  inst->priorModelPars[2] = (float)1.0;      //sgn_map par for spectral measure:
  // 1 for flatness measure
  inst->priorModelPars[3] = (float)0.5;      //threshold for template-difference feature:
  // determined on-line
  inst->priorModelPars[4] = (float)1.0;      //default weighting parameter for lrt feature
  inst->priorModelPars[5] = (float)0.0;      //default weighting parameter for
  // spectral flatness feature
  inst->priorModelPars[6] = (float)0.0;      //default weighting parameter for
  // spectral difference feature

  inst->modelUpdatePars[0] = 2;   //update flag for parameters:
  // 0 no update, 1=update once, 2=update every window
  inst->modelUpdatePars[1] = 500; //window for update
  inst->modelUpdatePars[2] = 0;   //counter for update of conservative noise spectrum
  //counter if the feature thresholds are updated during the sequence
  inst->modelUpdatePars[3] = inst->modelUpdatePars[1];

  inst->signalEnergy = 0.0;
  inst->sumMagn = 0.0;
  inst->whiteNoiseLevel = 0.0;
  inst->pinkNoiseNumerator = 0.0;
  inst->pinkNoiseExp = 0.0;

  WebRtcNs_set_feature_extraction_parameters(inst); // Set feature configuration

  //default mode
  WebRtcNs_set_policy_core(inst, 0);


  memset(inst->outBuf, 0, sizeof(float) * 3 * BLOCKL_MAX);

  inst->initFlag = 1;
  return 0;
}

int WebRtcNs_set_policy_core(NSinst_t* inst, int mode) {
  // allow for modes:0,1,2,3
  if (mode < 0 || mode > 3) {
    return (-1);
  }

  inst->aggrMode = mode;
  if (mode == 0) {
    inst->overdrive = (float)1.0;
    inst->denoiseBound = (float)0.5;
    inst->gainmap = 0;
  } else if (mode == 1) {
    //inst->overdrive = (float)1.25;
    inst->overdrive = (float)1.0;
    inst->denoiseBound = (float)0.25;
    inst->gainmap = 1;
  } else if (mode == 2) {
    //inst->overdrive = (float)1.25;
    inst->overdrive = (float)1.1;
    inst->denoiseBound = (float)0.125;
    inst->gainmap = 1;
  } else if (mode == 3) {
    //inst->overdrive = (float)1.30;
    inst->overdrive = (float)1.25;
    inst->denoiseBound = (float)0.09;
    inst->gainmap = 1;
  }
  return 0;
}

// Estimate noise
void WebRtcNs_NoiseEstimation(NSinst_t* inst, float* magn, float* noise) {
  int i, s, offset;
  float lmagn[HALF_ANAL_BLOCKL], delta;

  if (inst->updates < END_STARTUP_LONG) {
    inst->updates++;
  }

  for (i = 0; i < inst->magnLen; i++) {
    lmagn[i] = (float)log(magn[i]);
  }

  // loop over simultaneous estimates
  for (s = 0; s < SIMULT; s++) {
    offset = s * inst->magnLen;

    // newquantest(...)
    for (i = 0; i < inst->magnLen; i++) {
      // compute delta
      if (inst->density[offset + i] > 1.0) {
        delta = FACTOR * (float)1.0 / inst->density[offset + i];
      } else {
        delta = FACTOR;
      }

      // update log quantile estimate
      if (lmagn[i] > inst->lquantile[offset + i]) {
        inst->lquantile[offset + i] += QUANTILE * delta
                                       / (float)(inst->counter[s] + 1);
      } else {
        inst->lquantile[offset + i] -= ((float)1.0 - QUANTILE) * delta
                                       / (float)(inst->counter[s] + 1);
      }

      // update density estimate
      if (fabs(lmagn[i] - inst->lquantile[offset + i]) < WIDTH) {
        inst->density[offset + i] = ((float)inst->counter[s] * inst->density[offset
            + i] + (float)1.0 / ((float)2.0 * WIDTH)) / (float)(inst->counter[s] + 1);
      }
    } // end loop over magnitude spectrum

    if (inst->counter[s] >= END_STARTUP_LONG) {
      inst->counter[s] = 0;
      if (inst->updates >= END_STARTUP_LONG) {
        for (i = 0; i < inst->magnLen; i++) {
          inst->quantile[i] = (float)exp(inst->lquantile[offset + i]);
        }
      }
    }

    inst->counter[s]++;
  } // end loop over simultaneous estimates

  // Sequentially update the noise during startup
  if (inst->updates < END_STARTUP_LONG) {
    // Use the last "s" to get noise during startup that differ from zero.
    for (i = 0; i < inst->magnLen; i++) {
      inst->quantile[i] = (float)exp(inst->lquantile[offset + i]);
    }
  }

  for (i = 0; i < inst->magnLen; i++) {
    noise[i] = inst->quantile[i];
  }
}

// Extract thresholds for feature parameters
// histograms are computed over some window_size (given by inst->modelUpdatePars[1])
// thresholds and weights are extracted every window
// flag 0 means update histogram only, flag 1 means compute the thresholds/weights
// threshold and weights are returned in: inst->priorModelPars
void WebRtcNs_FeatureParameterExtraction(NSinst_t* inst, int flag) {
  int i, useFeatureSpecFlat, useFeatureSpecDiff, numHistLrt;
  int maxPeak1, maxPeak2;
  int weightPeak1SpecFlat, weightPeak2SpecFlat, weightPeak1SpecDiff, weightPeak2SpecDiff;

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

  //3 features: lrt, flatness, difference
  //lrt_feature = inst->featureData[3];
  //flat_feature = inst->featureData[0];
  //diff_feature = inst->featureData[4];

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

  // extract parameters for speech/noise probability
  if (flag == 1) {
    //lrt feature: compute the average over inst->featureExtractionParams.rangeAvgHistLrt
    avgHistLrt = 0.0;
    avgHistLrtCompl = 0.0;
    avgSquareHistLrt = 0.0;
    numHistLrt = 0;
    for (i = 0; i < HIST_PAR_EST; i++) {
      binMid = ((float)i + (float)0.5) * inst->featureExtractionParams.binSizeLrt;
      if (binMid <= inst->featureExtractionParams.rangeAvgHistLrt) {
        avgHistLrt += inst->histLrt[i] * binMid;
        numHistLrt += inst->histLrt[i];
      }
      avgSquareHistLrt += inst->histLrt[i] * binMid * binMid;
      avgHistLrtCompl += inst->histLrt[i] * binMid;
    }
    if (numHistLrt > 0) {
      avgHistLrt = avgHistLrt / ((float)numHistLrt);
    }
    avgHistLrtCompl = avgHistLrtCompl / ((float)inst->modelUpdatePars[1]);
    avgSquareHistLrt = avgSquareHistLrt / ((float)inst->modelUpdatePars[1]);
    fluctLrt = avgSquareHistLrt - avgHistLrt * avgHistLrtCompl;
    // get threshold for lrt feature:
    if (fluctLrt < inst->featureExtractionParams.thresFluctLrt) {
      //very low fluct, so likely noise
      inst->priorModelPars[0] = inst->featureExtractionParams.maxLrt;
    } else {
      inst->priorModelPars[0] = inst->featureExtractionParams.factor1ModelPars
                                * avgHistLrt;
      // check if value is within min/max range
      if (inst->priorModelPars[0] < inst->featureExtractionParams.minLrt) {
        inst->priorModelPars[0] = inst->featureExtractionParams.minLrt;
      }
      if (inst->priorModelPars[0] > inst->featureExtractionParams.maxLrt) {
        inst->priorModelPars[0] = inst->featureExtractionParams.maxLrt;
      }
    }
    // done with lrt feature

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

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

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

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

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

    // for spectrum flatness feature
    useFeatureSpecFlat = 1;
    // merge the two peaks if they are close
    if ((fabs(posPeak2SpecFlat - posPeak1SpecFlat)
         < inst->featureExtractionParams.limitPeakSpacingSpecFlat)
        && (weightPeak2SpecFlat
            > inst->featureExtractionParams.limitPeakWeightsSpecFlat
            * weightPeak1SpecFlat)) {
      weightPeak1SpecFlat += weightPeak2SpecFlat;
      posPeak1SpecFlat = (float)0.5 * (posPeak1SpecFlat + posPeak2SpecFlat);
    }
    //reject if weight of peaks is not large enough, or peak value too small
    if (weightPeak1SpecFlat < inst->featureExtractionParams.thresWeightSpecFlat
        || posPeak1SpecFlat < inst->featureExtractionParams.thresPosSpecFlat) {
      useFeatureSpecFlat = 0;
    }
    // if selected, get the threshold
    if (useFeatureSpecFlat == 1) {
      // compute the threshold
      inst->priorModelPars[1] = inst->featureExtractionParams.factor2ModelPars
                                * posPeak1SpecFlat;
      //check if value is within min/max range
      if (inst->priorModelPars[1] < inst->featureExtractionParams.minSpecFlat) {
        inst->priorModelPars[1] = inst->featureExtractionParams.minSpecFlat;
      }
      if (inst->priorModelPars[1] > inst->featureExtractionParams.maxSpecFlat) {
        inst->priorModelPars[1] = inst->featureExtractionParams.maxSpecFlat;
      }
    }
    // done with flatness feature

    // for template feature
    useFeatureSpecDiff = 1;
    // merge the two peaks if they are close
    if ((fabs(posPeak2SpecDiff - posPeak1SpecDiff)
         < inst->featureExtractionParams.limitPeakSpacingSpecDiff)
        && (weightPeak2SpecDiff
            > inst->featureExtractionParams.limitPeakWeightsSpecDiff
            * weightPeak1SpecDiff)) {
      weightPeak1SpecDiff += weightPeak2SpecDiff;
      posPeak1SpecDiff = (float)0.5 * (posPeak1SpecDiff + posPeak2SpecDiff);
    }
    // get the threshold value
    inst->priorModelPars[3] = inst->featureExtractionParams.factor1ModelPars
                              * posPeak1SpecDiff;
    //reject if weight of peaks is not large enough
    if (weightPeak1SpecDiff < inst->featureExtractionParams.thresWeightSpecDiff) {
      useFeatureSpecDiff = 0;
    }
    //check if value is within min/max range
    if (inst->priorModelPars[3] < inst->featureExtractionParams.minSpecDiff) {
      inst->priorModelPars[3] = inst->featureExtractionParams.minSpecDiff;
    }
    if (inst->priorModelPars[3] > inst->featureExtractionParams.maxSpecDiff) {
      inst->priorModelPars[3] = inst->featureExtractionParams.maxSpecDiff;
    }
    // done with spectral difference feature

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

    // select the weights between the features
    // inst->priorModelPars[4] is weight for lrt: always selected
    // inst->priorModelPars[5] is weight for spectral flatness
    // inst->priorModelPars[6] is weight for spectral difference
    featureSum = (float)(1 + useFeatureSpecFlat + useFeatureSpecDiff);
    inst->priorModelPars[4] = (float)1.0 / featureSum;
    inst->priorModelPars[5] = ((float)useFeatureSpecFlat) / featureSum;
    inst->priorModelPars[6] = ((float)useFeatureSpecDiff) / featureSum;

    // set hists to zero for next update
    if (inst->modelUpdatePars[0] >= 1) {
      for (i = 0; i < HIST_PAR_EST; i++) {
        inst->histLrt[i] = 0;
        inst->histSpecFlat[i] = 0;
        inst->histSpecDiff[i] = 0;
      }
    }
  } // end of flag == 1
}

// Compute spectral flatness on input spectrum
// magnIn is the magnitude spectrum
// spectral flatness is returned in inst->featureData[0]
void WebRtcNs_ComputeSpectralFlatness(NSinst_t* inst, float* magnIn) {
  int i;
  int shiftLP = 1; //option to remove first bin(s) from spectral measures
  float avgSpectralFlatnessNum, avgSpectralFlatnessDen, spectralTmp;

  // comute spectral measures
  // for flatness
  avgSpectralFlatnessNum = 0.0;
  avgSpectralFlatnessDen = inst->sumMagn;
  for (i = 0; i < shiftLP; i++) {
    avgSpectralFlatnessDen -= magnIn[i];
  }
  // compute log of ratio of the geometric to arithmetic mean: check for log(0) case
  for (i = shiftLP; i < inst->magnLen; i++) {
    if (magnIn[i] > 0.0) {
      avgSpectralFlatnessNum += (float)log(magnIn[i]);
    } else {
      inst->featureData[0] -= SPECT_FL_TAVG * inst->featureData[0];
      return;
    }
  }
  //normalize
  avgSpectralFlatnessDen = avgSpectralFlatnessDen / inst->magnLen;
  avgSpectralFlatnessNum = avgSpectralFlatnessNum / inst->magnLen;

  //ratio and inverse log: check for case of log(0)
  spectralTmp = (float)exp(avgSpectralFlatnessNum) / avgSpectralFlatnessDen;

  //time-avg update of spectral flatness feature
  inst->featureData[0] += SPECT_FL_TAVG * (spectralTmp - inst->featureData[0]);
  // done with flatness feature
}

// Compute the difference measure between input spectrum and a template/learned noise spectrum
// magnIn is the input spectrum
// the reference/template spectrum is inst->magnAvgPause[i]
// returns (normalized) spectral difference in inst->featureData[4]
void WebRtcNs_ComputeSpectralDifference(NSinst_t* inst, float* magnIn) {
  // avgDiffNormMagn = var(magnIn) - cov(magnIn, magnAvgPause)^2 / var(magnAvgPause)
  int i;
  float avgPause, avgMagn, covMagnPause, varPause, varMagn, avgDiffNormMagn;

  avgPause = 0.0;
  avgMagn = inst->sumMagn;
  // compute average quantities
  for (i = 0; i < inst->magnLen; i++) {
    //conservative smooth noise spectrum from pause frames
    avgPause += inst->magnAvgPause[i];
  }
  avgPause = avgPause / ((float)inst->magnLen);
  avgMagn = avgMagn / ((float)inst->magnLen);

  covMagnPause = 0.0;
  varPause = 0.0;
  varMagn = 0.0;
  // compute variance and covariance quantities
  for (i = 0; i < inst->magnLen; i++) {
    covMagnPause += (magnIn[i] - avgMagn) * (inst->magnAvgPause[i] - avgPause);
    varPause += (inst->magnAvgPause[i] - avgPause) * (inst->magnAvgPause[i] - avgPause);
    varMagn += (magnIn[i] - avgMagn) * (magnIn[i] - avgMagn);
  }
  covMagnPause = covMagnPause / ((float)inst->magnLen);
  varPause = varPause / ((float)inst->magnLen);
  varMagn = varMagn / ((float)inst->magnLen);
  // update of average magnitude spectrum
  inst->featureData[6] += inst->signalEnergy;

  avgDiffNormMagn = varMagn - (covMagnPause * covMagnPause) / (varPause + (float)0.0001);
  // normalize and compute time-avg update of difference feature
  avgDiffNormMagn = (float)(avgDiffNormMagn / (inst->featureData[5] + (float)0.0001));
  inst->featureData[4] += SPECT_DIFF_TAVG * (avgDiffNormMagn - inst->featureData[4]);
}

// 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 freq.
//snr loc_post is the post snr for each freq.
void WebRtcNs_SpeechNoiseProb(NSinst_t* inst, float* probSpeechFinal, float* snrLocPrior,
                              float* snrLocPost) {
  int i, sgnMap;
  float invLrt, gainPrior, indPrior;
  float logLrtTimeAvgKsum, besselTmp;
  float indicator0, indicator1, indicator2;
  float tmpFloat1, tmpFloat2;
  float weightIndPrior0, weightIndPrior1, weightIndPrior2;
  float threshPrior0, threshPrior1, threshPrior2;
  float widthPrior, widthPrior0, widthPrior1, widthPrior2;

  widthPrior0 = WIDTH_PR_MAP;
  widthPrior1 = (float)2.0 * WIDTH_PR_MAP; //width for pause region:
  // lower range, so increase width in tanh map
  widthPrior2 = (float)2.0 * WIDTH_PR_MAP; //for spectral-difference measure

  //threshold parameters for features
  threshPrior0 = inst->priorModelPars[0];
  threshPrior1 = inst->priorModelPars[1];
  threshPrior2 = inst->priorModelPars[3];

  //sign for flatness feature
  sgnMap = (int)(inst->priorModelPars[2]);

  //weight parameters for features
  weightIndPrior0 = inst->priorModelPars[4];
  weightIndPrior1 = inst->priorModelPars[5];
  weightIndPrior2 = inst->priorModelPars[6];

  // compute feature based on average LR factor
  // this is the average over all frequencies of the smooth log lrt
  logLrtTimeAvgKsum = 0.0;
  for (i = 0; i < inst->magnLen; i++) {
    tmpFloat1 = (float)1.0 + (float)2.0 * snrLocPrior[i];
    tmpFloat2 = (float)2.0 * snrLocPrior[i] / (tmpFloat1 + (float)0.0001);
    besselTmp = (snrLocPost[i] + (float)1.0) * tmpFloat2;
    inst->logLrtTimeAvg[i] += LRT_TAVG * (besselTmp - (float)log(tmpFloat1)
                                          - inst->logLrtTimeAvg[i]);
    logLrtTimeAvgKsum += inst->logLrtTimeAvg[i];
  }
  logLrtTimeAvgKsum = (float)logLrtTimeAvgKsum / (inst->magnLen);
  inst->featureData[3] = logLrtTimeAvgKsum;
  // 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 = (float)0.5 * ((float)tanh(widthPrior *
      (logLrtTimeAvgKsum - threshPrior0)) + (float)1.0);

  //spectral flatness feature
  tmpFloat1 = inst->featureData[0];
  widthPrior = widthPrior0;
  //use larger width in tanh map for pause regions
  if (sgnMap == 1 && (tmpFloat1 > threshPrior1)) {
    widthPrior = widthPrior1;
  }
  if (sgnMap == -1 && (tmpFloat1 < threshPrior1)) {
    widthPrior = widthPrior1;
  }
  // compute indicator function: sigmoid map
  indicator1 = (float)0.5 * ((float)tanh((float)sgnMap * 
      widthPrior * (threshPrior1 - tmpFloat1)) + (float)1.0);

  //for template spectrum-difference
  tmpFloat1 = inst->featureData[4];
  widthPrior = widthPrior0;
  //use larger width in tanh map for pause regions
  if (tmpFloat1 < threshPrior2) {
    widthPrior = widthPrior2;
  }
  // compute indicator function: sigmoid map
  indicator2 = (float)0.5 * ((float)tanh(widthPrior * (tmpFloat1 - threshPrior2))
                             + (float)1.0);

  //combine the indicator function with the feature weights
  indPrior = weightIndPrior0 * indicator0 + weightIndPrior1 * indicator1 + weightIndPrior2
             * indicator2;
  // done with computing indicator function

  //compute the prior probability
  inst->priorSpeechProb += PRIOR_UPDATE * (indPrior - inst->priorSpeechProb);
  // make sure probabilities are within range: keep floor to 0.01
  if (inst->priorSpeechProb > 1.0) {
    inst->priorSpeechProb = (float)1.0;
  }
  if (inst->priorSpeechProb < 0.01) {
    inst->priorSpeechProb = (float)0.01;
  }

  //final speech probability: combine prior model with LR factor:
  gainPrior = ((float)1.0 - inst->priorSpeechProb) / (inst->priorSpeechProb + (float)0.0001);
  for (i = 0; i < inst->magnLen; i++) {
    invLrt = (float)exp(-inst->logLrtTimeAvg[i]);
    invLrt = (float)gainPrior * invLrt;
    probSpeechFinal[i] = (float)1.0 / ((float)1.0 + invLrt);
  }
}

int WebRtcNs_ProcessCore(NSinst_t* inst,
                         short* speechFrame,
                         short* speechFrameHB,
                         short* outFrame,
                         short* outFrameHB) {
  // main routine for noise reduction

  int     flagHB = 0;
  int     i;
  const int kStartBand = 5; // Skip first frequency bins during estimation.
  int     updateParsFlag;

  float   energy1, energy2, gain, factor, factor1, factor2;
  float   signalEnergy, sumMagn;
  float   snrPrior, currentEstimateStsa;
  float   tmpFloat1, tmpFloat2, tmpFloat3, probSpeech, probNonSpeech;
  float   gammaNoiseTmp, gammaNoiseOld;
  float   noiseUpdateTmp, fTmp, dTmp;
  float   fin[BLOCKL_MAX], fout[BLOCKL_MAX];
  float   winData[ANAL_BLOCKL_MAX];
  float   magn[HALF_ANAL_BLOCKL], noise[HALF_ANAL_BLOCKL];
  float   theFilter[HALF_ANAL_BLOCKL], theFilterTmp[HALF_ANAL_BLOCKL];
  float   snrLocPost[HALF_ANAL_BLOCKL], snrLocPrior[HALF_ANAL_BLOCKL];
  float   probSpeechFinal[HALF_ANAL_BLOCKL], previousEstimateStsa[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_noise = 0.0;
  float   parametric_exp = 0.0;
  float   parametric_num = 0.0;

  // 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;

  // Check that initiation has been done
  if (inst->initFlag != 1) {
    return (-1);
  }
  // Check for valid pointers based on sampling rate
  if (inst->fs == 32000) {
    if (speechFrameHB == NULL) {
      return -1;
    }
    flagHB = 1;
    // range for averaging low band quantities for H band gain
    deltaBweHB = (int)inst->magnLen / 4;
    deltaGainHB = deltaBweHB;
  }
  //
  updateParsFlag = inst->modelUpdatePars[0];
  //

  //for LB do all processing
  // convert to float
  for (i = 0; i < inst->blockLen10ms; i++) {
    fin[i] = (float)speechFrame[i];
  }
  // update analysis buffer for L band
  memcpy(inst->dataBuf, inst->dataBuf + inst->blockLen10ms,
         sizeof(float) * (inst->anaLen - inst->blockLen10ms));
  memcpy(inst->dataBuf + inst->anaLen - inst->blockLen10ms, fin,
         sizeof(float) * inst->blockLen10ms);

  if (flagHB == 1) {
    // convert to float
    for (i = 0; i < inst->blockLen10ms; i++) {
      fin[i] = (float)speechFrameHB[i];
    }
    // update analysis buffer for H band
    memcpy(inst->dataBufHB, inst->dataBufHB + inst->blockLen10ms,
           sizeof(float) * (inst->anaLen - inst->blockLen10ms));
    memcpy(inst->dataBufHB + inst->anaLen - inst->blockLen10ms, fin,
           sizeof(float) * inst->blockLen10ms);
  }

  // check if processing needed
  if (inst->outLen == 0) {
    // windowing
    energy1 = 0.0;
    for (i = 0; i < inst->anaLen; i++) {
      winData[i] = inst->window[i] * inst->dataBuf[i];
      energy1 += winData[i] * winData[i];
    }
    if (energy1 == 0.0) {
      // synthesize the special case of zero input
      // 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.

      // read out fully processed segment
      for (i = inst->windShift; i < inst->blockLen + inst->windShift; i++) {
        fout[i - inst->windShift] = inst->syntBuf[i];
      }
      // update synthesis buffer
      memcpy(inst->syntBuf, inst->syntBuf + inst->blockLen,
             sizeof(float) * (inst->anaLen - inst->blockLen));
      memset(inst->syntBuf + inst->anaLen - inst->blockLen, 0,
             sizeof(float) * inst->blockLen);

      // out buffer
      inst->outLen = inst->blockLen - inst->blockLen10ms;
      if (inst->blockLen > inst->blockLen10ms) {
        for (i = 0; i < inst->outLen; i++) {
          inst->outBuf[i] = fout[i + inst->blockLen10ms];
        }
      }
      // convert to short
      for (i = 0; i < inst->blockLen10ms; i++) {
        dTmp = fout[i];
        if (dTmp < WEBRTC_SPL_WORD16_MIN) {
          dTmp = WEBRTC_SPL_WORD16_MIN;
        } else if (dTmp > WEBRTC_SPL_WORD16_MAX) {
          dTmp = WEBRTC_SPL_WORD16_MAX;
        }
        outFrame[i] = (short)dTmp;
      }

      // for time-domain gain of HB
      if (flagHB == 1) {
        for (i = 0; i < inst->blockLen10ms; i++) {
          dTmp = inst->dataBufHB[i];
          if (dTmp < WEBRTC_SPL_WORD16_MIN) {
            dTmp = WEBRTC_SPL_WORD16_MIN;
          } else if (dTmp > WEBRTC_SPL_WORD16_MAX) {
            dTmp = WEBRTC_SPL_WORD16_MAX;
          }
          outFrameHB[i] = (short)dTmp;
        }
      } // end of H band gain computation
      //
      return 0;
    }

    //
    inst->blockInd++; // Update the block index only when we process a block.
    // FFT
    WebRtc_rdft(inst->anaLen, 1, winData, inst->ip, inst->wfft);

    imag[0] = 0;
    real[0] = winData[0];
    magn[0] = (float)(fabs(real[0]) + 1.0f);
    imag[inst->magnLen - 1] = 0;
    real[inst->magnLen - 1] = winData[1];
    magn[inst->magnLen - 1] = (float)(fabs(real[inst->magnLen - 1]) + 1.0f);
    signalEnergy = (float)(real[0] * real[0]) + 
                   (float)(real[inst->magnLen - 1] * real[inst->magnLen - 1]);
    sumMagn = magn[0] + magn[inst->magnLen - 1];
    if (inst->blockInd < END_STARTUP_SHORT) {
      inst->initMagnEst[0] += magn[0];
      inst->initMagnEst[inst->magnLen - 1] += magn[inst->magnLen - 1];
      tmpFloat2 = log((float)(inst->magnLen - 1));
      sum_log_i = tmpFloat2;
      sum_log_i_square = tmpFloat2 * tmpFloat2;
      tmpFloat1 = log(magn[inst->magnLen - 1]);
      sum_log_magn = tmpFloat1;
      sum_log_i_log_magn = tmpFloat2 * tmpFloat1;
    }
    for (i = 1; i < inst->magnLen - 1; i++) {
      real[i] = winData[2 * i];
      imag[i] = winData[2 * i + 1];
      // magnitude spectrum
      fTmp = real[i] * real[i];
      fTmp += imag[i] * imag[i];
      signalEnergy += fTmp;
      magn[i] = ((float)sqrt(fTmp)) + 1.0f;
      sumMagn += magn[i];
      if (inst->blockInd < END_STARTUP_SHORT) {
        inst->initMagnEst[i] += magn[i];
        if (i >= kStartBand) {
          tmpFloat2 = log((float)i);
          sum_log_i += tmpFloat2;
          sum_log_i_square += tmpFloat2 * tmpFloat2;
          tmpFloat1 = log(magn[i]);
          sum_log_magn += tmpFloat1;
          sum_log_i_log_magn += tmpFloat2 * tmpFloat1;
        }
      }
    }
    signalEnergy = signalEnergy / ((float)inst->magnLen);
    inst->signalEnergy = signalEnergy;
    inst->sumMagn = sumMagn;

    //compute spectral flatness on input spectrum
    WebRtcNs_ComputeSpectralFlatness(inst, magn);
    // quantile noise estimate
    WebRtcNs_NoiseEstimation(inst, magn, noise);
    //compute simplified noise model during startup
    if (inst->blockInd < END_STARTUP_SHORT) {
      // Estimate White noise
      inst->whiteNoiseLevel += sumMagn / ((float)inst->magnLen) * inst->overdrive;
      // Estimate Pink noise parameters
      tmpFloat1 = sum_log_i_square * ((float)(inst->magnLen - kStartBand));
      tmpFloat1 -= (sum_log_i * sum_log_i);
      tmpFloat2 = (sum_log_i_square * sum_log_magn - sum_log_i * sum_log_i_log_magn);
      tmpFloat3 = tmpFloat2 / tmpFloat1;
      // Constrain the estimated spectrum to be positive
      if (tmpFloat3 < 0.0f) {
        tmpFloat3 = 0.0f;
      }
      inst->pinkNoiseNumerator += tmpFloat3;
      tmpFloat2 = (sum_log_i * sum_log_magn);
      tmpFloat2 -= ((float)(inst->magnLen - kStartBand)) * sum_log_i_log_magn;
      tmpFloat3 = tmpFloat2 / tmpFloat1;
      // Constrain the pink noise power to be in the interval [0, 1];
      if (tmpFloat3 < 0.0f) {
        tmpFloat3 = 0.0f;
      }
      if (tmpFloat3 > 1.0f) {
        tmpFloat3 = 1.0f;
      }
      inst->pinkNoiseExp += tmpFloat3;

      // Calculate frequency independent parts of parametric noise estimate.
      if (inst->pinkNoiseExp == 0.0f) {
        // Use white noise estimate
        parametric_noise = inst->whiteNoiseLevel;
      } else {
        // Use pink noise estimate
        parametric_num = exp(inst->pinkNoiseNumerator / (float)(inst->blockInd + 1));
        parametric_num *= (float)(inst->blockInd + 1);
        parametric_exp = inst->pinkNoiseExp / (float)(inst->blockInd + 1);
        parametric_noise = parametric_num / pow((float)kStartBand, parametric_exp);
      }
      for (i = 0; i < inst->magnLen; i++) {
        // Estimate the background noise using the white and pink noise parameters
        if ((inst->pinkNoiseExp > 0.0f) && (i >= kStartBand)) {
          // Use pink noise estimate
          parametric_noise = parametric_num / pow((float)i, parametric_exp);
        }
        theFilterTmp[i] = (inst->initMagnEst[i] - inst->overdrive * parametric_noise);
        theFilterTmp[i] /= (inst->initMagnEst[i] + (float)0.0001);
        // Weight quantile noise with modeled noise
        noise[i] *= (inst->blockInd);
        tmpFloat2 = parametric_noise * (END_STARTUP_SHORT - inst->blockInd);
        noise[i] += (tmpFloat2 / (float)(inst->blockInd + 1));
        noise[i] /= END_STARTUP_SHORT;
      }
    }
    //compute average signal during END_STARTUP_LONG time:
    // used to normalize spectral difference measure
    if (inst->blockInd < END_STARTUP_LONG) {
      inst->featureData[5] *= inst->blockInd;
      inst->featureData[5] += signalEnergy;
      inst->featureData[5] /= (inst->blockInd + 1);
    }

#ifdef PROCESS_FLOW_0
    if (inst->blockInd > END_STARTUP_LONG) {
      //option: average the quantile noise: for check with AEC2
      for (i = 0; i < inst->magnLen; i++) {
        noise[i] = (float)0.6 * inst->noisePrev[i] + (float)0.4 * noise[i];
      }
      for (i = 0; i < inst->magnLen; i++) {
        // Wiener with over sub-substraction:
        theFilter[i] = (magn[i] - inst->overdrive * noise[i]) / (magn[i] + (float)0.0001);
      }
    }
#else
    //start processing at frames == converged+1
    //
    // STEP 1: compute  prior and post snr based on quantile noise est
    //

    // compute DD estimate of prior SNR: needed for new method
    for (i = 0; i < inst->magnLen; i++) {
      // post snr
      snrLocPost[i] = (float)0.0;
      if (magn[i] > noise[i]) {
        snrLocPost[i] = magn[i] / (noise[i] + (float)0.0001) - (float)1.0;
      }
      // previous post snr
      // previous estimate: based on previous frame with gain filter
      previousEstimateStsa[i] = inst->magnPrev[i] / (inst->noisePrev[i] + (float)0.0001)
                                * (inst->smooth[i]);
      // DD estimate is sum of two terms: current estimate and previous estimate
      // directed decision update of snrPrior
      snrLocPrior[i] = DD_PR_SNR * previousEstimateStsa[i] + ((float)1.0 - DD_PR_SNR)
                       * snrLocPost[i];
      // post and prior snr needed for step 2
    } // end of loop over freqs
#ifdef PROCESS_FLOW_1
    for (i = 0; i < inst->magnLen; i++) {
      // gain filter
      tmpFloat1 = inst->overdrive + snrLocPrior[i];
      tmpFloat2 = (float)snrLocPrior[i] / tmpFloat1;
      theFilter[i] = (float)tmpFloat2;
    } // end of loop over freqs
#endif
    // done with step 1: dd computation of prior and post snr

    //
    //STEP 2: compute speech/noise likelihood
    //
#ifdef PROCESS_FLOW_2
    // compute difference of input spectrum with learned/estimated noise spectrum
    WebRtcNs_ComputeSpectralDifference(inst, magn);
    // compute histograms for parameter decisions (thresholds and weights for features)
    // parameters are extracted once every window time (=inst->modelUpdatePars[1])
    if (updateParsFlag >= 1) {
      // counter update
      inst->modelUpdatePars[3]--;
      // update histogram
      if (inst->modelUpdatePars[3] > 0) {
        WebRtcNs_FeatureParameterExtraction(inst, 0);
      }
      // compute model parameters
      if (inst->modelUpdatePars[3] == 0) {
        WebRtcNs_FeatureParameterExtraction(inst, 1);
        inst->modelUpdatePars[3] = inst->modelUpdatePars[1];
        // if wish to update only once, set flag to zero
        if (updateParsFlag == 1) {
          inst->modelUpdatePars[0] = 0;
        } else {
          // update every window:
          // get normalization for spectral difference for next window estimate
          inst->featureData[6] = inst->featureData[6]
                                 / ((float)inst->modelUpdatePars[1]);
          inst->featureData[5] = (float)0.5 * (inst->featureData[6]
                                               + inst->featureData[5]);
          inst->featureData[6] = (float)0.0;
        }
      }
    }
    // compute speech/noise probability
    WebRtcNs_SpeechNoiseProb(inst, probSpeechFinal, snrLocPrior, snrLocPost);
    // time-avg parameter for noise update
    gammaNoiseTmp = NOISE_UPDATE;
    for (i = 0; i < inst->magnLen; i++) {
      probSpeech = probSpeechFinal[i];
      probNonSpeech = (float)1.0 - probSpeech;
      // temporary noise update:
      // use it for speech frames if update value is less than previous
      noiseUpdateTmp = gammaNoiseTmp * inst->noisePrev[i] + ((float)1.0 - gammaNoiseTmp)
                       * (probNonSpeech * magn[i] + probSpeech * inst->noisePrev[i]);
      //
      // time-constant based on speech/noise state
      gammaNoiseOld = gammaNoiseTmp;
      gammaNoiseTmp = NOISE_UPDATE;
      // increase gamma (i.e., less noise update) for frame likely to be speech
      if (probSpeech > PROB_RANGE) {
        gammaNoiseTmp = SPEECH_UPDATE;
      }
      // conservative noise update
      if (probSpeech < PROB_RANGE) {
        inst->magnAvgPause[i] += GAMMA_PAUSE * (magn[i] - inst->magnAvgPause[i]);
      }
      // noise update
      if (gammaNoiseTmp == gammaNoiseOld) {
        noise[i] = noiseUpdateTmp;
      } else {
        noise[i] = gammaNoiseTmp * inst->noisePrev[i] + ((float)1.0 - gammaNoiseTmp)
                   * (probNonSpeech * magn[i] + probSpeech * inst->noisePrev[i]);
        // allow for noise update downwards:
        //  if noise update decreases the noise, it is safe, so allow it to happen
        if (noiseUpdateTmp < noise[i]) {
          noise[i] = noiseUpdateTmp;
        }
      }
    } // end of freq loop
    // done with step 2: noise update

    //
    // STEP 3: compute dd update of prior snr and post snr based on new noise estimate
    //
    for (i = 0; i < inst->magnLen; i++) {
      // post and prior snr
      currentEstimateStsa = (float)0.0;
      if (magn[i] > noise[i]) {
        currentEstimateStsa = magn[i] / (noise[i] + (float)0.0001) - (float)1.0;
      }
      // DD estimate is sume of two terms: current estimate and previous estimate
      // directed decision update of snrPrior
      snrPrior = DD_PR_SNR * previousEstimateStsa[i] + ((float)1.0 - DD_PR_SNR)
                 * currentEstimateStsa;
      // gain filter
      tmpFloat1 = inst->overdrive + snrPrior;
      tmpFloat2 = (float)snrPrior / tmpFloat1;
      theFilter[i] = (float)tmpFloat2;
    } // end of loop over freqs
    // done with step3
#endif
#endif

    for (i = 0; i < inst->magnLen; i++) {
      // flooring bottom
      if (theFilter[i] < inst->denoiseBound) {
        theFilter[i] = inst->denoiseBound;
      }
      // flooring top
      if (theFilter[i] > (float)1.0) {
        theFilter[i] = 1.0;
      }
      if (inst->blockInd < END_STARTUP_SHORT) {
        // flooring bottom
        if (theFilterTmp[i] < inst->denoiseBound) {
          theFilterTmp[i] = inst->denoiseBound;
        }
        // flooring top
        if (theFilterTmp[i] > (float)1.0) {
          theFilterTmp[i] = 1.0;
        }
        // Weight the two suppression filters
        theFilter[i] *= (inst->blockInd);
        theFilterTmp[i] *= (END_STARTUP_SHORT - inst->blockInd);
        theFilter[i] += theFilterTmp[i];
        theFilter[i] /= (END_STARTUP_SHORT);
      }
      // smoothing
#ifdef PROCESS_FLOW_0
      inst->smooth[i] *= SMOOTH; // value set to 0.7 in define.h file
      inst->smooth[i] += ((float)1.0 - SMOOTH) * theFilter[i];
#else
      inst->smooth[i] = theFilter[i];
#endif
      real[i] *= inst->smooth[i];
      imag[i] *= inst->smooth[i];
    }
    // keep track of noise and magn spectrum for next frame
    for (i = 0; i < inst->magnLen; i++) {
      inst->noisePrev[i] = noise[i];
      inst->magnPrev[i] = magn[i];
    }
    // back to time domain
    winData[0] = real[0];
    winData[1] = real[inst->magnLen - 1];
    for (i = 1; i < inst->magnLen - 1; i++) {
      winData[2 * i] = real[i];
      winData[2 * i + 1] = imag[i];
    }
    WebRtc_rdft(inst->anaLen, -1, winData, inst->ip, inst->wfft);

    for (i = 0; i < inst->anaLen; i++) {
      real[i] = 2.0f * winData[i] / inst->anaLen; // fft scaling
    }

    //scale factor: only do it after END_STARTUP_LONG time
    factor = (float)1.0;
    if (inst->gainmap == 1 && inst->blockInd > END_STARTUP_LONG) {
      factor1 = (float)1.0;
      factor2 = (float)1.0;

      energy2 = 0.0;
      for (i = 0; i < inst->anaLen; i++) {
        energy2 += (float)real[i] * (float)real[i];
      }
      gain = (float)sqrt(energy2 / (energy1 + (float)1.0));

#ifdef PROCESS_FLOW_2
      // scaling for new version
      if (gain > B_LIM) {
        factor1 = (float)1.0 + (float)1.3 * (gain - B_LIM);
        if (gain * factor1 > (float)1.0) {
          factor1 = (float)1.0 / gain;
        }
      }
      if (gain < B_LIM) {
        //don't reduce scale too much for pause regions:
        // attenuation here should be controlled by flooring
        if (gain <= inst->denoiseBound) {
          gain = inst->denoiseBound;
        }
        factor2 = (float)1.0 - (float)0.3 * (B_LIM - gain);
      }
      //combine both scales with speech/noise prob:
      // note prior (priorSpeechProb) is not frequency dependent
      factor = inst->priorSpeechProb * factor1 + ((float)1.0 - inst->priorSpeechProb)
               * factor2;
#else
      if (gain > B_LIM) {
        factor = (float)1.0 + (float)1.3 * (gain - B_LIM);
      } else {
        factor = (float)1.0 + (float)2.0 * (gain - B_LIM);
      }
      if (gain * factor > (float)1.0) {
        factor = (float)1.0 / gain;
      }
#endif
    } // out of inst->gainmap==1

    // synthesis
    for (i = 0; i < inst->anaLen; i++) {
      inst->syntBuf[i] += factor * inst->window[i] * (float)real[i];
    }
    // read out fully processed segment
    for (i = inst->windShift; i < inst->blockLen + inst->windShift; i++) {
      fout[i - inst->windShift] = inst->syntBuf[i];
    }
    // update synthesis buffer
    memcpy(inst->syntBuf, inst->syntBuf + inst->blockLen,
           sizeof(float) * (inst->anaLen - inst->blockLen));
    memset(inst->syntBuf + inst->anaLen - inst->blockLen, 0,
           sizeof(float) * inst->blockLen);

    // out buffer
    inst->outLen = inst->blockLen - inst->blockLen10ms;
    if (inst->blockLen > inst->blockLen10ms) {
      for (i = 0; i < inst->outLen; i++) {
        inst->outBuf[i] = fout[i + inst->blockLen10ms];
      }
    }
  } // end of if out.len==0
  else {
    for (i = 0; i < inst->blockLen10ms; i++) {
      fout[i] = inst->outBuf[i];
    }
    memcpy(inst->outBuf, inst->outBuf + inst->blockLen10ms,
           sizeof(float) * (inst->outLen - inst->blockLen10ms));
    memset(inst->outBuf + inst->outLen - inst->blockLen10ms, 0,
           sizeof(float) * inst->blockLen10ms);
    inst->outLen -= inst->blockLen10ms;
  }

  // convert to short
  for (i = 0; i < inst->blockLen10ms; i++) {
    dTmp = fout[i];
    if (dTmp < WEBRTC_SPL_WORD16_MIN) {
      dTmp = WEBRTC_SPL_WORD16_MIN;
    } else if (dTmp > WEBRTC_SPL_WORD16_MAX) {
      dTmp = WEBRTC_SPL_WORD16_MAX;
    }
    outFrame[i] = (short)dTmp;
  }

  // for time-domain gain of HB
  if (flagHB == 1) {
    for (i = 0; i < inst->magnLen; i++) {
      inst->speechProbHB[i] = probSpeechFinal[i];
    }
    // average speech prob from low band
    // avg over second half (i.e., 4->8kHz) of freq. spectrum
    avgProbSpeechHB = 0.0;
    for (i = inst->magnLen - deltaBweHB - 1; i < inst->magnLen - 1; i++) {
      avgProbSpeechHB += inst->speechProbHB[i];
    }
    avgProbSpeechHB = avgProbSpeechHB / ((float)deltaBweHB);
    // average filter gain from low band
    // average over second half (i.e., 4->8kHz) of freq. spectrum
    avgFilterGainHB = 0.0;
    for (i = inst->magnLen - deltaGainHB - 1; i < inst->magnLen - 1; i++) {
      avgFilterGainHB += inst->smooth[i];
    }
    avgFilterGainHB = avgFilterGainHB / ((float)(deltaGainHB));
    avgProbSpeechHBTmp = (float)2.0 * avgProbSpeechHB - (float)1.0;
    // gain based on speech prob:
    gainModHB = (float)0.5 * ((float)1.0 + (float)tanh(gainMapParHB * avgProbSpeechHBTmp));
    //combine gain with low band gain
    gainTimeDomainHB = (float)0.5 * gainModHB + (float)0.5 * avgFilterGainHB;
    if (avgProbSpeechHB >= (float)0.5) {
      gainTimeDomainHB = (float)0.25 * gainModHB + (float)0.75 * avgFilterGainHB;
    }
    gainTimeDomainHB = gainTimeDomainHB * decayBweHB;
    //make sure gain is within flooring range
    // flooring bottom
    if (gainTimeDomainHB < inst->denoiseBound) {
      gainTimeDomainHB = inst->denoiseBound;
    }
    // flooring top
    if (gainTimeDomainHB > (float)1.0) {
      gainTimeDomainHB = 1.0;
    }
    //apply gain
    for (i = 0; i < inst->blockLen10ms; i++) {
      dTmp = gainTimeDomainHB * inst->dataBufHB[i];
      if (dTmp < WEBRTC_SPL_WORD16_MIN) {
        dTmp = WEBRTC_SPL_WORD16_MIN;
      } else if (dTmp > WEBRTC_SPL_WORD16_MAX) {
        dTmp = WEBRTC_SPL_WORD16_MAX;
      }
      outFrameHB[i] = (short)dTmp;
    }
  } // end of H band gain computation
  //

  return 0;
}