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