src/third_party/webrtc/modules/audio_processing/aec/aec_core.c

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