AEC3: Increasing the accuracy of the detection for early reverb

This CL introduces an adaptive estimation of the early reverb
in the estimation for the room reverberation. The benefits of
this is that for room with long early reflections there is
a lower risk of underestimating the reverberation.

This CL is for a landing the code in
https://webrtc-review.googlesource.com/c/src/+/87420,
and the review of the code was done in that CL. The author of
code is devicentepena@webrtc.org

Bug: webrtc:9479, chromium:865397
Change-Id: Id6f57e2a684664aef96e8c502e66775f37da59da
Reviewed-on: https://webrtc-review.googlesource.com/91162
Commit-Queue: Per Åhgren <peah@webrtc.org>
Reviewed-by: Sam Zackrisson <saza@webrtc.org>
Cr-Commit-Position: refs/heads/master@{#24146}

modules/audio_processing/aec3/BUILD.gn

@@ -82,6 +82,10 @@ rtc_static_library("aec3") {
     "render_signal_analyzer.h",
     "residual_echo_estimator.cc",
     "residual_echo_estimator.h",
+    "reverb_decay_estimator.cc",
+    "reverb_decay_estimator.h",
+    "reverb_frequency_response.cc",
+    "reverb_frequency_response.h",
     "reverb_model.cc",
     "reverb_model.h",
     "reverb_model_estimator.cc",

modules/audio_processing/aec3/aec3_common.h

@@ -36,6 +36,7 @@ constexpr size_t kFftLengthBy2 = 64;
 constexpr size_t kFftLengthBy2Plus1 = kFftLengthBy2 + 1;
 constexpr size_t kFftLengthBy2Minus1 = kFftLengthBy2 - 1;
 constexpr size_t kFftLength = 2 * kFftLengthBy2;
+constexpr size_t kFftLengthBy2Log2 = 6;
 
 constexpr int kMaxAdaptiveFilterLength = 50;
 constexpr int kRenderTransferQueueSizeFrames = 100;
@@ -44,7 +45,7 @@ constexpr size_t kMaxNumBands = 3;
 constexpr size_t kSubFrameLength = 80;
 
 constexpr size_t kBlockSize = kFftLengthBy2;
-constexpr size_t kBlockSizeLog2 = 6;
+constexpr size_t kBlockSizeLog2 = kFftLengthBy2Log2;
 
 constexpr size_t kExtendedBlockSize = 2 * kFftLengthBy2;
 constexpr size_t kMatchedFilterWindowSizeSubBlocks = 32;
@@ -96,6 +97,9 @@ float Log2TodB(const float in_log2);
 static_assert(1 << kBlockSizeLog2 == kBlockSize,
               "Proper number of shifts for blocksize");
 
+static_assert(1 << kFftLengthBy2Log2 == kFftLengthBy2,
+              "Proper number of shifts for the fft length");
+
 static_assert(1 == NumBandsForRate(8000), "Number of bands for 8 kHz");
 static_assert(1 == NumBandsForRate(16000), "Number of bands for 16 kHz");
 static_assert(2 == NumBandsForRate(32000), "Number of bands for 32 kHz");

modules/audio_processing/aec3/aec_state.cc

@@ -301,14 +301,16 @@ void AecState::Update(
   use_linear_filter_output_ = usable_linear_estimate_ && !TransparentMode();
   diverged_linear_filter_ = diverged_filter;
 
+  const bool stationary_block =
+      use_stationary_properties_ && echo_audibility_.IsBlockStationary();
+
   reverb_model_estimator_.Update(
-      adaptive_filter_impulse_response, adaptive_filter_frequency_response,
+      filter_analyzer_.GetAdjustedFilter(), adaptive_filter_frequency_response,
       erle_estimator_.GetInstLinearQualityEstimate(), filter_delay_blocks_,
-      usable_linear_estimate_, config_.ep_strength.default_len,
+      usable_linear_estimate_, stationary_block);
-      IsBlockStationary());
 
   erle_estimator_.Dump(data_dumper_);
-  reverb_model_estimator_.Dump(data_dumper_);
+  reverb_model_estimator_.Dump(data_dumper_.get());
   data_dumper_->DumpRaw("aec3_erl", Erl());
   data_dumper_->DumpRaw("aec3_erl_time_domain", ErlTimeDomain());
   data_dumper_->DumpRaw("aec3_usable_linear_estimate", UsableLinearEstimate());
@@ -338,8 +340,8 @@ void AecState::Update(
                         recently_converged_filter);
   data_dumper_->DumpRaw("aec3_suppresion_gain_limiter_running",
                         IsSuppressionGainLimitActive());
-  data_dumper_->DumpRaw("aec3_filter_tail_freq_resp_est", GetFreqRespTail());
+  data_dumper_->DumpRaw("aec3_filter_tail_freq_resp_est",
-
+                        GetReverbFrequencyResponse());
 }
 
 bool AecState::DetectActiveRender(rtc::ArrayView<const float> x) const {

modules/audio_processing/aec3/aec_state.h

@@ -67,17 +67,6 @@ class AecState {
   // aec.
   bool UseStationaryProperties() const { return use_stationary_properties_; }
 
-  // Returns true if the current render block is estimated as stationary.
-  bool IsBlockStationary() const {
-    if (UseStationaryProperties()) {
-      return echo_audibility_.IsBlockStationary();
-    } else {
-      // Assume that a non stationary block when the use of
-      // stationary properties are not enabled.
-      return false;
-    }
-  }
-
   // Returns the ERLE.
   const std::array<float, kFftLengthBy2Plus1>& Erle() const {
     return erle_estimator_.Erle();
@@ -130,6 +119,11 @@ class AecState {
   // Returns the decay factor for the echo reverberation.
   float ReverbDecay() const { return reverb_model_estimator_.ReverbDecay(); }
 
+  // Return the frequency response of the reverberant echo.
+  rtc::ArrayView<const float> GetReverbFrequencyResponse() const {
+    return reverb_model_estimator_.GetReverbFrequencyResponse();
+  }
+
   // Returns the upper limit for the echo suppression gain.
   float SuppressionGainLimit() const {
     return suppression_gain_limiter_.Limit();
@@ -159,11 +153,6 @@ class AecState {
               const SubtractorOutput& subtractor_output,
               rtc::ArrayView<const float> y);
 
-  // Returns the tail freq. response of the linear filter.
-  rtc::ArrayView<const float> GetFreqRespTail() const {
-    return reverb_model_estimator_.GetFreqRespTail();
-  }
-
   // Returns filter length in blocks.
   int FilterLengthBlocks() const {
     return filter_analyzer_.FilterLengthBlocks();

modules/audio_processing/aec3/filter_analyzer.cc

@@ -95,10 +95,9 @@ void FilterAnalyzer::Update(
     const RenderBuffer& render_buffer) {
   // Preprocess the filter to avoid issues with low-frequency components in the
   // filter.
-  if (use_preprocessed_filter_) {
+  PreProcessFilter(filter_time_domain);
-    PreProcessFilter(filter_time_domain);
+  data_dumper_->DumpRaw("aec3_linear_filter_processed_td", h_highpass_);
-    data_dumper_->DumpRaw("aec3_linear_filter_processed_td", h_highpass_);
+
-  }
   const auto& filter_to_analyze =
       use_preprocessed_filter_ ? h_highpass_ : filter_time_domain;
   RTC_DCHECK_EQ(filter_to_analyze.size(), filter_time_domain.size());

modules/audio_processing/aec3/filter_analyzer.h

@@ -54,16 +54,14 @@ class FilterAnalyzer {
   // Returns the number of blocks for the current used filter.
   float FilterLengthBlocks() const { return filter_length_blocks_; }
 
+  // Returns the preprocessed filter.
+  rtc::ArrayView<const float> GetAdjustedFilter() const { return h_highpass_; }
+
  private:
   void UpdateFilterGain(rtc::ArrayView<const float> filter_time_domain,
                         size_t max_index);
   void PreProcessFilter(rtc::ArrayView<const float> filter_time_domain);
 
-  // Updates the estimation of the frequency response at the filter tails.
-  void UpdateFreqRespTail(
-      const std::vector<std::array<float, kFftLengthBy2Plus1>>&
-          filter_freq_response);
-
   static int instance_count_;
   std::unique_ptr<ApmDataDumper> data_dumper_;
   const bool use_preprocessed_filter_;

modules/audio_processing/aec3/residual_echo_estimator.cc

@@ -113,7 +113,7 @@ void ResidualEchoEstimator::Estimate(
     if (echo_reverb_) {
       echo_reverb_->AddReverb(
           render_buffer.Spectrum(aec_state.FilterLengthBlocks() + 1),
-          aec_state.GetFreqRespTail(), aec_state.ReverbDecay(), *R2);
+          aec_state.GetReverbFrequencyResponse(), aec_state.ReverbDecay(), *R2);
 
     } else {
       RTC_DCHECK(echo_reverb_fallback);

modules/audio_processing/aec3/reverb_decay_estimator.cc

@@ -0,0 +1,361 @@
+/*
+ *  Copyright (c) 2018 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.
+ */
+
+#include "modules/audio_processing/aec3/reverb_decay_estimator.h"
+
+#include <algorithm>
+#include <cmath>
+#include <numeric>
+
+#include "api/array_view.h"
+#include "modules/audio_processing/logging/apm_data_dumper.h"
+#include "rtc_base/checks.h"
+#include "system_wrappers/include/field_trial.h"
+
+namespace webrtc {
+
+namespace {
+
+bool EnforceAdaptiveEchoReverbEstimation() {
+  return field_trial::IsEnabled(
+      "WebRTC-Aec3EnableAdaptiveEchoReverbEstimation");
+}
+
+constexpr int kEarlyReverbMinSizeBlocks = 3;
+constexpr int kBlocksPerSection = 3;
+
+// Averages the values in a block of size kFftLengthBy2;
+float BlockAverage(rtc::ArrayView<const float> v, size_t block_index) {
+  constexpr float kOneByFftLengthBy2 = 1.f / kFftLengthBy2;
+  const int i = block_index * kFftLengthBy2;
+  RTC_DCHECK_GE(v.size(), i + kFftLengthBy2);
+  const float sum =
+      std::accumulate(v.begin() + i, v.begin() + i + kFftLengthBy2, 0.f);
+  return sum * kOneByFftLengthBy2;
+}
+
+// Analyzes the gain in a block.
+void AnalyzeBlockGain(const std::array<float, kFftLengthBy2>& h2,
+                      float floor_gain,
+                      float* previous_gain,
+                      bool* block_adapting,
+                      bool* decaying_gain) {
+  float gain = std::max(BlockAverage(h2, 0), 1e-32f);
+  *block_adapting =
+      *previous_gain > 1.1f * gain || *previous_gain < 0.9f * gain;
+  *decaying_gain = gain > floor_gain;
+  *previous_gain = gain;
+}
+
+// Arithmetic sum of $2 \sum_{i=0.5}^{(N-1)/2}i^2$ calculated directly.
+constexpr float SymmetricArithmetricSum(int N) {
+  return N * (N * N - 1.0f) * (1.f / 12.f);
+}
+
+}  // namespace
+
+ReverbDecayEstimator::ReverbDecayEstimator(const EchoCanceller3Config& config)
+    : filter_length_blocks_(config.filter.main.length_blocks),
+      filter_length_coefficients_(GetTimeDomainLength(filter_length_blocks_)),
+      use_adaptive_echo_decay_(config.ep_strength.default_len < 0.f ||
+                               EnforceAdaptiveEchoReverbEstimation()),
+      early_reverb_estimator_(config.filter.main.length_blocks -
+                              kEarlyReverbMinSizeBlocks),
+      late_reverb_start_(kEarlyReverbMinSizeBlocks),
+      late_reverb_end_(kEarlyReverbMinSizeBlocks),
+      decay_(std::fabs(config.ep_strength.default_len)) {
+  previous_gains_.fill(0.f);
+  RTC_DCHECK_GT(config.filter.main.length_blocks,
+                static_cast<size_t>(kEarlyReverbMinSizeBlocks));
+}
+
+ReverbDecayEstimator::~ReverbDecayEstimator() = default;
+
+void ReverbDecayEstimator::Update(rtc::ArrayView<const float> filter,
+                                  const absl::optional<float>& filter_quality,
+                                  int filter_delay_blocks,
+                                  bool usable_linear_filter,
+                                  bool stationary_signal) {
+  const int filter_size = static_cast<int>(filter.size());
+
+  if (stationary_signal) {
+    return;
+  }
+
+  // TODO(devicentepena): Verify that the below is correct.
+  bool estimation_feasible =
+      filter_delay_blocks <=
+      filter_length_blocks_ - kEarlyReverbMinSizeBlocks - 1;
+  estimation_feasible =
+      estimation_feasible && filter_size == filter_length_coefficients_;
+  estimation_feasible = estimation_feasible && filter_delay_blocks > 0;
+  estimation_feasible = estimation_feasible && usable_linear_filter;
+
+  if (!estimation_feasible) {
+    ResetDecayEstimation();
+    return;
+  }
+
+  if (!use_adaptive_echo_decay_) {
+    return;
+  }
+
+  const float new_smoothing = filter_quality ? *filter_quality * 0.2f : 0.f;
+  smoothing_constant_ = std::max(new_smoothing, smoothing_constant_);
+  if (smoothing_constant_ == 0.f) {
+    return;
+  }
+
+  if (block_to_analyze_ < filter_length_blocks_) {
+    // Analyze the filter and accumulate data for reverb estimation.
+    AnalyzeFilter(filter);
+    ++block_to_analyze_;
+  } else {
+    // When the filter is fully analyzed, estimate the reverb decay and reset
+    // the block_to_analyze_ counter.
+    EstimateDecay(filter);
+  }
+}
+
+void ReverbDecayEstimator::ResetDecayEstimation() {
+  early_reverb_estimator_.Reset();
+  late_reverb_decay_estimator_.Reset(0);
+  block_to_analyze_ = 0;
+  estimation_region_candidate_size_ = 0;
+  estimation_region_identified_ = false;
+  smoothing_constant_ = 0.f;
+  late_reverb_start_ = 0;
+  late_reverb_end_ = 0;
+}
+
+void ReverbDecayEstimator::EstimateDecay(rtc::ArrayView<const float> filter) {
+  auto& h = filter;
+
+  RTC_DCHECK_EQ(0, h.size() % kFftLengthBy2);
+
+  // Compute the squared filter coefficients.
+  std::array<float, GetTimeDomainLength(kMaxAdaptiveFilterLength)> h2_data;
+  RTC_DCHECK_GE(h2_data.size(), filter_length_coefficients_);
+  rtc::ArrayView<float> h2(h2_data.data(), filter_length_coefficients_);
+  std::transform(h.begin(), h.end(), h2.begin(), [](float a) { return a * a; });
+
+  // Identify the peak index of the filter.
+  const int peak_coefficient =
+      std::distance(h2.begin(), std::max_element(h2.begin(), h2.end()));
+  int peak_block = peak_coefficient >> kFftLengthBy2Log2;
+
+  // Reset the block analysis counter.
+  block_to_analyze_ =
+      std::min(peak_block + kEarlyReverbMinSizeBlocks, filter_length_blocks_);
+
+  // To estimate the reverb decay, the energy of the first filter section must
+  // be substantially larger than the last. Also, the first filter section
+  // energy must not deviate too much from the max peak.
+  const float first_reverb_gain = BlockAverage(h2, block_to_analyze_);
+  tail_gain_ = BlockAverage(h2, (h2.size() >> kFftLengthBy2Log2) - 1);
+  const bool sufficient_reverb_decay = first_reverb_gain > 4.f * tail_gain_;
+  const bool valid_filter =
+      first_reverb_gain > 2.f * tail_gain_ && h2[peak_coefficient] < 100.f;
+
+  // Estimate the size of the regions with early and late reflections.
+  const int size_early_reverb = early_reverb_estimator_.Estimate();
+  const int size_late_reverb =
+      std::max(estimation_region_candidate_size_ - size_early_reverb, 0);
+
+  // Only update the reverb decay estimate if the size of the identified late
+  // reverb is sufficiently large.
+  if (size_late_reverb >= 5) {
+    if (valid_filter && late_reverb_decay_estimator_.EstimateAvailable()) {
+      float decay = std::pow(
+          2.0f, late_reverb_decay_estimator_.Estimate() * kFftLengthBy2);
+      constexpr float kMaxDecay = 0.95f;  // ~1 sec min RT60.
+      constexpr float kMinDecay = 0.02f;  // ~15 ms max RT60.
+      decay = std::max(.97f * decay_, decay);
+      decay = std::min(decay, kMaxDecay);
+      decay = std::max(decay, kMinDecay);
+      decay_ += smoothing_constant_ * (decay - decay_);
+    }
+
+    // Update length of decay. Must have enough data (number of sections) in
+    // order to estimate decay rate.
+    late_reverb_decay_estimator_.Reset(size_late_reverb * kFftLengthBy2);
+    late_reverb_start_ =
+        peak_block + kEarlyReverbMinSizeBlocks + size_early_reverb;
+    late_reverb_end_ =
+        block_to_analyze_ + estimation_region_candidate_size_ - 1;
+  } else {
+    late_reverb_decay_estimator_.Reset(0);
+    late_reverb_start_ = 0;
+    late_reverb_end_ = 0;
+  }
+
+  // Reset variables for the identification of the region for reverb decay
+  // estimation.
+  estimation_region_identified_ = !(valid_filter && sufficient_reverb_decay);
+  estimation_region_candidate_size_ = 0;
+
+  // Stop estimation of the decay until another good filter is received.
+  smoothing_constant_ = 0.f;
+
+  // Reset early reflections detector.
+  early_reverb_estimator_.Reset();
+}
+
+void ReverbDecayEstimator::AnalyzeFilter(rtc::ArrayView<const float> filter) {
+  auto h = rtc::ArrayView<const float>(
+      filter.begin() + block_to_analyze_ * kFftLengthBy2, kFftLengthBy2);
+
+  // Compute squared filter coeffiecients for the block to analyze_;
+  std::array<float, kFftLengthBy2> h2;
+  std::transform(h.begin(), h.end(), h2.begin(), [](float a) { return a * a; });
+
+  // Map out the region for estimating the reverb decay.
+  bool adapting;
+  bool above_noise_floor;
+  AnalyzeBlockGain(h2, tail_gain_, &previous_gains_[block_to_analyze_],
+                   &adapting, &above_noise_floor);
+
+  // Count consecutive number of "good" filter sections, where "good" means:
+  // 1) energy is above noise floor.
+  // 2) energy of current section has not changed too much from last check.
+  estimation_region_identified_ =
+      estimation_region_identified_ || adapting || !above_noise_floor;
+  if (!estimation_region_identified_) {
+    ++estimation_region_candidate_size_;
+  }
+
+  // Accumulate data for reverb decay estimation and for the estimation of early
+  // reflections.
+  if (block_to_analyze_ <= late_reverb_end_) {
+    if (block_to_analyze_ >= late_reverb_start_) {
+      for (float h2_k : h2) {
+        float h2_log2 = FastApproxLog2f(h2_k + 1e-10);
+        late_reverb_decay_estimator_.Accumulate(h2_log2);
+        early_reverb_estimator_.Accumulate(h2_log2, smoothing_constant_);
+      }
+    } else {
+      for (float h2_k : h2) {
+        float h2_log2 = FastApproxLog2f(h2_k + 1e-10);
+        early_reverb_estimator_.Accumulate(h2_log2, smoothing_constant_);
+      }
+    }
+  }
+}
+
+void ReverbDecayEstimator::Dump(ApmDataDumper* data_dumper) const {
+  data_dumper->DumpRaw("aec3_reverb_decay", decay_);
+  data_dumper->DumpRaw("aec3_reverb_tail_energy", tail_gain_);
+  data_dumper->DumpRaw("aec3_reverb_alpha", smoothing_constant_);
+  data_dumper->DumpRaw("aec3_num_reverb_decay_blocks",
+                       late_reverb_end_ - late_reverb_start_);
+  data_dumper->DumpRaw("aec3_blocks_after_early_reflections",
+                       late_reverb_start_);
+  early_reverb_estimator_.Dump(data_dumper);
+}
+
+void ReverbDecayEstimator::LateReverbLinearRegressor::Reset(
+    int num_data_points) {
+  RTC_DCHECK_LE(0, num_data_points);
+  RTC_DCHECK_EQ(0, num_data_points % 2);
+  const int N = num_data_points;
+  nz_ = 0.f;
+  // Arithmetic sum of $2 \sum_{i=0.5}^{(N-1)/2}i^2$ calculated directly.
+  nn_ = SymmetricArithmetricSum(N);
+  // The linear regression approach assumes symmetric index around 0.
+  count_ = N > 0 ? count_ = -N * 0.5f + 0.5f : 0.f;
+  N_ = N;
+  n_ = 0;
+}
+
+void ReverbDecayEstimator::LateReverbLinearRegressor::Accumulate(float z) {
+  nz_ += count_ * z;
+  ++count_;
+  ++n_;
+}
+
+float ReverbDecayEstimator::LateReverbLinearRegressor::Estimate() {
+  RTC_DCHECK(EstimateAvailable());
+  if (nn_ == 0.f) {
+    RTC_NOTREACHED();
+    return 0.f;
+  }
+  return nz_ / nn_;
+}
+
+ReverbDecayEstimator::EarlyReverbLengthEstimator::EarlyReverbLengthEstimator(
+    int max_blocks)
+    : numerators_(1 + max_blocks / kBlocksPerSection, 0.f),
+      nz_(numerators_.size(), 0.f),
+      count_(numerators_.size(), 0.f) {
+  RTC_DCHECK_LE(0, max_blocks);
+}
+
+ReverbDecayEstimator::EarlyReverbLengthEstimator::
+    ~EarlyReverbLengthEstimator() = default;
+
+void ReverbDecayEstimator::EarlyReverbLengthEstimator::Reset() {
+  // Linear regression approach assumes symmetric index around 0.
+  constexpr float kCount = -0.5f * kBlocksPerSection * kFftLengthBy2 + 0.5f;
+  std::fill(count_.begin(), count_.end(), kCount);
+  std::fill(nz_.begin(), nz_.end(), 0.f);
+  section_ = 0;
+  section_update_counter_ = 0;
+}
+
+void ReverbDecayEstimator::EarlyReverbLengthEstimator::Accumulate(
+    float value,
+    float smoothing) {
+  nz_[section_] += count_[section_] * value;
+  ++count_[section_];
+
+  if (++section_update_counter_ == kBlocksPerSection * kFftLengthBy2) {
+    RTC_DCHECK_GT(nz_.size(), section_);
+    RTC_DCHECK_GT(numerators_.size(), section_);
+    numerators_[section_] +=
+        smoothing * (nz_[section_] - numerators_[section_]);
+    section_update_counter_ = 0;
+    ++section_;
+  }
+}
+
+int ReverbDecayEstimator::EarlyReverbLengthEstimator::Estimate() {
+  constexpr float N = kBlocksPerSection * kFftLengthBy2;
+  constexpr float nn = SymmetricArithmetricSum(N);
+  // numerator_11 refers to the quantity that the linear regressor needs in the
+  // numerator for getting a decay equal to 1.1 (which is not a decay).
+  // log2(1.1) * nn / kFftLengthBy2.
+  constexpr float numerator_11 = 0.13750352374993502f * nn / kFftLengthBy2;
+  // log2(0.8) *  nn / kFftLengthBy2.
+  constexpr float numerator_08 = -0.32192809488736229f * nn / kFftLengthBy2;
+  constexpr int kNumSectionsToAnalyze = 3;
+
+  // Analyze the first kNumSectionsToAnalyze regions.
+  // TODO(devicentepena): Add a more thorough comment for explaining the logic
+  // below.
+  const float min_stable_region = *std::min_element(
+      numerators_.begin() + kNumSectionsToAnalyze, numerators_.end());
+  int early_reverb_size = 0;
+  for (int k = 0; k < kNumSectionsToAnalyze; ++k) {
+    if ((numerators_[k] > numerator_11) ||
+        (numerators_[k] < numerator_08 &&
+         numerators_[k] < 0.9f * min_stable_region)) {
+      early_reverb_size = (k + 1) * kBlocksPerSection;
+    }
+  }
+
+  return early_reverb_size;
+}
+
+void ReverbDecayEstimator::EarlyReverbLengthEstimator::Dump(
+    ApmDataDumper* data_dumper) const {
+  data_dumper->DumpRaw("aec3_er_acum_numerator", numerators_);
+}
+
+}  // namespace webrtc

modules/audio_processing/aec3/reverb_decay_estimator.h

@@ -0,0 +1,110 @@
+/*
+ *  Copyright (c) 2018 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.
+ */
+
+#ifndef MODULES_AUDIO_PROCESSING_AEC3_REVERB_DECAY_ESTIMATOR_H_
+#define MODULES_AUDIO_PROCESSING_AEC3_REVERB_DECAY_ESTIMATOR_H_
+
+#include <array>
+#include <vector>
+
+#include "absl/types/optional.h"
+#include "api/array_view.h"
+#include "api/audio/echo_canceller3_config.h"
+#include "modules/audio_processing/aec3/aec3_common.h"
+
+namespace webrtc {
+
+class ApmDataDumper;
+
+// Class for estimating the decay of the late reverb.
+class ReverbDecayEstimator {
+ public:
+  explicit ReverbDecayEstimator(const EchoCanceller3Config& config);
+  ~ReverbDecayEstimator();
+  // Updates the decay estimate.
+  void Update(rtc::ArrayView<const float> filter,
+              const absl::optional<float>& filter_quality,
+              int filter_delay_blocks,
+              bool usable_linear_filter,
+              bool stationary_signal);
+  // Returns the decay for the exponential model.
+  float Decay() const { return decay_; }
+  // Dumps debug data.
+  void Dump(ApmDataDumper* data_dumper) const;
+
+ private:
+  void EstimateDecay(rtc::ArrayView<const float> filter);
+  void AnalyzeFilter(rtc::ArrayView<const float> filter);
+
+  void ResetDecayEstimation();
+
+  // Class for estimating the decay of the late reverb from the linear filter.
+  class LateReverbLinearRegressor {
+   public:
+    // Resets the estimator to receive a specified number of data points.
+    void Reset(int num_data_points);
+    // Accumulates estimation data.
+    void Accumulate(float z);
+    // Estimates the decay.
+    float Estimate();
+    // Returns whether an estimate is available.
+    bool EstimateAvailable() const { return n_ == N_ && N_ != 0; }
+
+   public:
+    float nz_ = 0.f;
+    float nn_ = 0.f;
+    float count_ = 0.f;
+    int N_ = 0;
+    int n_ = 0;
+  };
+
+  // Class for identifying the length of the early reverb from the linear
+  // filter.
+  class EarlyReverbLengthEstimator {
+   public:
+    explicit EarlyReverbLengthEstimator(int max_blocks);
+    ~EarlyReverbLengthEstimator();
+
+    // Resets the estimator.
+    void Reset();
+    // Accumulates estimation data.
+    void Accumulate(float value, float smoothing);
+    // Estimates the size in blocks of the early reverb.
+    int Estimate();
+    // Dumps debug data.
+    void Dump(ApmDataDumper* data_dumper) const;
+
+   private:
+    std::vector<float> numerators_;
+    std::vector<float> nz_;
+    std::vector<float> count_;
+    int section_ = 0;
+    int section_update_counter_ = 0;
+  };
+
+  const int filter_length_blocks_;
+  const int filter_length_coefficients_;
+  const bool use_adaptive_echo_decay_;
+  LateReverbLinearRegressor late_reverb_decay_estimator_;
+  EarlyReverbLengthEstimator early_reverb_estimator_;
+  int late_reverb_start_;
+  int late_reverb_end_;
+  int block_to_analyze_ = 0;
+  int estimation_region_candidate_size_ = 0;
+  bool estimation_region_identified_ = false;
+  std::array<float, kMaxAdaptiveFilterLength> previous_gains_;
+  float decay_;
+  float tail_gain_ = 0.f;
+  float smoothing_constant_ = 0.f;
+};
+
+}  // namespace webrtc
+
+#endif  // MODULES_AUDIO_PROCESSING_AEC3_REVERB_DECAY_ESTIMATOR_H_

modules/audio_processing/aec3/reverb_frequency_response.cc

@@ -0,0 +1,113 @@
+/*
+ *  Copyright (c) 2018 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.
+ */
+
+#include "modules/audio_processing/aec3/reverb_frequency_response.h"
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <memory>
+#include <numeric>
+
+#include "api/array_view.h"
+#include "api/audio/echo_canceller3_config.h"
+#include "modules/audio_processing/aec3/aec3_common.h"
+#include "rtc_base/checks.h"
+#include "system_wrappers/include/field_trial.h"
+
+namespace webrtc {
+
+namespace {
+
+bool EnableSmoothUpdatesTailFreqResp() {
+  return !field_trial::IsEnabled(
+      "WebRTC-Aec3SmoothUpdatesTailFreqRespKillSwitch");
+}
+
+// Computes the ratio of the energies between the direct path and the tail. The
+// energy is computed in the power spectrum domain discarding the DC
+// contributions.
+float AverageDecayWithinFilter(
+    rtc::ArrayView<const float> freq_resp_direct_path,
+    rtc::ArrayView<const float> freq_resp_tail) {
+  // Skipping the DC for the ratio computation
+  constexpr size_t kSkipBins = 1;
+  RTC_CHECK_EQ(freq_resp_direct_path.size(), freq_resp_tail.size());
+
+  float direct_path_energy =
+      std::accumulate(freq_resp_direct_path.begin() + kSkipBins,
+                      freq_resp_direct_path.end(), 0.f);
+
+  if (direct_path_energy == 0.f) {
+    return 0.f;
+  }
+
+  float tail_energy = std::accumulate(freq_resp_tail.begin() + kSkipBins,
+                                      freq_resp_tail.end(), 0.f);
+  return tail_energy / direct_path_energy;
+}
+
+}  // namespace
+
+ReverbFrequencyResponse::ReverbFrequencyResponse()
+    : enable_smooth_tail_response_updates_(EnableSmoothUpdatesTailFreqResp()) {
+  tail_response_.fill(0.f);
+}
+ReverbFrequencyResponse::~ReverbFrequencyResponse() = default;
+
+void ReverbFrequencyResponse::Update(
+    const std::vector<std::array<float, kFftLengthBy2Plus1>>&
+        frequency_response,
+    int filter_delay_blocks,
+    const absl::optional<float>& linear_filter_quality,
+    bool stationary_block) {
+  if (!enable_smooth_tail_response_updates_) {
+    Update(frequency_response, filter_delay_blocks, 0.1f);
+    return;
+  }
+
+  if (stationary_block || !linear_filter_quality) {
+    return;
+  }
+
+  Update(frequency_response, filter_delay_blocks, *linear_filter_quality);
+}
+
+void ReverbFrequencyResponse::Update(
+    const std::vector<std::array<float, kFftLengthBy2Plus1>>&
+        frequency_response,
+    int filter_delay_blocks,
+    float linear_filter_quality) {
+  rtc::ArrayView<const float> freq_resp_tail(
+      frequency_response[frequency_response.size() - 1]);
+
+  rtc::ArrayView<const float> freq_resp_direct_path(
+      frequency_response[filter_delay_blocks]);
+
+  float average_decay =
+      AverageDecayWithinFilter(freq_resp_direct_path, freq_resp_tail);
+
+  const float smoothing = 0.2f * linear_filter_quality;
+  average_decay_ += smoothing * (average_decay - average_decay_);
+
+  for (size_t k = 0; k < kFftLengthBy2Plus1; ++k) {
+    tail_response_[k] = freq_resp_direct_path[k] * average_decay_;
+  }
+
+  // TODO(devicentepena): Check if this should be done using a max that weights
+  // both the lower and upper bands equally.
+  for (size_t k = 1; k < kFftLengthBy2; ++k) {
+    const float avg_neighbour =
+        0.5f * (tail_response_[k - 1] + tail_response_[k + 1]);
+    tail_response_[k] = std::max(tail_response_[k], avg_neighbour);
+  }
+}
+
+}  // namespace webrtc

modules/audio_processing/aec3/reverb_frequency_response.h

@@ -0,0 +1,55 @@
+/*
+ *  Copyright (c) 2018 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.
+ */
+
+#ifndef MODULES_AUDIO_PROCESSING_AEC3_REVERB_FREQUENCY_RESPONSE_H_
+#define MODULES_AUDIO_PROCESSING_AEC3_REVERB_FREQUENCY_RESPONSE_H_
+
+#include <memory>
+#include <vector>
+
+#include "absl/types/optional.h"
+#include "api/array_view.h"
+#include "modules/audio_processing/aec3/aec3_common.h"
+#include "modules/audio_processing/logging/apm_data_dumper.h"
+
+namespace webrtc {
+
+// Class for updating the frequency response for the reverb.
+class ReverbFrequencyResponse {
+ public:
+  ReverbFrequencyResponse();
+  ~ReverbFrequencyResponse();
+
+  // Updates the frequency response estimate of the reverb.
+  void Update(const std::vector<std::array<float, kFftLengthBy2Plus1>>&
+                  frequency_response,
+              int filter_delay_blocks,
+              const absl::optional<float>& linear_filter_quality,
+              bool stationary_block);
+
+  // Returns the estimated frequency response for the reverb.
+  rtc::ArrayView<const float> FrequencyResponse() const {
+    return tail_response_;
+  }
+
+ private:
+  void Update(const std::vector<std::array<float, kFftLengthBy2Plus1>>&
+                  frequency_response,
+              int filter_delay_blocks,
+              float linear_filter_quality);
+
+  const bool enable_smooth_tail_response_updates_;
+  float average_decay_ = 0.f;
+  std::array<float, kFftLengthBy2Plus1> tail_response_;
+};
+
+}  // namespace webrtc
+
+#endif  // MODULES_AUDIO_PROCESSING_AEC3_REVERB_FREQUENCY_RESPONSE_H_

modules/audio_processing/aec3/reverb_model_estimator.cc

@@ -10,313 +10,28 @@
 
 #include "modules/audio_processing/aec3/reverb_model_estimator.h"
 
-#include <algorithm>
-#include <array>
-#include <cmath>
-#include <memory>
-#include <numeric>
-
-#include "api/array_view.h"
-#include "api/audio/echo_canceller3_config.h"
-#include "modules/audio_processing/aec3/aec3_common.h"
-#include "rtc_base/checks.h"
-#include "system_wrappers/include/field_trial.h"
-
 namespace webrtc {
 
-namespace {
-
-bool EnableSmoothUpdatesTailFreqResp() {
-  return !field_trial::IsEnabled(
-      "WebRTC-Aec3SmoothUpdatesTailFreqRespKillSwitch");
-}
-
-// Computes the ratio of the energies between the direct path and the tail. The
-// energy is computed in the power spectrum domain discarding the DC
-// contributions.
-float ComputeRatioEnergies(
-    const rtc::ArrayView<const float>& freq_resp_direct_path,
-    const rtc::ArrayView<const float>& freq_resp_tail) {
-  // Skipping the DC for the ratio computation
-  constexpr size_t n_skip_bins = 1;
-  RTC_CHECK_EQ(freq_resp_direct_path.size(), freq_resp_tail.size());
-
-  float direct_path_energy =
-      std::accumulate(freq_resp_direct_path.begin() + n_skip_bins,
-                      freq_resp_direct_path.end(), 0.f);
-
-  float tail_energy = std::accumulate(freq_resp_tail.begin() + n_skip_bins,
-                                      freq_resp_tail.end(), 0.f);
-
-  if (direct_path_energy > 0) {
-    return tail_energy / direct_path_energy;
-  } else {
-    return 0.f;
-  }
-}
-
-}  // namespace
-
 ReverbModelEstimator::ReverbModelEstimator(const EchoCanceller3Config& config)
-    : filter_main_length_blocks_(config.filter.main.length_blocks),
+    : reverb_decay_estimator_(config) {}
-      reverb_decay_(std::fabs(config.ep_strength.default_len)),
-      enable_smooth_freq_resp_tail_updates_(EnableSmoothUpdatesTailFreqResp()) {
-  block_energies_.fill(0.f);
-  freq_resp_tail_.fill(0.f);
-}
 
 ReverbModelEstimator::~ReverbModelEstimator() = default;
 
-bool ReverbModelEstimator::IsAGoodFilterForDecayEstimation(
-    int filter_delay_blocks,
-    bool usable_linear_estimate,
-    size_t length_filter) {
-  if ((filter_delay_blocks && usable_linear_estimate) &&
-      (filter_delay_blocks <=
-       static_cast<int>(filter_main_length_blocks_) - 4) &&
-      (length_filter >=
-       static_cast<size_t>(GetTimeDomainLength(filter_main_length_blocks_)))) {
-    return true;
-  } else {
-    return false;
-  }
-}
-
 void ReverbModelEstimator::Update(
-    const std::vector<float>& impulse_response,
+    rtc::ArrayView<const float> impulse_response,
     const std::vector<std::array<float, kFftLengthBy2Plus1>>&
-        filter_freq_response,
+        frequency_response,
-    const absl::optional<float>& quality_linear,
+    const absl::optional<float>& linear_filter_quality,
     int filter_delay_blocks,
     bool usable_linear_estimate,
-    float default_decay,
     bool stationary_block) {
-  if (enable_smooth_freq_resp_tail_updates_) {
+  // Estimate the frequency response for the reverb.
-    if (!stationary_block) {
+  reverb_frequency_response_.Update(frequency_response, filter_delay_blocks,
-      float alpha = 0;
+                                    linear_filter_quality, stationary_block);
-      if (quality_linear) {
+
-        alpha = 0.2f * quality_linear.value();
+  // Estimate the reverb decay,
-        UpdateFreqRespTail(filter_freq_response, filter_delay_blocks, alpha);
+  reverb_decay_estimator_.Update(impulse_response, linear_filter_quality,
-      }
+                                 filter_delay_blocks, usable_linear_estimate,
-      if (IsAGoodFilterForDecayEstimation(filter_delay_blocks,
+                                 stationary_block);
-                                          usable_linear_estimate,
-                                          impulse_response.size())) {
-        alpha_ = std::max(alpha, alpha_);
-        if ((alpha_ > 0.f) && (default_decay < 0.f)) {
-          // Echo tail decay estimation if default_decay is negative.
-          UpdateReverbDecay(impulse_response);
-        }
-      } else {
-        ResetDecayEstimation();
-      }
-    }
-  } else {
-    UpdateFreqRespTail(filter_freq_response, filter_delay_blocks, 0.1f);
-  }
 }
-
-void ReverbModelEstimator::ResetDecayEstimation() {
-  accumulated_nz_ = 0.f;
-  accumulated_nn_ = 0.f;
-  accumulated_count_ = 0.f;
-  current_reverb_decay_section_ = 0;
-  num_reverb_decay_sections_ = 0;
-  num_reverb_decay_sections_next_ = 0;
-  found_end_of_reverb_decay_ = false;
-  alpha_ = 0.f;
-}
-
-void ReverbModelEstimator::UpdateReverbDecay(
-    const std::vector<float>& impulse_response) {
-  constexpr float kOneByFftLengthBy2 = 1.f / kFftLengthBy2;
-
-  // Form the data to match against by squaring the impulse response
-  // coefficients.
-  std::array<float, GetTimeDomainLength(kMaxAdaptiveFilterLength)>
-      matching_data_data;
-  RTC_DCHECK_LE(GetTimeDomainLength(filter_main_length_blocks_),
-                matching_data_data.size());
-  rtc::ArrayView<float> matching_data(
-      matching_data_data.data(),
-      GetTimeDomainLength(filter_main_length_blocks_));
-  std::transform(
-      impulse_response.begin(), impulse_response.end(), matching_data.begin(),
-      [](float a) { return a * a; });  // TODO(devicentepena) check if focusing
-                                       // on one block would be enough.
-
-  if (current_reverb_decay_section_ < filter_main_length_blocks_) {
-    // Update accumulated variables for the current filter section.
-
-    const size_t start_index = current_reverb_decay_section_ * kFftLengthBy2;
-
-    RTC_DCHECK_GT(matching_data.size(), start_index);
-    RTC_DCHECK_GE(matching_data.size(), start_index + kFftLengthBy2);
-    float section_energy =
-        std::accumulate(matching_data.begin() + start_index,
-                        matching_data.begin() + start_index + kFftLengthBy2,
-                        0.f) *
-        kOneByFftLengthBy2;
-
-    section_energy = std::max(
-        section_energy, 1e-32f);  // Regularization to avoid division by 0.
-
-    RTC_DCHECK_LT(current_reverb_decay_section_, block_energies_.size());
-    const float energy_ratio =
-        block_energies_[current_reverb_decay_section_] / section_energy;
-
-    found_end_of_reverb_decay_ = found_end_of_reverb_decay_ ||
-                                 (energy_ratio > 1.1f || energy_ratio < 0.9f);
-
-    // Count consecutive number of "good" filter sections, where "good" means:
-    // 1) energy is above noise floor.
-    // 2) energy of current section has not changed too much from last check.
-    if (!found_end_of_reverb_decay_ && section_energy > tail_energy_) {
-      ++num_reverb_decay_sections_next_;
-    } else {
-      found_end_of_reverb_decay_ = true;
-    }
-
-    block_energies_[current_reverb_decay_section_] = section_energy;
-
-    if (num_reverb_decay_sections_ > 0) {
-      // Linear regression of log squared magnitude of impulse response.
-      for (size_t i = 0; i < kFftLengthBy2; i++) {
-        RTC_DCHECK_GT(matching_data.size(), start_index + i);
-        float z = FastApproxLog2f(matching_data[start_index + i] + 1e-10);
-        accumulated_nz_ += accumulated_count_ * z;
-        ++accumulated_count_;
-      }
-    }
-
-    num_reverb_decay_sections_ =
-        num_reverb_decay_sections_ > 0 ? num_reverb_decay_sections_ - 1 : 0;
-    ++current_reverb_decay_section_;
-
-  } else {
-    constexpr float kMaxDecay = 0.95f;  // ~1 sec min RT60.
-    constexpr float kMinDecay = 0.02f;  // ~15 ms max RT60.
-
-    // Accumulated variables throughout whole filter.
-
-    // Solve for decay rate.
-
-    float decay = reverb_decay_;
-
-    if (accumulated_nn_ != 0.f) {
-      const float exp_candidate = -accumulated_nz_ / accumulated_nn_;
-      decay = std::pow(2.0f, -exp_candidate * kFftLengthBy2);
-      decay = std::min(decay, kMaxDecay);
-      decay = std::max(decay, kMinDecay);
-    }
-
-    // Filter tail energy (assumed to be noise).
-    constexpr size_t kTailLength = kFftLengthBy2;
-
-    constexpr float k1ByTailLength = 1.f / kTailLength;
-    const size_t tail_index =
-        GetTimeDomainLength(filter_main_length_blocks_) - kTailLength;
-
-    RTC_DCHECK_GT(matching_data.size(), tail_index);
-
-    tail_energy_ = std::accumulate(matching_data.begin() + tail_index,
-                                   matching_data.end(), 0.f) *
-                   k1ByTailLength;
-
-    // Update length of decay.
-    num_reverb_decay_sections_ = num_reverb_decay_sections_next_;
-    num_reverb_decay_sections_next_ = 0;
-    // Must have enough data (number of sections) in order
-    // to estimate decay rate.
-    if (num_reverb_decay_sections_ < 5) {
-      num_reverb_decay_sections_ = 0;
-    }
-
-    const float N = num_reverb_decay_sections_ * kFftLengthBy2;
-    accumulated_nz_ = 0.f;
-    const float k1By12 = 1.f / 12.f;
-    // Arithmetic sum $2 \sum_{i=0.5}^{(N-1)/2}i^2$ calculated directly.
-    accumulated_nn_ = N * (N * N - 1.0f) * k1By12;
-    accumulated_count_ = -N * 0.5f;
-    // Linear regression approach assumes symmetric index around 0.
-    accumulated_count_ += 0.5f;
-
-    // Identify the peak index of the impulse response.
-    const size_t peak_index = std::distance(
-        matching_data.begin(),
-        std::max_element(matching_data.begin(), matching_data.end()));
-
-    current_reverb_decay_section_ = peak_index * kOneByFftLengthBy2 + 3;
-    // Make sure we're not out of bounds.
-    if (current_reverb_decay_section_ + 1 >= filter_main_length_blocks_) {
-      current_reverb_decay_section_ = filter_main_length_blocks_;
-    }
-    size_t start_index = current_reverb_decay_section_ * kFftLengthBy2;
-    float first_section_energy =
-        std::accumulate(matching_data.begin() + start_index,
-                        matching_data.begin() + start_index + kFftLengthBy2,
-                        0.f) *
-        kOneByFftLengthBy2;
-
-    // To estimate the reverb decay, the energy of the first filter section
-    // must be substantially larger than the last.
-    // Also, the first filter section energy must not deviate too much
-    // from the max peak.
-    bool main_filter_has_reverb = first_section_energy > 4.f * tail_energy_;
-    bool main_filter_is_sane = first_section_energy > 2.f * tail_energy_ &&
-                               matching_data[peak_index] < 100.f;
-
-    // Not detecting any decay, but tail is over noise - assume max decay.
-    if (num_reverb_decay_sections_ == 0 && main_filter_is_sane &&
-        main_filter_has_reverb) {
-      decay = kMaxDecay;
-    }
-
-    if (main_filter_is_sane && num_reverb_decay_sections_ > 0) {
-      decay = std::max(.97f * reverb_decay_, decay);
-      reverb_decay_ -= alpha_ * (reverb_decay_ - decay);
-    }
-
-    found_end_of_reverb_decay_ =
-        !(main_filter_is_sane && main_filter_has_reverb);
-    alpha_ = 0.f;  // Stop estimation of the decay until another good filter is
-                   // received
-  }
-}
-
-// Updates the estimation of the frequency response at the filter tail.
-void ReverbModelEstimator::UpdateFreqRespTail(
-    const std::vector<std::array<float, kFftLengthBy2Plus1>>&
-        filter_freq_response,
-    int filter_delay_blocks,
-    float alpha) {
-  size_t num_blocks = filter_freq_response.size();
-  rtc::ArrayView<const float> freq_resp_tail(
-      filter_freq_response[num_blocks - 1]);
-  rtc::ArrayView<const float> freq_resp_direct_path(
-      filter_freq_response[filter_delay_blocks]);
-  float ratio_energies =
-      ComputeRatioEnergies(freq_resp_direct_path, freq_resp_tail);
-  ratio_tail_to_direct_path_ +=
-      alpha * (ratio_energies - ratio_tail_to_direct_path_);
-
-  for (size_t k = 0; k < kFftLengthBy2Plus1; ++k) {
-    freq_resp_tail_[k] = freq_resp_direct_path[k] * ratio_tail_to_direct_path_;
-  }
-
-  for (size_t k = 1; k < kFftLengthBy2; ++k) {
-    float avg_neighbour =
-        0.5f * (freq_resp_tail_[k - 1] + freq_resp_tail_[k + 1]);
-    freq_resp_tail_[k] = std::max(freq_resp_tail_[k], avg_neighbour);
-  }
-}
-
-void ReverbModelEstimator::Dump(
-    const std::unique_ptr<ApmDataDumper>& data_dumper) {
-  data_dumper->DumpRaw("aec3_reverb_decay", reverb_decay_);
-  data_dumper->DumpRaw("aec3_reverb_tail_energy", tail_energy_);
-  data_dumper->DumpRaw("aec3_reverb_alpha", alpha_);
-  data_dumper->DumpRaw("aec3_num_reverb_decay_sections",
-                       static_cast<int>(num_reverb_decay_sections_));
-}
-
 }  // namespace webrtc

modules/audio_processing/aec3/reverb_model_estimator.h

@@ -11,72 +11,49 @@
 #ifndef MODULES_AUDIO_PROCESSING_AEC3_REVERB_MODEL_ESTIMATOR_H_
 #define MODULES_AUDIO_PROCESSING_AEC3_REVERB_MODEL_ESTIMATOR_H_
 
-#include <memory>
 #include <vector>
 
 #include "absl/types/optional.h"
 #include "api/array_view.h"
 #include "api/audio/echo_canceller3_config.h"
-#include "modules/audio_processing/aec3/aec3_common.h"
+#include "modules/audio_processing/aec3/reverb_decay_estimator.h"
-#include "modules/audio_processing/logging/apm_data_dumper.h"
+#include "modules/audio_processing/aec3/reverb_frequency_response.h"
 
 namespace webrtc {
 
-// The ReverbModelEstimator class describes an estimator of the parameters
+class ApmDataDumper;
-// that are used for the reverberant model.
+
+// Class for estimating the model parameters for the reverberant echo.
 class ReverbModelEstimator {
  public:
   explicit ReverbModelEstimator(const EchoCanceller3Config& config);
   ~ReverbModelEstimator();
-  // Updates the model.
+
-  void Update(const std::vector<float>& impulse_response,
+  // Updates the estimates based on new data.
+  void Update(rtc::ArrayView<const float> impulse_response,
               const std::vector<std::array<float, kFftLengthBy2Plus1>>&
-                  filter_freq_response,
+                  frequency_response,
-              const absl::optional<float>& quality_linear,
+              const absl::optional<float>& linear_filter_quality,
               int filter_delay_blocks,
               bool usable_linear_estimate,
-              float default_decay,
               bool stationary_block);
-  // Returns the decay for the exponential model.
-  float ReverbDecay() const { return reverb_decay_; }
 
-  void Dump(const std::unique_ptr<ApmDataDumper>& data_dumper);
+  // Returns the exponential decay of the reverberant echo.
+  float ReverbDecay() const { return reverb_decay_estimator_.Decay(); }
 
-  // Return the estimated freq. response of the tail of the filter.
+  // Return the frequency response of the reverberant echo.
-  rtc::ArrayView<const float> GetFreqRespTail() const {
+  rtc::ArrayView<const float> GetReverbFrequencyResponse() const {
-    return freq_resp_tail_;
+    return reverb_frequency_response_.FrequencyResponse();
+  }
+
+  // Dumps debug data.
+  void Dump(ApmDataDumper* data_dumper) const {
+    reverb_decay_estimator_.Dump(data_dumper);
   }
 
  private:
-  bool IsAGoodFilterForDecayEstimation(int filter_delay_blocks,
+  ReverbDecayEstimator reverb_decay_estimator_;
-                                       bool usable_linear_estimate,
+  ReverbFrequencyResponse reverb_frequency_response_;
-                                       size_t length_filter);
-  void UpdateReverbDecay(const std::vector<float>& impulse_response);
-
-  void UpdateFreqRespTail(
-      const std::vector<std::array<float, kFftLengthBy2Plus1>>&
-          filter_freq_response,
-      int filter_delay_blocks,
-      float alpha);
-
-  void ResetDecayEstimation();
-
-  const size_t filter_main_length_blocks_;
-
-  float accumulated_nz_ = 0.f;
-  float accumulated_nn_ = 0.f;
-  float accumulated_count_ = 0.f;
-  size_t current_reverb_decay_section_ = 0;
-  size_t num_reverb_decay_sections_ = 0;
-  size_t num_reverb_decay_sections_next_ = 0;
-  bool found_end_of_reverb_decay_ = false;
-  std::array<float, kMaxAdaptiveFilterLength> block_energies_;
-  float reverb_decay_;
-  float tail_energy_ = 0.f;
-  float alpha_ = 0.f;
-  std::array<float, kFftLengthBy2Plus1> freq_resp_tail_;
-  float ratio_tail_to_direct_path_ = 0.f;
-  bool enable_smooth_freq_resp_tail_updates_;
 };
 
 }  // namespace webrtc

modules/audio_processing/aec3/reverb_model_estimator_unittest.cc

@@ -29,16 +29,16 @@ namespace webrtc {
 class ReverbModelEstimatorTest {
  public:
   explicit ReverbModelEstimatorTest(float default_decay)
-      : default_decay_(default_decay),
+      : default_decay_(default_decay), estimated_decay_(default_decay) {
-        estimated_decay_(default_decay),
-        h_(aec3_config_.filter.main.length_blocks * kBlockSize, 0.f),
-        H2_(aec3_config_.filter.main.length_blocks) {
     aec3_config_.ep_strength.default_len = default_decay_;
+    aec3_config_.filter.main.length_blocks = 40;
+    h_.resize(aec3_config_.filter.main.length_blocks * kBlockSize);
+    H2_.resize(aec3_config_.filter.main.length_blocks);
     CreateImpulseResponseWithDecay();
   }
   void RunEstimator();
   float GetDecay() { return estimated_decay_; }
-  float GetTrueDecay() { return true_power_decay_; }
+  float GetTrueDecay() { return kTruePowerDecay; }
   float GetPowerTailDb() { return 10.f * log10(estimated_power_tail_); }
   float GetTruePowerTailDb() { return 10.f * log10(true_power_tail_); }
 
@@ -46,10 +46,10 @@ class ReverbModelEstimatorTest {
   void CreateImpulseResponseWithDecay();
 
   absl::optional<float> quality_linear_ = 1.0f;
-  static constexpr int filter_delay_blocks_ = 2;
+  static constexpr int kFilterDelayBlocks = 2;
-  static constexpr bool usable_linear_estimate_ = true;
+  static constexpr bool kUsableLinearEstimate = true;
-  static constexpr bool stationary_block_ = false;
+  static constexpr bool kStationaryBlock = false;
-  static constexpr float true_power_decay_ = 0.5f;
+  static constexpr float kTruePowerDecay = 0.5f;
   EchoCanceller3Config aec3_config_;
   float default_decay_;
   float estimated_decay_;
@@ -63,43 +63,35 @@ void ReverbModelEstimatorTest::CreateImpulseResponseWithDecay() {
   const Aec3Fft fft;
   RTC_DCHECK_EQ(h_.size(), aec3_config_.filter.main.length_blocks * kBlockSize);
   RTC_DCHECK_EQ(H2_.size(), aec3_config_.filter.main.length_blocks);
-  RTC_DCHECK_EQ(filter_delay_blocks_, 2);
+  RTC_DCHECK_EQ(kFilterDelayBlocks, 2);
-  const float peak = 1.0f;
+
-  float decay_power_sample = std::sqrt(true_power_decay_);
+  float decay_sample = std::sqrt(powf(kTruePowerDecay, 1.f / kBlockSize));
-  for (size_t k = 1; k < kBlockSizeLog2; k++) {
+  const size_t filter_delay_coefficients = kFilterDelayBlocks * kBlockSize;
-    decay_power_sample = std::sqrt(decay_power_sample);
+  std::fill(h_.begin(), h_.end(), 0.f);
-  }
+  h_[filter_delay_coefficients] = 1.f;
-  h_[filter_delay_blocks_ * kBlockSize] = peak;
+  for (size_t k = filter_delay_coefficients + 1; k < h_.size(); ++k) {
-  for (size_t k = filter_delay_blocks_ * kBlockSize + 1; k < h_.size(); ++k) {
+    h_[k] = h_[k - 1] * decay_sample;
-    h_[k] = h_[k - 1] * std::sqrt(decay_power_sample);
   }
 
-  for (size_t block = 0; block < H2_.size(); ++block) {
+  std::array<float, kFftLength> fft_data;
-    std::array<float, kFftLength> h_block;
+  FftData H_j;
-    h_block.fill(0.f);
+  for (size_t j = 0, k = 0; j < H2_.size(); ++j, k += kBlockSize) {
-    FftData H_block;
+    fft_data.fill(0.f);
-    rtc::ArrayView<float> H2_block(H2_[block]);
+    std::copy(h_.begin() + k, h_.begin() + k + kBlockSize, fft_data.begin());
-    std::copy(h_.begin() + block * kBlockSize,
+    fft.Fft(&fft_data, &H_j);
-              h_.begin() + block * (kBlockSize + 1), h_block.begin());
+    H_j.Spectrum(Aec3Optimization::kNone, H2_[j]);
-
-    fft.Fft(&h_block, &H_block);
-    for (size_t k = 0; k < H2_block.size(); ++k) {
-      H2_block[k] =
-          H_block.re[k] * H_block.re[k] + H_block.im[k] * H_block.im[k];
-    }
   }
   rtc::ArrayView<float> H2_tail(H2_[H2_.size() - 1]);
   true_power_tail_ = std::accumulate(H2_tail.begin(), H2_tail.end(), 0.f);
 }
 void ReverbModelEstimatorTest::RunEstimator() {
   ReverbModelEstimator estimator(aec3_config_);
-  for (size_t k = 0; k < 1000; ++k) {
+  for (size_t k = 0; k < 3000; ++k) {
-    estimator.Update(h_, H2_, quality_linear_, filter_delay_blocks_,
+    estimator.Update(h_, H2_, quality_linear_, kFilterDelayBlocks,
-                     usable_linear_estimate_, default_decay_,
+                     kUsableLinearEstimate, kStationaryBlock);
-                     stationary_block_);
   }
   estimated_decay_ = estimator.ReverbDecay();
-  rtc::ArrayView<const float> freq_resp_tail = estimator.GetFreqRespTail();
+  auto freq_resp_tail = estimator.GetReverbFrequencyResponse();
   estimated_power_tail_ =
       std::accumulate(freq_resp_tail.begin(), freq_resp_tail.end(), 0.f);
 }

test/fuzzers/audio_processing_configs_fuzzer.cc

@@ -50,7 +50,8 @@ const std::string kFieldTrialNames[] = {
     "WebRTC-Aec3ShadowFilterJumpstartKillSwitch",
     "WebRTC-Aec3EarlyLinearFilterUsageKillSwitch",
     "WebRTC-Aec3ShortInitialStateKillSwitch",
-    "WebRTC-Aec3StandardNonlinearReverbModelKillSwitch"};
+    "WebRTC-Aec3StandardNonlinearReverbModelKillSwitch",
+    "WebRTC-Aec3EnableAdaptiveEchoReverbEstimation"};
 
 std::unique_ptr<AudioProcessing> CreateApm(test::FuzzDataHelper* fuzz_data,
                                            std::string* field_trial_string) {