Added support in the AEC for refined filter adaptation.

webrtc/modules/audio_processing/aec/aec_core.cc

 /*
  *  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.
  */
 
 /*
  * The core AEC algorithm, which is presented with time-aligned signals.
  */
 
 #include "webrtc/modules/audio_processing/aec/aec_core.h"
 
 #ifdef WEBRTC_AEC_DEBUG_DUMP
 #include <stdio.h>
 #endif
 
 #include <algorithm>
 #include <assert.h>
 #include <math.h>
 #include <stddef.h>  // size_t
 #include <stdlib.h>
 #include <string.h>
 
+#include "webrtc/base/checks.h"
 extern "C" {
 #include "webrtc/common_audio/ring_buffer.h"
 }
 #include "webrtc/common_audio/signal_processing/include/signal_processing_library.h"
 #include "webrtc/modules/audio_processing/aec/aec_common.h"
 #include "webrtc/modules/audio_processing/aec/aec_core_internal.h"
 #include "webrtc/modules/audio_processing/aec/aec_rdft.h"
 #include "webrtc/modules/audio_processing/logging/aec_logging.h"
 #include "webrtc/modules/audio_processing/utility/delay_estimator_wrapper.h"
 #include "webrtc/system_wrappers/include/cpu_features_wrapper.h"
@@ skipping 193 matching lines @@
 
     for (j = 0; j < PART_LEN1; j++) {
       y_fft[0][j] += MulRe(x_fft_buf[0][xPos + j], x_fft_buf[1][xPos + j],
                            h_fft_buf[0][pos + j], h_fft_buf[1][pos + j]);
       y_fft[1][j] += MulIm(x_fft_buf[0][xPos + j], x_fft_buf[1][xPos + j],
                            h_fft_buf[0][pos + j], h_fft_buf[1][pos + j]);
     }
   }
 }
 
 static void ScaleErrorSignal(int extended_filter_enabled,

[tlegrand-webrtc, 2016/04/14 14:32:42]

This function don't really need the first parameter, if you instead call it with
normal_error_threshold = kExtendedErrorThreshold (and rename
normal_error_threshold to simply error_threshold).

[peah-webrtc, 2016/04/14 22:12:41]

That is definitely a good change. I did not want to do it in this CL, as it is
not really related to the changes being done. But I now grab the opportunity and
do the change :-)
Done.

[tlegrand-webrtc, 2016/04/15 09:07:04]

Great! I think it's fine, since the CL isn't that big, and it makes the code
easier to follow.

-                             float normal_mu,
+                             float mu,

[tlegrand-webrtc, 2016/04/14 14:32:42]

With removal of mu = kExtendedMu if extended filter is enabled, is the function
now called in a different way? Otherwise, you have changed behavior.

[hlundin-webrtc, 2016/04/14 14:47:46]

It is now called (through EchoSubtraction) with aec->filter_step_size instead of
aec->norma_mu. 

[peah-webrtc, 2016/04/14 22:12:41]

Fully true.

[peah-webrtc, 2016/04/14 22:12:41]

It is now called in a different way but with a mu parameter that is updated in
another place. As it was, mu was selected in several places which did not make
sense. Now the decision is done in one place, in such a way that it is identical
to what it was before (if not the new functionality in this CL is active as the
modification of mu is part of that).

[tlegrand-webrtc, 2016/04/15 09:07:04]

Acknowledged.

                              float normal_error_threshold,
                              float x_pow[PART_LEN1],
                              float ef[2][PART_LEN1]) {
-  const float mu = extended_filter_enabled ? kExtendedMu : normal_mu;
   const float error_threshold = extended_filter_enabled
                                     ? kExtendedErrorThreshold
                                     : normal_error_threshold;
   int i;
   float abs_ef;
   for (i = 0; i < (PART_LEN1); i++) {
     ef[0][i] /= (x_pow[i] + 1e-10f);
     ef[1][i] /= (x_pow[i] + 1e-10f);
     abs_ef = sqrtf(ef[0][i] * ef[0][i] + ef[1][i] * ef[1][i]);
 
@@ skipping 672 matching lines @@
         (delay_quality > kDelayQualityThresholdMax ? kDelayQualityThresholdMax
                                                    : delay_quality);
     self->delay_quality_threshold =
         (delay_quality > self->delay_quality_threshold
              ? delay_quality
              : self->delay_quality_threshold);
   }
   return delay_correction;
 }
 
