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