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 /*
12 * The core AEC algorithm, which is presented with time-aligned signals.
13 */
14
15 #include "webrtc/modules/audio_processing/aec/aec_core.h"
16
17 #include <assert.h>
18 #include <math.h>
19 #include <stddef.h> // size_t
20 #include <stdlib.h>
21 #include <string.h>
22
23 #include "webrtc/common_audio/signal_processing/include/signal_processing_library.h"
24 #include "webrtc/modules/audio_processing/aec/aec_core_internal.h"
25 #include "webrtc/modules/audio_processing/aec/aec_rdft.h"
26 #include "webrtc/modules/audio_processing/utility/delay_estimator_wrapper.h"
27 #include "webrtc/modules/audio_processing/utility/ring_buffer.h"
28 #include "webrtc/system_wrappers/interface/cpu_features_wrapper.h"
29 #include "webrtc/typedefs.h"
30
31 // Buffer size (samples)
32 static const size_t kBufSizePartitions = 250; // 1 second of audio in 16 kHz.
33
34 // Metrics
35 static const int subCountLen = 4;
36 static const int countLen = 50;
37
38 // Quantities to control H band scaling for SWB input
39 static const int flagHbandCn = 1; // flag for adding comfort noise in H band
40 static const float cnScaleHband =
41 (float)0.4; // scale for comfort noise in H band
42 // Initial bin for averaging nlp gain in low band
43 static const int freqAvgIc = PART_LEN / 2;
44
45 // Matlab code to produce table:
46 // win = sqrt(hanning(63)); win = [0 ; win(1:32)];
47 // fprintf(1, '\t%.14f, %.14f, %.14f,\n', win);
48 static const float sqrtHanning[65] = {
49 0.00000000000000f, 0.02454122852291f, 0.04906767432742f, 0.07356456359967f,
50 0.09801714032956f, 0.12241067519922f, 0.14673047445536f, 0.17096188876030f,
51 0.19509032201613f, 0.21910124015687f, 0.24298017990326f, 0.26671275747490f,
52 0.29028467725446f, 0.31368174039889f, 0.33688985339222f, 0.35989503653499f,
53 0.38268343236509f, 0.40524131400499f, 0.42755509343028f, 0.44961132965461f,
54 0.47139673682600f, 0.49289819222978f, 0.51410274419322f, 0.53499761988710f,
55 0.55557023301960f, 0.57580819141785f, 0.59569930449243f, 0.61523159058063f,
56 0.63439328416365f, 0.65317284295378f, 0.67155895484702f, 0.68954054473707f,
57 0.70710678118655f, 0.72424708295147f, 0.74095112535496f, 0.75720884650648f,
58 0.77301045336274f, 0.78834642762661f, 0.80320753148064f, 0.81758481315158f,
59 0.83146961230255f, 0.84485356524971f, 0.85772861000027f, 0.87008699110871f,
60 0.88192126434835f, 0.89322430119552f, 0.90398929312344f, 0.91420975570353f,
61 0.92387953251129f, 0.93299279883474f, 0.94154406518302f, 0.94952818059304f,
62 0.95694033573221f, 0.96377606579544f, 0.97003125319454f, 0.97570213003853f,
63 0.98078528040323f, 0.98527764238894f, 0.98917650996478f, 0.99247953459871f,
64 0.99518472667220f, 0.99729045667869f, 0.99879545620517f, 0.99969881869620f,
65 1.00000000000000f};
66
67 // Matlab code to produce table:
68 // weightCurve = [0 ; 0.3 * sqrt(linspace(0,1,64))' + 0.1];
69 // fprintf(1, '\t%.4f, %.4f, %.4f, %.4f, %.4f, %.4f,\n', weightCurve);
70 ALIGN16_BEG const float ALIGN16_END WebRtcAec_weightCurve[65] = {
71 0.0000f, 0.1000f, 0.1378f, 0.1535f, 0.1655f, 0.1756f, 0.1845f, 0.1926f,
72 0.2000f, 0.2069f, 0.2134f, 0.2195f, 0.2254f, 0.2309f, 0.2363f, 0.2414f,
73 0.2464f, 0.2512f, 0.2558f, 0.2604f, 0.2648f, 0.2690f, 0.2732f, 0.2773f,
74 0.2813f, 0.2852f, 0.2890f, 0.2927f, 0.2964f, 0.3000f, 0.3035f, 0.3070f,
75 0.3104f, 0.3138f, 0.3171f, 0.3204f, 0.3236f, 0.3268f, 0.3299f, 0.3330f,
76 0.3360f, 0.3390f, 0.3420f, 0.3449f, 0.3478f, 0.3507f, 0.3535f, 0.3563f,
77 0.3591f, 0.3619f, 0.3646f, 0.3673f, 0.3699f, 0.3726f, 0.3752f, 0.3777f,
78 0.3803f, 0.3828f, 0.3854f, 0.3878f, 0.3903f, 0.3928f, 0.3952f, 0.3976f,
79 0.4000f};
80
81 // Matlab code to produce table:
82 // overDriveCurve = [sqrt(linspace(0,1,65))' + 1];
83 // fprintf(1, '\t%.4f, %.4f, %.4f, %.4f, %.4f, %.4f,\n', overDriveCurve);
84 ALIGN16_BEG const float ALIGN16_END WebRtcAec_overDriveCurve[65] = {
85 1.0000f, 1.1250f, 1.1768f, 1.2165f, 1.2500f, 1.2795f, 1.3062f, 1.3307f,
86 1.3536f, 1.3750f, 1.3953f, 1.4146f, 1.4330f, 1.4507f, 1.4677f, 1.4841f,
87 1.5000f, 1.5154f, 1.5303f, 1.5449f, 1.5590f, 1.5728f, 1.5863f, 1.5995f,
88 1.6124f, 1.6250f, 1.6374f, 1.6495f, 1.6614f, 1.6731f, 1.6847f, 1.6960f,
89 1.7071f, 1.7181f, 1.7289f, 1.7395f, 1.7500f, 1.7603f, 1.7706f, 1.7806f,
90 1.7906f, 1.8004f, 1.8101f, 1.8197f, 1.8292f, 1.8385f, 1.8478f, 1.8570f,
91 1.8660f, 1.8750f, 1.8839f, 1.8927f, 1.9014f, 1.9100f, 1.9186f, 1.9270f,
92 1.9354f, 1.9437f, 1.9520f, 1.9601f, 1.9682f, 1.9763f, 1.9843f, 1.9922f,
93 2.0000f};
94
95 // Target suppression levels for nlp modes.
96 // log{0.001, 0.00001, 0.00000001}
97 static const float kTargetSupp[3] = {-6.9f, -11.5f, -18.4f};
98
99 // Two sets of parameters, one for the extended filter mode.
100 static const float kExtendedMinOverDrive[3] = {3.0f, 6.0f, 15.0f};
101 static const float kNormalMinOverDrive[3] = {1.0f, 2.0f, 5.0f};
102 static const float kExtendedSmoothingCoefficients[2][2] = {{0.9f, 0.1f},
103 {0.92f, 0.08f}};
104 static const float kNormalSmoothingCoefficients[2][2] = {{0.9f, 0.1f},
105 {0.93f, 0.07f}};
106
107 // Number of partitions forming the NLP's "preferred" bands.
108 enum {
109 kPrefBandSize = 24
110 };
111
112 #ifdef WEBRTC_AEC_DEBUG_DUMP
113 extern int webrtc_aec_instance_count;
114 #endif
115
116 // "Private" function prototypes.
117 static void ProcessBlock(AecCore* aec);
118
119 static void NonLinearProcessing(AecCore* aec, float* output, float* outputH);
120
121 static void GetHighbandGain(const float* lambda, float* nlpGainHband);
122
123 // Comfort_noise also computes noise for H band returned in comfortNoiseHband
124 static void ComfortNoise(AecCore* aec,
125 float efw[2][PART_LEN1],
126 complex_t* comfortNoiseHband,
127 const float* noisePow,
128 const float* lambda);
129
130 static void InitLevel(PowerLevel* level);
131 static void InitStats(Stats* stats);
132 static void InitMetrics(AecCore* aec);
133 static void UpdateLevel(PowerLevel* level, float in[2][PART_LEN1]);
134 static void UpdateMetrics(AecCore* aec);
135 // Convert from time domain to frequency domain. Note that |time_data| are
136 // overwritten.
137 static void TimeToFrequency(float time_data[PART_LEN2],
138 float freq_data[2][PART_LEN1],
139 int window);
140
MulRe(float aRe,float aIm,float bRe,float bIm)141 __inline static float MulRe(float aRe, float aIm, float bRe, float bIm) {
142 return aRe * bRe - aIm * bIm;
143 }
144
MulIm(float aRe,float aIm,float bRe,float bIm)145 __inline static float MulIm(float aRe, float aIm, float bRe, float bIm) {
146 return aRe * bIm + aIm * bRe;
147 }
148
CmpFloat(const void * a,const void * b)149 static int CmpFloat(const void* a, const void* b) {
150 const float* da = (const float*)a;
151 const float* db = (const float*)b;
152
153 return (*da > *db) - (*da < *db);
154 }
155
WebRtcAec_CreateAec(AecCore ** aecInst)156 int WebRtcAec_CreateAec(AecCore** aecInst) {
157 AecCore* aec = malloc(sizeof(AecCore));
158 *aecInst = aec;
159 if (aec == NULL) {
160 return -1;
161 }
162
163 aec->nearFrBuf = WebRtc_CreateBuffer(FRAME_LEN + PART_LEN, sizeof(float));
164 if (!aec->nearFrBuf) {
165 WebRtcAec_FreeAec(aec);
166 aec = NULL;
167 return -1;
168 }
169
170 aec->outFrBuf = WebRtc_CreateBuffer(FRAME_LEN + PART_LEN, sizeof(float));
171 if (!aec->outFrBuf) {
172 WebRtcAec_FreeAec(aec);
173 aec = NULL;
174 return -1;
175 }
176
177 aec->nearFrBufH = WebRtc_CreateBuffer(FRAME_LEN + PART_LEN, sizeof(float));
178 if (!aec->nearFrBufH) {
179 WebRtcAec_FreeAec(aec);
180 aec = NULL;
181 return -1;
182 }
183
184 aec->outFrBufH = WebRtc_CreateBuffer(FRAME_LEN + PART_LEN, sizeof(float));
185 if (!aec->outFrBufH) {
186 WebRtcAec_FreeAec(aec);
187 aec = NULL;
188 return -1;
189 }
190
191 // Create far-end buffers.
192 aec->far_buf =
193 WebRtc_CreateBuffer(kBufSizePartitions, sizeof(float) * 2 * PART_LEN1);
194 if (!aec->far_buf) {
195 WebRtcAec_FreeAec(aec);
196 aec = NULL;
197 return -1;
198 }
199 aec->far_buf_windowed =
200 WebRtc_CreateBuffer(kBufSizePartitions, sizeof(float) * 2 * PART_LEN1);
201 if (!aec->far_buf_windowed) {
202 WebRtcAec_FreeAec(aec);
203 aec = NULL;
204 return -1;
205 }
206 #ifdef WEBRTC_AEC_DEBUG_DUMP
207 aec->far_time_buf =
208 WebRtc_CreateBuffer(kBufSizePartitions, sizeof(int16_t) * PART_LEN);
209 if (!aec->far_time_buf) {
210 WebRtcAec_FreeAec(aec);
211 aec = NULL;
212 return -1;
213 }
214 {
215 char filename[64];
216 sprintf(filename, "aec_far%d.pcm", webrtc_aec_instance_count);
217 aec->farFile = fopen(filename, "wb");
218 sprintf(filename, "aec_near%d.pcm", webrtc_aec_instance_count);
219 aec->nearFile = fopen(filename, "wb");
220 sprintf(filename, "aec_out%d.pcm", webrtc_aec_instance_count);
221 aec->outFile = fopen(filename, "wb");
222 sprintf(filename, "aec_out_linear%d.pcm", webrtc_aec_instance_count);
223 aec->outLinearFile = fopen(filename, "wb");
224 }
225 #endif
226 aec->delay_estimator_farend =
227 WebRtc_CreateDelayEstimatorFarend(PART_LEN1, kHistorySizeBlocks);
228 if (aec->delay_estimator_farend == NULL) {
229 WebRtcAec_FreeAec(aec);
230 aec = NULL;
231 return -1;
232 }
233 aec->delay_estimator = WebRtc_CreateDelayEstimator(
234 aec->delay_estimator_farend, kLookaheadBlocks);
235 if (aec->delay_estimator == NULL) {
236 WebRtcAec_FreeAec(aec);
237 aec = NULL;
238 return -1;
239 }
240
241 return 0;
242 }
243
WebRtcAec_FreeAec(AecCore * aec)244 int WebRtcAec_FreeAec(AecCore* aec) {
245 if (aec == NULL) {
246 return -1;
247 }
248
249 WebRtc_FreeBuffer(aec->nearFrBuf);
250 WebRtc_FreeBuffer(aec->outFrBuf);
251
252 WebRtc_FreeBuffer(aec->nearFrBufH);
253 WebRtc_FreeBuffer(aec->outFrBufH);
254
255 WebRtc_FreeBuffer(aec->far_buf);
256 WebRtc_FreeBuffer(aec->far_buf_windowed);
257 #ifdef WEBRTC_AEC_DEBUG_DUMP
258 WebRtc_FreeBuffer(aec->far_time_buf);
259 fclose(aec->farFile);
260 fclose(aec->nearFile);
261 fclose(aec->outFile);
262 fclose(aec->outLinearFile);
263 #endif
264 WebRtc_FreeDelayEstimator(aec->delay_estimator);
265 WebRtc_FreeDelayEstimatorFarend(aec->delay_estimator_farend);
266
267 free(aec);
268 return 0;
269 }
270
FilterFar(AecCore * aec,float yf[2][PART_LEN1])271 static void FilterFar(AecCore* aec, float yf[2][PART_LEN1]) {
272 int i;
273 for (i = 0; i < aec->num_partitions; i++) {
274 int j;
275 int xPos = (i + aec->xfBufBlockPos) * PART_LEN1;
276 int pos = i * PART_LEN1;
277 // Check for wrap
278 if (i + aec->xfBufBlockPos >= aec->num_partitions) {
279 xPos -= aec->num_partitions * (PART_LEN1);
280 }
281
282 for (j = 0; j < PART_LEN1; j++) {
283 yf[0][j] += MulRe(aec->xfBuf[0][xPos + j],
284 aec->xfBuf[1][xPos + j],
285 aec->wfBuf[0][pos + j],
286 aec->wfBuf[1][pos + j]);
287 yf[1][j] += MulIm(aec->xfBuf[0][xPos + j],
288 aec->xfBuf[1][xPos + j],
289 aec->wfBuf[0][pos + j],
290 aec->wfBuf[1][pos + j]);
291 }
292 }
293 }
294
ScaleErrorSignal(AecCore * aec,float ef[2][PART_LEN1])295 static void ScaleErrorSignal(AecCore* aec, float ef[2][PART_LEN1]) {
296 const float mu = aec->extended_filter_enabled ? kExtendedMu : aec->normal_mu;
297 const float error_threshold = aec->extended_filter_enabled
298 ? kExtendedErrorThreshold
299 : aec->normal_error_threshold;
300 int i;
301 float abs_ef;
302 for (i = 0; i < (PART_LEN1); i++) {
303 ef[0][i] /= (aec->xPow[i] + 1e-10f);
304 ef[1][i] /= (aec->xPow[i] + 1e-10f);
305 abs_ef = sqrtf(ef[0][i] * ef[0][i] + ef[1][i] * ef[1][i]);
306
307 if (abs_ef > error_threshold) {
308 abs_ef = error_threshold / (abs_ef + 1e-10f);
309 ef[0][i] *= abs_ef;
310 ef[1][i] *= abs_ef;
311 }
312
313 // Stepsize factor
314 ef[0][i] *= mu;
315 ef[1][i] *= mu;
316 }
317 }
318
319 // Time-unconstrined filter adaptation.
320 // TODO(andrew): consider for a low-complexity mode.
321 // static void FilterAdaptationUnconstrained(AecCore* aec, float *fft,
322 // float ef[2][PART_LEN1]) {
323 // int i, j;
324 // for (i = 0; i < aec->num_partitions; i++) {
325 // int xPos = (i + aec->xfBufBlockPos)*(PART_LEN1);
326 // int pos;
327 // // Check for wrap
328 // if (i + aec->xfBufBlockPos >= aec->num_partitions) {
329 // xPos -= aec->num_partitions * PART_LEN1;
330 // }
331 //
332 // pos = i * PART_LEN1;
333 //
334 // for (j = 0; j < PART_LEN1; j++) {
335 // aec->wfBuf[0][pos + j] += MulRe(aec->xfBuf[0][xPos + j],
336 // -aec->xfBuf[1][xPos + j],
337 // ef[0][j], ef[1][j]);
338 // aec->wfBuf[1][pos + j] += MulIm(aec->xfBuf[0][xPos + j],
339 // -aec->xfBuf[1][xPos + j],
340 // ef[0][j], ef[1][j]);
341 // }
342 // }
343 //}
344
FilterAdaptation(AecCore * aec,float * fft,float ef[2][PART_LEN1])345 static void FilterAdaptation(AecCore* aec, float* fft, float ef[2][PART_LEN1]) {
346 int i, j;
347 for (i = 0; i < aec->num_partitions; i++) {
348 int xPos = (i + aec->xfBufBlockPos) * (PART_LEN1);
349 int pos;
350 // Check for wrap
351 if (i + aec->xfBufBlockPos >= aec->num_partitions) {
352 xPos -= aec->num_partitions * PART_LEN1;
353 }
354
355 pos = i * PART_LEN1;
356
357 for (j = 0; j < PART_LEN; j++) {
358
359 fft[2 * j] = MulRe(aec->xfBuf[0][xPos + j],
360 -aec->xfBuf[1][xPos + j],
361 ef[0][j],
362 ef[1][j]);
363 fft[2 * j + 1] = MulIm(aec->xfBuf[0][xPos + j],
364 -aec->xfBuf[1][xPos + j],
365 ef[0][j],
366 ef[1][j]);
367 }
368 fft[1] = MulRe(aec->xfBuf[0][xPos + PART_LEN],
369 -aec->xfBuf[1][xPos + PART_LEN],
370 ef[0][PART_LEN],
371 ef[1][PART_LEN]);
372
373 aec_rdft_inverse_128(fft);
374 memset(fft + PART_LEN, 0, sizeof(float) * PART_LEN);
375
376 // fft scaling
377 {
378 float scale = 2.0f / PART_LEN2;
379 for (j = 0; j < PART_LEN; j++) {
380 fft[j] *= scale;
381 }
382 }
383 aec_rdft_forward_128(fft);
384
385 aec->wfBuf[0][pos] += fft[0];
386 aec->wfBuf[0][pos + PART_LEN] += fft[1];
387
388 for (j = 1; j < PART_LEN; j++) {
389 aec->wfBuf[0][pos + j] += fft[2 * j];
390 aec->wfBuf[1][pos + j] += fft[2 * j + 1];
391 }
392 }
393 }
394
OverdriveAndSuppress(AecCore * aec,float hNl[PART_LEN1],const float hNlFb,float efw[2][PART_LEN1])395 static void OverdriveAndSuppress(AecCore* aec,
396 float hNl[PART_LEN1],
397 const float hNlFb,
398 float efw[2][PART_LEN1]) {
399 int i;
400 for (i = 0; i < PART_LEN1; i++) {
401 // Weight subbands
402 if (hNl[i] > hNlFb) {
403 hNl[i] = WebRtcAec_weightCurve[i] * hNlFb +
404 (1 - WebRtcAec_weightCurve[i]) * hNl[i];
405 }
406 hNl[i] = powf(hNl[i], aec->overDriveSm * WebRtcAec_overDriveCurve[i]);
407
408 // Suppress error signal
409 efw[0][i] *= hNl[i];
410 efw[1][i] *= hNl[i];
411
412 // Ooura fft returns incorrect sign on imaginary component. It matters here
413 // because we are making an additive change with comfort noise.
414 efw[1][i] *= -1;
415 }
416 }
417
418 WebRtcAec_FilterFar_t WebRtcAec_FilterFar;
419 WebRtcAec_ScaleErrorSignal_t WebRtcAec_ScaleErrorSignal;
420 WebRtcAec_FilterAdaptation_t WebRtcAec_FilterAdaptation;
421 WebRtcAec_OverdriveAndSuppress_t WebRtcAec_OverdriveAndSuppress;
422 WebRtcAec_ComfortNoise_t WebRtcAec_ComfortNoise;
423
WebRtcAec_InitAec(AecCore * aec,int sampFreq)424 int WebRtcAec_InitAec(AecCore* aec, int sampFreq) {
425 int i;
426
427 aec->sampFreq = sampFreq;
428
429 if (sampFreq == 8000) {
430 aec->normal_mu = 0.6f;
431 aec->normal_error_threshold = 2e-6f;
432 } else {
433 aec->normal_mu = 0.5f;
434 aec->normal_error_threshold = 1.5e-6f;
435 }
436
437 if (WebRtc_InitBuffer(aec->nearFrBuf) == -1) {
438 return -1;
439 }
440
441 if (WebRtc_InitBuffer(aec->outFrBuf) == -1) {
442 return -1;
443 }
444
445 if (WebRtc_InitBuffer(aec->nearFrBufH) == -1) {
446 return -1;
447 }
448
449 if (WebRtc_InitBuffer(aec->outFrBufH) == -1) {
450 return -1;
451 }
452
453 // Initialize far-end buffers.
454 if (WebRtc_InitBuffer(aec->far_buf) == -1) {
455 return -1;
456 }
457 if (WebRtc_InitBuffer(aec->far_buf_windowed) == -1) {
458 return -1;
459 }
460 #ifdef WEBRTC_AEC_DEBUG_DUMP
461 if (WebRtc_InitBuffer(aec->far_time_buf) == -1) {
462 return -1;
463 }
464 #endif
465 aec->system_delay = 0;
466
467 if (WebRtc_InitDelayEstimatorFarend(aec->delay_estimator_farend) != 0) {
468 return -1;
469 }
470 if (WebRtc_InitDelayEstimator(aec->delay_estimator) != 0) {
471 return -1;
472 }
473 aec->delay_logging_enabled = 0;
474 memset(aec->delay_histogram, 0, sizeof(aec->delay_histogram));
475
476 aec->reported_delay_enabled = 1;
477 aec->extended_filter_enabled = 0;
478 aec->num_partitions = kNormalNumPartitions;
479
480 // Update the delay estimator with filter length. We use half the
481 // |num_partitions| to take the echo path into account. In practice we say
482 // that the echo has a duration of maximum half |num_partitions|, which is not
483 // true, but serves as a crude measure.
484 WebRtc_set_allowed_offset(aec->delay_estimator, aec->num_partitions / 2);
485 // TODO(bjornv): I currently hard coded the enable. Once we've established
486 // that AECM has no performance regression, robust_validation will be enabled
487 // all the time and the APIs to turn it on/off will be removed. Hence, remove
488 // this line then.
489 WebRtc_enable_robust_validation(aec->delay_estimator, 1);
490
491 // Default target suppression mode.
492 aec->nlp_mode = 1;
493
494 // Sampling frequency multiplier
495 // SWB is processed as 160 frame size
496 if (aec->sampFreq == 32000) {
497 aec->mult = (short)aec->sampFreq / 16000;
498 } else {
499 aec->mult = (short)aec->sampFreq / 8000;
500 }
501
502 aec->farBufWritePos = 0;
503 aec->farBufReadPos = 0;
504
505 aec->inSamples = 0;
506 aec->outSamples = 0;
507 aec->knownDelay = 0;
508
509 // Initialize buffers
510 memset(aec->dBuf, 0, sizeof(aec->dBuf));
511 memset(aec->eBuf, 0, sizeof(aec->eBuf));
512 // For H band
513 memset(aec->dBufH, 0, sizeof(aec->dBufH));
514
515 memset(aec->xPow, 0, sizeof(aec->xPow));
516 memset(aec->dPow, 0, sizeof(aec->dPow));
517 memset(aec->dInitMinPow, 0, sizeof(aec->dInitMinPow));
518 aec->noisePow = aec->dInitMinPow;
519 aec->noiseEstCtr = 0;
520
521 // Initial comfort noise power
522 for (i = 0; i < PART_LEN1; i++) {
523 aec->dMinPow[i] = 1.0e6f;
524 }
525
526 // Holds the last block written to
527 aec->xfBufBlockPos = 0;
528 // TODO: Investigate need for these initializations. Deleting them doesn't
529 // change the output at all and yields 0.4% overall speedup.
530 memset(aec->xfBuf, 0, sizeof(complex_t) * kExtendedNumPartitions * PART_LEN1);
531 memset(aec->wfBuf, 0, sizeof(complex_t) * kExtendedNumPartitions * PART_LEN1);
532 memset(aec->sde, 0, sizeof(complex_t) * PART_LEN1);
533 memset(aec->sxd, 0, sizeof(complex_t) * PART_LEN1);
534 memset(
535 aec->xfwBuf, 0, sizeof(complex_t) * kExtendedNumPartitions * PART_LEN1);
536 memset(aec->se, 0, sizeof(float) * PART_LEN1);
537
538 // To prevent numerical instability in the first block.
539 for (i = 0; i < PART_LEN1; i++) {
540 aec->sd[i] = 1;
541 }
542 for (i = 0; i < PART_LEN1; i++) {
543 aec->sx[i] = 1;
544 }
545
546 memset(aec->hNs, 0, sizeof(aec->hNs));
547 memset(aec->outBuf, 0, sizeof(float) * PART_LEN);
548
549 aec->hNlFbMin = 1;
550 aec->hNlFbLocalMin = 1;
551 aec->hNlXdAvgMin = 1;
552 aec->hNlNewMin = 0;
553 aec->hNlMinCtr = 0;
554 aec->overDrive = 2;
555 aec->overDriveSm = 2;
556 aec->delayIdx = 0;
557 aec->stNearState = 0;
558 aec->echoState = 0;
559 aec->divergeState = 0;
560
561 aec->seed = 777;
562 aec->delayEstCtr = 0;
563
564 // Metrics disabled by default
565 aec->metricsMode = 0;
566 InitMetrics(aec);
567
568 // Assembly optimization
569 WebRtcAec_FilterFar = FilterFar;
570 WebRtcAec_ScaleErrorSignal = ScaleErrorSignal;
571 WebRtcAec_FilterAdaptation = FilterAdaptation;
572 WebRtcAec_OverdriveAndSuppress = OverdriveAndSuppress;
573 WebRtcAec_ComfortNoise = ComfortNoise;
574
575 #if defined(WEBRTC_ARCH_X86_FAMILY)
576 if (WebRtc_GetCPUInfo(kSSE2)) {
577 WebRtcAec_InitAec_SSE2();
578 }
579 #endif
580
581 #if defined(MIPS_FPU_LE)
582 WebRtcAec_InitAec_mips();
583 #endif
584
585 #if defined(WEBRTC_DETECT_ARM_NEON) || defined(WEBRTC_ARCH_ARM_NEON)
586 WebRtcAec_InitAec_neon();
587 #endif
588
589 aec_rdft_init();
590
591 return 0;
592 }
593
WebRtcAec_BufferFarendPartition(AecCore * aec,const float * farend)594 void WebRtcAec_BufferFarendPartition(AecCore* aec, const float* farend) {
595 float fft[PART_LEN2];
596 float xf[2][PART_LEN1];
597
598 // Check if the buffer is full, and in that case flush the oldest data.
599 if (WebRtc_available_write(aec->far_buf) < 1) {
600 WebRtcAec_MoveFarReadPtr(aec, 1);
601 }
602 // Convert far-end partition to the frequency domain without windowing.
603 memcpy(fft, farend, sizeof(float) * PART_LEN2);
604 TimeToFrequency(fft, xf, 0);
605 WebRtc_WriteBuffer(aec->far_buf, &xf[0][0], 1);
606
607 // Convert far-end partition to the frequency domain with windowing.
608 memcpy(fft, farend, sizeof(float) * PART_LEN2);
609 TimeToFrequency(fft, xf, 1);
610 WebRtc_WriteBuffer(aec->far_buf_windowed, &xf[0][0], 1);
611 }
612
WebRtcAec_MoveFarReadPtr(AecCore * aec,int elements)613 int WebRtcAec_MoveFarReadPtr(AecCore* aec, int elements) {
614 int elements_moved = WebRtc_MoveReadPtr(aec->far_buf_windowed, elements);
615 WebRtc_MoveReadPtr(aec->far_buf, elements);
616 #ifdef WEBRTC_AEC_DEBUG_DUMP
617 WebRtc_MoveReadPtr(aec->far_time_buf, elements);
618 #endif
619 aec->system_delay -= elements_moved * PART_LEN;
620 return elements_moved;
621 }
622
WebRtcAec_ProcessFrame(AecCore * aec,const float * nearend,const float * nearendH,int knownDelay,float * out,float * outH)623 void WebRtcAec_ProcessFrame(AecCore* aec,
624 const float* nearend,
625 const float* nearendH,
626 int knownDelay,
627 float* out,
628 float* outH) {
629 int out_elements = 0;
630
631 // For each frame the process is as follows:
632 // 1) If the system_delay indicates on being too small for processing a
633 // frame we stuff the buffer with enough data for 10 ms.
634 // 2) Adjust the buffer to the system delay, by moving the read pointer.
635 // 3) TODO(bjornv): Investigate if we need to add this:
636 // If we can't move read pointer due to buffer size limitations we
637 // flush/stuff the buffer.
638 // 4) Process as many partitions as possible.
639 // 5) Update the |system_delay| with respect to a full frame of FRAME_LEN
640 // samples. Even though we will have data left to process (we work with
641 // partitions) we consider updating a whole frame, since that's the
642 // amount of data we input and output in audio_processing.
643 // 6) Update the outputs.
644
645 // TODO(bjornv): Investigate how we should round the delay difference; right
646 // now we know that incoming |knownDelay| is underestimated when it's less
647 // than |aec->knownDelay|. We therefore, round (-32) in that direction. In
648 // the other direction, we don't have this situation, but might flush one
649 // partition too little. This can cause non-causality, which should be
650 // investigated. Maybe, allow for a non-symmetric rounding, like -16.
651 int move_elements = (aec->knownDelay - knownDelay - 32) / PART_LEN;
652 int moved_elements = 0;
653
654 // TODO(bjornv): Change the near-end buffer handling to be the same as for
655 // far-end, that is, with a near_pre_buf.
656 // Buffer the near-end frame.
657 WebRtc_WriteBuffer(aec->nearFrBuf, nearend, FRAME_LEN);
658 // For H band
659 if (aec->sampFreq == 32000) {
660 WebRtc_WriteBuffer(aec->nearFrBufH, nearendH, FRAME_LEN);
661 }
662
663 // 1) At most we process |aec->mult|+1 partitions in 10 ms. Make sure we
664 // have enough far-end data for that by stuffing the buffer if the
665 // |system_delay| indicates others.
666 if (aec->system_delay < FRAME_LEN) {
667 // We don't have enough data so we rewind 10 ms.
668 WebRtcAec_MoveFarReadPtr(aec, -(aec->mult + 1));
669 }
670
671 // 2) Compensate for a possible change in the system delay.
672 WebRtc_MoveReadPtr(aec->far_buf_windowed, move_elements);
673 moved_elements = WebRtc_MoveReadPtr(aec->far_buf, move_elements);
674 aec->knownDelay -= moved_elements * PART_LEN;
675 #ifdef WEBRTC_AEC_DEBUG_DUMP
676 WebRtc_MoveReadPtr(aec->far_time_buf, move_elements);
677 #endif
678
679 // 4) Process as many blocks as possible.
680 while (WebRtc_available_read(aec->nearFrBuf) >= PART_LEN) {
681 ProcessBlock(aec);
682 }
683
684 // 5) Update system delay with respect to the entire frame.
685 aec->system_delay -= FRAME_LEN;
686
687 // 6) Update output frame.
688 // Stuff the out buffer if we have less than a frame to output.
689 // This should only happen for the first frame.
690 out_elements = (int)WebRtc_available_read(aec->outFrBuf);
691 if (out_elements < FRAME_LEN) {
692 WebRtc_MoveReadPtr(aec->outFrBuf, out_elements - FRAME_LEN);
693 if (aec->sampFreq == 32000) {
694 WebRtc_MoveReadPtr(aec->outFrBufH, out_elements - FRAME_LEN);
695 }
696 }
697 // Obtain an output frame.
698 WebRtc_ReadBuffer(aec->outFrBuf, NULL, out, FRAME_LEN);
699 // For H band.
700 if (aec->sampFreq == 32000) {
701 WebRtc_ReadBuffer(aec->outFrBufH, NULL, outH, FRAME_LEN);
702 }
703 }
704
WebRtcAec_GetDelayMetricsCore(AecCore * self,int * median,int * std)705 int WebRtcAec_GetDelayMetricsCore(AecCore* self, int* median, int* std) {
706 int i = 0;
707 int delay_values = 0;
708 int num_delay_values = 0;
709 int my_median = 0;
710 const int kMsPerBlock = PART_LEN / (self->mult * 8);
711 float l1_norm = 0;
712
713 assert(self != NULL);
714 assert(median != NULL);
715 assert(std != NULL);
716
717 if (self->delay_logging_enabled == 0) {
718 // Logging disabled.
719 return -1;
720 }
721
722 // Get number of delay values since last update.
723 for (i = 0; i < kHistorySizeBlocks; i++) {
724 num_delay_values += self->delay_histogram[i];
725 }
726 if (num_delay_values == 0) {
727 // We have no new delay value data. Even though -1 is a valid estimate, it
728 // will practically never be used since multiples of |kMsPerBlock| will
729 // always be returned.
730 *median = -1;
731 *std = -1;
732 return 0;
733 }
734
735 delay_values = num_delay_values >> 1; // Start value for median count down.
736 // Get median of delay values since last update.
737 for (i = 0; i < kHistorySizeBlocks; i++) {
738 delay_values -= self->delay_histogram[i];
739 if (delay_values < 0) {
740 my_median = i;
741 break;
742 }
743 }
744 // Account for lookahead.
745 *median = (my_median - kLookaheadBlocks) * kMsPerBlock;
746
747 // Calculate the L1 norm, with median value as central moment.
748 for (i = 0; i < kHistorySizeBlocks; i++) {
749 l1_norm += (float)abs(i - my_median) * self->delay_histogram[i];
750 }
751 *std = (int)(l1_norm / (float)num_delay_values + 0.5f) * kMsPerBlock;
752
753 // Reset histogram.
754 memset(self->delay_histogram, 0, sizeof(self->delay_histogram));
755
756 return 0;
757 }
758
WebRtcAec_echo_state(AecCore * self)759 int WebRtcAec_echo_state(AecCore* self) { return self->echoState; }
760
WebRtcAec_GetEchoStats(AecCore * self,Stats * erl,Stats * erle,Stats * a_nlp)761 void WebRtcAec_GetEchoStats(AecCore* self,
762 Stats* erl,
763 Stats* erle,
764 Stats* a_nlp) {
765 assert(erl != NULL);
766 assert(erle != NULL);
767 assert(a_nlp != NULL);
768 *erl = self->erl;
769 *erle = self->erle;
770 *a_nlp = self->aNlp;
771 }
772
773 #ifdef WEBRTC_AEC_DEBUG_DUMP
WebRtcAec_far_time_buf(AecCore * self)774 void* WebRtcAec_far_time_buf(AecCore* self) { return self->far_time_buf; }
775 #endif
776
WebRtcAec_SetConfigCore(AecCore * self,int nlp_mode,int metrics_mode,int delay_logging)777 void WebRtcAec_SetConfigCore(AecCore* self,
778 int nlp_mode,
779 int metrics_mode,
780 int delay_logging) {
781 assert(nlp_mode >= 0 && nlp_mode < 3);
782 self->nlp_mode = nlp_mode;
783 self->metricsMode = metrics_mode;
784 if (self->metricsMode) {
785 InitMetrics(self);
786 }
787 self->delay_logging_enabled = delay_logging;
788 if (self->delay_logging_enabled) {
789 memset(self->delay_histogram, 0, sizeof(self->delay_histogram));
790 }
791 }
792
WebRtcAec_enable_reported_delay(AecCore * self,int enable)793 void WebRtcAec_enable_reported_delay(AecCore* self, int enable) {
794 self->reported_delay_enabled = enable;
795 }
796
WebRtcAec_reported_delay_enabled(AecCore * self)797 int WebRtcAec_reported_delay_enabled(AecCore* self) {
798 return self->reported_delay_enabled;
799 }
800
WebRtcAec_enable_delay_correction(AecCore * self,int enable)801 void WebRtcAec_enable_delay_correction(AecCore* self, int enable) {
802 self->extended_filter_enabled = enable;
803 self->num_partitions = enable ? kExtendedNumPartitions : kNormalNumPartitions;
804 // Update the delay estimator with filter length. See InitAEC() for details.
805 WebRtc_set_allowed_offset(self->delay_estimator, self->num_partitions / 2);
806 }
807
WebRtcAec_delay_correction_enabled(AecCore * self)808 int WebRtcAec_delay_correction_enabled(AecCore* self) {
809 return self->extended_filter_enabled;
810 }
811
WebRtcAec_system_delay(AecCore * self)812 int WebRtcAec_system_delay(AecCore* self) { return self->system_delay; }
813
WebRtcAec_SetSystemDelay(AecCore * self,int delay)814 void WebRtcAec_SetSystemDelay(AecCore* self, int delay) {
815 assert(delay >= 0);
816 self->system_delay = delay;
817 }
818
ProcessBlock(AecCore * aec)819 static void ProcessBlock(AecCore* aec) {
820 int i;
821 float y[PART_LEN], e[PART_LEN];
822 float scale;
823
824 float fft[PART_LEN2];
825 float xf[2][PART_LEN1], yf[2][PART_LEN1], ef[2][PART_LEN1];
826 float df[2][PART_LEN1];
827 float far_spectrum = 0.0f;
828 float near_spectrum = 0.0f;
829 float abs_far_spectrum[PART_LEN1];
830 float abs_near_spectrum[PART_LEN1];
831
832 const float gPow[2] = {0.9f, 0.1f};
833
834 // Noise estimate constants.
835 const int noiseInitBlocks = 500 * aec->mult;
836 const float step = 0.1f;
837 const float ramp = 1.0002f;
838 const float gInitNoise[2] = {0.999f, 0.001f};
839
840 float nearend[PART_LEN];
841 float* nearend_ptr = NULL;
842 float output[PART_LEN];
843 float outputH[PART_LEN];
844
845 float* xf_ptr = NULL;
846
847 // Concatenate old and new nearend blocks.
848 if (aec->sampFreq == 32000) {
849 WebRtc_ReadBuffer(aec->nearFrBufH, (void**)&nearend_ptr, nearend, PART_LEN);
850 memcpy(aec->dBufH + PART_LEN, nearend_ptr, sizeof(nearend));
851 }
852 WebRtc_ReadBuffer(aec->nearFrBuf, (void**)&nearend_ptr, nearend, PART_LEN);
853 memcpy(aec->dBuf + PART_LEN, nearend_ptr, sizeof(nearend));
854
855 // ---------- Ooura fft ----------
856
857 #ifdef WEBRTC_AEC_DEBUG_DUMP
858 {
859 int16_t farend[PART_LEN];
860 int16_t* farend_ptr = NULL;
861 WebRtc_ReadBuffer(aec->far_time_buf, (void**)&farend_ptr, farend, 1);
862 (void)fwrite(farend_ptr, sizeof(int16_t), PART_LEN, aec->farFile);
863 (void)fwrite(nearend_ptr, sizeof(int16_t), PART_LEN, aec->nearFile);
864 }
865 #endif
866
867 // We should always have at least one element stored in |far_buf|.
868 assert(WebRtc_available_read(aec->far_buf) > 0);
869 WebRtc_ReadBuffer(aec->far_buf, (void**)&xf_ptr, &xf[0][0], 1);
870
871 // Near fft
872 memcpy(fft, aec->dBuf, sizeof(float) * PART_LEN2);
873 TimeToFrequency(fft, df, 0);
874
875 // Power smoothing
876 for (i = 0; i < PART_LEN1; i++) {
877 far_spectrum = (xf_ptr[i] * xf_ptr[i]) +
878 (xf_ptr[PART_LEN1 + i] * xf_ptr[PART_LEN1 + i]);
879 aec->xPow[i] =
880 gPow[0] * aec->xPow[i] + gPow[1] * aec->num_partitions * far_spectrum;
881 // Calculate absolute spectra
882 abs_far_spectrum[i] = sqrtf(far_spectrum);
883
884 near_spectrum = df[0][i] * df[0][i] + df[1][i] * df[1][i];
885 aec->dPow[i] = gPow[0] * aec->dPow[i] + gPow[1] * near_spectrum;
886 // Calculate absolute spectra
887 abs_near_spectrum[i] = sqrtf(near_spectrum);
888 }
889
890 // Estimate noise power. Wait until dPow is more stable.
891 if (aec->noiseEstCtr > 50) {
892 for (i = 0; i < PART_LEN1; i++) {
893 if (aec->dPow[i] < aec->dMinPow[i]) {
894 aec->dMinPow[i] =
895 (aec->dPow[i] + step * (aec->dMinPow[i] - aec->dPow[i])) * ramp;
896 } else {
897 aec->dMinPow[i] *= ramp;
898 }
899 }
900 }
901
902 // Smooth increasing noise power from zero at the start,
903 // to avoid a sudden burst of comfort noise.
904 if (aec->noiseEstCtr < noiseInitBlocks) {
905 aec->noiseEstCtr++;
906 for (i = 0; i < PART_LEN1; i++) {
907 if (aec->dMinPow[i] > aec->dInitMinPow[i]) {
908 aec->dInitMinPow[i] = gInitNoise[0] * aec->dInitMinPow[i] +
909 gInitNoise[1] * aec->dMinPow[i];
910 } else {
911 aec->dInitMinPow[i] = aec->dMinPow[i];
912 }
913 }
914 aec->noisePow = aec->dInitMinPow;
915 } else {
916 aec->noisePow = aec->dMinPow;
917 }
918
919 // Block wise delay estimation used for logging
920 if (aec->delay_logging_enabled) {
921 int delay_estimate = 0;
922 if (WebRtc_AddFarSpectrumFloat(
923 aec->delay_estimator_farend, abs_far_spectrum, PART_LEN1) == 0) {
924 delay_estimate = WebRtc_DelayEstimatorProcessFloat(
925 aec->delay_estimator, abs_near_spectrum, PART_LEN1);
926 if (delay_estimate >= 0) {
927 // Update delay estimate buffer.
928 aec->delay_histogram[delay_estimate]++;
929 }
930 }
931 }
932
933 // Update the xfBuf block position.
934 aec->xfBufBlockPos--;
935 if (aec->xfBufBlockPos == -1) {
936 aec->xfBufBlockPos = aec->num_partitions - 1;
937 }
938
939 // Buffer xf
940 memcpy(aec->xfBuf[0] + aec->xfBufBlockPos * PART_LEN1,
941 xf_ptr,
942 sizeof(float) * PART_LEN1);
943 memcpy(aec->xfBuf[1] + aec->xfBufBlockPos * PART_LEN1,
944 &xf_ptr[PART_LEN1],
945 sizeof(float) * PART_LEN1);
946
947 memset(yf, 0, sizeof(yf));
948
949 // Filter far
950 WebRtcAec_FilterFar(aec, yf);
951
952 // Inverse fft to obtain echo estimate and error.
953 fft[0] = yf[0][0];
954 fft[1] = yf[0][PART_LEN];
955 for (i = 1; i < PART_LEN; i++) {
956 fft[2 * i] = yf[0][i];
957 fft[2 * i + 1] = yf[1][i];
958 }
959 aec_rdft_inverse_128(fft);
960
961 scale = 2.0f / PART_LEN2;
962 for (i = 0; i < PART_LEN; i++) {
963 y[i] = fft[PART_LEN + i] * scale; // fft scaling
964 }
965
966 for (i = 0; i < PART_LEN; i++) {
967 e[i] = nearend_ptr[i] - y[i];
968 }
969
970 // Error fft
971 memcpy(aec->eBuf + PART_LEN, e, sizeof(float) * PART_LEN);
972 memset(fft, 0, sizeof(float) * PART_LEN);
973 memcpy(fft + PART_LEN, e, sizeof(float) * PART_LEN);
974 // TODO(bjornv): Change to use TimeToFrequency().
975 aec_rdft_forward_128(fft);
976
977 ef[1][0] = 0;
978 ef[1][PART_LEN] = 0;
979 ef[0][0] = fft[0];
980 ef[0][PART_LEN] = fft[1];
981 for (i = 1; i < PART_LEN; i++) {
982 ef[0][i] = fft[2 * i];
983 ef[1][i] = fft[2 * i + 1];
984 }
985
986 if (aec->metricsMode == 1) {
987 // Note that the first PART_LEN samples in fft (before transformation) are
988 // zero. Hence, the scaling by two in UpdateLevel() should not be
989 // performed. That scaling is taken care of in UpdateMetrics() instead.
990 UpdateLevel(&aec->linoutlevel, ef);
991 }
992
993 // Scale error signal inversely with far power.
994 WebRtcAec_ScaleErrorSignal(aec, ef);
995 WebRtcAec_FilterAdaptation(aec, fft, ef);
996 NonLinearProcessing(aec, output, outputH);
997
998 if (aec->metricsMode == 1) {
999 // Update power levels and echo metrics
1000 UpdateLevel(&aec->farlevel, (float(*)[PART_LEN1])xf_ptr);
1001 UpdateLevel(&aec->nearlevel, df);
1002 UpdateMetrics(aec);
1003 }
1004
1005 // Store the output block.
1006 WebRtc_WriteBuffer(aec->outFrBuf, output, PART_LEN);
1007 // For H band
1008 if (aec->sampFreq == 32000) {
1009 WebRtc_WriteBuffer(aec->outFrBufH, outputH, PART_LEN);
1010 }
1011
1012 #ifdef WEBRTC_AEC_DEBUG_DUMP
1013 {
1014 int16_t eInt16[PART_LEN];
1015 for (i = 0; i < PART_LEN; i++) {
1016 eInt16[i] = (int16_t)WEBRTC_SPL_SAT(
1017 WEBRTC_SPL_WORD16_MAX, e[i], WEBRTC_SPL_WORD16_MIN);
1018 }
1019
1020 (void)fwrite(eInt16, sizeof(int16_t), PART_LEN, aec->outLinearFile);
1021 (void)fwrite(output, sizeof(int16_t), PART_LEN, aec->outFile);
1022 }
1023 #endif
1024 }
1025
NonLinearProcessing(AecCore * aec,float * output,float * outputH)1026 static void NonLinearProcessing(AecCore* aec, float* output, float* outputH) {
1027 float efw[2][PART_LEN1], dfw[2][PART_LEN1], xfw[2][PART_LEN1];
1028 complex_t comfortNoiseHband[PART_LEN1];
1029 float fft[PART_LEN2];
1030 float scale, dtmp;
1031 float nlpGainHband;
1032 int i, j, pos;
1033
1034 // Coherence and non-linear filter
1035 float cohde[PART_LEN1], cohxd[PART_LEN1];
1036 float hNlDeAvg, hNlXdAvg;
1037 float hNl[PART_LEN1];
1038 float hNlPref[kPrefBandSize];
1039 float hNlFb = 0, hNlFbLow = 0;
1040 const float prefBandQuant = 0.75f, prefBandQuantLow = 0.5f;
1041 const int prefBandSize = kPrefBandSize / aec->mult;
1042 const int minPrefBand = 4 / aec->mult;
1043
1044 // Near and error power sums
1045 float sdSum = 0, seSum = 0;
1046
1047 // Power estimate smoothing coefficients.
1048 const float* ptrGCoh = aec->extended_filter_enabled
1049 ? kExtendedSmoothingCoefficients[aec->mult - 1]
1050 : kNormalSmoothingCoefficients[aec->mult - 1];
1051 const float* min_overdrive = aec->extended_filter_enabled
1052 ? kExtendedMinOverDrive
1053 : kNormalMinOverDrive;
1054
1055 // Filter energy
1056 float wfEnMax = 0, wfEn = 0;
1057 const int delayEstInterval = 10 * aec->mult;
1058
1059 float* xfw_ptr = NULL;
1060
1061 aec->delayEstCtr++;
1062 if (aec->delayEstCtr == delayEstInterval) {
1063 aec->delayEstCtr = 0;
1064 }
1065
1066 // initialize comfort noise for H band
1067 memset(comfortNoiseHband, 0, sizeof(comfortNoiseHband));
1068 nlpGainHband = (float)0.0;
1069 dtmp = (float)0.0;
1070
1071 // Measure energy in each filter partition to determine delay.
1072 // TODO: Spread by computing one partition per block?
1073 if (aec->delayEstCtr == 0) {
1074 wfEnMax = 0;
1075 aec->delayIdx = 0;
1076 for (i = 0; i < aec->num_partitions; i++) {
1077 pos = i * PART_LEN1;
1078 wfEn = 0;
1079 for (j = 0; j < PART_LEN1; j++) {
1080 wfEn += aec->wfBuf[0][pos + j] * aec->wfBuf[0][pos + j] +
1081 aec->wfBuf[1][pos + j] * aec->wfBuf[1][pos + j];
1082 }
1083
1084 if (wfEn > wfEnMax) {
1085 wfEnMax = wfEn;
1086 aec->delayIdx = i;
1087 }
1088 }
1089 }
1090
1091 // We should always have at least one element stored in |far_buf|.
1092 assert(WebRtc_available_read(aec->far_buf_windowed) > 0);
1093 // NLP
1094 WebRtc_ReadBuffer(aec->far_buf_windowed, (void**)&xfw_ptr, &xfw[0][0], 1);
1095
1096 // TODO(bjornv): Investigate if we can reuse |far_buf_windowed| instead of
1097 // |xfwBuf|.
1098 // Buffer far.
1099 memcpy(aec->xfwBuf, xfw_ptr, sizeof(float) * 2 * PART_LEN1);
1100
1101 // Use delayed far.
1102 memcpy(xfw, aec->xfwBuf + aec->delayIdx * PART_LEN1, sizeof(xfw));
1103
1104 // Windowed near fft
1105 for (i = 0; i < PART_LEN; i++) {
1106 fft[i] = aec->dBuf[i] * sqrtHanning[i];
1107 fft[PART_LEN + i] = aec->dBuf[PART_LEN + i] * sqrtHanning[PART_LEN - i];
1108 }
1109 aec_rdft_forward_128(fft);
1110
1111 dfw[1][0] = 0;
1112 dfw[1][PART_LEN] = 0;
1113 dfw[0][0] = fft[0];
1114 dfw[0][PART_LEN] = fft[1];
1115 for (i = 1; i < PART_LEN; i++) {
1116 dfw[0][i] = fft[2 * i];
1117 dfw[1][i] = fft[2 * i + 1];
1118 }
1119
1120 // Windowed error fft
1121 for (i = 0; i < PART_LEN; i++) {
1122 fft[i] = aec->eBuf[i] * sqrtHanning[i];
1123 fft[PART_LEN + i] = aec->eBuf[PART_LEN + i] * sqrtHanning[PART_LEN - i];
1124 }
1125 aec_rdft_forward_128(fft);
1126 efw[1][0] = 0;
1127 efw[1][PART_LEN] = 0;
1128 efw[0][0] = fft[0];
1129 efw[0][PART_LEN] = fft[1];
1130 for (i = 1; i < PART_LEN; i++) {
1131 efw[0][i] = fft[2 * i];
1132 efw[1][i] = fft[2 * i + 1];
1133 }
1134
1135 // Smoothed PSD
1136 for (i = 0; i < PART_LEN1; i++) {
1137 aec->sd[i] = ptrGCoh[0] * aec->sd[i] +
1138 ptrGCoh[1] * (dfw[0][i] * dfw[0][i] + dfw[1][i] * dfw[1][i]);
1139 aec->se[i] = ptrGCoh[0] * aec->se[i] +
1140 ptrGCoh[1] * (efw[0][i] * efw[0][i] + efw[1][i] * efw[1][i]);
1141 // We threshold here to protect against the ill-effects of a zero farend.
1142 // The threshold is not arbitrarily chosen, but balances protection and
1143 // adverse interaction with the algorithm's tuning.
1144 // TODO: investigate further why this is so sensitive.
1145 aec->sx[i] =
1146 ptrGCoh[0] * aec->sx[i] +
1147 ptrGCoh[1] *
1148 WEBRTC_SPL_MAX(xfw[0][i] * xfw[0][i] + xfw[1][i] * xfw[1][i], 15);
1149
1150 aec->sde[i][0] =
1151 ptrGCoh[0] * aec->sde[i][0] +
1152 ptrGCoh[1] * (dfw[0][i] * efw[0][i] + dfw[1][i] * efw[1][i]);
1153 aec->sde[i][1] =
1154 ptrGCoh[0] * aec->sde[i][1] +
1155 ptrGCoh[1] * (dfw[0][i] * efw[1][i] - dfw[1][i] * efw[0][i]);
1156
1157 aec->sxd[i][0] =
1158 ptrGCoh[0] * aec->sxd[i][0] +
1159 ptrGCoh[1] * (dfw[0][i] * xfw[0][i] + dfw[1][i] * xfw[1][i]);
1160 aec->sxd[i][1] =
1161 ptrGCoh[0] * aec->sxd[i][1] +
1162 ptrGCoh[1] * (dfw[0][i] * xfw[1][i] - dfw[1][i] * xfw[0][i]);
1163
1164 sdSum += aec->sd[i];
1165 seSum += aec->se[i];
1166 }
1167
1168 // Divergent filter safeguard.
1169 if (aec->divergeState == 0) {
1170 if (seSum > sdSum) {
1171 aec->divergeState = 1;
1172 }
1173 } else {
1174 if (seSum * 1.05f < sdSum) {
1175 aec->divergeState = 0;
1176 }
1177 }
1178
1179 if (aec->divergeState == 1) {
1180 memcpy(efw, dfw, sizeof(efw));
1181 }
1182
1183 if (!aec->extended_filter_enabled) {
1184 // Reset if error is significantly larger than nearend (13 dB).
1185 if (seSum > (19.95f * sdSum)) {
1186 memset(aec->wfBuf, 0, sizeof(aec->wfBuf));
1187 }
1188 }
1189
1190 // Subband coherence
1191 for (i = 0; i < PART_LEN1; i++) {
1192 cohde[i] =
1193 (aec->sde[i][0] * aec->sde[i][0] + aec->sde[i][1] * aec->sde[i][1]) /
1194 (aec->sd[i] * aec->se[i] + 1e-10f);
1195 cohxd[i] =
1196 (aec->sxd[i][0] * aec->sxd[i][0] + aec->sxd[i][1] * aec->sxd[i][1]) /
1197 (aec->sx[i] * aec->sd[i] + 1e-10f);
1198 }
1199
1200 hNlXdAvg = 0;
1201 for (i = minPrefBand; i < prefBandSize + minPrefBand; i++) {
1202 hNlXdAvg += cohxd[i];
1203 }
1204 hNlXdAvg /= prefBandSize;
1205 hNlXdAvg = 1 - hNlXdAvg;
1206
1207 hNlDeAvg = 0;
1208 for (i = minPrefBand; i < prefBandSize + minPrefBand; i++) {
1209 hNlDeAvg += cohde[i];
1210 }
1211 hNlDeAvg /= prefBandSize;
1212
1213 if (hNlXdAvg < 0.75f && hNlXdAvg < aec->hNlXdAvgMin) {
1214 aec->hNlXdAvgMin = hNlXdAvg;
1215 }
1216
1217 if (hNlDeAvg > 0.98f && hNlXdAvg > 0.9f) {
1218 aec->stNearState = 1;
1219 } else if (hNlDeAvg < 0.95f || hNlXdAvg < 0.8f) {
1220 aec->stNearState = 0;
1221 }
1222
1223 if (aec->hNlXdAvgMin == 1) {
1224 aec->echoState = 0;
1225 aec->overDrive = min_overdrive[aec->nlp_mode];
1226
1227 if (aec->stNearState == 1) {
1228 memcpy(hNl, cohde, sizeof(hNl));
1229 hNlFb = hNlDeAvg;
1230 hNlFbLow = hNlDeAvg;
1231 } else {
1232 for (i = 0; i < PART_LEN1; i++) {
1233 hNl[i] = 1 - cohxd[i];
1234 }
1235 hNlFb = hNlXdAvg;
1236 hNlFbLow = hNlXdAvg;
1237 }
1238 } else {
1239
1240 if (aec->stNearState == 1) {
1241 aec->echoState = 0;
1242 memcpy(hNl, cohde, sizeof(hNl));
1243 hNlFb = hNlDeAvg;
1244 hNlFbLow = hNlDeAvg;
1245 } else {
1246 aec->echoState = 1;
1247 for (i = 0; i < PART_LEN1; i++) {
1248 hNl[i] = WEBRTC_SPL_MIN(cohde[i], 1 - cohxd[i]);
1249 }
1250
1251 // Select an order statistic from the preferred bands.
1252 // TODO: Using quicksort now, but a selection algorithm may be preferred.
1253 memcpy(hNlPref, &hNl[minPrefBand], sizeof(float) * prefBandSize);
1254 qsort(hNlPref, prefBandSize, sizeof(float), CmpFloat);
1255 hNlFb = hNlPref[(int)floor(prefBandQuant * (prefBandSize - 1))];
1256 hNlFbLow = hNlPref[(int)floor(prefBandQuantLow * (prefBandSize - 1))];
1257 }
1258 }
1259
1260 // Track the local filter minimum to determine suppression overdrive.
1261 if (hNlFbLow < 0.6f && hNlFbLow < aec->hNlFbLocalMin) {
1262 aec->hNlFbLocalMin = hNlFbLow;
1263 aec->hNlFbMin = hNlFbLow;
1264 aec->hNlNewMin = 1;
1265 aec->hNlMinCtr = 0;
1266 }
1267 aec->hNlFbLocalMin =
1268 WEBRTC_SPL_MIN(aec->hNlFbLocalMin + 0.0008f / aec->mult, 1);
1269 aec->hNlXdAvgMin = WEBRTC_SPL_MIN(aec->hNlXdAvgMin + 0.0006f / aec->mult, 1);
1270
1271 if (aec->hNlNewMin == 1) {
1272 aec->hNlMinCtr++;
1273 }
1274 if (aec->hNlMinCtr == 2) {
1275 aec->hNlNewMin = 0;
1276 aec->hNlMinCtr = 0;
1277 aec->overDrive =
1278 WEBRTC_SPL_MAX(kTargetSupp[aec->nlp_mode] /
1279 ((float)log(aec->hNlFbMin + 1e-10f) + 1e-10f),
1280 min_overdrive[aec->nlp_mode]);
1281 }
1282
1283 // Smooth the overdrive.
1284 if (aec->overDrive < aec->overDriveSm) {
1285 aec->overDriveSm = 0.99f * aec->overDriveSm + 0.01f * aec->overDrive;
1286 } else {
1287 aec->overDriveSm = 0.9f * aec->overDriveSm + 0.1f * aec->overDrive;
1288 }
1289
1290 WebRtcAec_OverdriveAndSuppress(aec, hNl, hNlFb, efw);
1291
1292 // Add comfort noise.
1293 WebRtcAec_ComfortNoise(aec, efw, comfortNoiseHband, aec->noisePow, hNl);
1294
1295 // TODO(bjornv): Investigate how to take the windowing below into account if
1296 // needed.
1297 if (aec->metricsMode == 1) {
1298 // Note that we have a scaling by two in the time domain |eBuf|.
1299 // In addition the time domain signal is windowed before transformation,
1300 // losing half the energy on the average. We take care of the first
1301 // scaling only in UpdateMetrics().
1302 UpdateLevel(&aec->nlpoutlevel, efw);
1303 }
1304 // Inverse error fft.
1305 fft[0] = efw[0][0];
1306 fft[1] = efw[0][PART_LEN];
1307 for (i = 1; i < PART_LEN; i++) {
1308 fft[2 * i] = efw[0][i];
1309 // Sign change required by Ooura fft.
1310 fft[2 * i + 1] = -efw[1][i];
1311 }
1312 aec_rdft_inverse_128(fft);
1313
1314 // Overlap and add to obtain output.
1315 scale = 2.0f / PART_LEN2;
1316 for (i = 0; i < PART_LEN; i++) {
1317 fft[i] *= scale; // fft scaling
1318 fft[i] = fft[i] * sqrtHanning[i] + aec->outBuf[i];
1319
1320 fft[PART_LEN + i] *= scale; // fft scaling
1321 aec->outBuf[i] = fft[PART_LEN + i] * sqrtHanning[PART_LEN - i];
1322
1323 // Saturate output to keep it in the allowed range.
1324 output[i] = WEBRTC_SPL_SAT(
1325 WEBRTC_SPL_WORD16_MAX, fft[i], WEBRTC_SPL_WORD16_MIN);
1326 }
1327
1328 // For H band
1329 if (aec->sampFreq == 32000) {
1330
1331 // H band gain
1332 // average nlp over low band: average over second half of freq spectrum
1333 // (4->8khz)
1334 GetHighbandGain(hNl, &nlpGainHband);
1335
1336 // Inverse comfort_noise
1337 if (flagHbandCn == 1) {
1338 fft[0] = comfortNoiseHband[0][0];
1339 fft[1] = comfortNoiseHband[PART_LEN][0];
1340 for (i = 1; i < PART_LEN; i++) {
1341 fft[2 * i] = comfortNoiseHband[i][0];
1342 fft[2 * i + 1] = comfortNoiseHband[i][1];
1343 }
1344 aec_rdft_inverse_128(fft);
1345 scale = 2.0f / PART_LEN2;
1346 }
1347
1348 // compute gain factor
1349 for (i = 0; i < PART_LEN; i++) {
1350 dtmp = aec->dBufH[i];
1351 dtmp = dtmp * nlpGainHband; // for variable gain
1352
1353 // add some comfort noise where Hband is attenuated
1354 if (flagHbandCn == 1) {
1355 fft[i] *= scale; // fft scaling
1356 dtmp += cnScaleHband * fft[i];
1357 }
1358
1359 // Saturate output to keep it in the allowed range.
1360 outputH[i] = WEBRTC_SPL_SAT(
1361 WEBRTC_SPL_WORD16_MAX, dtmp, WEBRTC_SPL_WORD16_MIN);
1362 }
1363 }
1364
1365 // Copy the current block to the old position.
1366 memcpy(aec->dBuf, aec->dBuf + PART_LEN, sizeof(float) * PART_LEN);
1367 memcpy(aec->eBuf, aec->eBuf + PART_LEN, sizeof(float) * PART_LEN);
1368
1369 // Copy the current block to the old position for H band
1370 if (aec->sampFreq == 32000) {
1371 memcpy(aec->dBufH, aec->dBufH + PART_LEN, sizeof(float) * PART_LEN);
1372 }
1373
1374 memmove(aec->xfwBuf + PART_LEN1,
1375 aec->xfwBuf,
1376 sizeof(aec->xfwBuf) - sizeof(complex_t) * PART_LEN1);
1377 }
1378
GetHighbandGain(const float * lambda,float * nlpGainHband)1379 static void GetHighbandGain(const float* lambda, float* nlpGainHband) {
1380 int i;
1381
1382 nlpGainHband[0] = (float)0.0;
1383 for (i = freqAvgIc; i < PART_LEN1 - 1; i++) {
1384 nlpGainHband[0] += lambda[i];
1385 }
1386 nlpGainHband[0] /= (float)(PART_LEN1 - 1 - freqAvgIc);
1387 }
1388
ComfortNoise(AecCore * aec,float efw[2][PART_LEN1],complex_t * comfortNoiseHband,const float * noisePow,const float * lambda)1389 static void ComfortNoise(AecCore* aec,
1390 float efw[2][PART_LEN1],
1391 complex_t* comfortNoiseHband,
1392 const float* noisePow,
1393 const float* lambda) {
1394 int i, num;
1395 float rand[PART_LEN];
1396 float noise, noiseAvg, tmp, tmpAvg;
1397 int16_t randW16[PART_LEN];
1398 complex_t u[PART_LEN1];
1399
1400 const float pi2 = 6.28318530717959f;
1401
1402 // Generate a uniform random array on [0 1]
1403 WebRtcSpl_RandUArray(randW16, PART_LEN, &aec->seed);
1404 for (i = 0; i < PART_LEN; i++) {
1405 rand[i] = ((float)randW16[i]) / 32768;
1406 }
1407
1408 // Reject LF noise
1409 u[0][0] = 0;
1410 u[0][1] = 0;
1411 for (i = 1; i < PART_LEN1; i++) {
1412 tmp = pi2 * rand[i - 1];
1413
1414 noise = sqrtf(noisePow[i]);
1415 u[i][0] = noise * cosf(tmp);
1416 u[i][1] = -noise * sinf(tmp);
1417 }
1418 u[PART_LEN][1] = 0;
1419
1420 for (i = 0; i < PART_LEN1; i++) {
1421 // This is the proper weighting to match the background noise power
1422 tmp = sqrtf(WEBRTC_SPL_MAX(1 - lambda[i] * lambda[i], 0));
1423 // tmp = 1 - lambda[i];
1424 efw[0][i] += tmp * u[i][0];
1425 efw[1][i] += tmp * u[i][1];
1426 }
1427
1428 // For H band comfort noise
1429 // TODO: don't compute noise and "tmp" twice. Use the previous results.
1430 noiseAvg = 0.0;
1431 tmpAvg = 0.0;
1432 num = 0;
1433 if (aec->sampFreq == 32000 && flagHbandCn == 1) {
1434
1435 // average noise scale
1436 // average over second half of freq spectrum (i.e., 4->8khz)
1437 // TODO: we shouldn't need num. We know how many elements we're summing.
1438 for (i = PART_LEN1 >> 1; i < PART_LEN1; i++) {
1439 num++;
1440 noiseAvg += sqrtf(noisePow[i]);
1441 }
1442 noiseAvg /= (float)num;
1443
1444 // average nlp scale
1445 // average over second half of freq spectrum (i.e., 4->8khz)
1446 // TODO: we shouldn't need num. We know how many elements we're summing.
1447 num = 0;
1448 for (i = PART_LEN1 >> 1; i < PART_LEN1; i++) {
1449 num++;
1450 tmpAvg += sqrtf(WEBRTC_SPL_MAX(1 - lambda[i] * lambda[i], 0));
1451 }
1452 tmpAvg /= (float)num;
1453
1454 // Use average noise for H band
1455 // TODO: we should probably have a new random vector here.
1456 // Reject LF noise
1457 u[0][0] = 0;
1458 u[0][1] = 0;
1459 for (i = 1; i < PART_LEN1; i++) {
1460 tmp = pi2 * rand[i - 1];
1461
1462 // Use average noise for H band
1463 u[i][0] = noiseAvg * (float)cos(tmp);
1464 u[i][1] = -noiseAvg * (float)sin(tmp);
1465 }
1466 u[PART_LEN][1] = 0;
1467
1468 for (i = 0; i < PART_LEN1; i++) {
1469 // Use average NLP weight for H band
1470 comfortNoiseHband[i][0] = tmpAvg * u[i][0];
1471 comfortNoiseHband[i][1] = tmpAvg * u[i][1];
1472 }
1473 }
1474 }
1475
InitLevel(PowerLevel * level)1476 static void InitLevel(PowerLevel* level) {
1477 const float kBigFloat = 1E17f;
1478
1479 level->averagelevel = 0;
1480 level->framelevel = 0;
1481 level->minlevel = kBigFloat;
1482 level->frsum = 0;
1483 level->sfrsum = 0;
1484 level->frcounter = 0;
1485 level->sfrcounter = 0;
1486 }
1487
InitStats(Stats * stats)1488 static void InitStats(Stats* stats) {
1489 stats->instant = kOffsetLevel;
1490 stats->average = kOffsetLevel;
1491 stats->max = kOffsetLevel;
1492 stats->min = kOffsetLevel * (-1);
1493 stats->sum = 0;
1494 stats->hisum = 0;
1495 stats->himean = kOffsetLevel;
1496 stats->counter = 0;
1497 stats->hicounter = 0;
1498 }
1499
InitMetrics(AecCore * self)1500 static void InitMetrics(AecCore* self) {
1501 self->stateCounter = 0;
1502 InitLevel(&self->farlevel);
1503 InitLevel(&self->nearlevel);
1504 InitLevel(&self->linoutlevel);
1505 InitLevel(&self->nlpoutlevel);
1506
1507 InitStats(&self->erl);
1508 InitStats(&self->erle);
1509 InitStats(&self->aNlp);
1510 InitStats(&self->rerl);
1511 }
1512
UpdateLevel(PowerLevel * level,float in[2][PART_LEN1])1513 static void UpdateLevel(PowerLevel* level, float in[2][PART_LEN1]) {
1514 // Do the energy calculation in the frequency domain. The FFT is performed on
1515 // a segment of PART_LEN2 samples due to overlap, but we only want the energy
1516 // of half that data (the last PART_LEN samples). Parseval's relation states
1517 // that the energy is preserved according to
1518 //
1519 // \sum_{n=0}^{N-1} |x(n)|^2 = 1/N * \sum_{n=0}^{N-1} |X(n)|^2
1520 // = ENERGY,
1521 //
1522 // where N = PART_LEN2. Since we are only interested in calculating the energy
1523 // for the last PART_LEN samples we approximate by calculating ENERGY and
1524 // divide by 2,
1525 //
1526 // \sum_{n=N/2}^{N-1} |x(n)|^2 ~= ENERGY / 2
1527 //
1528 // Since we deal with real valued time domain signals we only store frequency
1529 // bins [0, PART_LEN], which is what |in| consists of. To calculate ENERGY we
1530 // need to add the contribution from the missing part in
1531 // [PART_LEN+1, PART_LEN2-1]. These values are, up to a phase shift, identical
1532 // with the values in [1, PART_LEN-1], hence multiply those values by 2. This
1533 // is the values in the for loop below, but multiplication by 2 and division
1534 // by 2 cancel.
1535
1536 // TODO(bjornv): Investigate reusing energy calculations performed at other
1537 // places in the code.
1538 int k = 1;
1539 // Imaginary parts are zero at end points and left out of the calculation.
1540 float energy = (in[0][0] * in[0][0]) / 2;
1541 energy += (in[0][PART_LEN] * in[0][PART_LEN]) / 2;
1542
1543 for (k = 1; k < PART_LEN; k++) {
1544 energy += (in[0][k] * in[0][k] + in[1][k] * in[1][k]);
1545 }
1546 energy /= PART_LEN2;
1547
1548 level->sfrsum += energy;
1549 level->sfrcounter++;
1550
1551 if (level->sfrcounter > subCountLen) {
1552 level->framelevel = level->sfrsum / (subCountLen * PART_LEN);
1553 level->sfrsum = 0;
1554 level->sfrcounter = 0;
1555 if (level->framelevel > 0) {
1556 if (level->framelevel < level->minlevel) {
1557 level->minlevel = level->framelevel; // New minimum.
1558 } else {
1559 level->minlevel *= (1 + 0.001f); // Small increase.
1560 }
1561 }
1562 level->frcounter++;
1563 level->frsum += level->framelevel;
1564 if (level->frcounter > countLen) {
1565 level->averagelevel = level->frsum / countLen;
1566 level->frsum = 0;
1567 level->frcounter = 0;
1568 }
1569 }
1570 }
1571
UpdateMetrics(AecCore * aec)1572 static void UpdateMetrics(AecCore* aec) {
1573 float dtmp, dtmp2;
1574
1575 const float actThresholdNoisy = 8.0f;
1576 const float actThresholdClean = 40.0f;
1577 const float safety = 0.99995f;
1578 const float noisyPower = 300000.0f;
1579
1580 float actThreshold;
1581 float echo, suppressedEcho;
1582
1583 if (aec->echoState) { // Check if echo is likely present
1584 aec->stateCounter++;
1585 }
1586
1587 if (aec->farlevel.frcounter == 0) {
1588
1589 if (aec->farlevel.minlevel < noisyPower) {
1590 actThreshold = actThresholdClean;
1591 } else {
1592 actThreshold = actThresholdNoisy;
1593 }
1594
1595 if ((aec->stateCounter > (0.5f * countLen * subCountLen)) &&
1596 (aec->farlevel.sfrcounter == 0)
1597
1598 // Estimate in active far-end segments only
1599 &&
1600 (aec->farlevel.averagelevel >
1601 (actThreshold * aec->farlevel.minlevel))) {
1602
1603 // Subtract noise power
1604 echo = aec->nearlevel.averagelevel - safety * aec->nearlevel.minlevel;
1605
1606 // ERL
1607 dtmp = 10 * (float)log10(aec->farlevel.averagelevel /
1608 aec->nearlevel.averagelevel +
1609 1e-10f);
1610 dtmp2 = 10 * (float)log10(aec->farlevel.averagelevel / echo + 1e-10f);
1611
1612 aec->erl.instant = dtmp;
1613 if (dtmp > aec->erl.max) {
1614 aec->erl.max = dtmp;
1615 }
1616
1617 if (dtmp < aec->erl.min) {
1618 aec->erl.min = dtmp;
1619 }
1620
1621 aec->erl.counter++;
1622 aec->erl.sum += dtmp;
1623 aec->erl.average = aec->erl.sum / aec->erl.counter;
1624
1625 // Upper mean
1626 if (dtmp > aec->erl.average) {
1627 aec->erl.hicounter++;
1628 aec->erl.hisum += dtmp;
1629 aec->erl.himean = aec->erl.hisum / aec->erl.hicounter;
1630 }
1631
1632 // A_NLP
1633 dtmp = 10 * (float)log10(aec->nearlevel.averagelevel /
1634 (2 * aec->linoutlevel.averagelevel) +
1635 1e-10f);
1636
1637 // subtract noise power
1638 suppressedEcho = 2 * (aec->linoutlevel.averagelevel -
1639 safety * aec->linoutlevel.minlevel);
1640
1641 dtmp2 = 10 * (float)log10(echo / suppressedEcho + 1e-10f);
1642
1643 aec->aNlp.instant = dtmp2;
1644 if (dtmp > aec->aNlp.max) {
1645 aec->aNlp.max = dtmp;
1646 }
1647
1648 if (dtmp < aec->aNlp.min) {
1649 aec->aNlp.min = dtmp;
1650 }
1651
1652 aec->aNlp.counter++;
1653 aec->aNlp.sum += dtmp;
1654 aec->aNlp.average = aec->aNlp.sum / aec->aNlp.counter;
1655
1656 // Upper mean
1657 if (dtmp > aec->aNlp.average) {
1658 aec->aNlp.hicounter++;
1659 aec->aNlp.hisum += dtmp;
1660 aec->aNlp.himean = aec->aNlp.hisum / aec->aNlp.hicounter;
1661 }
1662
1663 // ERLE
1664
1665 // subtract noise power
1666 suppressedEcho = 2 * (aec->nlpoutlevel.averagelevel -
1667 safety * aec->nlpoutlevel.minlevel);
1668
1669 dtmp = 10 * (float)log10(aec->nearlevel.averagelevel /
1670 (2 * aec->nlpoutlevel.averagelevel) +
1671 1e-10f);
1672 dtmp2 = 10 * (float)log10(echo / suppressedEcho + 1e-10f);
1673
1674 dtmp = dtmp2;
1675 aec->erle.instant = dtmp;
1676 if (dtmp > aec->erle.max) {
1677 aec->erle.max = dtmp;
1678 }
1679
1680 if (dtmp < aec->erle.min) {
1681 aec->erle.min = dtmp;
1682 }
1683
1684 aec->erle.counter++;
1685 aec->erle.sum += dtmp;
1686 aec->erle.average = aec->erle.sum / aec->erle.counter;
1687
1688 // Upper mean
1689 if (dtmp > aec->erle.average) {
1690 aec->erle.hicounter++;
1691 aec->erle.hisum += dtmp;
1692 aec->erle.himean = aec->erle.hisum / aec->erle.hicounter;
1693 }
1694 }
1695
1696 aec->stateCounter = 0;
1697 }
1698 }
1699
TimeToFrequency(float time_data[PART_LEN2],float freq_data[2][PART_LEN1],int window)1700 static void TimeToFrequency(float time_data[PART_LEN2],
1701 float freq_data[2][PART_LEN1],
1702 int window) {
1703 int i = 0;
1704
1705 // TODO(bjornv): Should we have a different function/wrapper for windowed FFT?
1706 if (window) {
1707 for (i = 0; i < PART_LEN; i++) {
1708 time_data[i] *= sqrtHanning[i];
1709 time_data[PART_LEN + i] *= sqrtHanning[PART_LEN - i];
1710 }
1711 }
1712
1713 aec_rdft_forward_128(time_data);
1714 // Reorder.
1715 freq_data[1][0] = 0;
1716 freq_data[1][PART_LEN] = 0;
1717 freq_data[0][0] = time_data[0];
1718 freq_data[0][PART_LEN] = time_data[1];
1719 for (i = 1; i < PART_LEN; i++) {
1720 freq_data[0][i] = time_data[2 * i];
1721 freq_data[1][i] = time_data[2 * i + 1];
1722 }
1723 }
1724