+static void RegressorPower(int num_partitions,
+                           int latest_added_partition,
+                           float x_fft_buf[2]
+                                          [kExtendedNumPartitions * PART_LEN1],
+                           float x_pow[PART_LEN1]) {
+  RTC_DCHECK_LT(latest_added_partition, num_partitions);
+  memset(x_pow, 0, PART_LEN1 * sizeof(x_pow[0]));
+
+  int partition = latest_added_partition;
+  int x_fft_buf_position = partition * PART_LEN1;
+  for (int i = 0; i < num_partitions; ++i) {
+    for (int bin = 0; bin < PART_LEN1; ++bin) {
+      float re = x_fft_buf[0][x_fft_buf_position];
+      float im = x_fft_buf[1][x_fft_buf_position];
+      x_pow[bin] += re * re + im * im;
+      ++x_fft_buf_position;
+    }
+
+    ++partition;
+    if (partition == num_partitions) {
+      partition = 0;
+      x_fft_buf_position = 0;

[hlundin-webrtc, 2016/04/14 14:02:26]

Can you DCHECK what you expect x_fft_buf_position to be before setting it to 0?

[peah-webrtc, 2016/04/14 22:12:41]

Good point! 
Done.

+    }
+  }
+}
+
 static void EchoSubtraction(AecCore* aec,
                             int num_partitions,
                             int extended_filter_enabled,
-                            float normal_mu,
+                            float filter_step_size,
                             float normal_error_threshold,
                             float* x_fft,
                             int* x_fft_buf_block_pos,
                             float x_fft_buf[2]
                                            [kExtendedNumPartitions * PART_LEN1],
                             float* const y,
                             float x_pow[PART_LEN1],
                             float h_fft_buf[2]
                                            [kExtendedNumPartitions * PART_LEN1],
                             float echo_subtractor_output[PART_LEN]) {
@@ skipping 41 matching lines @@
 
   // Compute the frequency domain echo prediction error.
   memset(e_extended, 0, sizeof(float) * PART_LEN);
   memcpy(e_extended + PART_LEN, e, sizeof(float) * PART_LEN);
   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);
 
   // Scale error signal inversely with far power.
