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