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"
 extern "C" {
 #include "webrtc/modules/audio_processing/aec/aec_rdft.h"
 }
 #include "webrtc/modules/audio_processing/logging/aec_logging.h"
@@ skipping 197 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,
-                             float normal_mu,
-                             float normal_error_threshold,
+static void ScaleErrorSignal(float mu,
+                             float 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]);
 
     if (abs_ef > error_threshold) {
       abs_ef = error_threshold / (abs_ef + 1e-10f);
       ef[0][i] *= abs_ef;
@@ skipping 669 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;
+      RTC_DCHECK_EQ(num_partitions * PART_LEN1, x_fft_buf_position);
+      x_fft_buf_position = 0;
+    }
+  }
+}
+
 static void EchoSubtraction(AecCore* aec,
                             int num_partitions,
                             int extended_filter_enabled,
-                            float normal_mu,
-                            float normal_error_threshold,
+                            float filter_step_size,
+                            float 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]) {
   float s_fft[2][PART_LEN1];
@@ skipping 40 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,
-                             normal_error_threshold, x_pow, e_fft);
+  WebRtcAec_ScaleErrorSignal(filter_step_size, 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,
                             float* const* outputH) {
@@ skipping 292 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);
+  // Power smoothing.
+  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 the magnitude spectrum.
+      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 the magnitude spectrum.
+      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
+    // Calculate the magnitude spectrum.
     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]) {
         aec->dMinPow[i] =
             (aec->dPow[i] + step * (aec->dMinPow[i] - aec->dPow[i])) * ramp;
       } else {
@@ skipping 32 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->filter_step_size, aec->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));
   }
 
@@ skipping 85 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;
+  } 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;
+      }
+    }
+  }
+}
+
+static void SetErrorThreshold(AecCore* aec) {
+  // Extended filter adaptation parameter.
+  // TODO(ajm): No narrowband tuning yet.
+  static const float kExtendedErrorThreshold = 1.0e-6f;
+
+  if (aec->extended_filter_enabled) {
+    aec->error_threshold = kExtendedErrorThreshold;
+  } else {
+    if (aec->sampFreq == 8000) {
+      aec->error_threshold = 2e-6f;
+    } else {
+      aec->error_threshold = 1.5e-6f;
+    }
+  }
+}
+
 int WebRtcAec_InitAec(AecCore* aec, int sampFreq) {
   int i;
 
   aec->sampFreq = sampFreq;
 
+  SetAdaptiveFilterStepSize(aec);
+  SetErrorThreshold(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 350 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);
+  SetErrorThreshold(self);
+}
+
+bool WebRtcAec_refined_adaptive_filter_enabled(const AecCore* self) {
+  return self->refined_adaptive_filter_enabled;
+}
 
 void WebRtcAec_enable_extended_filter(AecCore* self, int enable) {
   self->extended_filter_enabled = enable;
+  SetAdaptiveFilterStepSize(self);
+  SetErrorThreshold(self);
   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);
 }
 
 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