1/*
2 *  Copyright (c) 2011 The WebRTC project authors. All Rights Reserved.
3 *
4 *  Use of this source code is governed by a BSD-style license
5 *  that can be found in the LICENSE file in the root of the source
6 *  tree. An additional intellectual property rights grant can be found
7 *  in the file PATENTS.  All contributing project authors may
8 *  be found in the AUTHORS file in the root of the source tree.
9 */
10
11#include <string.h>
12#include <math.h>
13//#include <stdio.h>
14#include <stdlib.h>
15#include "noise_suppression.h"
16#include "ns_core.h"
17#include "windows_private.h"
18#include "fft4g.h"
19#include "signal_processing_library.h"
20
21// Set Feature Extraction Parameters
22void WebRtcNs_set_feature_extraction_parameters(NSinst_t* inst) {
23  //bin size of histogram
24  inst->featureExtractionParams.binSizeLrt      = (float)0.1;
25  inst->featureExtractionParams.binSizeSpecFlat = (float)0.05;
26  inst->featureExtractionParams.binSizeSpecDiff = (float)0.1;
27
28  //range of histogram over which lrt threshold is computed
29  inst->featureExtractionParams.rangeAvgHistLrt = (float)1.0;
30
31  //scale parameters: multiply dominant peaks of the histograms by scale factor to obtain
32  // thresholds for prior model
33  inst->featureExtractionParams.factor1ModelPars = (float)1.20; //for lrt and spectral diff
34  inst->featureExtractionParams.factor2ModelPars = (float)0.9;  //for spectral_flatness:
35  // used when noise is flatter than speech
36
37  //peak limit for spectral flatness (varies between 0 and 1)
38  inst->featureExtractionParams.thresPosSpecFlat = (float)0.6;
39
40  //limit on spacing of two highest peaks in histogram: spacing determined by bin size
41  inst->featureExtractionParams.limitPeakSpacingSpecFlat =
42      2 * inst->featureExtractionParams.binSizeSpecFlat;
43  inst->featureExtractionParams.limitPeakSpacingSpecDiff =
44      2 * inst->featureExtractionParams.binSizeSpecDiff;
45
46  //limit on relevance of second peak:
47  inst->featureExtractionParams.limitPeakWeightsSpecFlat = (float)0.5;
48  inst->featureExtractionParams.limitPeakWeightsSpecDiff = (float)0.5;
49
50  // fluctuation limit of lrt feature
51  inst->featureExtractionParams.thresFluctLrt = (float)0.05;
52
53  //limit on the max and min values for the feature thresholds
54  inst->featureExtractionParams.maxLrt = (float)1.0;
55  inst->featureExtractionParams.minLrt = (float)0.20;
56
57  inst->featureExtractionParams.maxSpecFlat = (float)0.95;
58  inst->featureExtractionParams.minSpecFlat = (float)0.10;
59
60  inst->featureExtractionParams.maxSpecDiff = (float)1.0;
61  inst->featureExtractionParams.minSpecDiff = (float)0.16;
62
63  //criteria of weight of histogram peak  to accept/reject feature
64  inst->featureExtractionParams.thresWeightSpecFlat = (int)(0.3
65      * (inst->modelUpdatePars[1])); //for spectral flatness
66  inst->featureExtractionParams.thresWeightSpecDiff = (int)(0.3
67      * (inst->modelUpdatePars[1])); //for spectral difference
68}
69
70// Initialize state
71int WebRtcNs_InitCore(NSinst_t* inst, WebRtc_UWord32 fs) {
72  int i;
73  //We only support 10ms frames
74
75  //check for valid pointer
76  if (inst == NULL) {
77    return -1;
78  }
79
80  // Initialization of struct
81  if (fs == 8000 || fs == 16000 || fs == 32000) {
82    inst->fs = fs;
83  } else {
84    return -1;
85  }
86  inst->windShift = 0;
87  if (fs == 8000) {
88    // We only support 10ms frames
89    inst->blockLen = 80;
90    inst->blockLen10ms = 80;
91    inst->anaLen = 128;
92    inst->window = kBlocks80w128;
93    inst->outLen = 0;
94  } else if (fs == 16000) {
95    // We only support 10ms frames
96    inst->blockLen = 160;
97    inst->blockLen10ms = 160;
98    inst->anaLen = 256;
99    inst->window = kBlocks160w256;
100    inst->outLen = 0;
101  } else if (fs == 32000) {
102    // We only support 10ms frames
103    inst->blockLen = 160;
104    inst->blockLen10ms = 160;
105    inst->anaLen = 256;
106    inst->window = kBlocks160w256;
107    inst->outLen = 0;
108  }
109  inst->magnLen = inst->anaLen / 2 + 1; // Number of frequency bins
110
111  // Initialize fft work arrays.
112  inst->ip[0] = 0; // Setting this triggers initialization.
113  memset(inst->dataBuf, 0, sizeof(float) * ANAL_BLOCKL_MAX);
114  WebRtc_rdft(inst->anaLen, 1, inst->dataBuf, inst->ip, inst->wfft);
115
116  memset(inst->dataBuf, 0, sizeof(float) * ANAL_BLOCKL_MAX);
117  memset(inst->syntBuf, 0, sizeof(float) * ANAL_BLOCKL_MAX);
118
119  //for HB processing
120  memset(inst->dataBufHB, 0, sizeof(float) * ANAL_BLOCKL_MAX);
121
122  //for quantile noise estimation
123  memset(inst->quantile, 0, sizeof(float) * HALF_ANAL_BLOCKL);
124  for (i = 0; i < SIMULT * HALF_ANAL_BLOCKL; i++) {
125    inst->lquantile[i] = (float)8.0;
126    inst->density[i] = (float)0.3;
127  }
128
129  for (i = 0; i < SIMULT; i++) {
130    inst->counter[i] = (int)floor((float)(END_STARTUP_LONG * (i + 1)) / (float)SIMULT);
131  }
132
133  inst->updates = 0;
134
135  // Wiener filter initialization
136  for (i = 0; i < HALF_ANAL_BLOCKL; i++) {
137    inst->smooth[i] = (float)1.0;
138  }
139
140  // Set the aggressiveness: default
141  inst->aggrMode = 0;
142
143  //initialize variables for new method
144  inst->priorSpeechProb = (float)0.5; //prior prob for speech/noise
145  for (i = 0; i < HALF_ANAL_BLOCKL; i++) {
146    inst->magnPrev[i]      = (float)0.0; //previous mag spectrum
147    inst->noisePrev[i]     = (float)0.0; //previous noise-spectrum
148    inst->logLrtTimeAvg[i] = LRT_FEATURE_THR; //smooth LR ratio (same as threshold)
149    inst->magnAvgPause[i]  = (float)0.0; //conservative noise spectrum estimate
150    inst->speechProbHB[i]  = (float)0.0; //for estimation of HB in second pass
151    inst->initMagnEst[i]   = (float)0.0; //initial average mag spectrum
152  }
153
154  //feature quantities
155  inst->featureData[0] = SF_FEATURE_THR;  //spectral flatness (start on threshold)
156  inst->featureData[1] = (float)0.0;      //spectral entropy: not used in this version
157  inst->featureData[2] = (float)0.0;      //spectral variance: not used in this version
158  inst->featureData[3] = LRT_FEATURE_THR; //average lrt factor (start on threshold)
159  inst->featureData[4] = SF_FEATURE_THR;  //spectral template diff (start on threshold)
160  inst->featureData[5] = (float)0.0;      //normalization for spectral-diff
161  inst->featureData[6] = (float)0.0;      //window time-average of input magnitude spectrum
162
163  //histogram quantities: used to estimate/update thresholds for features
164  for (i = 0; i < HIST_PAR_EST; i++) {
165    inst->histLrt[i] = 0;
166    inst->histSpecFlat[i] = 0;
167    inst->histSpecDiff[i] = 0;
168  }
169
170  inst->blockInd = -1; //frame counter
171  inst->priorModelPars[0] = LRT_FEATURE_THR; //default threshold for lrt feature
172  inst->priorModelPars[1] = (float)0.5;      //threshold for spectral flatness:
173  // determined on-line
174  inst->priorModelPars[2] = (float)1.0;      //sgn_map par for spectral measure:
175  // 1 for flatness measure
176  inst->priorModelPars[3] = (float)0.5;      //threshold for template-difference feature:
177  // determined on-line
178  inst->priorModelPars[4] = (float)1.0;      //default weighting parameter for lrt feature
179  inst->priorModelPars[5] = (float)0.0;      //default weighting parameter for
180  // spectral flatness feature
181  inst->priorModelPars[6] = (float)0.0;      //default weighting parameter for
182  // spectral difference feature
183
184  inst->modelUpdatePars[0] = 2;   //update flag for parameters:
185  // 0 no update, 1=update once, 2=update every window
186  inst->modelUpdatePars[1] = 500; //window for update
187  inst->modelUpdatePars[2] = 0;   //counter for update of conservative noise spectrum
188  //counter if the feature thresholds are updated during the sequence
189  inst->modelUpdatePars[3] = inst->modelUpdatePars[1];
190
191  inst->signalEnergy = 0.0;
192  inst->sumMagn = 0.0;
193  inst->whiteNoiseLevel = 0.0;
194  inst->pinkNoiseNumerator = 0.0;
195  inst->pinkNoiseExp = 0.0;
196
197  WebRtcNs_set_feature_extraction_parameters(inst); // Set feature configuration
198
199  //default mode
200  WebRtcNs_set_policy_core(inst, 0);
201
202
203  memset(inst->outBuf, 0, sizeof(float) * 3 * BLOCKL_MAX);
204
205  inst->initFlag = 1;
206  return 0;
207}
208
209int WebRtcNs_set_policy_core(NSinst_t* inst, int mode) {
210  // allow for modes:0,1,2,3
211  if (mode < 0 || mode > 3) {
212    return (-1);
213  }
214
215  inst->aggrMode = mode;
216  if (mode == 0) {
217    inst->overdrive = (float)1.0;
218    inst->denoiseBound = (float)0.5;
219    inst->gainmap = 0;
220  } else if (mode == 1) {
221    //inst->overdrive = (float)1.25;
222    inst->overdrive = (float)1.0;
223    inst->denoiseBound = (float)0.25;
224    inst->gainmap = 1;
225  } else if (mode == 2) {
226    //inst->overdrive = (float)1.25;
227    inst->overdrive = (float)1.1;
228    inst->denoiseBound = (float)0.125;
229    inst->gainmap = 1;
230  } else if (mode == 3) {
231    //inst->overdrive = (float)1.30;
232    inst->overdrive = (float)1.25;
233    inst->denoiseBound = (float)0.09;
234    inst->gainmap = 1;
235  }
236  return 0;
237}
238
239// Estimate noise
240void WebRtcNs_NoiseEstimation(NSinst_t* inst, float* magn, float* noise) {
241  int i, s, offset;
242  float lmagn[HALF_ANAL_BLOCKL], delta;
243
244  if (inst->updates < END_STARTUP_LONG) {
245    inst->updates++;
246  }
247
248  for (i = 0; i < inst->magnLen; i++) {
249    lmagn[i] = (float)log(magn[i]);
250  }
251
252  // loop over simultaneous estimates
253  for (s = 0; s < SIMULT; s++) {
254    offset = s * inst->magnLen;
255
256    // newquantest(...)
257    for (i = 0; i < inst->magnLen; i++) {
258      // compute delta
259      if (inst->density[offset + i] > 1.0) {
260        delta = FACTOR * (float)1.0 / inst->density[offset + i];
261      } else {
262        delta = FACTOR;
263      }
264
265      // update log quantile estimate
266      if (lmagn[i] > inst->lquantile[offset + i]) {
267        inst->lquantile[offset + i] += QUANTILE * delta
268                                       / (float)(inst->counter[s] + 1);
269      } else {
270        inst->lquantile[offset + i] -= ((float)1.0 - QUANTILE) * delta
271                                       / (float)(inst->counter[s] + 1);
272      }
273
274      // update density estimate
275      if (fabs(lmagn[i] - inst->lquantile[offset + i]) < WIDTH) {
276        inst->density[offset + i] = ((float)inst->counter[s] * inst->density[offset
277            + i] + (float)1.0 / ((float)2.0 * WIDTH)) / (float)(inst->counter[s] + 1);
278      }
279    } // end loop over magnitude spectrum
280
281    if (inst->counter[s] >= END_STARTUP_LONG) {
282      inst->counter[s] = 0;
283      if (inst->updates >= END_STARTUP_LONG) {
284        for (i = 0; i < inst->magnLen; i++) {
285          inst->quantile[i] = (float)exp(inst->lquantile[offset + i]);
286        }
287      }
288    }
289
290    inst->counter[s]++;
291  } // end loop over simultaneous estimates
292
293  // Sequentially update the noise during startup
294  if (inst->updates < END_STARTUP_LONG) {
295    // Use the last "s" to get noise during startup that differ from zero.
296    for (i = 0; i < inst->magnLen; i++) {
297      inst->quantile[i] = (float)exp(inst->lquantile[offset + i]);
298    }
299  }
300
301  for (i = 0; i < inst->magnLen; i++) {
302    noise[i] = inst->quantile[i];
303  }
304}
305
306// Extract thresholds for feature parameters
307// histograms are computed over some window_size (given by inst->modelUpdatePars[1])
308// thresholds and weights are extracted every window
309// flag 0 means update histogram only, flag 1 means compute the thresholds/weights
310// threshold and weights are returned in: inst->priorModelPars
311void WebRtcNs_FeatureParameterExtraction(NSinst_t* inst, int flag) {
312  int i, useFeatureSpecFlat, useFeatureSpecDiff, numHistLrt;
313  int maxPeak1, maxPeak2;
314  int weightPeak1SpecFlat, weightPeak2SpecFlat, weightPeak1SpecDiff, weightPeak2SpecDiff;
315
316  float binMid, featureSum;
317  float posPeak1SpecFlat, posPeak2SpecFlat, posPeak1SpecDiff, posPeak2SpecDiff;
318  float fluctLrt, avgHistLrt, avgSquareHistLrt, avgHistLrtCompl;
319
320  //3 features: lrt, flatness, difference
321  //lrt_feature = inst->featureData[3];
322  //flat_feature = inst->featureData[0];
323  //diff_feature = inst->featureData[4];
324
325  //update histograms
326  if (flag == 0) {
327    // LRT
328    if ((inst->featureData[3] < HIST_PAR_EST * inst->featureExtractionParams.binSizeLrt)
329        && (inst->featureData[3] >= 0.0)) {
330      i = (int)(inst->featureData[3] / inst->featureExtractionParams.binSizeLrt);
331      inst->histLrt[i]++;
332    }
333    // Spectral flatness
334    if ((inst->featureData[0] < HIST_PAR_EST
335         * inst->featureExtractionParams.binSizeSpecFlat)
336        && (inst->featureData[0] >= 0.0)) {
337      i = (int)(inst->featureData[0] / inst->featureExtractionParams.binSizeSpecFlat);
338      inst->histSpecFlat[i]++;
339    }
340    // Spectral difference
341    if ((inst->featureData[4] < HIST_PAR_EST
342         * inst->featureExtractionParams.binSizeSpecDiff)
343        && (inst->featureData[4] >= 0.0)) {
344      i = (int)(inst->featureData[4] / inst->featureExtractionParams.binSizeSpecDiff);
345      inst->histSpecDiff[i]++;
346    }
347  }
348
349  // extract parameters for speech/noise probability
350  if (flag == 1) {
351    //lrt feature: compute the average over inst->featureExtractionParams.rangeAvgHistLrt
352    avgHistLrt = 0.0;
353    avgHistLrtCompl = 0.0;
354    avgSquareHistLrt = 0.0;
355    numHistLrt = 0;
356    for (i = 0; i < HIST_PAR_EST; i++) {
357      binMid = ((float)i + (float)0.5) * inst->featureExtractionParams.binSizeLrt;
358      if (binMid <= inst->featureExtractionParams.rangeAvgHistLrt) {
359        avgHistLrt += inst->histLrt[i] * binMid;
360        numHistLrt += inst->histLrt[i];
361      }
362      avgSquareHistLrt += inst->histLrt[i] * binMid * binMid;
363      avgHistLrtCompl += inst->histLrt[i] * binMid;
364    }
365    if (numHistLrt > 0) {
366      avgHistLrt = avgHistLrt / ((float)numHistLrt);
367    }
368    avgHistLrtCompl = avgHistLrtCompl / ((float)inst->modelUpdatePars[1]);
369    avgSquareHistLrt = avgSquareHistLrt / ((float)inst->modelUpdatePars[1]);
370    fluctLrt = avgSquareHistLrt - avgHistLrt * avgHistLrtCompl;
371    // get threshold for lrt feature:
372    if (fluctLrt < inst->featureExtractionParams.thresFluctLrt) {
373      //very low fluct, so likely noise
374      inst->priorModelPars[0] = inst->featureExtractionParams.maxLrt;
375    } else {
376      inst->priorModelPars[0] = inst->featureExtractionParams.factor1ModelPars
377                                * avgHistLrt;
378      // check if value is within min/max range
379      if (inst->priorModelPars[0] < inst->featureExtractionParams.minLrt) {
380        inst->priorModelPars[0] = inst->featureExtractionParams.minLrt;
381      }
382      if (inst->priorModelPars[0] > inst->featureExtractionParams.maxLrt) {
383        inst->priorModelPars[0] = inst->featureExtractionParams.maxLrt;
384      }
385    }
386    // done with lrt feature
387
388    //
389    // for spectral flatness and spectral difference: compute the main peaks of histogram
390    maxPeak1 = 0;
391    maxPeak2 = 0;
392    posPeak1SpecFlat = 0.0;
393    posPeak2SpecFlat = 0.0;
394    weightPeak1SpecFlat = 0;
395    weightPeak2SpecFlat = 0;
396
397    // peaks for flatness
398    for (i = 0; i < HIST_PAR_EST; i++) {
399      binMid = ((float)i + (float)0.5) * inst->featureExtractionParams.binSizeSpecFlat;
400      if (inst->histSpecFlat[i] > maxPeak1) {
401        // Found new "first" peak
402        maxPeak2 = maxPeak1;
403        weightPeak2SpecFlat = weightPeak1SpecFlat;
404        posPeak2SpecFlat = posPeak1SpecFlat;
405
406        maxPeak1 = inst->histSpecFlat[i];
407        weightPeak1SpecFlat = inst->histSpecFlat[i];
408        posPeak1SpecFlat = binMid;
409      } else if (inst->histSpecFlat[i] > maxPeak2) {
410        // Found new "second" peak
411        maxPeak2 = inst->histSpecFlat[i];
412        weightPeak2SpecFlat = inst->histSpecFlat[i];
413        posPeak2SpecFlat = binMid;
414      }
415    }
416
417    //compute two peaks for spectral difference
418    maxPeak1 = 0;
419    maxPeak2 = 0;
420    posPeak1SpecDiff = 0.0;
421    posPeak2SpecDiff = 0.0;
422    weightPeak1SpecDiff = 0;
423    weightPeak2SpecDiff = 0;
424    // peaks for spectral difference
425    for (i = 0; i < HIST_PAR_EST; i++) {
426      binMid = ((float)i + (float)0.5) * inst->featureExtractionParams.binSizeSpecDiff;
427      if (inst->histSpecDiff[i] > maxPeak1) {
428        // Found new "first" peak
429        maxPeak2 = maxPeak1;
430        weightPeak2SpecDiff = weightPeak1SpecDiff;
431        posPeak2SpecDiff = posPeak1SpecDiff;
432
433        maxPeak1 = inst->histSpecDiff[i];
434        weightPeak1SpecDiff = inst->histSpecDiff[i];
435        posPeak1SpecDiff = binMid;
436      } else if (inst->histSpecDiff[i] > maxPeak2) {
437        // Found new "second" peak
438        maxPeak2 = inst->histSpecDiff[i];
439        weightPeak2SpecDiff = inst->histSpecDiff[i];
440        posPeak2SpecDiff = binMid;
441      }
442    }
443
444    // for spectrum flatness feature
445    useFeatureSpecFlat = 1;
446    // merge the two peaks if they are close
447    if ((fabs(posPeak2SpecFlat - posPeak1SpecFlat)
448         < inst->featureExtractionParams.limitPeakSpacingSpecFlat)
449        && (weightPeak2SpecFlat
450            > inst->featureExtractionParams.limitPeakWeightsSpecFlat
451            * weightPeak1SpecFlat)) {
452      weightPeak1SpecFlat += weightPeak2SpecFlat;
453      posPeak1SpecFlat = (float)0.5 * (posPeak1SpecFlat + posPeak2SpecFlat);
454    }
455    //reject if weight of peaks is not large enough, or peak value too small
456    if (weightPeak1SpecFlat < inst->featureExtractionParams.thresWeightSpecFlat
457        || posPeak1SpecFlat < inst->featureExtractionParams.thresPosSpecFlat) {
458      useFeatureSpecFlat = 0;
459    }
460    // if selected, get the threshold
461    if (useFeatureSpecFlat == 1) {
462      // compute the threshold
463      inst->priorModelPars[1] = inst->featureExtractionParams.factor2ModelPars
464                                * posPeak1SpecFlat;
465      //check if value is within min/max range
466      if (inst->priorModelPars[1] < inst->featureExtractionParams.minSpecFlat) {
467        inst->priorModelPars[1] = inst->featureExtractionParams.minSpecFlat;
468      }
469      if (inst->priorModelPars[1] > inst->featureExtractionParams.maxSpecFlat) {
470        inst->priorModelPars[1] = inst->featureExtractionParams.maxSpecFlat;
471      }
472    }
473    // done with flatness feature
474
475    // for template feature
476    useFeatureSpecDiff = 1;
477    // merge the two peaks if they are close
478    if ((fabs(posPeak2SpecDiff - posPeak1SpecDiff)
479         < inst->featureExtractionParams.limitPeakSpacingSpecDiff)
480        && (weightPeak2SpecDiff
481            > inst->featureExtractionParams.limitPeakWeightsSpecDiff
482            * weightPeak1SpecDiff)) {
483      weightPeak1SpecDiff += weightPeak2SpecDiff;
484      posPeak1SpecDiff = (float)0.5 * (posPeak1SpecDiff + posPeak2SpecDiff);
485    }
486    // get the threshold value
487    inst->priorModelPars[3] = inst->featureExtractionParams.factor1ModelPars
488                              * posPeak1SpecDiff;
489    //reject if weight of peaks is not large enough
490    if (weightPeak1SpecDiff < inst->featureExtractionParams.thresWeightSpecDiff) {
491      useFeatureSpecDiff = 0;
492    }
493    //check if value is within min/max range
494    if (inst->priorModelPars[3] < inst->featureExtractionParams.minSpecDiff) {
495      inst->priorModelPars[3] = inst->featureExtractionParams.minSpecDiff;
496    }
497    if (inst->priorModelPars[3] > inst->featureExtractionParams.maxSpecDiff) {
498      inst->priorModelPars[3] = inst->featureExtractionParams.maxSpecDiff;
499    }
500    // done with spectral difference feature
501
502    // don't use template feature if fluctuation of lrt feature is very low:
503    //  most likely just noise state
504    if (fluctLrt < inst->featureExtractionParams.thresFluctLrt) {
505      useFeatureSpecDiff = 0;
506    }
507
508    // select the weights between the features
509    // inst->priorModelPars[4] is weight for lrt: always selected
510    // inst->priorModelPars[5] is weight for spectral flatness
511    // inst->priorModelPars[6] is weight for spectral difference
512    featureSum = (float)(1 + useFeatureSpecFlat + useFeatureSpecDiff);
513    inst->priorModelPars[4] = (float)1.0 / featureSum;
514    inst->priorModelPars[5] = ((float)useFeatureSpecFlat) / featureSum;
515    inst->priorModelPars[6] = ((float)useFeatureSpecDiff) / featureSum;
516
517    // set hists to zero for next update
518    if (inst->modelUpdatePars[0] >= 1) {
519      for (i = 0; i < HIST_PAR_EST; i++) {
520        inst->histLrt[i] = 0;
521        inst->histSpecFlat[i] = 0;
522        inst->histSpecDiff[i] = 0;
523      }
524    }
525  } // end of flag == 1
526}
527
528// Compute spectral flatness on input spectrum
529// magnIn is the magnitude spectrum
530// spectral flatness is returned in inst->featureData[0]
531void WebRtcNs_ComputeSpectralFlatness(NSinst_t* inst, float* magnIn) {
532  int i;
533  int shiftLP = 1; //option to remove first bin(s) from spectral measures
534  float avgSpectralFlatnessNum, avgSpectralFlatnessDen, spectralTmp;
535
536  // comute spectral measures
537  // for flatness
538  avgSpectralFlatnessNum = 0.0;
539  avgSpectralFlatnessDen = inst->sumMagn;
540  for (i = 0; i < shiftLP; i++) {
541    avgSpectralFlatnessDen -= magnIn[i];
542  }
543  // compute log of ratio of the geometric to arithmetic mean: check for log(0) case
544  for (i = shiftLP; i < inst->magnLen; i++) {
545    if (magnIn[i] > 0.0) {
546      avgSpectralFlatnessNum += (float)log(magnIn[i]);
547    } else {
548      inst->featureData[0] -= SPECT_FL_TAVG * inst->featureData[0];
549      return;
550    }
551  }
552  //normalize
553  avgSpectralFlatnessDen = avgSpectralFlatnessDen / inst->magnLen;
554  avgSpectralFlatnessNum = avgSpectralFlatnessNum / inst->magnLen;
555
556  //ratio and inverse log: check for case of log(0)
557  spectralTmp = (float)exp(avgSpectralFlatnessNum) / avgSpectralFlatnessDen;
558
559  //time-avg update of spectral flatness feature
560  inst->featureData[0] += SPECT_FL_TAVG * (spectralTmp - inst->featureData[0]);
561  // done with flatness feature
562}
563
564// Compute the difference measure between input spectrum and a template/learned noise spectrum
565// magnIn is the input spectrum
566// the reference/template spectrum is inst->magnAvgPause[i]
567// returns (normalized) spectral difference in inst->featureData[4]
568void WebRtcNs_ComputeSpectralDifference(NSinst_t* inst, float* magnIn) {
569  // avgDiffNormMagn = var(magnIn) - cov(magnIn, magnAvgPause)^2 / var(magnAvgPause)
570  int i;
571  float avgPause, avgMagn, covMagnPause, varPause, varMagn, avgDiffNormMagn;
572
573  avgPause = 0.0;
574  avgMagn = inst->sumMagn;
575  // compute average quantities
576  for (i = 0; i < inst->magnLen; i++) {
577    //conservative smooth noise spectrum from pause frames
578    avgPause += inst->magnAvgPause[i];
579  }
580  avgPause = avgPause / ((float)inst->magnLen);
581  avgMagn = avgMagn / ((float)inst->magnLen);
582
583  covMagnPause = 0.0;
584  varPause = 0.0;
585  varMagn = 0.0;
586  // compute variance and covariance quantities
587  for (i = 0; i < inst->magnLen; i++) {
588    covMagnPause += (magnIn[i] - avgMagn) * (inst->magnAvgPause[i] - avgPause);
589    varPause += (inst->magnAvgPause[i] - avgPause) * (inst->magnAvgPause[i] - avgPause);
590    varMagn += (magnIn[i] - avgMagn) * (magnIn[i] - avgMagn);
591  }
592  covMagnPause = covMagnPause / ((float)inst->magnLen);
593  varPause = varPause / ((float)inst->magnLen);
594  varMagn = varMagn / ((float)inst->magnLen);
595  // update of average magnitude spectrum
596  inst->featureData[6] += inst->signalEnergy;
597
598  avgDiffNormMagn = varMagn - (covMagnPause * covMagnPause) / (varPause + (float)0.0001);
599  // normalize and compute time-avg update of difference feature
600  avgDiffNormMagn = (float)(avgDiffNormMagn / (inst->featureData[5] + (float)0.0001));
601  inst->featureData[4] += SPECT_DIFF_TAVG * (avgDiffNormMagn - inst->featureData[4]);
602}
603
604// Compute speech/noise probability
605// speech/noise probability is returned in: probSpeechFinal
606//magn is the input magnitude spectrum
607//noise is the noise spectrum
608//snrLocPrior is the prior snr for each freq.
609//snr loc_post is the post snr for each freq.
610void WebRtcNs_SpeechNoiseProb(NSinst_t* inst, float* probSpeechFinal, float* snrLocPrior,
611                              float* snrLocPost) {
612  int i, sgnMap;
613  float invLrt, gainPrior, indPrior;
614  float logLrtTimeAvgKsum, besselTmp;
615  float indicator0, indicator1, indicator2;
616  float tmpFloat1, tmpFloat2;
617  float weightIndPrior0, weightIndPrior1, weightIndPrior2;
618  float threshPrior0, threshPrior1, threshPrior2;
619  float widthPrior, widthPrior0, widthPrior1, widthPrior2;
620
621  widthPrior0 = WIDTH_PR_MAP;
622  widthPrior1 = (float)2.0 * WIDTH_PR_MAP; //width for pause region:
623  // lower range, so increase width in tanh map
624  widthPrior2 = (float)2.0 * WIDTH_PR_MAP; //for spectral-difference measure
625
626  //threshold parameters for features
627  threshPrior0 = inst->priorModelPars[0];
628  threshPrior1 = inst->priorModelPars[1];
629  threshPrior2 = inst->priorModelPars[3];
630
631  //sign for flatness feature
632  sgnMap = (int)(inst->priorModelPars[2]);
633
634  //weight parameters for features
635  weightIndPrior0 = inst->priorModelPars[4];
636  weightIndPrior1 = inst->priorModelPars[5];
637  weightIndPrior2 = inst->priorModelPars[6];
638
639  // compute feature based on average LR factor
640  // this is the average over all frequencies of the smooth log lrt
641  logLrtTimeAvgKsum = 0.0;
642  for (i = 0; i < inst->magnLen; i++) {
643    tmpFloat1 = (float)1.0 + (float)2.0 * snrLocPrior[i];
644    tmpFloat2 = (float)2.0 * snrLocPrior[i] / (tmpFloat1 + (float)0.0001);
645    besselTmp = (snrLocPost[i] + (float)1.0) * tmpFloat2;
646    inst->logLrtTimeAvg[i] += LRT_TAVG * (besselTmp - (float)log(tmpFloat1)
647                                          - inst->logLrtTimeAvg[i]);
648    logLrtTimeAvgKsum += inst->logLrtTimeAvg[i];
649  }
650  logLrtTimeAvgKsum = (float)logLrtTimeAvgKsum / (inst->magnLen);
651  inst->featureData[3] = logLrtTimeAvgKsum;
652  // done with computation of LR factor
653
654  //
655  //compute the indicator functions
656  //
657
658  // average lrt feature
659  widthPrior = widthPrior0;
660  //use larger width in tanh map for pause regions
661  if (logLrtTimeAvgKsum < threshPrior0) {
662    widthPrior = widthPrior1;
663  }
664  // compute indicator function: sigmoid map
665  indicator0 = (float)0.5 * ((float)tanh(widthPrior *
666      (logLrtTimeAvgKsum - threshPrior0)) + (float)1.0);
667
668  //spectral flatness feature
669  tmpFloat1 = inst->featureData[0];
670  widthPrior = widthPrior0;
671  //use larger width in tanh map for pause regions
672  if (sgnMap == 1 && (tmpFloat1 > threshPrior1)) {
673    widthPrior = widthPrior1;
674  }
675  if (sgnMap == -1 && (tmpFloat1 < threshPrior1)) {
676    widthPrior = widthPrior1;
677  }
678  // compute indicator function: sigmoid map
679  indicator1 = (float)0.5 * ((float)tanh((float)sgnMap *
680      widthPrior * (threshPrior1 - tmpFloat1)) + (float)1.0);
681
682  //for template spectrum-difference
683  tmpFloat1 = inst->featureData[4];
684  widthPrior = widthPrior0;
685  //use larger width in tanh map for pause regions
686  if (tmpFloat1 < threshPrior2) {
687    widthPrior = widthPrior2;
688  }
689  // compute indicator function: sigmoid map
690  indicator2 = (float)0.5 * ((float)tanh(widthPrior * (tmpFloat1 - threshPrior2))
691                             + (float)1.0);
692
693  //combine the indicator function with the feature weights
694  indPrior = weightIndPrior0 * indicator0 + weightIndPrior1 * indicator1 + weightIndPrior2
695             * indicator2;
696  // done with computing indicator function
697
698  //compute the prior probability
699  inst->priorSpeechProb += PRIOR_UPDATE * (indPrior - inst->priorSpeechProb);
700  // make sure probabilities are within range: keep floor to 0.01
701  if (inst->priorSpeechProb > 1.0) {
702    inst->priorSpeechProb = (float)1.0;
703  }
704  if (inst->priorSpeechProb < 0.01) {
705    inst->priorSpeechProb = (float)0.01;
706  }
707
708  //final speech probability: combine prior model with LR factor:
709  gainPrior = ((float)1.0 - inst->priorSpeechProb) / (inst->priorSpeechProb + (float)0.0001);
710  for (i = 0; i < inst->magnLen; i++) {
711    invLrt = (float)exp(-inst->logLrtTimeAvg[i]);
712    invLrt = (float)gainPrior * invLrt;
713    probSpeechFinal[i] = (float)1.0 / ((float)1.0 + invLrt);
714  }
715}
716
717int WebRtcNs_ProcessCore(NSinst_t* inst,
718                         short* speechFrame,
719                         short* speechFrameHB,
720                         short* outFrame,
721                         short* outFrameHB) {
722  // main routine for noise reduction
723
724  int     flagHB = 0;
725  int     i;
726  const int kStartBand = 5; // Skip first frequency bins during estimation.
727  int     updateParsFlag;
728
729  float   energy1, energy2, gain, factor, factor1, factor2;
730  float   signalEnergy, sumMagn;
731  float   snrPrior, currentEstimateStsa;
732  float   tmpFloat1, tmpFloat2, tmpFloat3, probSpeech, probNonSpeech;
733  float   gammaNoiseTmp, gammaNoiseOld;
734  float   noiseUpdateTmp, fTmp, dTmp;
735  float   fin[BLOCKL_MAX], fout[BLOCKL_MAX];
736  float   winData[ANAL_BLOCKL_MAX];
737  float   magn[HALF_ANAL_BLOCKL], noise[HALF_ANAL_BLOCKL];
738  float   theFilter[HALF_ANAL_BLOCKL], theFilterTmp[HALF_ANAL_BLOCKL];
739  float   snrLocPost[HALF_ANAL_BLOCKL], snrLocPrior[HALF_ANAL_BLOCKL];
740  float   probSpeechFinal[HALF_ANAL_BLOCKL], previousEstimateStsa[HALF_ANAL_BLOCKL];
741  float   real[ANAL_BLOCKL_MAX], imag[HALF_ANAL_BLOCKL];
742  // Variables during startup
743  float   sum_log_i = 0.0;
744  float   sum_log_i_square = 0.0;
745  float   sum_log_magn = 0.0;
746  float   sum_log_i_log_magn = 0.0;
747  float   parametric_noise = 0.0;
748  float   parametric_exp = 0.0;
749  float   parametric_num = 0.0;
750
751  // SWB variables
752  int     deltaBweHB = 1;
753  int     deltaGainHB = 1;
754  float   decayBweHB = 1.0;
755  float   gainMapParHB = 1.0;
756  float   gainTimeDomainHB = 1.0;
757  float   avgProbSpeechHB, avgProbSpeechHBTmp, avgFilterGainHB, gainModHB;
758
759  // Check that initiation has been done
760  if (inst->initFlag != 1) {
761    return (-1);
762  }
763  // Check for valid pointers based on sampling rate
764  if (inst->fs == 32000) {
765    if (speechFrameHB == NULL) {
766      return -1;
767    }
768    flagHB = 1;
769    // range for averaging low band quantities for H band gain
770    deltaBweHB = (int)inst->magnLen / 4;
771    deltaGainHB = deltaBweHB;
772  }
773  //
774  updateParsFlag = inst->modelUpdatePars[0];
775  //
776
777  //for LB do all processing
778  // convert to float
779  for (i = 0; i < inst->blockLen10ms; i++) {
780    fin[i] = (float)speechFrame[i];
781  }
782  // update analysis buffer for L band
783  memcpy(inst->dataBuf, inst->dataBuf + inst->blockLen10ms,
784         sizeof(float) * (inst->anaLen - inst->blockLen10ms));
785  memcpy(inst->dataBuf + inst->anaLen - inst->blockLen10ms, fin,
786         sizeof(float) * inst->blockLen10ms);
787
788  if (flagHB == 1) {
789    // convert to float
790    for (i = 0; i < inst->blockLen10ms; i++) {
791      fin[i] = (float)speechFrameHB[i];
792    }
793    // update analysis buffer for H band
794    memcpy(inst->dataBufHB, inst->dataBufHB + inst->blockLen10ms,
795           sizeof(float) * (inst->anaLen - inst->blockLen10ms));
796    memcpy(inst->dataBufHB + inst->anaLen - inst->blockLen10ms, fin,
797           sizeof(float) * inst->blockLen10ms);
798  }
799
800  // check if processing needed
801  if (inst->outLen == 0) {
802    // windowing
803    energy1 = 0.0;
804    for (i = 0; i < inst->anaLen; i++) {
805      winData[i] = inst->window[i] * inst->dataBuf[i];
806      energy1 += winData[i] * winData[i];
807    }
808    if (energy1 == 0.0) {
809      // synthesize the special case of zero input
810      // we want to avoid updating statistics in this case:
811      // Updating feature statistics when we have zeros only will cause thresholds to
812      // move towards zero signal situations. This in turn has the effect that once the
813      // signal is "turned on" (non-zero values) everything will be treated as speech
814      // and there is no noise suppression effect. Depending on the duration of the
815      // inactive signal it takes a considerable amount of time for the system to learn
816      // what is noise and what is speech.
817
818      // read out fully processed segment
819      for (i = inst->windShift; i < inst->blockLen + inst->windShift; i++) {
820        fout[i - inst->windShift] = inst->syntBuf[i];
821      }
822      // update synthesis buffer
823      memcpy(inst->syntBuf, inst->syntBuf + inst->blockLen,
824             sizeof(float) * (inst->anaLen - inst->blockLen));
825      memset(inst->syntBuf + inst->anaLen - inst->blockLen, 0,
826             sizeof(float) * inst->blockLen);
827
828      // out buffer
829      inst->outLen = inst->blockLen - inst->blockLen10ms;
830      if (inst->blockLen > inst->blockLen10ms) {
831        for (i = 0; i < inst->outLen; i++) {
832          inst->outBuf[i] = fout[i + inst->blockLen10ms];
833        }
834      }
835      // convert to short
836      for (i = 0; i < inst->blockLen10ms; i++) {
837        dTmp = fout[i];
838        if (dTmp < WEBRTC_SPL_WORD16_MIN) {
839          dTmp = WEBRTC_SPL_WORD16_MIN;
840        } else if (dTmp > WEBRTC_SPL_WORD16_MAX) {
841          dTmp = WEBRTC_SPL_WORD16_MAX;
842        }
843        outFrame[i] = (short)dTmp;
844      }
845
846      // for time-domain gain of HB
847      if (flagHB == 1) {
848        for (i = 0; i < inst->blockLen10ms; i++) {
849          dTmp = inst->dataBufHB[i];
850          if (dTmp < WEBRTC_SPL_WORD16_MIN) {
851            dTmp = WEBRTC_SPL_WORD16_MIN;
852          } else if (dTmp > WEBRTC_SPL_WORD16_MAX) {
853            dTmp = WEBRTC_SPL_WORD16_MAX;
854          }
855          outFrameHB[i] = (short)dTmp;
856        }
857      } // end of H band gain computation
858      //
859      return 0;
860    }
861
862    //
863    inst->blockInd++; // Update the block index only when we process a block.
864    // FFT
865    WebRtc_rdft(inst->anaLen, 1, winData, inst->ip, inst->wfft);
866
867    imag[0] = 0;
868    real[0] = winData[0];
869    magn[0] = (float)(fabs(real[0]) + 1.0f);
870    imag[inst->magnLen - 1] = 0;
871    real[inst->magnLen - 1] = winData[1];
872    magn[inst->magnLen - 1] = (float)(fabs(real[inst->magnLen - 1]) + 1.0f);
873    signalEnergy = (float)(real[0] * real[0]) +
874                   (float)(real[inst->magnLen - 1] * real[inst->magnLen - 1]);
875    sumMagn = magn[0] + magn[inst->magnLen - 1];
876    if (inst->blockInd < END_STARTUP_SHORT) {
877      inst->initMagnEst[0] += magn[0];
878      inst->initMagnEst[inst->magnLen - 1] += magn[inst->magnLen - 1];
879      tmpFloat2 = log((float)(inst->magnLen - 1));
880      sum_log_i = tmpFloat2;
881      sum_log_i_square = tmpFloat2 * tmpFloat2;
882      tmpFloat1 = log(magn[inst->magnLen - 1]);
883      sum_log_magn = tmpFloat1;
884      sum_log_i_log_magn = tmpFloat2 * tmpFloat1;
885    }
886    for (i = 1; i < inst->magnLen - 1; i++) {
887      real[i] = winData[2 * i];
888      imag[i] = winData[2 * i + 1];
889      // magnitude spectrum
890      fTmp = real[i] * real[i];
891      fTmp += imag[i] * imag[i];
892      signalEnergy += fTmp;
893      magn[i] = ((float)sqrt(fTmp)) + 1.0f;
894      sumMagn += magn[i];
895      if (inst->blockInd < END_STARTUP_SHORT) {
896        inst->initMagnEst[i] += magn[i];
897        if (i >= kStartBand) {
898          tmpFloat2 = log((float)i);
899          sum_log_i += tmpFloat2;
900          sum_log_i_square += tmpFloat2 * tmpFloat2;
901          tmpFloat1 = log(magn[i]);
902          sum_log_magn += tmpFloat1;
903          sum_log_i_log_magn += tmpFloat2 * tmpFloat1;
904        }
905      }
906    }
907    signalEnergy = signalEnergy / ((float)inst->magnLen);
908    inst->signalEnergy = signalEnergy;
909    inst->sumMagn = sumMagn;
910
911    //compute spectral flatness on input spectrum
912    WebRtcNs_ComputeSpectralFlatness(inst, magn);
913    // quantile noise estimate
914    WebRtcNs_NoiseEstimation(inst, magn, noise);
915    //compute simplified noise model during startup
916    if (inst->blockInd < END_STARTUP_SHORT) {
917      // Estimate White noise
918      inst->whiteNoiseLevel += sumMagn / ((float)inst->magnLen) * inst->overdrive;
919      // Estimate Pink noise parameters
920      tmpFloat1 = sum_log_i_square * ((float)(inst->magnLen - kStartBand));
921      tmpFloat1 -= (sum_log_i * sum_log_i);
922      tmpFloat2 = (sum_log_i_square * sum_log_magn - sum_log_i * sum_log_i_log_magn);
923      tmpFloat3 = tmpFloat2 / tmpFloat1;
924      // Constrain the estimated spectrum to be positive
925      if (tmpFloat3 < 0.0f) {
926        tmpFloat3 = 0.0f;
927      }
928      inst->pinkNoiseNumerator += tmpFloat3;
929      tmpFloat2 = (sum_log_i * sum_log_magn);
930      tmpFloat2 -= ((float)(inst->magnLen - kStartBand)) * sum_log_i_log_magn;
931      tmpFloat3 = tmpFloat2 / tmpFloat1;
932      // Constrain the pink noise power to be in the interval [0, 1];
933      if (tmpFloat3 < 0.0f) {
934        tmpFloat3 = 0.0f;
935      }
936      if (tmpFloat3 > 1.0f) {
937        tmpFloat3 = 1.0f;
938      }
939      inst->pinkNoiseExp += tmpFloat3;
940
941      // Calculate frequency independent parts of parametric noise estimate.
942      if (inst->pinkNoiseExp == 0.0f) {
943        // Use white noise estimate
944        parametric_noise = inst->whiteNoiseLevel;
945      } else {
946        // Use pink noise estimate
947        parametric_num = exp(inst->pinkNoiseNumerator / (float)(inst->blockInd + 1));
948        parametric_num *= (float)(inst->blockInd + 1);
949        parametric_exp = inst->pinkNoiseExp / (float)(inst->blockInd + 1);
950        parametric_noise = parametric_num / pow((float)kStartBand, parametric_exp);
951      }
952      for (i = 0; i < inst->magnLen; i++) {
953        // Estimate the background noise using the white and pink noise parameters
954        if ((inst->pinkNoiseExp > 0.0f) && (i >= kStartBand)) {
955          // Use pink noise estimate
956          parametric_noise = parametric_num / pow((float)i, parametric_exp);
957        }
958        theFilterTmp[i] = (inst->initMagnEst[i] - inst->overdrive * parametric_noise);
959        theFilterTmp[i] /= (inst->initMagnEst[i] + (float)0.0001);
960        // Weight quantile noise with modeled noise
961        noise[i] *= (inst->blockInd);
962        tmpFloat2 = parametric_noise * (END_STARTUP_SHORT - inst->blockInd);
963        noise[i] += (tmpFloat2 / (float)(inst->blockInd + 1));
964        noise[i] /= END_STARTUP_SHORT;
965      }
966    }
967    //compute average signal during END_STARTUP_LONG time:
968    // used to normalize spectral difference measure
969    if (inst->blockInd < END_STARTUP_LONG) {
970      inst->featureData[5] *= inst->blockInd;
971      inst->featureData[5] += signalEnergy;
972      inst->featureData[5] /= (inst->blockInd + 1);
973    }
974
975#ifdef PROCESS_FLOW_0
976    if (inst->blockInd > END_STARTUP_LONG) {
977      //option: average the quantile noise: for check with AEC2
978      for (i = 0; i < inst->magnLen; i++) {
979        noise[i] = (float)0.6 * inst->noisePrev[i] + (float)0.4 * noise[i];
980      }
981      for (i = 0; i < inst->magnLen; i++) {
982        // Wiener with over sub-substraction:
983        theFilter[i] = (magn[i] - inst->overdrive * noise[i]) / (magn[i] + (float)0.0001);
984      }
985    }
986#else
987    //start processing at frames == converged+1
988    //
989    // STEP 1: compute  prior and post snr based on quantile noise est
990    //
991
992    // compute DD estimate of prior SNR: needed for new method
993    for (i = 0; i < inst->magnLen; i++) {
994      // post snr
995      snrLocPost[i] = (float)0.0;
996      if (magn[i] > noise[i]) {
997        snrLocPost[i] = magn[i] / (noise[i] + (float)0.0001) - (float)1.0;
998      }
999      // previous post snr
1000      // previous estimate: based on previous frame with gain filter
1001      previousEstimateStsa[i] = inst->magnPrev[i] / (inst->noisePrev[i] + (float)0.0001)
1002                                * (inst->smooth[i]);
1003      // DD estimate is sum of two terms: current estimate and previous estimate
1004      // directed decision update of snrPrior
1005      snrLocPrior[i] = DD_PR_SNR * previousEstimateStsa[i] + ((float)1.0 - DD_PR_SNR)
1006                       * snrLocPost[i];
1007      // post and prior snr needed for step 2
1008    } // end of loop over freqs
1009#ifdef PROCESS_FLOW_1
1010    for (i = 0; i < inst->magnLen; i++) {
1011      // gain filter
1012      tmpFloat1 = inst->overdrive + snrLocPrior[i];
1013      tmpFloat2 = (float)snrLocPrior[i] / tmpFloat1;
1014      theFilter[i] = (float)tmpFloat2;
1015    } // end of loop over freqs
1016#endif
1017    // done with step 1: dd computation of prior and post snr
1018
1019    //
1020    //STEP 2: compute speech/noise likelihood
1021    //
1022#ifdef PROCESS_FLOW_2
1023    // compute difference of input spectrum with learned/estimated noise spectrum
1024    WebRtcNs_ComputeSpectralDifference(inst, magn);
1025    // compute histograms for parameter decisions (thresholds and weights for features)
1026    // parameters are extracted once every window time (=inst->modelUpdatePars[1])
1027    if (updateParsFlag >= 1) {
1028      // counter update
1029      inst->modelUpdatePars[3]--;
1030      // update histogram
1031      if (inst->modelUpdatePars[3] > 0) {
1032        WebRtcNs_FeatureParameterExtraction(inst, 0);
1033      }
1034      // compute model parameters
1035      if (inst->modelUpdatePars[3] == 0) {
1036        WebRtcNs_FeatureParameterExtraction(inst, 1);
1037        inst->modelUpdatePars[3] = inst->modelUpdatePars[1];
1038        // if wish to update only once, set flag to zero
1039        if (updateParsFlag == 1) {
1040          inst->modelUpdatePars[0] = 0;
1041        } else {
1042          // update every window:
1043          // get normalization for spectral difference for next window estimate
1044          inst->featureData[6] = inst->featureData[6]
1045                                 / ((float)inst->modelUpdatePars[1]);
1046          inst->featureData[5] = (float)0.5 * (inst->featureData[6]
1047                                               + inst->featureData[5]);
1048          inst->featureData[6] = (float)0.0;
1049        }
1050      }
1051    }
1052    // compute speech/noise probability
1053    WebRtcNs_SpeechNoiseProb(inst, probSpeechFinal, snrLocPrior, snrLocPost);
1054    // time-avg parameter for noise update
1055    gammaNoiseTmp = NOISE_UPDATE;
1056    for (i = 0; i < inst->magnLen; i++) {
1057      probSpeech = probSpeechFinal[i];
1058      probNonSpeech = (float)1.0 - probSpeech;
1059      // temporary noise update:
1060      // use it for speech frames if update value is less than previous
1061      noiseUpdateTmp = gammaNoiseTmp * inst->noisePrev[i] + ((float)1.0 - gammaNoiseTmp)
1062                       * (probNonSpeech * magn[i] + probSpeech * inst->noisePrev[i]);
1063      //
1064      // time-constant based on speech/noise state
1065      gammaNoiseOld = gammaNoiseTmp;
1066      gammaNoiseTmp = NOISE_UPDATE;
1067      // increase gamma (i.e., less noise update) for frame likely to be speech
1068      if (probSpeech > PROB_RANGE) {
1069        gammaNoiseTmp = SPEECH_UPDATE;
1070      }
1071      // conservative noise update
1072      if (probSpeech < PROB_RANGE) {
1073        inst->magnAvgPause[i] += GAMMA_PAUSE * (magn[i] - inst->magnAvgPause[i]);
1074      }
1075      // noise update
1076      if (gammaNoiseTmp == gammaNoiseOld) {
1077        noise[i] = noiseUpdateTmp;
1078      } else {
1079        noise[i] = gammaNoiseTmp * inst->noisePrev[i] + ((float)1.0 - gammaNoiseTmp)
1080                   * (probNonSpeech * magn[i] + probSpeech * inst->noisePrev[i]);
1081        // allow for noise update downwards:
1082        //  if noise update decreases the noise, it is safe, so allow it to happen
1083        if (noiseUpdateTmp < noise[i]) {
1084          noise[i] = noiseUpdateTmp;
1085        }
1086      }
1087    } // end of freq loop
1088    // done with step 2: noise update
1089
1090    //
1091    // STEP 3: compute dd update of prior snr and post snr based on new noise estimate
1092    //
1093    for (i = 0; i < inst->magnLen; i++) {
1094      // post and prior snr
1095      currentEstimateStsa = (float)0.0;
1096      if (magn[i] > noise[i]) {
1097        currentEstimateStsa = magn[i] / (noise[i] + (float)0.0001) - (float)1.0;
1098      }
1099      // DD estimate is sume of two terms: current estimate and previous estimate
1100      // directed decision update of snrPrior
1101      snrPrior = DD_PR_SNR * previousEstimateStsa[i] + ((float)1.0 - DD_PR_SNR)
1102                 * currentEstimateStsa;
1103      // gain filter
1104      tmpFloat1 = inst->overdrive + snrPrior;
1105      tmpFloat2 = (float)snrPrior / tmpFloat1;
1106      theFilter[i] = (float)tmpFloat2;
1107    } // end of loop over freqs
1108    // done with step3
1109#endif
1110#endif
1111
1112    for (i = 0; i < inst->magnLen; i++) {
1113      // flooring bottom
1114      if (theFilter[i] < inst->denoiseBound) {
1115        theFilter[i] = inst->denoiseBound;
1116      }
1117      // flooring top
1118      if (theFilter[i] > (float)1.0) {
1119        theFilter[i] = 1.0;
1120      }
1121      if (inst->blockInd < END_STARTUP_SHORT) {
1122        // flooring bottom
1123        if (theFilterTmp[i] < inst->denoiseBound) {
1124          theFilterTmp[i] = inst->denoiseBound;
1125        }
1126        // flooring top
1127        if (theFilterTmp[i] > (float)1.0) {
1128          theFilterTmp[i] = 1.0;
1129        }
1130        // Weight the two suppression filters
1131        theFilter[i] *= (inst->blockInd);
1132        theFilterTmp[i] *= (END_STARTUP_SHORT - inst->blockInd);
1133        theFilter[i] += theFilterTmp[i];
1134        theFilter[i] /= (END_STARTUP_SHORT);
1135      }
1136      // smoothing
1137#ifdef PROCESS_FLOW_0
1138      inst->smooth[i] *= SMOOTH; // value set to 0.7 in define.h file
1139      inst->smooth[i] += ((float)1.0 - SMOOTH) * theFilter[i];
1140#else
1141      inst->smooth[i] = theFilter[i];
1142#endif
1143      real[i] *= inst->smooth[i];
1144      imag[i] *= inst->smooth[i];
1145    }
1146    // keep track of noise and magn spectrum for next frame
1147    for (i = 0; i < inst->magnLen; i++) {
1148      inst->noisePrev[i] = noise[i];
1149      inst->magnPrev[i] = magn[i];
1150    }
1151    // back to time domain
1152    winData[0] = real[0];
1153    winData[1] = real[inst->magnLen - 1];
1154    for (i = 1; i < inst->magnLen - 1; i++) {
1155      winData[2 * i] = real[i];
1156      winData[2 * i + 1] = imag[i];
1157    }
1158    WebRtc_rdft(inst->anaLen, -1, winData, inst->ip, inst->wfft);
1159
1160    for (i = 0; i < inst->anaLen; i++) {
1161      real[i] = 2.0f * winData[i] / inst->anaLen; // fft scaling
1162    }
1163
1164    //scale factor: only do it after END_STARTUP_LONG time
1165    factor = (float)1.0;
1166    if (inst->gainmap == 1 && inst->blockInd > END_STARTUP_LONG) {
1167      factor1 = (float)1.0;
1168      factor2 = (float)1.0;
1169
1170      energy2 = 0.0;
1171      for (i = 0; i < inst->anaLen; i++) {
1172        energy2 += (float)real[i] * (float)real[i];
1173      }
1174      gain = (float)sqrt(energy2 / (energy1 + (float)1.0));
1175
1176#ifdef PROCESS_FLOW_2
1177      // scaling for new version
1178      if (gain > B_LIM) {
1179        factor1 = (float)1.0 + (float)1.3 * (gain - B_LIM);
1180        if (gain * factor1 > (float)1.0) {
1181          factor1 = (float)1.0 / gain;
1182        }
1183      }
1184      if (gain < B_LIM) {
1185        //don't reduce scale too much for pause regions:
1186        // attenuation here should be controlled by flooring
1187        if (gain <= inst->denoiseBound) {
1188          gain = inst->denoiseBound;
1189        }
1190        factor2 = (float)1.0 - (float)0.3 * (B_LIM - gain);
1191      }
1192      //combine both scales with speech/noise prob:
1193      // note prior (priorSpeechProb) is not frequency dependent
1194      factor = inst->priorSpeechProb * factor1 + ((float)1.0 - inst->priorSpeechProb)
1195               * factor2;
1196#else
1197      if (gain > B_LIM) {
1198        factor = (float)1.0 + (float)1.3 * (gain - B_LIM);
1199      } else {
1200        factor = (float)1.0 + (float)2.0 * (gain - B_LIM);
1201      }
1202      if (gain * factor > (float)1.0) {
1203        factor = (float)1.0 / gain;
1204      }
1205#endif
1206    } // out of inst->gainmap==1
1207
1208    // synthesis
1209    for (i = 0; i < inst->anaLen; i++) {
1210      inst->syntBuf[i] += factor * inst->window[i] * (float)real[i];
1211    }
1212    // read out fully processed segment
1213    for (i = inst->windShift; i < inst->blockLen + inst->windShift; i++) {
1214      fout[i - inst->windShift] = inst->syntBuf[i];
1215    }
1216    // update synthesis buffer
1217    memcpy(inst->syntBuf, inst->syntBuf + inst->blockLen,
1218           sizeof(float) * (inst->anaLen - inst->blockLen));
1219    memset(inst->syntBuf + inst->anaLen - inst->blockLen, 0,
1220           sizeof(float) * inst->blockLen);
1221
1222    // out buffer
1223    inst->outLen = inst->blockLen - inst->blockLen10ms;
1224    if (inst->blockLen > inst->blockLen10ms) {
1225      for (i = 0; i < inst->outLen; i++) {
1226        inst->outBuf[i] = fout[i + inst->blockLen10ms];
1227      }
1228    }
1229  } // end of if out.len==0
1230  else {
1231    for (i = 0; i < inst->blockLen10ms; i++) {
1232      fout[i] = inst->outBuf[i];
1233    }
1234    memcpy(inst->outBuf, inst->outBuf + inst->blockLen10ms,
1235           sizeof(float) * (inst->outLen - inst->blockLen10ms));
1236    memset(inst->outBuf + inst->outLen - inst->blockLen10ms, 0,
1237           sizeof(float) * inst->blockLen10ms);
1238    inst->outLen -= inst->blockLen10ms;
1239  }
1240
1241  // convert to short
1242  for (i = 0; i < inst->blockLen10ms; i++) {
1243    dTmp = fout[i];
1244    if (dTmp < WEBRTC_SPL_WORD16_MIN) {
1245      dTmp = WEBRTC_SPL_WORD16_MIN;
1246    } else if (dTmp > WEBRTC_SPL_WORD16_MAX) {
1247      dTmp = WEBRTC_SPL_WORD16_MAX;
1248    }
1249    outFrame[i] = (short)dTmp;
1250  }
1251
1252  // for time-domain gain of HB
1253  if (flagHB == 1) {
1254    for (i = 0; i < inst->magnLen; i++) {
1255      inst->speechProbHB[i] = probSpeechFinal[i];
1256    }
1257    if (inst->blockInd > END_STARTUP_LONG) {
1258      // average speech prob from low band
1259      // avg over second half (i.e., 4->8kHz) of freq. spectrum
1260      avgProbSpeechHB = 0.0;
1261      for (i = inst->magnLen - deltaBweHB - 1; i < inst->magnLen - 1; i++) {
1262        avgProbSpeechHB += inst->speechProbHB[i];
1263      }
1264      avgProbSpeechHB = avgProbSpeechHB / ((float)deltaBweHB);
1265      // average filter gain from low band
1266      // average over second half (i.e., 4->8kHz) of freq. spectrum
1267      avgFilterGainHB = 0.0;
1268      for (i = inst->magnLen - deltaGainHB - 1; i < inst->magnLen - 1; i++) {
1269        avgFilterGainHB += inst->smooth[i];
1270      }
1271      avgFilterGainHB = avgFilterGainHB / ((float)(deltaGainHB));
1272      avgProbSpeechHBTmp = (float)2.0 * avgProbSpeechHB - (float)1.0;
1273      // gain based on speech prob:
1274      gainModHB = (float)0.5 * ((float)1.0 + (float)tanh(gainMapParHB * avgProbSpeechHBTmp));
1275      //combine gain with low band gain
1276      gainTimeDomainHB = (float)0.5 * gainModHB + (float)0.5 * avgFilterGainHB;
1277      if (avgProbSpeechHB >= (float)0.5) {
1278        gainTimeDomainHB = (float)0.25 * gainModHB + (float)0.75 * avgFilterGainHB;
1279      }
1280      gainTimeDomainHB = gainTimeDomainHB * decayBweHB;
1281    } // end of converged
1282    //make sure gain is within flooring range
1283    // flooring bottom
1284    if (gainTimeDomainHB < inst->denoiseBound) {
1285      gainTimeDomainHB = inst->denoiseBound;
1286    }
1287    // flooring top
1288    if (gainTimeDomainHB > (float)1.0) {
1289      gainTimeDomainHB = 1.0;
1290    }
1291    //apply gain
1292    for (i = 0; i < inst->blockLen10ms; i++) {
1293      dTmp = gainTimeDomainHB * inst->dataBufHB[i];
1294      if (dTmp < WEBRTC_SPL_WORD16_MIN) {
1295        dTmp = WEBRTC_SPL_WORD16_MIN;
1296      } else if (dTmp > WEBRTC_SPL_WORD16_MAX) {
1297        dTmp = WEBRTC_SPL_WORD16_MAX;
1298      }
1299      outFrameHB[i] = (short)dTmp;
1300    }
1301  } // end of H band gain computation
1302  //
1303
1304  return 0;
1305}
1306