Calculate audio levels in AEC in time domain.

webrtc/modules/audio_processing/aec/aec_core.c

 /*
  *  Copyright (c) 2012 The WebRTC project authors. All Rights Reserved.
  *
  *  Use of this source code is governed by a BSD-style license
  *  that can be found in the LICENSE file in the root of the source
  *  tree. An additional intellectual property rights grant can be found
  *  in the file PATENTS.  All contributing project authors may
  *  be found in the AUTHORS file in the root of the source tree.
  */
 
@@ skipping 547 matching lines @@
   InitLevel(&self->nearlevel);
   InitLevel(&self->linoutlevel);
   InitLevel(&self->nlpoutlevel);
 
   InitStats(&self->erl);
   InitStats(&self->erle);
   InitStats(&self->aNlp);
   InitStats(&self->rerl);
 }
 
-static void UpdateLevel(PowerLevel* level, float in[2][PART_LEN1]) {
-  // Do the energy calculation in the frequency domain. The FFT is performed on
-  // a segment of PART_LEN2 samples due to overlap, but we only want the energy
-  // of half that data (the last PART_LEN samples). Parseval's relation states
-  // that the energy is preserved according to
-  //
-  // \sum_{n=0}^{N-1} |x(n)|^2 = 1/N * \sum_{n=0}^{N-1} |X(n)|^2
-  //                           = ENERGY,
-  //
-  // where N = PART_LEN2. Since we are only interested in calculating the energy
-  // for the last PART_LEN samples we approximate by calculating ENERGY and
-  // divide by 2,
-  //
-  // \sum_{n=N/2}^{N-1} |x(n)|^2 ~= ENERGY / 2
-  //
-  // Since we deal with real valued time domain signals we only store frequency
-  // bins [0, PART_LEN], which is what |in| consists of. To calculate ENERGY we
-  // need to add the contribution from the missing part in
-  // [PART_LEN+1, PART_LEN2-1]. These values are, up to a phase shift, identical
-  // with the values in [1, PART_LEN-1], hence multiply those values by 2. This
-  // is the values in the for loop below, but multiplication by 2 and division
-  // by 2 cancel.
+static float CalculatePower(const float* in, size_t num_samples) {
+  size_t k;
+  float energy = 0.0f;
 
-  // TODO(bjornv): Investigate reusing energy calculations performed at other
-  // places in the code.
-  int k = 1;
-  // Imaginary parts are zero at end points and left out of the calculation.
-  float energy = (in[0][0] * in[0][0]) / 2;
-  energy += (in[0][PART_LEN] * in[0][PART_LEN]) / 2;
+  for (k = 0; k < num_samples; ++k) {
+    energy += in[k] * in[k];
+  }
+  return energy / num_samples;
+}
 
-  for (k = 1; k < PART_LEN; k++) {
-    energy += (in[0][k] * in[0][k] + in[1][k] * in[1][k]);
-  }
-  energy /= PART_LEN2;
-
+static void UpdateLevel(PowerLevel* level, float energy) {
   level->sfrsum += energy;
   level->sfrcounter++;
 
   if (level->sfrcounter > subCountLen) {
     level->framelevel = level->sfrsum / (subCountLen * PART_LEN);
     level->sfrsum = 0;
     level->sfrcounter = 0;
     if (level->framelevel > 0) {
       if (level->framelevel < level->minlevel) {
         level->minlevel = level->framelevel;  // New minimum.
@@ skipping 10 matching lines @@
     }
   }
 }
 
 static void UpdateMetrics(AecCore* aec) {
   float dtmp, dtmp2;
 
   const float actThresholdNoisy = 8.0f;
   const float actThresholdClean = 40.0f;
   const float safety = 0.99995f;
-  const float noisyPower = 300000.0f;
+
+  // To make noisePower consistent with the legacy code, a factor of
+  // 2.0f / PART_LEN2 is applied to noisyPower, since the legacy code uses
+  // the energy of a frame as the audio levels, while the new code uses a
+  // a per-sample energy (i.e., power).
+  const float noisyPower = 300000.0f * 2.0f / PART_LEN2;
 
   float actThreshold;
   float echo, suppressedEcho;
 
   if (aec->echoState) {  // Check if echo is likely present
     aec->stateCounter++;
   }
 
   if (aec->farlevel.frcounter == 0) {
 
@@ skipping 195 matching lines @@
   freq_data[1][0] = 0;
   freq_data[1][PART_LEN] = 0;
   freq_data[0][0] = time_data[0];
   freq_data[0][PART_LEN] = time_data[1];
   for (i = 1; i < PART_LEN; i++) {
     freq_data[0][i] = time_data[2 * i];
     freq_data[1][i] = time_data[2 * i + 1];
   }
 }
 
-
 static int SignalBasedDelayCorrection(AecCore* self) {
   int delay_correction = 0;
   int last_delay = -2;
   assert(self != NULL);
 #if !defined(WEBRTC_ANDROID)
   // On desktops, turn on correction after |kDelayCorrectionStart| frames.  This
   // is to let the delay estimation get a chance to converge.  Also, if the
   // playout audio volume is low (or even muted) the delay estimation can return
   // a very large delay, which will break the AEC if it is applied.
   if (self->frame_count < kDelayCorrectionStart) {
@@ skipping 112 matching lines @@
   Fft(e_extended, e_fft);
 
   RTC_AEC_DEBUG_RAW_WRITE(aec->e_fft_file,
                           &e_fft[0][0],
                           sizeof(e_fft[0][0]) * PART_LEN1 * 2);
 
   if (metrics_mode == 1) {
     // Note that the first PART_LEN samples in fft (before transformation) are
     // zero. Hence, the scaling by two in UpdateLevel() should not be
     // performed. That scaling is taken care of in UpdateMetrics() instead.
-    UpdateLevel(linout_level, e_fft);
+    UpdateLevel(linout_level, CalculatePower(e, PART_LEN) / 2.0f);
   }
 
   // Scale error signal inversely with far power.
   WebRtcAec_ScaleErrorSignal(extended_filter_enabled,
                              normal_mu,
                              normal_error_threshold,
                              x_pow,
                              e_fft);
   WebRtcAec_FilterAdaptation(num_partitions,
                              x_fft_buf_block_pos,
@@ skipping 171 matching lines @@
     aec->overDriveSm = 0.99f * aec->overDriveSm + 0.01f * aec->overDrive;
   } else {
     aec->overDriveSm = 0.9f * aec->overDriveSm + 0.1f * aec->overDrive;
   }
 
   WebRtcAec_OverdriveAndSuppress(aec, hNl, hNlFb, efw);
 
   // Add comfort noise.
   WebRtcAec_ComfortNoise(aec, efw, comfortNoiseHband, aec->noisePow, hNl);
 
+  // Inverse error fft.
+  ScaledInverseFft(efw, fft, 2.0f, 1);
+
   // TODO(bjornv): Investigate how to take the windowing below into account if
   // needed.
   if (aec->metricsMode == 1) {
     // Note that we have a scaling by two in the time domain |eBuf|.
     // In addition the time domain signal is windowed before transformation,
     // losing half the energy on the average. We take care of the first
     // scaling only in UpdateMetrics().
-    UpdateLevel(&aec->nlpoutlevel, efw);
+    UpdateLevel(&aec->nlpoutlevel, CalculatePower(fft, PART_LEN2));
   }
 
-  // Inverse error fft.
-  ScaledInverseFft(efw, fft, 2.0f, 1);
-
   // Overlap and add to obtain output.
   for (i = 0; i < PART_LEN; i++) {
     output[i] = (fft[i] * WebRtcAec_sqrtHanning[i] +
                  aec->outBuf[i] * WebRtcAec_sqrtHanning[PART_LEN - i]);
 
     // Saturate output to keep it in the allowed range.
     output[i] = WEBRTC_SPL_SAT(
         WEBRTC_SPL_WORD16_MAX, output[i], WEBRTC_SPL_WORD16_MIN);
   }
   memcpy(aec->outBuf, &fft[PART_LEN], PART_LEN * sizeof(aec->outBuf[0]));
@@ skipping 104 matching lines @@
 
   // Convert far-end signal to the frequency domain.
   memcpy(fft, farend_ptr, sizeof(float) * PART_LEN2);
   Fft(fft, xf);
   xf_ptr = &xf[0][0];
 
   // Near fft
   memcpy(fft, aec->dBuf, sizeof(float) * PART_LEN2);
   Fft(fft, df);
 
+  if (aec->metricsMode == 1) {
+    // Update power levels
+    UpdateLevel(&aec->farlevel, CalculatePower(farend_ptr, PART_LEN2));
+    UpdateLevel(&aec->nearlevel, CalculatePower(aec->dBuf, PART_LEN2));
+  }
+
   // Power smoothing
   for (i = 0; i < PART_LEN1; i++) {
     far_spectrum = (xf_ptr[i] * xf_ptr[i]) +
                    (xf_ptr[PART_LEN1 + i] * xf_ptr[PART_LEN1 + i]);
     aec->xPow[i] =
         gPow[0] * aec->xPow[i] + gPow[1] * aec->num_partitions * far_spectrum;
     // Calculate absolute spectra
     abs_far_spectrum[i] = sqrtf(far_spectrum);
 
     near_spectrum = df[0][i] * df[0][i] + df[1][i] * df[1][i];
@@ skipping 77 matching lines @@
                   aec->wfBuf,
                   &aec->linoutlevel,
                   echo_subtractor_output);
 
   RTC_AEC_DEBUG_WAV_WRITE(aec->outLinearFile, echo_subtractor_output, PART_LEN);
 
   // Perform echo suppression.
   EchoSuppression(aec, farend_ptr, echo_subtractor_output, output, outputH_ptr);
 
   if (aec->metricsMode == 1) {
-    // Update power levels and echo metrics
-    UpdateLevel(&aec->farlevel, (float(*)[PART_LEN1])xf_ptr);
-    UpdateLevel(&aec->nearlevel, df);
     UpdateMetrics(aec);
   }
 
   // Store the output block.
   WebRtc_WriteBuffer(aec->outFrBuf, output, PART_LEN);
   // For high bands
   for (i = 0; i < aec->num_bands - 1; ++i) {
     WebRtc_WriteBuffer(aec->outFrBufH[i], outputH[i], PART_LEN);
   }
 
@@ skipping 520 matching lines @@
 int WebRtcAec_extended_filter_enabled(AecCore* self) {
   return self->extended_filter_enabled;
 }
 
 int WebRtcAec_system_delay(AecCore* self) { return self->system_delay; }
 
 void WebRtcAec_SetSystemDelay(AecCore* self, int delay) {
   assert(delay >= 0);
   self->system_delay = delay;
 }