Bug Summary

File:media/webrtc/trunk/src/modules/audio_processing/ns/ns_core.c
Location:line 1124, column 13
Description:The left operand of '<' is a garbage value

Annotated Source Code

1/*
2 * Copyright (c) 2012 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((void*)0)) {
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_MAX256);
114 WebRtc_rdft(inst->anaLen, 1, inst->dataBuf, inst->ip, inst->wfft);
115
116 memset(inst->dataBuf, 0, sizeof(float) * ANAL_BLOCKL_MAX256);
117 memset(inst->syntBuf, 0, sizeof(float) * ANAL_BLOCKL_MAX256);
118
119 //for HB processing
120 memset(inst->dataBufHB, 0, sizeof(float) * ANAL_BLOCKL_MAX256);
121
122 //for quantile noise estimation
123 memset(inst->quantile, 0, sizeof(float) * HALF_ANAL_BLOCKL129);
124 for (i = 0; i < SIMULT3 * HALF_ANAL_BLOCKL129; i++) {
125 inst->lquantile[i] = (float)8.0;
126 inst->density[i] = (float)0.3;
127 }
128
129 for (i = 0; i < SIMULT3; i++) {
130 inst->counter[i] = (int)floor((float)(END_STARTUP_LONG200 * (i + 1)) / (float)SIMULT3);
131 }
132
133 inst->updates = 0;
134
135 // Wiener filter initialization
136 for (i = 0; i < HALF_ANAL_BLOCKL129; 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_BLOCKL129; 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(float)0.5; //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(float)0.5; //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(float)0.5; //average lrt factor (start on threshold)
159 inst->featureData[4] = SF_FEATURE_THR(float)0.5; //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_EST1000; 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(float)0.5; //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_MAX160);
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_BLOCKL129], delta;
243
244 if (inst->updates < END_STARTUP_LONG200) {
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 < SIMULT3; 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)40.0 * (float)1.0 / inst->density[offset + i];
261 } else {
262 delta = FACTOR(float)40.0;
263 }
264
265 // update log quantile estimate
266 if (lmagn[i] > inst->lquantile[offset + i]) {
267 inst->lquantile[offset + i] += QUANTILE(float)0.25 * delta
268 / (float)(inst->counter[s] + 1);
269 } else {
270 inst->lquantile[offset + i] -= ((float)1.0 - QUANTILE(float)0.25) * delta
271 / (float)(inst->counter[s] + 1);
272 }
273
274 // update density estimate
275 if (fabs(lmagn[i] - inst->lquantile[offset + i]) < WIDTH(float)0.01) {
276 inst->density[offset + i] = ((float)inst->counter[s] * inst->density[offset
277 + i] + (float)1.0 / ((float)2.0 * WIDTH(float)0.01)) / (float)(inst->counter[s] + 1);
278 }
279 } // end loop over magnitude spectrum
280
281 if (inst->counter[s] >= END_STARTUP_LONG200) {
282 inst->counter[s] = 0;
283 if (inst->updates >= END_STARTUP_LONG200) {
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_LONG200) {
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_EST1000 * 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_EST1000
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_EST1000
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_EST1000; 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_EST1000; 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_EST1000; 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_EST1000; 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(float)0.30 * 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(float)0.30 * (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(float)0.30 * (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(float)4.0;
622 widthPrior1 = (float)2.0 * WIDTH_PR_MAP(float)4.0; //width for pause region:
623 // lower range, so increase width in tanh map
624 widthPrior2 = (float)2.0 * WIDTH_PR_MAP(float)4.0; //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(float)0.50 * (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(float)0.10 * (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_MAX160], fout[BLOCKL_MAX160];
736 float winData[ANAL_BLOCKL_MAX256];
737 float magn[HALF_ANAL_BLOCKL129], noise[HALF_ANAL_BLOCKL129];
738 float theFilter[HALF_ANAL_BLOCKL129], theFilterTmp[HALF_ANAL_BLOCKL129];
739 float snrLocPost[HALF_ANAL_BLOCKL129], snrLocPrior[HALF_ANAL_BLOCKL129];
740 float probSpeechFinal[HALF_ANAL_BLOCKL129] = { 0 };
741 float previousEstimateStsa[HALF_ANAL_BLOCKL129];
742 float real[ANAL_BLOCKL_MAX256], imag[HALF_ANAL_BLOCKL129];
743 // Variables during startup
744 float sum_log_i = 0.0;
745 float sum_log_i_square = 0.0;
746 float sum_log_magn = 0.0;
747 float sum_log_i_log_magn = 0.0;
748 float parametric_noise = 0.0;
749 float parametric_exp = 0.0;
750 float parametric_num = 0.0;
751
752 // SWB variables
753 int deltaBweHB = 1;
754 int deltaGainHB = 1;
755 float decayBweHB = 1.0;
756 float gainMapParHB = 1.0;
757 float gainTimeDomainHB = 1.0;
758 float avgProbSpeechHB, avgProbSpeechHBTmp, avgFilterGainHB, gainModHB;
759
760 // Check that initiation has been done
761 if (inst->initFlag != 1) {
1
Taking false branch
→
762 return (-1);
763 }
764 // Check for valid pointers based on sampling rate
765 if (inst->fs == 32000) {
2
←
Taking false branch
→
766 if (speechFrameHB == NULL((void*)0)) {
767 return -1;
768 }
769 flagHB = 1;
770 // range for averaging low band quantities for H band gain
771 deltaBweHB = (int)inst->magnLen / 4;
772 deltaGainHB = deltaBweHB;
773 }
774 //
775 updateParsFlag = inst->modelUpdatePars[0];
776 //
777
778 //for LB do all processing
779 // convert to float
780 for (i = 0; i < inst->blockLen10ms; i++) {
3
←
Loop condition is false. Execution continues on line 785
→
781 fin[i] = (float)speechFrame[i];
782 }
783 // update analysis buffer for L band
784 memcpy(inst->dataBuf, inst->dataBuf + inst->blockLen10ms,
785 sizeof(float) * (inst->anaLen - inst->blockLen10ms));
786 memcpy(inst->dataBuf + inst->anaLen - inst->blockLen10ms, fin,
787 sizeof(float) * inst->blockLen10ms);
788
789 if (flagHB == 1) {
4
←
Taking false branch
→
790 // convert to float
791 for (i = 0; i < inst->blockLen10ms; i++) {
792 fin[i] = (float)speechFrameHB[i];
793 }
794 // update analysis buffer for H band
795 memcpy(inst->dataBufHB, inst->dataBufHB + inst->blockLen10ms,
796 sizeof(float) * (inst->anaLen - inst->blockLen10ms));
797 memcpy(inst->dataBufHB + inst->anaLen - inst->blockLen10ms, fin,
798 sizeof(float) * inst->blockLen10ms);
799 }
800
801 // check if processing needed
802 if (inst->outLen == 0) {
5
←
Taking true branch
→
803 // windowing
804 energy1 = 0.0;
805 for (i = 0; i < inst->anaLen; i++) {
6
←
Loop condition is false. Execution continues on line 809
→
806 winData[i] = inst->window[i] * inst->dataBuf[i];
807 energy1 += winData[i] * winData[i];
808 }
809 if (energy1 == 0.0) {
7
←
Taking false branch
→
810 // synthesize the special case of zero input
811 // we want to avoid updating statistics in this case:
812 // Updating feature statistics when we have zeros only will cause thresholds to
813 // move towards zero signal situations. This in turn has the effect that once the
814 // signal is "turned on" (non-zero values) everything will be treated as speech
815 // and there is no noise suppression effect. Depending on the duration of the
816 // inactive signal it takes a considerable amount of time for the system to learn
817 // what is noise and what is speech.
818
819 // read out fully processed segment
820 for (i = inst->windShift; i < inst->blockLen + inst->windShift; i++) {
821 fout[i - inst->windShift] = inst->syntBuf[i];
822 }
823 // update synthesis buffer
824 memcpy(inst->syntBuf, inst->syntBuf + inst->blockLen,
825 sizeof(float) * (inst->anaLen - inst->blockLen));
826 memset(inst->syntBuf + inst->anaLen - inst->blockLen, 0,
827 sizeof(float) * inst->blockLen);
828
829 // out buffer
830 inst->outLen = inst->blockLen - inst->blockLen10ms;
831 if (inst->blockLen > inst->blockLen10ms) {
832 for (i = 0; i < inst->outLen; i++) {
833 inst->outBuf[i] = fout[i + inst->blockLen10ms];
834 }
835 }
836 // convert to short
837 for (i = 0; i < inst->blockLen10ms; i++) {
838 dTmp = fout[i];
839 if (dTmp < WEBRTC_SPL_WORD16_MIN-32768) {
840 dTmp = WEBRTC_SPL_WORD16_MIN-32768;
841 } else if (dTmp > WEBRTC_SPL_WORD16_MAX32767) {
842 dTmp = WEBRTC_SPL_WORD16_MAX32767;
843 }
844 outFrame[i] = (short)dTmp;
845 }
846
847 // for time-domain gain of HB
848 if (flagHB == 1) {
849 for (i = 0; i < inst->blockLen10ms; i++) {
850 dTmp = inst->dataBufHB[i];
851 if (dTmp < WEBRTC_SPL_WORD16_MIN-32768) {
852 dTmp = WEBRTC_SPL_WORD16_MIN-32768;
853 } else if (dTmp > WEBRTC_SPL_WORD16_MAX32767) {
854 dTmp = WEBRTC_SPL_WORD16_MAX32767;
855 }
856 outFrameHB[i] = (short)dTmp;
857 }
858 } // end of H band gain computation
859 //
860 return 0;
861 }
862
863 //
864 inst->blockInd++; // Update the block index only when we process a block.
865 // FFT
866 WebRtc_rdft(inst->anaLen, 1, winData, inst->ip, inst->wfft);
867
868 imag[0] = 0;
869 real[0] = winData[0];
870 magn[0] = (float)(fabs(real[0]) + 1.0f);
871 imag[inst->magnLen - 1] = 0;
872 real[inst->magnLen - 1] = winData[1];
873 magn[inst->magnLen - 1] = (float)(fabs(real[inst->magnLen - 1]) + 1.0f);
874 signalEnergy = (float)(real[0] * real[0]) +
875 (float)(real[inst->magnLen - 1] * real[inst->magnLen - 1]);
876 sumMagn = magn[0] + magn[inst->magnLen - 1];
877 if (inst->blockInd < END_STARTUP_SHORT50) {
8
←
Taking false branch
→
878 inst->initMagnEst[0] += magn[0];
879 inst->initMagnEst[inst->magnLen - 1] += magn[inst->magnLen - 1];
880 tmpFloat2 = log((float)(inst->magnLen - 1));
881 sum_log_i = tmpFloat2;
882 sum_log_i_square = tmpFloat2 * tmpFloat2;
883 tmpFloat1 = log(magn[inst->magnLen - 1]);
884 sum_log_magn = tmpFloat1;
885 sum_log_i_log_magn = tmpFloat2 * tmpFloat1;
886 }
887 for (i = 1; i < inst->magnLen - 1; i++) {
9
←
Loop condition is false. Execution continues on line 908
→
888 real[i] = winData[2 * i];
889 imag[i] = winData[2 * i + 1];
890 // magnitude spectrum
891 fTmp = real[i] * real[i];
892 fTmp += imag[i] * imag[i];
893 signalEnergy += fTmp;
894 magn[i] = ((float)sqrt(fTmp)) + 1.0f;
895 sumMagn += magn[i];
896 if (inst->blockInd < END_STARTUP_SHORT50) {
897 inst->initMagnEst[i] += magn[i];
898 if (i >= kStartBand) {
899 tmpFloat2 = log((float)i);
900 sum_log_i += tmpFloat2;
901 sum_log_i_square += tmpFloat2 * tmpFloat2;
902 tmpFloat1 = log(magn[i]);
903 sum_log_magn += tmpFloat1;
904 sum_log_i_log_magn += tmpFloat2 * tmpFloat1;
905 }
906 }
907 }
908 signalEnergy = signalEnergy / ((float)inst->magnLen);
909 inst->signalEnergy = signalEnergy;
910 inst->sumMagn = sumMagn;
911
912 //compute spectral flatness on input spectrum
913 WebRtcNs_ComputeSpectralFlatness(inst, magn);
914 // quantile noise estimate
915 WebRtcNs_NoiseEstimation(inst, magn, noise);
916 //compute simplified noise model during startup
917 if (inst->blockInd < END_STARTUP_SHORT50) {
10
←
Taking false branch
→
918 // Estimate White noise
919 inst->whiteNoiseLevel += sumMagn / ((float)inst->magnLen) * inst->overdrive;
920 // Estimate Pink noise parameters
921 tmpFloat1 = sum_log_i_square * ((float)(inst->magnLen - kStartBand));
922 tmpFloat1 -= (sum_log_i * sum_log_i);
923 tmpFloat2 = (sum_log_i_square * sum_log_magn - sum_log_i * sum_log_i_log_magn);
924 tmpFloat3 = tmpFloat2 / tmpFloat1;
925 // Constrain the estimated spectrum to be positive
926 if (tmpFloat3 < 0.0f) {
927 tmpFloat3 = 0.0f;
928 }
929 inst->pinkNoiseNumerator += tmpFloat3;
930 tmpFloat2 = (sum_log_i * sum_log_magn);
931 tmpFloat2 -= ((float)(inst->magnLen - kStartBand)) * sum_log_i_log_magn;
932 tmpFloat3 = tmpFloat2 / tmpFloat1;
933 // Constrain the pink noise power to be in the interval [0, 1];
934 if (tmpFloat3 < 0.0f) {
935 tmpFloat3 = 0.0f;
936 }
937 if (tmpFloat3 > 1.0f) {
938 tmpFloat3 = 1.0f;
939 }
940 inst->pinkNoiseExp += tmpFloat3;
941
942 // Calculate frequency independent parts of parametric noise estimate.
943 if (inst->pinkNoiseExp == 0.0f) {
944 // Use white noise estimate
945 parametric_noise = inst->whiteNoiseLevel;
946 } else {
947 // Use pink noise estimate
948 parametric_num = exp(inst->pinkNoiseNumerator / (float)(inst->blockInd + 1));
949 parametric_num *= (float)(inst->blockInd + 1);
950 parametric_exp = inst->pinkNoiseExp / (float)(inst->blockInd + 1);
951 parametric_noise = parametric_num / pow((float)kStartBand, parametric_exp);
952 }
953 for (i = 0; i < inst->magnLen; i++) {
954 // Estimate the background noise using the white and pink noise parameters
955 if ((inst->pinkNoiseExp > 0.0f) && (i >= kStartBand)) {
956 // Use pink noise estimate
957 parametric_noise = parametric_num / pow((float)i, parametric_exp);
958 }
959 theFilterTmp[i] = (inst->initMagnEst[i] - inst->overdrive * parametric_noise);
960 theFilterTmp[i] /= (inst->initMagnEst[i] + (float)0.0001);
961 // Weight quantile noise with modeled noise
962 noise[i] *= (inst->blockInd);
963 tmpFloat2 = parametric_noise * (END_STARTUP_SHORT50 - inst->blockInd);
964 noise[i] += (tmpFloat2 / (float)(inst->blockInd + 1));
965 noise[i] /= END_STARTUP_SHORT50;
966 }
967 }
968 //compute average signal during END_STARTUP_LONG time:
969 // used to normalize spectral difference measure
970 if (inst->blockInd < END_STARTUP_LONG200) {
11
←
Taking false branch
→
971 inst->featureData[5] *= inst->blockInd;
972 inst->featureData[5] += signalEnergy;
973 inst->featureData[5] /= (inst->blockInd + 1);
974 }
975
976#ifdef PROCESS_FLOW_0
977 if (inst->blockInd > END_STARTUP_LONG200) {
978 //option: average the quantile noise: for check with AEC2
979 for (i = 0; i < inst->magnLen; i++) {
980 noise[i] = (float)0.6 * inst->noisePrev[i] + (float)0.4 * noise[i];
981 }
982 for (i = 0; i < inst->magnLen; i++) {
983 // Wiener with over sub-substraction:
984 theFilter[i] = (magn[i] - inst->overdrive * noise[i]) / (magn[i] + (float)0.0001);
985 }
986 }
987#else
988 //start processing at frames == converged+1
989 //
990 // STEP 1: compute prior and post snr based on quantile noise est
991 //
992
993 // compute DD estimate of prior SNR: needed for new method
994 for (i = 0; i < inst->magnLen; i++) {
12
←
Loop condition is true. Entering loop body
→
14
←
Loop condition is false. Execution continues on line 1025
→
995 // post snr
996 snrLocPost[i] = (float)0.0;
997 if (magn[i] > noise[i]) {
13
←
Taking false branch
→
998 snrLocPost[i] = magn[i] / (noise[i] + (float)0.0001) - (float)1.0;
999 }
1000 // previous post snr
1001 // previous estimate: based on previous frame with gain filter
1002 previousEstimateStsa[i] = inst->magnPrev[i] / (inst->noisePrev[i] + (float)0.0001)
1003 * (inst->smooth[i]);
1004 // DD estimate is sum of two terms: current estimate and previous estimate
1005 // directed decision update of snrPrior
1006 snrLocPrior[i] = DD_PR_SNR(float)0.98 * previousEstimateStsa[i] + ((float)1.0 - DD_PR_SNR(float)0.98)
1007 * snrLocPost[i];
1008 // post and prior snr needed for step 2
1009 } // end of loop over freqs
1010#ifdef PROCESS_FLOW_1
1011 for (i = 0; i < inst->magnLen; i++) {
1012 // gain filter
1013 tmpFloat1 = inst->overdrive + snrLocPrior[i];
1014 tmpFloat2 = (float)snrLocPrior[i] / tmpFloat1;
1015 theFilter[i] = (float)tmpFloat2;
1016 } // end of loop over freqs
1017#endif
1018 // done with step 1: dd computation of prior and post snr
1019
1020 //
1021 //STEP 2: compute speech/noise likelihood
1022 //
1023#ifdef PROCESS_FLOW_2
1024 // compute difference of input spectrum with learned/estimated noise spectrum
1025 WebRtcNs_ComputeSpectralDifference(inst, magn);
1026 // compute histograms for parameter decisions (thresholds and weights for features)
1027 // parameters are extracted once every window time (=inst->modelUpdatePars[1])
1028 if (updateParsFlag >= 1) {
15
←
Assuming 'updateParsFlag' is >= 1
→
16
←
Taking true branch
→
1029 // counter update
1030 inst->modelUpdatePars[3]--;
1031 // update histogram
1032 if (inst->modelUpdatePars[3] > 0) {
17
←
Taking true branch
→
1033 WebRtcNs_FeatureParameterExtraction(inst, 0);
1034 }
1035 // compute model parameters
1036 if (inst->modelUpdatePars[3] == 0) {
18
←
Taking false branch
→
1037 WebRtcNs_FeatureParameterExtraction(inst, 1);
1038 inst->modelUpdatePars[3] = inst->modelUpdatePars[1];
1039 // if wish to update only once, set flag to zero
1040 if (updateParsFlag == 1) {
1041 inst->modelUpdatePars[0] = 0;
1042 } else {
1043 // update every window:
1044 // get normalization for spectral difference for next window estimate
1045 inst->featureData[6] = inst->featureData[6]
1046 / ((float)inst->modelUpdatePars[1]);
1047 inst->featureData[5] = (float)0.5 * (inst->featureData[6]
1048 + inst->featureData[5]);
1049 inst->featureData[6] = (float)0.0;
1050 }
1051 }
1052 }
1053 // compute speech/noise probability
1054 WebRtcNs_SpeechNoiseProb(inst, probSpeechFinal, snrLocPrior, snrLocPost);
1055 // time-avg parameter for noise update
1056 gammaNoiseTmp = NOISE_UPDATE(float)0.90;
1057 for (i = 0; i < inst->magnLen; i++) {
19
←
Loop condition is true. Entering loop body
→
23
←
Loop condition is false. Execution continues on line 1094
→
1058 probSpeech = probSpeechFinal[i];
1059 probNonSpeech = (float)1.0 - probSpeech;
1060 // temporary noise update:
1061 // use it for speech frames if update value is less than previous
1062 noiseUpdateTmp = gammaNoiseTmp * inst->noisePrev[i] + ((float)1.0 - gammaNoiseTmp)
1063 * (probNonSpeech * magn[i] + probSpeech * inst->noisePrev[i]);
1064 //
1065 // time-constant based on speech/noise state
1066 gammaNoiseOld = gammaNoiseTmp;
1067 gammaNoiseTmp = NOISE_UPDATE(float)0.90;
1068 // increase gamma (i.e., less noise update) for frame likely to be speech
1069 if (probSpeech > PROB_RANGE(float)0.20) {
20
←
Taking false branch
→
1070 gammaNoiseTmp = SPEECH_UPDATE(float)0.99;
1071 }
1072 // conservative noise update
1073 if (probSpeech < PROB_RANGE(float)0.20) {
21
←
Taking false branch
→
1074 inst->magnAvgPause[i] += GAMMA_PAUSE(float)0.05 * (magn[i] - inst->magnAvgPause[i]);
1075 }
1076 // noise update
1077 if (gammaNoiseTmp == gammaNoiseOld) {
22
←
Taking true branch
→
1078 noise[i] = noiseUpdateTmp;
1079 } else {
1080 noise[i] = gammaNoiseTmp * inst->noisePrev[i] + ((float)1.0 - gammaNoiseTmp)
1081 * (probNonSpeech * magn[i] + probSpeech * inst->noisePrev[i]);
1082 // allow for noise update downwards:
1083 // if noise update decreases the noise, it is safe, so allow it to happen
1084 if (noiseUpdateTmp < noise[i]) {
1085 noise[i] = noiseUpdateTmp;
1086 }
1087 }
1088 } // end of freq loop
1089 // done with step 2: noise update
1090
1091 //
1092 // STEP 3: compute dd update of prior snr and post snr based on new noise estimate
1093 //
1094 for (i = 0; i < inst->magnLen; i++) {
24
←
Loop condition is true. Entering loop body
→
26
←
Loop condition is false. Execution continues on line 1113
→
1095 // post and prior snr
1096 currentEstimateStsa = (float)0.0;
1097 if (magn[i] > noise[i]) {
25
←
Taking false branch
→
1098 currentEstimateStsa = magn[i] / (noise[i] + (float)0.0001) - (float)1.0;
1099 }
1100 // DD estimate is sume of two terms: current estimate and previous estimate
1101 // directed decision update of snrPrior
1102 snrPrior = DD_PR_SNR(float)0.98 * previousEstimateStsa[i] + ((float)1.0 - DD_PR_SNR(float)0.98)
1103 * currentEstimateStsa;
1104 // gain filter
1105 tmpFloat1 = inst->overdrive + snrPrior;
1106 tmpFloat2 = (float)snrPrior / tmpFloat1;
1107 theFilter[i] = (float)tmpFloat2;
1108 } // end of loop over freqs
1109 // done with step3
1110#endif
1111#endif
1112
1113 for (i = 0; i < inst->magnLen; i++) {
27
←
Loop condition is true. Entering loop body
→
1114 // flooring bottom
1115 if (theFilter[i] < inst->denoiseBound) {
28
←
Taking false branch
→
1116 theFilter[i] = inst->denoiseBound;
1117 }
1118 // flooring top
1119 if (theFilter[i] > (float)1.0) {
29
←
Taking false branch
→
1120 theFilter[i] = 1.0;
1121 }
1122 if (inst->blockInd < END_STARTUP_SHORT50) {
30
←
Taking true branch
→
1123 // flooring bottom
1124 if (theFilterTmp[i] < inst->denoiseBound) {
31
←
The left operand of '<' is a garbage value
1125 theFilterTmp[i] = inst->denoiseBound;
1126 }
1127 // flooring top
1128 if (theFilterTmp[i] > (float)1.0) {
1129 theFilterTmp[i] = 1.0;
1130 }
1131 // Weight the two suppression filters
1132 theFilter[i] *= (inst->blockInd);
1133 theFilterTmp[i] *= (END_STARTUP_SHORT50 - inst->blockInd);
1134 theFilter[i] += theFilterTmp[i];
1135 theFilter[i] /= (END_STARTUP_SHORT50);
1136 }
1137 // smoothing
1138#ifdef PROCESS_FLOW_0
1139 inst->smooth[i] *= SMOOTH(float)0.75; // value set to 0.7 in define.h file
1140 inst->smooth[i] += ((float)1.0 - SMOOTH(float)0.75) * theFilter[i];
1141#else
1142 inst->smooth[i] = theFilter[i];
1143#endif
1144 real[i] *= inst->smooth[i];
1145 imag[i] *= inst->smooth[i];
1146 }
1147 // keep track of noise and magn spectrum for next frame
1148 for (i = 0; i < inst->magnLen; i++) {
1149 inst->noisePrev[i] = noise[i];
1150 inst->magnPrev[i] = magn[i];
1151 }
1152 // back to time domain
1153 winData[0] = real[0];
1154 winData[1] = real[inst->magnLen - 1];
1155 for (i = 1; i < inst->magnLen - 1; i++) {
1156 winData[2 * i] = real[i];
1157 winData[2 * i + 1] = imag[i];
1158 }
1159 WebRtc_rdft(inst->anaLen, -1, winData, inst->ip, inst->wfft);
1160
1161 for (i = 0; i < inst->anaLen; i++) {
1162 real[i] = 2.0f * winData[i] / inst->anaLen; // fft scaling
1163 }
1164
1165 //scale factor: only do it after END_STARTUP_LONG time
1166 factor = (float)1.0;
1167 if (inst->gainmap == 1 && inst->blockInd > END_STARTUP_LONG200) {
1168 factor1 = (float)1.0;
1169 factor2 = (float)1.0;
1170
1171 energy2 = 0.0;
1172 for (i = 0; i < inst->anaLen; i++) {
1173 energy2 += (float)real[i] * (float)real[i];
1174 }
1175 gain = (float)sqrt(energy2 / (energy1 + (float)1.0));
1176
1177#ifdef PROCESS_FLOW_2
1178 // scaling for new version
1179 if (gain > B_LIM(float)0.5) {
1180 factor1 = (float)1.0 + (float)1.3 * (gain - B_LIM(float)0.5);
1181 if (gain * factor1 > (float)1.0) {
1182 factor1 = (float)1.0 / gain;
1183 }
1184 }
1185 if (gain < B_LIM(float)0.5) {
1186 //don't reduce scale too much for pause regions:
1187 // attenuation here should be controlled by flooring
1188 if (gain <= inst->denoiseBound) {
1189 gain = inst->denoiseBound;
1190 }
1191 factor2 = (float)1.0 - (float)0.3 * (B_LIM(float)0.5 - gain);
1192 }
1193 //combine both scales with speech/noise prob:
1194 // note prior (priorSpeechProb) is not frequency dependent
1195 factor = inst->priorSpeechProb * factor1 + ((float)1.0 - inst->priorSpeechProb)
1196 * factor2;
1197#else
1198 if (gain > B_LIM(float)0.5) {
1199 factor = (float)1.0 + (float)1.3 * (gain - B_LIM(float)0.5);
1200 } else {
1201 factor = (float)1.0 + (float)2.0 * (gain - B_LIM(float)0.5);
1202 }
1203 if (gain * factor > (float)1.0) {
1204 factor = (float)1.0 / gain;
1205 }
1206#endif
1207 } // out of inst->gainmap==1
1208
1209 // synthesis
1210 for (i = 0; i < inst->anaLen; i++) {
1211 inst->syntBuf[i] += factor * inst->window[i] * (float)real[i];
1212 }
1213 // read out fully processed segment
1214 for (i = inst->windShift; i < inst->blockLen + inst->windShift; i++) {
1215 fout[i - inst->windShift] = inst->syntBuf[i];
1216 }
1217 // update synthesis buffer
1218 memcpy(inst->syntBuf, inst->syntBuf + inst->blockLen,
1219 sizeof(float) * (inst->anaLen - inst->blockLen));
1220 memset(inst->syntBuf + inst->anaLen - inst->blockLen, 0,
1221 sizeof(float) * inst->blockLen);
1222
1223 // out buffer
1224 inst->outLen = inst->blockLen - inst->blockLen10ms;
1225 if (inst->blockLen > inst->blockLen10ms) {
1226 for (i = 0; i < inst->outLen; i++) {
1227 inst->outBuf[i] = fout[i + inst->blockLen10ms];
1228 }
1229 }
1230 } // end of if out.len==0
1231 else {
1232 for (i = 0; i < inst->blockLen10ms; i++) {
1233 fout[i] = inst->outBuf[i];
1234 }
1235 memcpy(inst->outBuf, inst->outBuf + inst->blockLen10ms,
1236 sizeof(float) * (inst->outLen - inst->blockLen10ms));
1237 memset(inst->outBuf + inst->outLen - inst->blockLen10ms, 0,
1238 sizeof(float) * inst->blockLen10ms);
1239 inst->outLen -= inst->blockLen10ms;
1240 }
1241
1242 // convert to short
1243 for (i = 0; i < inst->blockLen10ms; i++) {
1244 dTmp = fout[i];
1245 if (dTmp < WEBRTC_SPL_WORD16_MIN-32768) {
1246 dTmp = WEBRTC_SPL_WORD16_MIN-32768;
1247 } else if (dTmp > WEBRTC_SPL_WORD16_MAX32767) {
1248 dTmp = WEBRTC_SPL_WORD16_MAX32767;
1249 }
1250 outFrame[i] = (short)dTmp;
1251 }
1252
1253 // for time-domain gain of HB
1254 if (flagHB == 1) {
1255 for (i = 0; i < inst->magnLen; i++) {
1256 inst->speechProbHB[i] = probSpeechFinal[i];
1257 }
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 //make sure gain is within flooring range
1282 // flooring bottom
1283 if (gainTimeDomainHB < inst->denoiseBound) {
1284 gainTimeDomainHB = inst->denoiseBound;
1285 }
1286 // flooring top
1287 if (gainTimeDomainHB > (float)1.0) {
1288 gainTimeDomainHB = 1.0;
1289 }
1290 //apply gain
1291 for (i = 0; i < inst->blockLen10ms; i++) {
1292 dTmp = gainTimeDomainHB * inst->dataBufHB[i];
1293 if (dTmp < WEBRTC_SPL_WORD16_MIN-32768) {
1294 dTmp = WEBRTC_SPL_WORD16_MIN-32768;
1295 } else if (dTmp > WEBRTC_SPL_WORD16_MAX32767) {
1296 dTmp = WEBRTC_SPL_WORD16_MAX32767;
1297 }
1298 outFrameHB[i] = (short)dTmp;
1299 }
1300 } // end of H band gain computation
1301 //
1302
1303 return 0;
1304}