-  WebRtcAec_ScaleErrorSignal(extended_filter_enabled, normal_mu,
+  WebRtcAec_ScaleErrorSignal(extended_filter_enabled, filter_step_size,
                              normal_error_threshold, x_pow, e_fft);
   WebRtcAec_FilterAdaptation(num_partitions, *x_fft_buf_block_pos, x_fft_buf,
                              e_fft, h_fft_buf);
   memcpy(echo_subtractor_output, e, sizeof(float) * PART_LEN);
 }
 
 static void EchoSuppression(AecCore* aec,
                             float farend[PART_LEN2],
                             float* echo_subtractor_output,
                             float* output,
@@ skipping 294 matching lines @@
   // Convert far-end signal to the frequency domain.
   memcpy(fft, farend_ptr, sizeof(float) * PART_LEN2);
   Fft(fft, x_fft);
   x_fft_ptr = &x_fft[0][0];
 
   // Near fft
   memcpy(fft, aec->dBuf, sizeof(float) * PART_LEN2);
   Fft(fft, df);
 
   // Power smoothing
-  for (i = 0; i < PART_LEN1; i++) {
-    far_spectrum = (x_fft_ptr[i] * x_fft_ptr[i]) +
-                   (x_fft_ptr[PART_LEN1 + i] * x_fft_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);
+  if (aec->refined_adaptive_filter_enabled) {
+    for (i = 0; i < PART_LEN1; ++i) {
+      far_spectrum = (x_fft_ptr[i] * x_fft_ptr[i]) +
+                     (x_fft_ptr[PART_LEN1 + i] * x_fft_ptr[PART_LEN1 + i]);
+      // Calculate absolute spectra
+      abs_far_spectrum[i] = sqrtf(far_spectrum);
+    }
+    RegressorPower(aec->num_partitions, aec->xfBufBlockPos, aec->xfBuf,
+                   aec->xPow);
+  } else {
+    for (i = 0; i < PART_LEN1; ++i) {
+      far_spectrum = (x_fft_ptr[i] * x_fft_ptr[i]) +
+                     (x_fft_ptr[PART_LEN1 + i] * x_fft_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);
+    }
+  }
 
+  for (i = 0; i < PART_LEN1; ++i) {
     near_spectrum = df[0][i] * df[0][i] + df[1][i] * df[1][i];
     aec->dPow[i] = gPow[0] * aec->dPow[i] + gPow[1] * near_spectrum;
     // Calculate absolute spectra
     abs_near_spectrum[i] = sqrtf(near_spectrum);
   }
 
   // Estimate noise power. Wait until dPow is more stable.
   if (aec->noiseEstCtr > 50) {
     for (i = 0; i < PART_LEN1; i++) {
       if (aec->dPow[i] < aec->dMinPow[i]) {
@@ skipping 35 matching lines @@
       }
       if (aec->delay_metrics_delivered == 1 &&
           aec->num_delay_values >= kDelayMetricsAggregationWindow) {
         UpdateDelayMetrics(aec);
       }
     }
   }
 
   // Perform echo subtraction.
   EchoSubtraction(aec, aec->num_partitions, aec->extended_filter_enabled,
-                  aec->normal_mu, aec->normal_error_threshold, &x_fft[0][0],
-                  &aec->xfBufBlockPos, aec->xfBuf, nearend_ptr, aec->xPow,
-                  aec->wfBuf, echo_subtractor_output);
+                  aec->filter_step_size, aec->normal_error_threshold,
+                  &x_fft[0][0], &aec->xfBufBlockPos, aec->xfBuf, nearend_ptr,
+                  aec->xPow, aec->wfBuf, echo_subtractor_output);
 
   RTC_AEC_DEBUG_WAV_WRITE(aec->outLinearFile, echo_subtractor_output, PART_LEN);
 
   if (aec->metricsMode == 1) {
     UpdateLevel(&aec->linoutlevel,
                 CalculatePower(echo_subtractor_output, PART_LEN));
   }
 
   // Perform echo suppression.
   EchoSuppression(aec, farend_ptr, echo_subtractor_output, output, outputH_ptr);
@@ skipping 83 matching lines @@
   aec->delay_agnostic_enabled = 1;  // DA-AEC enabled by default.
   // DA-AEC assumes the system is causal from the beginning and will self adjust
   // the lookahead when shifting is required.
   WebRtc_set_lookahead(aec->delay_estimator, 0);
 #else
   aec->delay_agnostic_enabled = 0;
   WebRtc_set_lookahead(aec->delay_estimator, kLookaheadBlocks);
 #endif
   aec->extended_filter_enabled = 0;
   aec->aec3_enabled = 0;
+  aec->refined_adaptive_filter_enabled = false;
 
   // Assembly optimization
   WebRtcAec_FilterFar = FilterFar;
   WebRtcAec_ScaleErrorSignal = ScaleErrorSignal;
   WebRtcAec_FilterAdaptation = FilterAdaptation;
   WebRtcAec_OverdriveAndSuppress = OverdriveAndSuppress;
   WebRtcAec_ComfortNoise = ComfortNoise;
   WebRtcAec_SubbandCoherence = SubbandCoherence;
   WebRtcAec_StoreAsComplex = StoreAsComplex;
   WebRtcAec_PartitionDelay = PartitionDelay;
@@ skipping 43 matching lines @@
   RTC_AEC_DEBUG_WAV_CLOSE(aec->outFile);
   RTC_AEC_DEBUG_WAV_CLOSE(aec->outLinearFile);
   RTC_AEC_DEBUG_RAW_CLOSE(aec->e_fft_file);
 
   WebRtc_FreeDelayEstimator(aec->delay_estimator);
   WebRtc_FreeDelayEstimatorFarend(aec->delay_estimator_farend);
 
   delete aec;
 }
 
