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