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 "../other/signal_processing_library.h"
17 #include "include/noise_suppression.h"
18 #include "ns_core.h"
19 #include "windows_private.h"
20 #include "../other/fft4g.h"
21
22 // Set Feature Extraction Parameters.
set_feature_extraction_parameters(NSinst_t * self)23 static void set_feature_extraction_parameters(NSinst_t* 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(NSinst_t * self,uint32_t fs)74 int WebRtcNs_InitCore(NSinst_t* 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) {
83 self->fs = fs;
84 } else {
85 return -1;
86 }
87 self->windShift = 0;
88 if (fs == 8000) {
89 // We only support 10ms frames.
90 self->blockLen = 80;
91 self->anaLen = 128;
92 self->window = kBlocks80w128;
93 } else if (fs == 16000) {
94 // We only support 10ms frames.
95 self->blockLen = 160;
96 self->anaLen = 256;
97 self->window = kBlocks160w256;
98 } else if (fs == 32000) {
99 // We only support 10ms frames.
100 self->blockLen = 160;
101 self->anaLen = 256;
102 self->window = kBlocks160w256;
103 }
104 self->magnLen = self->anaLen / 2 + 1; // Number of frequency bins.
105
106 // Initialize FFT work arrays.
107 self->ip[0] = 0; // Setting this triggers initialization.
108 memset(self->dataBuf, 0, sizeof(float) * ANAL_BLOCKL_MAX);
109 WebRtc_rdft(self->anaLen, 1, self->dataBuf, self->ip, self->wfft);
110
111 memset(self->analyzeBuf, 0, sizeof(float) * ANAL_BLOCKL_MAX);
112 memset(self->dataBuf, 0, sizeof(float) * ANAL_BLOCKL_MAX);
113 memset(self->syntBuf, 0, sizeof(float) * ANAL_BLOCKL_MAX);
114
115 // For HB processing.
116 memset(self->dataBufHB, 0, sizeof(float) * ANAL_BLOCKL_MAX);
117
118 // For quantile noise estimation.
119 memset(self->quantile, 0, sizeof(float) * HALF_ANAL_BLOCKL);
120 for (i = 0; i < SIMULT * HALF_ANAL_BLOCKL; i++) {
121 self->lquantile[i] = 8.f;
122 self->density[i] = 0.3f;
123 }
124
125 for (i = 0; i < SIMULT; i++) {
126 self->counter[i] =
127 (int)floor((float)(END_STARTUP_LONG * (i + 1)) / (float)SIMULT);
128 }
129
130 self->updates = 0;
131
132 // Wiener filter initialization.
133 for (i = 0; i < HALF_ANAL_BLOCKL; i++) {
134 self->smooth[i] = 1.f;
135 }
136
137 // Set the aggressiveness: default.
138 self->aggrMode = 0;
139
140 // Initialize variables for new method.
141 self->priorSpeechProb = 0.5f; // Prior prob for speech/noise.
142 // Previous analyze mag spectrum.
143 memset(self->magnPrevAnalyze, 0, sizeof(float) * HALF_ANAL_BLOCKL);
144 // Previous process mag spectrum.
145 memset(self->magnPrevProcess, 0, sizeof(float) * HALF_ANAL_BLOCKL);
146 // Current noise-spectrum.
147 memset(self->noise, 0, sizeof(float) * HALF_ANAL_BLOCKL);
148 // Previous noise-spectrum.
149 memset(self->noisePrev, 0, sizeof(float) * HALF_ANAL_BLOCKL);
150 // Conservative noise spectrum estimate.
151 memset(self->magnAvgPause, 0, sizeof(float) * HALF_ANAL_BLOCKL);
152 // For estimation of HB in second pass.
153 memset(self->speechProb, 0, sizeof(float) * HALF_ANAL_BLOCKL);
154 // Initial average magnitude spectrum.
155 memset(self->initMagnEst, 0, sizeof(float) * HALF_ANAL_BLOCKL);
156 for (i = 0; i < HALF_ANAL_BLOCKL; i++) {
157 // Smooth LR (same as threshold).
158 self->logLrtTimeAvg[i] = LRT_FEATURE_THR;
159 }
160
161 // Feature quantities.
162 // Spectral flatness (start on threshold).
163 self->featureData[0] = SF_FEATURE_THR;
164 self->featureData[1] = 0.f; // Spectral entropy: not used in this version.
165 self->featureData[2] = 0.f; // Spectral variance: not used in this version.
166 // Average LRT factor (start on threshold).
167 self->featureData[3] = LRT_FEATURE_THR;
168 // Spectral template diff (start on threshold).
169 self->featureData[4] = SF_FEATURE_THR;
170 self->featureData[5] = 0.f; // Normalization for spectral difference.
171 // Window time-average of input magnitude spectrum.
172 self->featureData[6] = 0.f;
173
174 // Histogram quantities: used to estimate/update thresholds for features.
175 memset(self->histLrt, 0, sizeof(int) * HIST_PAR_EST);
176 memset(self->histSpecFlat, 0, sizeof(int) * HIST_PAR_EST);
177 memset(self->histSpecDiff, 0, sizeof(int) * HIST_PAR_EST);
178
179
180 self->blockInd = -1; // Frame counter.
181 // Default threshold for LRT feature.
182 self->priorModelPars[0] = LRT_FEATURE_THR;
183 // Threshold for spectral flatness: determined on-line.
184 self->priorModelPars[1] = 0.5f;
185 // sgn_map par for spectral measure: 1 for flatness measure.
186 self->priorModelPars[2] = 1.f;
187 // Threshold for template-difference feature: determined on-line.
188 self->priorModelPars[3] = 0.5f;
189 // Default weighting parameter for LRT feature.
190 self->priorModelPars[4] = 1.f;
191 // Default weighting parameter for spectral flatness feature.
192 self->priorModelPars[5] = 0.f;
193 // Default weighting parameter for spectral difference feature.
194 self->priorModelPars[6] = 0.f;
195
196 // Update flag for parameters:
197 // 0 no update, 1 = update once, 2 = update every window.
198 self->modelUpdatePars[0] = 2;
199 self->modelUpdatePars[1] = 500; // Window for update.
200 // Counter for update of conservative noise spectrum.
201 self->modelUpdatePars[2] = 0;
202 // Counter if the feature thresholds are updated during the sequence.
203 self->modelUpdatePars[3] = self->modelUpdatePars[1];
204
205 self->signalEnergy = 0.0;
206 self->sumMagn = 0.0;
207 self->whiteNoiseLevel = 0.0;
208 self->pinkNoiseNumerator = 0.0;
209 self->pinkNoiseExp = 0.0;
210
211 set_feature_extraction_parameters(self);
212
213 // Default mode.
214 WebRtcNs_set_policy_core(self, 0);
215
216 self->initFlag = 1;
217 return 0;
218 }
219
220 // Estimate noise.
NoiseEstimation(NSinst_t * self,float * magn,float * noise)221 static void NoiseEstimation(NSinst_t* self, float* magn, float* noise) {
222 int i, s, offset;
223 float lmagn[HALF_ANAL_BLOCKL], delta;
224
225 if (self->updates < END_STARTUP_LONG) {
226 self->updates++;
227 }
228
229 for (i = 0; i < self->magnLen; i++) {
230 lmagn[i] = (float)log(magn[i]);
231 }
232
233 // Loop over simultaneous estimates.
234 for (s = 0; s < SIMULT; s++) {
235 offset = s * self->magnLen;
236
237 // newquantest(...)
238 for (i = 0; i < self->magnLen; i++) {
239 // Compute delta.
240 if (self->density[offset + i] > 1.0) {
241 delta = FACTOR * 1.f / self->density[offset + i];
242 } else {
243 delta = FACTOR;
244 }
245
246 // Update log quantile estimate.
247 if (lmagn[i] > self->lquantile[offset + i]) {
248 self->lquantile[offset + i] +=
249 QUANTILE * delta / (float)(self->counter[s] + 1);
250 } else {
251 self->lquantile[offset + i] -=
252 (1.f - QUANTILE) * delta / (float)(self->counter[s] + 1);
253 }
254
255 // Update density estimate.
256 if (fabs(lmagn[i] - self->lquantile[offset + i]) < WIDTH) {
257 self->density[offset + i] =
258 ((float)self->counter[s] * self->density[offset + i] +
259 1.f / (2.f * WIDTH)) /
260 (float)(self->counter[s] + 1);
261 }
262 } // End loop over magnitude spectrum.
263
264 if (self->counter[s] >= END_STARTUP_LONG) {
265 self->counter[s] = 0;
266 if (self->updates >= END_STARTUP_LONG) {
267 for (i = 0; i < self->magnLen; i++) {
268 self->quantile[i] = (float)exp(self->lquantile[offset + i]);
269 }
270 }
271 }
272
273 self->counter[s]++;
274 } // End loop over simultaneous estimates.
275
276 // Sequentially update the noise during startup.
277 if (self->updates < END_STARTUP_LONG) {
278 // Use the last "s" to get noise during startup that differ from zero.
279 for (i = 0; i < self->magnLen; i++) {
280 self->quantile[i] = (float)exp(self->lquantile[offset + i]);
281 }
282 }
283
284 for (i = 0; i < self->magnLen; i++) {
285 noise[i] = self->quantile[i];
286 }
287 }
288
289 // Extract thresholds for feature parameters.
290 // Histograms are computed over some window size (given by
291 // self->modelUpdatePars[1]).
292 // Thresholds and weights are extracted every window.
293 // |flag| = 0 updates histogram only, |flag| = 1 computes the threshold/weights.
294 // Threshold and weights are returned in: self->priorModelPars.
FeatureParameterExtraction(NSinst_t * self,int flag)295 static void FeatureParameterExtraction(NSinst_t* self, int flag) {
296 int i, useFeatureSpecFlat, useFeatureSpecDiff, numHistLrt;
297 int maxPeak1, maxPeak2;
298 int weightPeak1SpecFlat, weightPeak2SpecFlat, weightPeak1SpecDiff,
299 weightPeak2SpecDiff;
300
301 float binMid, featureSum;
302 float posPeak1SpecFlat, posPeak2SpecFlat, posPeak1SpecDiff, posPeak2SpecDiff;
303 float fluctLrt, avgHistLrt, avgSquareHistLrt, avgHistLrtCompl;
304
305 // 3 features: LRT, flatness, difference.
306 // lrt_feature = self->featureData[3];
307 // flat_feature = self->featureData[0];
308 // diff_feature = self->featureData[4];
309
310 // Update histograms.
311 if (flag == 0) {
312 // LRT
313 if ((self->featureData[3] <
314 HIST_PAR_EST * self->featureExtractionParams.binSizeLrt) &&
315 (self->featureData[3] >= 0.0)) {
316 i = (int)(self->featureData[3] /
317 self->featureExtractionParams.binSizeLrt);
318 self->histLrt[i]++;
319 }
320 // Spectral flatness.
321 if ((self->featureData[0] <
322 HIST_PAR_EST * self->featureExtractionParams.binSizeSpecFlat) &&
323 (self->featureData[0] >= 0.0)) {
324 i = (int)(self->featureData[0] /
325 self->featureExtractionParams.binSizeSpecFlat);
326 self->histSpecFlat[i]++;
327 }
328 // Spectral difference.
329 if ((self->featureData[4] <
330 HIST_PAR_EST * self->featureExtractionParams.binSizeSpecDiff) &&
331 (self->featureData[4] >= 0.0)) {
332 i = (int)(self->featureData[4] /
333 self->featureExtractionParams.binSizeSpecDiff);
334 self->histSpecDiff[i]++;
335 }
336 }
337
338 // Extract parameters for speech/noise probability.
339 if (flag == 1) {
340 // LRT feature: compute the average over
341 // self->featureExtractionParams.rangeAvgHistLrt.
342 avgHistLrt = 0.0;
343 avgHistLrtCompl = 0.0;
344 avgSquareHistLrt = 0.0;
345 numHistLrt = 0;
346 for (i = 0; i < HIST_PAR_EST; i++) {
347 binMid = ((float)i + 0.5f) * self->featureExtractionParams.binSizeLrt;
348 if (binMid <= self->featureExtractionParams.rangeAvgHistLrt) {
349 avgHistLrt += self->histLrt[i] * binMid;
350 numHistLrt += self->histLrt[i];
351 }
352 avgSquareHistLrt += self->histLrt[i] * binMid * binMid;
353 avgHistLrtCompl += self->histLrt[i] * binMid;
354 }
355 if (numHistLrt > 0) {
356 avgHistLrt = avgHistLrt / ((float)numHistLrt);
357 }
358 avgHistLrtCompl = avgHistLrtCompl / ((float)self->modelUpdatePars[1]);
359 avgSquareHistLrt = avgSquareHistLrt / ((float)self->modelUpdatePars[1]);
360 fluctLrt = avgSquareHistLrt - avgHistLrt * avgHistLrtCompl;
361 // Get threshold for LRT feature.
362 if (fluctLrt < self->featureExtractionParams.thresFluctLrt) {
363 // Very low fluctuation, so likely noise.
364 self->priorModelPars[0] = self->featureExtractionParams.maxLrt;
365 } else {
366 self->priorModelPars[0] =
367 self->featureExtractionParams.factor1ModelPars * avgHistLrt;
368 // Check if value is within min/max range.
369 if (self->priorModelPars[0] < self->featureExtractionParams.minLrt) {
370 self->priorModelPars[0] = self->featureExtractionParams.minLrt;
371 }
372 if (self->priorModelPars[0] > self->featureExtractionParams.maxLrt) {
373 self->priorModelPars[0] = self->featureExtractionParams.maxLrt;
374 }
375 }
376 // Done with LRT feature.
377
378 // For spectral flatness and spectral difference: compute the main peaks of
379 // histogram.
380 maxPeak1 = 0;
381 maxPeak2 = 0;
382 posPeak1SpecFlat = 0.0;
383 posPeak2SpecFlat = 0.0;
384 weightPeak1SpecFlat = 0;
385 weightPeak2SpecFlat = 0;
386
387 // Peaks for flatness.
388 for (i = 0; i < HIST_PAR_EST; i++) {
389 binMid =
390 (i + 0.5f) * self->featureExtractionParams.binSizeSpecFlat;
391 if (self->histSpecFlat[i] > maxPeak1) {
392 // Found new "first" peak.
393 maxPeak2 = maxPeak1;
394 weightPeak2SpecFlat = weightPeak1SpecFlat;
395 posPeak2SpecFlat = posPeak1SpecFlat;
396
397 maxPeak1 = self->histSpecFlat[i];
398 weightPeak1SpecFlat = self->histSpecFlat[i];
399 posPeak1SpecFlat = binMid;
400 } else if (self->histSpecFlat[i] > maxPeak2) {
401 // Found new "second" peak.
402 maxPeak2 = self->histSpecFlat[i];
403 weightPeak2SpecFlat = self->histSpecFlat[i];
404 posPeak2SpecFlat = binMid;
405 }
406 }
407
408 // Compute two peaks for spectral difference.
409 maxPeak1 = 0;
410 maxPeak2 = 0;
411 posPeak1SpecDiff = 0.0;
412 posPeak2SpecDiff = 0.0;
413 weightPeak1SpecDiff = 0;
414 weightPeak2SpecDiff = 0;
415 // Peaks for spectral difference.
416 for (i = 0; i < HIST_PAR_EST; i++) {
417 binMid =
418 ((float)i + 0.5f) * self->featureExtractionParams.binSizeSpecDiff;
419 if (self->histSpecDiff[i] > maxPeak1) {
420 // Found new "first" peak.
421 maxPeak2 = maxPeak1;
422 weightPeak2SpecDiff = weightPeak1SpecDiff;
423 posPeak2SpecDiff = posPeak1SpecDiff;
424
425 maxPeak1 = self->histSpecDiff[i];
426 weightPeak1SpecDiff = self->histSpecDiff[i];
427 posPeak1SpecDiff = binMid;
428 } else if (self->histSpecDiff[i] > maxPeak2) {
429 // Found new "second" peak.
430 maxPeak2 = self->histSpecDiff[i];
431 weightPeak2SpecDiff = self->histSpecDiff[i];
432 posPeak2SpecDiff = binMid;
433 }
434 }
435
436 // For spectrum flatness feature.
437 useFeatureSpecFlat = 1;
438 // Merge the two peaks if they are close.
439 if ((fabs(posPeak2SpecFlat - posPeak1SpecFlat) <
440 self->featureExtractionParams.limitPeakSpacingSpecFlat) &&
441 (weightPeak2SpecFlat >
442 self->featureExtractionParams.limitPeakWeightsSpecFlat *
443 weightPeak1SpecFlat)) {
444 weightPeak1SpecFlat += weightPeak2SpecFlat;
445 posPeak1SpecFlat = 0.5f * (posPeak1SpecFlat + posPeak2SpecFlat);
446 }
447 // Reject if weight of peaks is not large enough, or peak value too small.
448 if (weightPeak1SpecFlat <
449 self->featureExtractionParams.thresWeightSpecFlat ||
450 posPeak1SpecFlat < self->featureExtractionParams.thresPosSpecFlat) {
451 useFeatureSpecFlat = 0;
452 }
453 // If selected, get the threshold.
454 if (useFeatureSpecFlat == 1) {
455 // Compute the threshold.
456 self->priorModelPars[1] =
457 self->featureExtractionParams.factor2ModelPars * posPeak1SpecFlat;
458 // Check if value is within min/max range.
459 if (self->priorModelPars[1] < self->featureExtractionParams.minSpecFlat) {
460 self->priorModelPars[1] = self->featureExtractionParams.minSpecFlat;
461 }
462 if (self->priorModelPars[1] > self->featureExtractionParams.maxSpecFlat) {
463 self->priorModelPars[1] = self->featureExtractionParams.maxSpecFlat;
464 }
465 }
466 // Done with flatness feature.
467
468 // For template feature.
469 useFeatureSpecDiff = 1;
470 // Merge the two peaks if they are close.
471 if ((fabs(posPeak2SpecDiff - posPeak1SpecDiff) <
472 self->featureExtractionParams.limitPeakSpacingSpecDiff) &&
473 (weightPeak2SpecDiff >
474 self->featureExtractionParams.limitPeakWeightsSpecDiff *
475 weightPeak1SpecDiff)) {
476 weightPeak1SpecDiff += weightPeak2SpecDiff;
477 posPeak1SpecDiff = 0.5f * (posPeak1SpecDiff + posPeak2SpecDiff);
478 }
479 // Get the threshold value.
480 self->priorModelPars[3] =
481 self->featureExtractionParams.factor1ModelPars * posPeak1SpecDiff;
482 // Reject if weight of peaks is not large enough.
483 if (weightPeak1SpecDiff <
484 self->featureExtractionParams.thresWeightSpecDiff) {
485 useFeatureSpecDiff = 0;
486 }
487 // Check if value is within min/max range.
488 if (self->priorModelPars[3] < self->featureExtractionParams.minSpecDiff) {
489 self->priorModelPars[3] = self->featureExtractionParams.minSpecDiff;
490 }
491 if (self->priorModelPars[3] > self->featureExtractionParams.maxSpecDiff) {
492 self->priorModelPars[3] = self->featureExtractionParams.maxSpecDiff;
493 }
494 // Done with spectral difference feature.
495
496 // Don't use template feature if fluctuation of LRT feature is very low:
497 // most likely just noise state.
498 if (fluctLrt < self->featureExtractionParams.thresFluctLrt) {
499 useFeatureSpecDiff = 0;
500 }
501
502 // Select the weights between the features.
503 // self->priorModelPars[4] is weight for LRT: always selected.
504 // self->priorModelPars[5] is weight for spectral flatness.
505 // self->priorModelPars[6] is weight for spectral difference.
506 featureSum = (float)(1 + useFeatureSpecFlat + useFeatureSpecDiff);
507 self->priorModelPars[4] = 1.f / featureSum;
508 self->priorModelPars[5] = ((float)useFeatureSpecFlat) / featureSum;
509 self->priorModelPars[6] = ((float)useFeatureSpecDiff) / featureSum;
510
511 // Set hists to zero for next update.
512 if (self->modelUpdatePars[0] >= 1) {
513 for (i = 0; i < HIST_PAR_EST; i++) {
514 self->histLrt[i] = 0;
515 self->histSpecFlat[i] = 0;
516 self->histSpecDiff[i] = 0;
517 }
518 }
519 } // End of flag == 1.
520 }
521
522 // Compute spectral flatness on input spectrum.
523 // |magnIn| is the magnitude spectrum.
524 // Spectral flatness is returned in self->featureData[0].
ComputeSpectralFlatness(NSinst_t * self,const float * magnIn)525 static void ComputeSpectralFlatness(NSinst_t* self, const float* magnIn) {
526 int i;
527 int shiftLP = 1; // Option to remove first bin(s) from spectral measures.
528 float avgSpectralFlatnessNum, avgSpectralFlatnessDen, spectralTmp;
529
530 // Compute spectral measures.
531 // For flatness.
532 avgSpectralFlatnessNum = 0.0;
533 avgSpectralFlatnessDen = self->sumMagn;
534 for (i = 0; i < shiftLP; i++) {
535 avgSpectralFlatnessDen -= magnIn[i];
536 }
537 // Compute log of ratio of the geometric to arithmetic mean: check for log(0)
538 // case.
539 for (i = shiftLP; i < self->magnLen; i++) {
540 if (magnIn[i] > 0.0) {
541 avgSpectralFlatnessNum += (float)log(magnIn[i]);
542 } else {
543 self->featureData[0] -= SPECT_FL_TAVG * self->featureData[0];
544 return;
545 }
546 }
547 // Normalize.
548 avgSpectralFlatnessDen = avgSpectralFlatnessDen / self->magnLen;
549 avgSpectralFlatnessNum = avgSpectralFlatnessNum / self->magnLen;
550
551 // Ratio and inverse log: check for case of log(0).
552 spectralTmp = (float)exp(avgSpectralFlatnessNum) / avgSpectralFlatnessDen;
553
554 // Time-avg update of spectral flatness feature.
555 self->featureData[0] += SPECT_FL_TAVG * (spectralTmp - self->featureData[0]);
556 // Done with flatness feature.
557 }
558
559 // Compute prior and post SNR based on quantile noise estimation.
560 // Compute DD estimate of prior SNR.
561 // Inputs:
562 // * |magn| is the signal magnitude spectrum estimate.
563 // * |noise| is the magnitude noise spectrum estimate.
564 // Outputs:
565 // * |snrLocPrior| is the computed prior SNR.
566 // * |snrLocPost| is the computed post SNR.
ComputeSnr(const NSinst_t * self,const float * magn,const float * noise,float * snrLocPrior,float * snrLocPost)567 static void ComputeSnr(const NSinst_t* self,
568 const float* magn,
569 const float* noise,
570 float* snrLocPrior,
571 float* snrLocPost) {
572 int i;
573
574 for (i = 0; i < self->magnLen; i++) {
575 // Previous post SNR.
576 // Previous estimate: based on previous frame with gain filter.
577 float previousEstimateStsa = self->magnPrevAnalyze[i] /
578 (self->noisePrev[i] + 0.0001f) * self->smooth[i];
579 // Post SNR.
580 snrLocPost[i] = 0.f;
581 if (magn[i] > noise[i]) {
582 snrLocPost[i] = magn[i] / (noise[i] + 0.0001f) - 1.f;
583 }
584 // DD estimate is sum of two terms: current estimate and previous estimate.
585 // Directed decision update of snrPrior.
586 snrLocPrior[i] =
587 DD_PR_SNR * previousEstimateStsa + (1.f - DD_PR_SNR) * snrLocPost[i];
588 } // End of loop over frequencies.
589 }
590
591 // Compute the difference measure between input spectrum and a template/learned
592 // noise spectrum.
593 // |magnIn| is the input spectrum.
594 // The reference/template spectrum is self->magnAvgPause[i].
595 // Returns (normalized) spectral difference in self->featureData[4].
ComputeSpectralDifference(NSinst_t * self,const float * magnIn)596 static void ComputeSpectralDifference(NSinst_t* self,
597 const float* magnIn) {
598 // avgDiffNormMagn = var(magnIn) - cov(magnIn, magnAvgPause)^2 /
599 // var(magnAvgPause)
600 int i;
601 float avgPause, avgMagn, covMagnPause, varPause, varMagn, avgDiffNormMagn;
602
603 avgPause = 0.0;
604 avgMagn = self->sumMagn;
605 // Compute average quantities.
606 for (i = 0; i < self->magnLen; i++) {
607 // Conservative smooth noise spectrum from pause frames.
608 avgPause += self->magnAvgPause[i];
609 }
610 avgPause = avgPause / ((float)self->magnLen);
611 avgMagn = avgMagn / ((float)self->magnLen);
612
613 covMagnPause = 0.0;
614 varPause = 0.0;
615 varMagn = 0.0;
616 // Compute variance and covariance quantities.
617 for (i = 0; i < self->magnLen; i++) {
618 covMagnPause += (magnIn[i] - avgMagn) * (self->magnAvgPause[i] - avgPause);
619 varPause +=
620 (self->magnAvgPause[i] - avgPause) * (self->magnAvgPause[i] - avgPause);
621 varMagn += (magnIn[i] - avgMagn) * (magnIn[i] - avgMagn);
622 }
623 covMagnPause = covMagnPause / ((float)self->magnLen);
624 varPause = varPause / ((float)self->magnLen);
625 varMagn = varMagn / ((float)self->magnLen);
626 // Update of average magnitude spectrum.
627 self->featureData[6] += self->signalEnergy;
628
629 avgDiffNormMagn =
630 varMagn - (covMagnPause * covMagnPause) / (varPause + 0.0001f);
631 // Normalize and compute time-avg update of difference feature.
632 avgDiffNormMagn = (float)(avgDiffNormMagn / (self->featureData[5] + 0.0001f));
633 self->featureData[4] +=
634 SPECT_DIFF_TAVG * (avgDiffNormMagn - self->featureData[4]);
635 }
636
637 // Compute speech/noise probability.
638 // Speech/noise probability is returned in |probSpeechFinal|.
639 // |magn| is the input magnitude spectrum.
640 // |noise| is the noise spectrum.
641 // |snrLocPrior| is the prior SNR for each frequency.
642 // |snrLocPost| is the post SNR for each frequency.
SpeechNoiseProb(NSinst_t * self,float * probSpeechFinal,const float * snrLocPrior,const float * snrLocPost)643 static void SpeechNoiseProb(NSinst_t* self,
644 float* probSpeechFinal,
645 const float* snrLocPrior,
646 const float* snrLocPost) {
647 int i, 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(NSinst_t * self,const float * magn,int updateParsFlag)756 static void FeatureUpdate(NSinst_t* 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(NSinst_t * self,const float * magn,const float * snrLocPrior,const float * snrLocPost,float * noise)801 static void UpdateNoiseEstimate(NSinst_t* self,
802 const float* magn,
803 const float* snrLocPrior,
804 const float* snrLocPost,
805 float* noise) {
806 int 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,int frame_length,int buffer_length,float * buffer)856 static void UpdateBuffer(const float* frame,
857 int frame_length,
858 int 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(NSinst_t * self,float * time_data,int time_data_length,int magnitude_length,float * real,float * imag,float * magn)887 static void FFT(NSinst_t* self,
888 float* time_data,
889 int time_data_length,
890 int magnitude_length,
891 float* real,
892 float* imag,
893 float* magn) {
894 int 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] = fabs(real[0]) + 1.f;
903 imag[magnitude_length - 1] = 0;
904 real[magnitude_length - 1] = time_data[1];
905 magn[magnitude_length - 1] = fabs(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(NSinst_t * self,const float * real,const float * imag,int magnitude_length,int time_data_length,float * time_data)924 static void IFFT(NSinst_t* self,
925 const float* real,
926 const float* imag,
927 int magnitude_length,
928 int time_data_length,
929 float* time_data) {
930 int 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,int length)952 static float Energy(const float* buffer, int length) {
953 int 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,int length,float * data_windowed)970 static void Windowing(const float* window,
971 const float* data,
972 int length,
973 float* data_windowed) {
974 int 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 NSinst_t * self,const float * magn,float * theFilter)986 static void ComputeDdBasedWienerFilter(const NSinst_t* self,
987 const float* magn,
988 float* theFilter) {
989 int 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(NSinst_t * self,int mode)1014 int WebRtcNs_set_policy_core(NSinst_t* 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(NSinst_t * self,float * speechFrame)1044 int WebRtcNs_AnalyzeCore(NSinst_t* self, float* speechFrame) {
1045 int i;
1046 const int 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 if (self->initFlag != 1) {
1066 return (-1);
1067 }
1068 updateParsFlag = self->modelUpdatePars[0];
1069
1070 // Update analysis buffer for L band.
1071 UpdateBuffer(speechFrame, self->blockLen, self->anaLen, self->analyzeBuf);
1072
1073 Windowing(self->window, self->analyzeBuf, self->anaLen, winData);
1074 energy = Energy(winData, self->anaLen);
1075 if (energy == 0.0) {
1076 // We want to avoid updating statistics in this case:
1077 // Updating feature statistics when we have zeros only will cause
1078 // thresholds to move towards zero signal situations. This in turn has the
1079 // effect that once the signal is "turned on" (non-zero values) everything
1080 // will be treated as speech and there is no noise suppression effect.
1081 // Depending on the duration of the inactive signal it takes a
1082 // considerable amount of time for the system to learn what is noise and
1083 // what is speech.
1084 return 0;
1085 }
1086
1087 self->blockInd++; // Update the block index only when we process a block.
1088
1089 FFT(self, winData, self->anaLen, self->magnLen, real, imag, magn);
1090
1091 for (i = 0; i < self->magnLen; i++) {
1092 signalEnergy += real[i] * real[i] + imag[i] * imag[i];
1093 sumMagn += magn[i];
1094 if (self->blockInd < END_STARTUP_SHORT) {
1095 if (i >= kStartBand) {
1096 tmpFloat2 = log((float)i);
1097 sum_log_i += tmpFloat2;
1098 sum_log_i_square += tmpFloat2 * tmpFloat2;
1099 tmpFloat1 = log(magn[i]);
1100 sum_log_magn += tmpFloat1;
1101 sum_log_i_log_magn += tmpFloat2 * tmpFloat1;
1102 }
1103 }
1104 }
1105 signalEnergy = signalEnergy / ((float)self->magnLen);
1106 self->signalEnergy = signalEnergy;
1107 self->sumMagn = sumMagn;
1108
1109 // Quantile noise estimate.
1110 NoiseEstimation(self, magn, noise);
1111 // Compute simplified noise model during startup.
1112 if (self->blockInd < END_STARTUP_SHORT) {
1113 // Estimate White noise.
1114 self->whiteNoiseLevel += sumMagn / ((float)self->magnLen) * self->overdrive;
1115 // Estimate Pink noise parameters.
1116 tmpFloat1 = sum_log_i_square * ((float)(self->magnLen - kStartBand));
1117 tmpFloat1 -= (sum_log_i * sum_log_i);
1118 tmpFloat2 =
1119 (sum_log_i_square * sum_log_magn - sum_log_i * sum_log_i_log_magn);
1120 tmpFloat3 = tmpFloat2 / tmpFloat1;
1121 // Constrain the estimated spectrum to be positive.
1122 if (tmpFloat3 < 0.f) {
1123 tmpFloat3 = 0.f;
1124 }
1125 self->pinkNoiseNumerator += tmpFloat3;
1126 tmpFloat2 = (sum_log_i * sum_log_magn);
1127 tmpFloat2 -= ((float)(self->magnLen - kStartBand)) * sum_log_i_log_magn;
1128 tmpFloat3 = tmpFloat2 / tmpFloat1;
1129 // Constrain the pink noise power to be in the interval [0, 1].
1130 if (tmpFloat3 < 0.f) {
1131 tmpFloat3 = 0.f;
1132 }
1133 if (tmpFloat3 > 1.f) {
1134 tmpFloat3 = 1.f;
1135 }
1136 self->pinkNoiseExp += tmpFloat3;
1137
1138 // Calculate frequency independent parts of parametric noise estimate.
1139 if (self->pinkNoiseExp > 0.f) {
1140 // Use pink noise estimate.
1141 parametric_num =
1142 exp(self->pinkNoiseNumerator / (float)(self->blockInd + 1));
1143 parametric_num *= (float)(self->blockInd + 1);
1144 parametric_exp = self->pinkNoiseExp / (float)(self->blockInd + 1);
1145 }
1146 for (i = 0; i < self->magnLen; i++) {
1147 // Estimate the background noise using the white and pink noise
1148 // parameters.
1149 if (self->pinkNoiseExp == 0.f) {
1150 // Use white noise estimate.
1151 self->parametricNoise[i] = self->whiteNoiseLevel;
1152 } else {
1153 // Use pink noise estimate.
1154 float use_band = (float)(i < kStartBand ? kStartBand : i);
1155 self->parametricNoise[i] =
1156 parametric_num / pow(use_band, parametric_exp);
1157 }
1158 // Weight quantile noise with modeled noise.
1159 noise[i] *= (self->blockInd);
1160 tmpFloat2 =
1161 self->parametricNoise[i] * (END_STARTUP_SHORT - self->blockInd);
1162 noise[i] += (tmpFloat2 / (float)(self->blockInd + 1));
1163 noise[i] /= END_STARTUP_SHORT;
1164 }
1165 }
1166 // Compute average signal during END_STARTUP_LONG time:
1167 // used to normalize spectral difference measure.
1168 if (self->blockInd < END_STARTUP_LONG) {
1169 self->featureData[5] *= self->blockInd;
1170 self->featureData[5] += signalEnergy;
1171 self->featureData[5] /= (self->blockInd + 1);
1172 }
1173
1174 // Post and prior SNR needed for SpeechNoiseProb.
1175 ComputeSnr(self, magn, noise, snrLocPrior, snrLocPost);
1176
1177 FeatureUpdate(self, magn, updateParsFlag);
1178 SpeechNoiseProb(self, self->speechProb, snrLocPrior, snrLocPost);
1179 UpdateNoiseEstimate(self, magn, snrLocPrior, snrLocPost, noise);
1180
1181 // Keep track of noise spectrum for next frame.
1182 memcpy(self->noise, noise, sizeof(*noise) * self->magnLen);
1183 memcpy(self->magnPrevAnalyze, magn, sizeof(*magn) * self->magnLen);
1184
1185 return 0;
1186 }
1187
WebRtcNs_ProcessCore(NSinst_t * self,float * speechFrame,float * speechFrameHB,float * outFrame,float * outFrameHB)1188 int WebRtcNs_ProcessCore(NSinst_t* self,
1189 float* speechFrame,
1190 float* speechFrameHB,
1191 float* outFrame,
1192 float* outFrameHB) {
1193 // Main routine for noise reduction.
1194 int flagHB = 0;
1195 int i;
1196
1197 float energy1, energy2, gain, factor, factor1, factor2;
1198 float fout[BLOCKL_MAX];
1199 float winData[ANAL_BLOCKL_MAX];
1200 float magn[HALF_ANAL_BLOCKL];
1201 float theFilter[HALF_ANAL_BLOCKL], theFilterTmp[HALF_ANAL_BLOCKL];
1202 float real[ANAL_BLOCKL_MAX], imag[HALF_ANAL_BLOCKL];
1203
1204 // SWB variables.
1205 int deltaBweHB = 1;
1206 int deltaGainHB = 1;
1207 float decayBweHB = 1.0;
1208 float gainMapParHB = 1.0;
1209 float gainTimeDomainHB = 1.0;
1210 float avgProbSpeechHB, avgProbSpeechHBTmp, avgFilterGainHB, gainModHB;
1211 float sumMagnAnalyze, sumMagnProcess;
1212
1213 // Check that initiation has been done.
1214 if (self->initFlag != 1) {
1215 return (-1);
1216 }
1217 // Check for valid pointers based on sampling rate.
1218 if (self->fs == 32000) {
1219 if (speechFrameHB == NULL) {
1220 return -1;
1221 }
1222 flagHB = 1;
1223 // Range for averaging low band quantities for H band gain.
1224 deltaBweHB = (int)self->magnLen / 4;
1225 deltaGainHB = deltaBweHB;
1226 }
1227
1228 // Update analysis buffer for L band.
1229 UpdateBuffer(speechFrame, self->blockLen, self->anaLen, self->dataBuf);
1230
1231 if (flagHB == 1) {
1232 // Update analysis buffer for H band.
1233 UpdateBuffer(speechFrameHB, self->blockLen, self->anaLen, self->dataBufHB);
1234 }
1235
1236 Windowing(self->window, self->dataBuf, self->anaLen, winData);
1237 energy1 = Energy(winData, self->anaLen);
1238 if (energy1 == 0.0) {
1239 // Synthesize the special case of zero input.
1240 // Read out fully processed segment.
1241 for (i = self->windShift; i < self->blockLen + self->windShift; i++) {
1242 fout[i - self->windShift] = self->syntBuf[i];
1243 }
1244 // Update synthesis buffer.
1245 UpdateBuffer(NULL, self->blockLen, self->anaLen, self->syntBuf);
1246
1247 for (i = 0; i < self->blockLen; ++i)
1248 outFrame[i] =
1249 WEBRTC_SPL_SAT(WEBRTC_SPL_WORD16_MAX, fout[i], WEBRTC_SPL_WORD16_MIN);
1250
1251 // For time-domain gain of HB.
1252 if (flagHB == 1)
1253 for (i = 0; i < self->blockLen; ++i)
1254 outFrameHB[i] = WEBRTC_SPL_SAT(
1255 WEBRTC_SPL_WORD16_MAX, self->dataBufHB[i], WEBRTC_SPL_WORD16_MIN);
1256
1257 return 0;
1258 }
1259
1260 FFT(self, winData, self->anaLen, self->magnLen, real, imag, magn);
1261
1262 if (self->blockInd < END_STARTUP_SHORT) {
1263 for (i = 0; i < self->magnLen; i++) {
1264 self->initMagnEst[i] += magn[i];
1265 }
1266 }
1267
1268 ComputeDdBasedWienerFilter(self, magn, theFilter);
1269
1270 for (i = 0; i < self->magnLen; i++) {
1271 // Flooring bottom.
1272 if (theFilter[i] < self->denoiseBound) {
1273 theFilter[i] = self->denoiseBound;
1274 }
1275 // Flooring top.
1276 if (theFilter[i] > 1.f) {
1277 theFilter[i] = 1.f;
1278 }
1279 if (self->blockInd < END_STARTUP_SHORT) {
1280 theFilterTmp[i] =
1281 (self->initMagnEst[i] - self->overdrive * self->parametricNoise[i]);
1282 theFilterTmp[i] /= (self->initMagnEst[i] + 0.0001f);
1283 // Flooring bottom.
1284 if (theFilterTmp[i] < self->denoiseBound) {
1285 theFilterTmp[i] = self->denoiseBound;
1286 }
1287 // Flooring top.
1288 if (theFilterTmp[i] > 1.f) {
1289 theFilterTmp[i] = 1.f;
1290 }
1291 // Weight the two suppression filters.
1292 theFilter[i] *= (self->blockInd);
1293 theFilterTmp[i] *= (END_STARTUP_SHORT - self->blockInd);
1294 theFilter[i] += theFilterTmp[i];
1295 theFilter[i] /= (END_STARTUP_SHORT);
1296 }
1297
1298 self->smooth[i] = theFilter[i];
1299 real[i] *= self->smooth[i];
1300 imag[i] *= self->smooth[i];
1301 }
1302 // Keep track of |magn| spectrum for next frame.
1303 memcpy(self->magnPrevProcess, magn, sizeof(*magn) * self->magnLen);
1304 memcpy(self->noisePrev, self->noise, sizeof(self->noise[0]) * self->magnLen);
1305 // Back to time domain.
1306 IFFT(self, real, imag, self->magnLen, self->anaLen, winData);
1307
1308 // Scale factor: only do it after END_STARTUP_LONG time.
1309 factor = 1.f;
1310 if (self->gainmap == 1 && self->blockInd > END_STARTUP_LONG) {
1311 factor1 = 1.f;
1312 factor2 = 1.f;
1313
1314 energy2 = Energy(winData, self->anaLen);
1315 gain = (float)sqrt(energy2 / (energy1 + 1.f));
1316
1317 // Scaling for new version.
1318 if (gain > B_LIM) {
1319 factor1 = 1.f + 1.3f * (gain - B_LIM);
1320 if (gain * factor1 > 1.f) {
1321 factor1 = 1.f / gain;
1322 }
1323 }
1324 if (gain < B_LIM) {
1325 // Don't reduce scale too much for pause regions:
1326 // attenuation here should be controlled by flooring.
1327 if (gain <= self->denoiseBound) {
1328 gain = self->denoiseBound;
1329 }
1330 factor2 = 1.f - 0.3f * (B_LIM - gain);
1331 }
1332 // Combine both scales with speech/noise prob:
1333 // note prior (priorSpeechProb) is not frequency dependent.
1334 factor = self->priorSpeechProb * factor1 +
1335 (1.f - self->priorSpeechProb) * factor2;
1336 } // Out of self->gainmap == 1.
1337
1338 Windowing(self->window, winData, self->anaLen, winData);
1339
1340 // Synthesis.
1341 for (i = 0; i < self->anaLen; i++) {
1342 self->syntBuf[i] += factor * winData[i];
1343 }
1344 // Read out fully processed segment.
1345 for (i = self->windShift; i < self->blockLen + self->windShift; i++) {
1346 fout[i - self->windShift] = self->syntBuf[i];
1347 }
1348 // Update synthesis buffer.
1349 UpdateBuffer(NULL, self->blockLen, self->anaLen, self->syntBuf);
1350
1351 for (i = 0; i < self->blockLen; ++i)
1352 outFrame[i] =
1353 WEBRTC_SPL_SAT(WEBRTC_SPL_WORD16_MAX, fout[i], WEBRTC_SPL_WORD16_MIN);
1354
1355 // For time-domain gain of HB.
1356 if (flagHB == 1) {
1357 // Average speech prob from low band.
1358 // Average over second half (i.e., 4->8kHz) of frequencies spectrum.
1359 avgProbSpeechHB = 0.0;
1360 for (i = self->magnLen - deltaBweHB - 1; i < self->magnLen - 1; i++) {
1361 avgProbSpeechHB += self->speechProb[i];
1362 }
1363 avgProbSpeechHB = avgProbSpeechHB / ((float)deltaBweHB);
1364 // If the speech was suppressed by a component between Analyze and
1365 // Process, for example the AEC, then it should not be considered speech
1366 // for high band suppression purposes.
1367 sumMagnAnalyze = 0;
1368 sumMagnProcess = 0;
1369 for (i = 0; i < self->magnLen; ++i) {
1370 sumMagnAnalyze += self->magnPrevAnalyze[i];
1371 sumMagnProcess += self->magnPrevProcess[i];
1372 }
1373 avgProbSpeechHB *= sumMagnProcess / sumMagnAnalyze;
1374 // Average filter gain from low band.
1375 // Average over second half (i.e., 4->8kHz) of frequencies spectrum.
1376 avgFilterGainHB = 0.0;
1377 for (i = self->magnLen - deltaGainHB - 1; i < self->magnLen - 1; i++) {
1378 avgFilterGainHB += self->smooth[i];
1379 }
1380 avgFilterGainHB = avgFilterGainHB / ((float)(deltaGainHB));
1381 avgProbSpeechHBTmp = 2.f * avgProbSpeechHB - 1.f;
1382 // Gain based on speech probability.
1383 gainModHB = 0.5f * (1.f + (float)tanh(gainMapParHB * avgProbSpeechHBTmp));
1384 // Combine gain with low band gain.
1385 gainTimeDomainHB = 0.5f * gainModHB + 0.5f * avgFilterGainHB;
1386 if (avgProbSpeechHB >= 0.5f) {
1387 gainTimeDomainHB = 0.25f * gainModHB + 0.75f * avgFilterGainHB;
1388 }
1389 gainTimeDomainHB = gainTimeDomainHB * decayBweHB;
1390 // Make sure gain is within flooring range.
1391 // Flooring bottom.
1392 if (gainTimeDomainHB < self->denoiseBound) {
1393 gainTimeDomainHB = self->denoiseBound;
1394 }
1395 // Flooring top.
1396 if (gainTimeDomainHB > 1.f) {
1397 gainTimeDomainHB = 1.f;
1398 }
1399 // Apply gain.
1400 for (i = 0; i < self->blockLen; i++) {
1401 float o = gainTimeDomainHB * self->dataBufHB[i];
1402 outFrameHB[i] =
1403 WEBRTC_SPL_SAT(WEBRTC_SPL_WORD16_MAX, o, WEBRTC_SPL_WORD16_MIN);
1404 }
1405 } // End of H band gain computation.
1406
1407 return 0;
1408 }
1409