+static void SetAdaptiveFilterStepSize(AecCore* aec) {
+  // Extended filter adaptation parameter.
+  // TODO(ajm): No narrowband tuning yet.
+  const float kExtendedMu = 0.4f;
+
+  if (aec->refined_adaptive_filter_enabled) {
+    aec->filter_step_size = 0.05f;

[hlundin-webrtc, 2016/04/14 14:02:26]

The old code seems to vary the step size depending on sample rate. Is this
motivated for the new code too?

[peah-webrtc, 2016/04/14 22:12:41]

Maybe, 
The only way the sample rate affects the adaptive filter is the way that the
echo path varies. And that variability could motivate both using a higher and
lower rate.
To me the sample-rate dependency would be better handled by having the filter
length (measured in samples) be proportional to the sample rate.

I have no idea how the decision to have the different adaptation step sizes for
the different sample rates were made. Both step sizes are very high.

Probably the step-size could be tuned differently for different sample rates,
but that tuning would require a lot of work and be very specific to specific AEC
implementation at hand. I'd rather limit the use to the new single step-size
value.

+  } else {
+    if (aec->extended_filter_enabled) {
+      aec->filter_step_size = kExtendedMu;
+    } else {
+      if (aec->sampFreq == 8000) {
+        aec->filter_step_size = 0.6f;
+      } else {
+        aec->filter_step_size = 0.5f;
+      }
+    }
+  }
+}
+
 int WebRtcAec_InitAec(AecCore* aec, int sampFreq) {
   int i;
 
   aec->sampFreq = sampFreq;
 
+  SetAdaptiveFilterStepSize(aec);
+
   if (sampFreq == 8000) {
-    aec->normal_mu = 0.6f;
     aec->normal_error_threshold = 2e-6f;
     aec->num_bands = 1;
   } else {
-    aec->normal_mu = 0.5f;
     aec->normal_error_threshold = 1.5e-6f;
     aec->num_bands = (size_t)(sampFreq / 16000);
   }
 
   WebRtc_InitBuffer(aec->nearFrBuf);
   WebRtc_InitBuffer(aec->outFrBuf);
   for (i = 0; i < NUM_HIGH_BANDS_MAX; ++i) {
     WebRtc_InitBuffer(aec->nearFrBufH[i]);
     WebRtc_InitBuffer(aec->outFrBufH[i]);
   }
@@ skipping 351 matching lines @@
 
 void WebRtcAec_enable_aec3(AecCore* self, int enable) {
   self->aec3_enabled = (enable != 0);
 }
 
 int WebRtcAec_aec3_enabled(AecCore* self) {
   assert(self->aec3_enabled == 0 || self->aec3_enabled == 1);
   return self->aec3_enabled;
 }
 
+void WebRtcAec_enable_refined_adaptive_filter(AecCore* self, bool enable) {
+  self->refined_adaptive_filter_enabled = enable;
+  SetAdaptiveFilterStepSize(self);
+}
+
+bool WebRtcAec_refined_adaptive_filter_enabled(AecCore* self) {

[hlundin-webrtc, 2016/04/14 14:02:27]

const AecCore*

[peah-webrtc, 2016/04/14 22:12:41]

Done.

+  return self->refined_adaptive_filter_enabled;
+}
 
 void WebRtcAec_enable_extended_filter(AecCore* self, int enable) {
   self->extended_filter_enabled = enable;
   self->num_partitions = enable ? kExtendedNumPartitions : kNormalNumPartitions;
   // Update the delay estimator with filter length.  See InitAEC() for details.
   WebRtc_set_allowed_offset(self->delay_estimator, self->num_partitions / 2);
 }

[hlundin-webrtc, 2016/04/14 14:47:46]

You must call SetAdaptiveFilterStepSize here too.

[peah-webrtc, 2016/04/14 22:12:41]

Great find!
Done.

 
 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;
 }
 }  // namespace webrtc