modules/audio_processing/aec3/reverb_decay_estimator.cc - src - Git at Google

 /*
  *  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 <stddef.h>

 #include <algorithm>
 #include <cmath>
 #include <numeric>

 #include "api/array_view.h"
 #include "api/audio/echo_canceller3_config.h"
 #include "modules/audio_processing/logging/apm_data_dumper.h"
 #include "rtc_base/checks.h"

 namespace webrtc {

 namespace {

 constexpr int kEarlyReverbMinSizeBlocks = 3;
 constexpr int kBlocksPerSection = 6;
 // Linear regression approach assumes symmetric index around 0.
 constexpr float kEarlyReverbFirstPointAtLinearRegressors =
     -0.5f * kBlocksPerSection * kFftLengthBy2 + 0.5f;

 // 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);
 }

 // Returns the peak energy of an impulse response.
 float BlockEnergyPeak(rtc::ArrayView<const float> h, int peak_block) {
   RTC_DCHECK_LE((peak_block + 1) * kFftLengthBy2, h.size());
   RTC_DCHECK_GE(peak_block, 0);
   float peak_value =
       *std::max_element(h.begin() + peak_block * kFftLengthBy2,
                         h.begin() + (peak_block + 1) * kFftLengthBy2,
                         [](float a, float b) { return a * a < b * b; });
   return peak_value * peak_value;
 }

 // Returns the average energy of an impulse response block.
 float BlockEnergyAverage(rtc::ArrayView<const float> h, int block_index) {
   RTC_DCHECK_LE((block_index + 1) * kFftLengthBy2, h.size());
   RTC_DCHECK_GE(block_index, 0);
   constexpr float kOneByFftLengthBy2 = 1.f / kFftLengthBy2;
   const auto sum_of_squares = [](float a, float b) { return a + b * b; };
   return std::accumulate(h.begin() + block_index * kFftLengthBy2,
                          h.begin() + (block_index + 1) * kFftLengthBy2, 0.f,
                          sum_of_squares) *
          kOneByFftLengthBy2;
 }

 }  // namespace

 ReverbDecayEstimator::ReverbDecayEstimator(const EchoCanceller3Config& config)
     : filter_length_blocks_(config.filter.refined.length_blocks),
       filter_length_coefficients_(GetTimeDomainLength(filter_length_blocks_)),
       use_adaptive_echo_decay_(config.ep_strength.default_len < 0.f),
       early_reverb_estimator_(config.filter.refined.length_blocks -
                               kEarlyReverbMinSizeBlocks),
       late_reverb_start_(kEarlyReverbMinSizeBlocks),
       late_reverb_end_(kEarlyReverbMinSizeBlocks),
       previous_gains_(config.filter.refined.length_blocks, 0.f),
       decay_(std::fabs(config.ep_strength.default_len)),
       mild_decay_(std::fabs(config.ep_strength.nearend_len)) {
   RTC_DCHECK_GT(config.filter.refined.length_blocks,
                 static_cast<size_t>(kEarlyReverbMinSizeBlocks));
 }

 ReverbDecayEstimator::~ReverbDecayEstimator() = default;

 void ReverbDecayEstimator::Update(rtc::ArrayView<const float> filter,
                                   const std::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;
   }

   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, filter_delay_blocks);
   }
 }

 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,
                                          int peak_block) {
   auto& h = filter;
   RTC_DCHECK_EQ(0, h.size() % kFftLengthBy2);

   // 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 = BlockEnergyAverage(h, block_to_analyze_);
   const size_t h_size_blocks = h.size() >> kFftLengthBy2Log2;
   tail_gain_ = BlockEnergyAverage(h, h_size_blocks - 1);
   float peak_energy = BlockEnergyPeak(h, peak_block);
   const bool sufficient_reverb_decay = first_reverb_gain > 4.f * tail_gain_;
   const bool valid_filter =
       first_reverb_gain > 2.f * tail_gain_ && peak_energy < 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_late_reverb_start", late_reverb_start_);
   data_dumper->DumpRaw("aec3_late_reverb_end", late_reverb_end_);
   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 ? -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_DCHECK_NOTREACHED();
     return 0.f;
   }
   return nz_ / nn_;
 }

 ReverbDecayEstimator::EarlyReverbLengthEstimator::EarlyReverbLengthEstimator(
     int max_blocks)
     : numerators_smooth_(max_blocks - kBlocksPerSection, 0.f),
       numerators_(numerators_smooth_.size(), 0.f),
       coefficients_counter_(0) {
   RTC_DCHECK_LE(0, max_blocks);
 }

 ReverbDecayEstimator::EarlyReverbLengthEstimator::
     ~EarlyReverbLengthEstimator() = default;

 void ReverbDecayEstimator::EarlyReverbLengthEstimator::Reset() {
   coefficients_counter_ = 0;
   std::fill(numerators_.begin(), numerators_.end(), 0.f);
   block_counter_ = 0;
 }

 void ReverbDecayEstimator::EarlyReverbLengthEstimator::Accumulate(
     float value,
     float smoothing) {
   // Each section is composed by kBlocksPerSection blocks and each section
   // overlaps with the next one in (kBlocksPerSection - 1) blocks. For example,
   // the first section covers the blocks [0:5], the second covers the blocks
   // [1:6] and so on. As a result, for each value, kBlocksPerSection sections
   // need to be updated.
   int first_section_index = std::max(block_counter_ - kBlocksPerSection + 1, 0);
   int last_section_index =
       std::min(block_counter_, static_cast<int>(numerators_.size() - 1));
   float x_value = static_cast<float>(coefficients_counter_) +
                   kEarlyReverbFirstPointAtLinearRegressors;
   const float value_to_inc = kFftLengthBy2 * value;
   float value_to_add =
       x_value * value + (block_counter_ - last_section_index) * value_to_inc;
   for (int section = last_section_index; section >= first_section_index;
        --section, value_to_add += value_to_inc) {
     numerators_[section] += value_to_add;
   }

   // Check if this update was the last coefficient of the current block. In that
   // case, check if we are at the end of one of the sections and update the
   // numerator of the linear regressor that is computed in such section.
   if (++coefficients_counter_ == kFftLengthBy2) {
     if (block_counter_ >= (kBlocksPerSection - 1)) {
       size_t section = block_counter_ - (kBlocksPerSection - 1);
       RTC_DCHECK_GT(numerators_.size(), section);
       RTC_DCHECK_GT(numerators_smooth_.size(), section);
       numerators_smooth_[section] +=
           smoothing * (numerators_[section] - numerators_smooth_[section]);
       n_sections_ = section + 1;
     }
     ++block_counter_;
     coefficients_counter_ = 0;
   }
 }

 // Estimates the size in blocks of the early reverb. The estimation is done by
 // comparing the tilt that is estimated in each section. As an optimization
 // detail and due to the fact that all the linear regressors that are computed
 // shared the same denominator, the comparison of the tilts is done by a
 // comparison of the numerator of the linear regressors.
 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 = 9;

   if (n_sections_ < kNumSectionsToAnalyze) {
     return 0;
   }

   // Estimation of the blocks that correspond to early reverberations. The
   // estimation is done by analyzing the impulse response. The portions of the
   // impulse response whose energy is not decreasing over its coefficients are
   // considered to be part of the early reverberations. Furthermore, the blocks
   // where the energy is decreasing faster than what it does at the end of the
   // impulse response are also considered to be part of the early
   // reverberations. The estimation is limited to the first
   // kNumSectionsToAnalyze sections.

   RTC_DCHECK_LE(n_sections_, numerators_smooth_.size());
   const float min_numerator_tail =
       *std::min_element(numerators_smooth_.begin() + kNumSectionsToAnalyze,
                         numerators_smooth_.begin() + n_sections_);
   int early_reverb_size_minus_1 = 0;
   for (int k = 0; k < kNumSectionsToAnalyze; ++k) {
     if ((numerators_smooth_[k] > numerator_11) ||
         (numerators_smooth_[k] < numerator_08 &&
          numerators_smooth_[k] < 0.9f * min_numerator_tail)) {
       early_reverb_size_minus_1 = k;
     }
   }

   return early_reverb_size_minus_1 == 0 ? 0 : early_reverb_size_minus_1 + 1;
 }

 void ReverbDecayEstimator::EarlyReverbLengthEstimator::Dump(
     ApmDataDumper* data_dumper) const {
   data_dumper->DumpRaw("aec3_er_acum_numerator", numerators_smooth_);
 }

 }  // namespace webrtc
	/*
	* 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 <stddef.h>

	#include <algorithm>
	#include <cmath>
	#include <numeric>

	#include "api/array_view.h"
	#include "api/audio/echo_canceller3_config.h"
	#include "modules/audio_processing/logging/apm_data_dumper.h"
	#include "rtc_base/checks.h"

	namespace webrtc {

	namespace {

	constexpr int kEarlyReverbMinSizeBlocks = 3;
	constexpr int kBlocksPerSection = 6;
	// Linear regression approach assumes symmetric index around 0.
	constexpr float kEarlyReverbFirstPointAtLinearRegressors =
	-0.5f * kBlocksPerSection * kFftLengthBy2 + 0.5f;

	// 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);
	}

	// Returns the peak energy of an impulse response.
	float BlockEnergyPeak(rtc::ArrayView<const float> h, int peak_block) {
	RTC_DCHECK_LE((peak_block + 1) * kFftLengthBy2, h.size());
	RTC_DCHECK_GE(peak_block, 0);
	float peak_value =
	std::max_element(h.begin() + peak_block kFftLengthBy2,
	h.begin() + (peak_block + 1) * kFftLengthBy2,
	[](float a, float b) { return a * a < b * b; });
	return peak_value * peak_value;
	}

	// Returns the average energy of an impulse response block.
	float BlockEnergyAverage(rtc::ArrayView<const float> h, int block_index) {
	RTC_DCHECK_LE((block_index + 1) * kFftLengthBy2, h.size());
	RTC_DCHECK_GE(block_index, 0);
	constexpr float kOneByFftLengthBy2 = 1.f / kFftLengthBy2;
	const auto sum_of_squares = [](float a, float b) { return a + b * b; };
	return std::accumulate(h.begin() + block_index * kFftLengthBy2,
	h.begin() + (block_index + 1) * kFftLengthBy2, 0.f,
	sum_of_squares) *
	kOneByFftLengthBy2;
	}

	} // namespace

	ReverbDecayEstimator::ReverbDecayEstimator(const EchoCanceller3Config& config)
	: filter_length_blocks_(config.filter.refined.length_blocks),
	filter_length_coefficients_(GetTimeDomainLength(filter_length_blocks_)),
	use_adaptive_echo_decay_(config.ep_strength.default_len < 0.f),
	early_reverb_estimator_(config.filter.refined.length_blocks -
	kEarlyReverbMinSizeBlocks),
	late_reverb_start_(kEarlyReverbMinSizeBlocks),
	late_reverb_end_(kEarlyReverbMinSizeBlocks),
	previous_gains_(config.filter.refined.length_blocks, 0.f),
	decay_(std::fabs(config.ep_strength.default_len)),
	mild_decay_(std::fabs(config.ep_strength.nearend_len)) {
	RTC_DCHECK_GT(config.filter.refined.length_blocks,
	static_cast<size_t>(kEarlyReverbMinSizeBlocks));
	}

	ReverbDecayEstimator::~ReverbDecayEstimator() = default;

	void ReverbDecayEstimator::Update(rtc::ArrayView<const float> filter,
	const std::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;
	}

	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, filter_delay_blocks);
	}
	}

	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,
	int peak_block) {
	auto& h = filter;
	RTC_DCHECK_EQ(0, h.size() % kFftLengthBy2);

	// 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 = BlockEnergyAverage(h, block_to_analyze_);
	const size_t h_size_blocks = h.size() >> kFftLengthBy2Log2;
	tail_gain_ = BlockEnergyAverage(h, h_size_blocks - 1);
	float peak_energy = BlockEnergyPeak(h, peak_block);
	const bool sufficient_reverb_decay = first_reverb_gain > 4.f * tail_gain_;
	const bool valid_filter =
	first_reverb_gain > 2.f * tail_gain_ && peak_energy < 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_late_reverb_start", late_reverb_start_);
	data_dumper->DumpRaw("aec3_late_reverb_end", late_reverb_end_);
	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 ? -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_DCHECK_NOTREACHED();
	return 0.f;
	}
	return nz_ / nn_;
	}

	ReverbDecayEstimator::EarlyReverbLengthEstimator::EarlyReverbLengthEstimator(
	int max_blocks)
	: numerators_smooth_(max_blocks - kBlocksPerSection, 0.f),
	numerators_(numerators_smooth_.size(), 0.f),
	coefficients_counter_(0) {
	RTC_DCHECK_LE(0, max_blocks);
	}

	ReverbDecayEstimator::EarlyReverbLengthEstimator::
	~EarlyReverbLengthEstimator() = default;

	void ReverbDecayEstimator::EarlyReverbLengthEstimator::Reset() {
	coefficients_counter_ = 0;
	std::fill(numerators_.begin(), numerators_.end(), 0.f);
	block_counter_ = 0;
	}

	void ReverbDecayEstimator::EarlyReverbLengthEstimator::Accumulate(
	float value,
	float smoothing) {
	// Each section is composed by kBlocksPerSection blocks and each section
	// overlaps with the next one in (kBlocksPerSection - 1) blocks. For example,
	// the first section covers the blocks [0:5], the second covers the blocks
	// [1:6] and so on. As a result, for each value, kBlocksPerSection sections
	// need to be updated.
	int first_section_index = std::max(block_counter_ - kBlocksPerSection + 1, 0);
	int last_section_index =
	std::min(block_counter_, static_cast<int>(numerators_.size() - 1));
	float x_value = static_cast<float>(coefficients_counter_) +
	kEarlyReverbFirstPointAtLinearRegressors;
	const float value_to_inc = kFftLengthBy2 * value;
	float value_to_add =
	x_value * value + (block_counter_ - last_section_index) * value_to_inc;
	for (int section = last_section_index; section >= first_section_index;
	--section, value_to_add += value_to_inc) {
	numerators_[section] += value_to_add;
	}

	// Check if this update was the last coefficient of the current block. In that
	// case, check if we are at the end of one of the sections and update the
	// numerator of the linear regressor that is computed in such section.
	if (++coefficients_counter_ == kFftLengthBy2) {
	if (block_counter_ >= (kBlocksPerSection - 1)) {
	size_t section = block_counter_ - (kBlocksPerSection - 1);
	RTC_DCHECK_GT(numerators_.size(), section);
	RTC_DCHECK_GT(numerators_smooth_.size(), section);
	numerators_smooth_[section] +=
	smoothing * (numerators_[section] - numerators_smooth_[section]);
	n_sections_ = section + 1;
	}
	++block_counter_;
	coefficients_counter_ = 0;
	}
	}

	// Estimates the size in blocks of the early reverb. The estimation is done by
	// comparing the tilt that is estimated in each section. As an optimization
	// detail and due to the fact that all the linear regressors that are computed
	// shared the same denominator, the comparison of the tilts is done by a
	// comparison of the numerator of the linear regressors.
	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 = 9;

	if (n_sections_ < kNumSectionsToAnalyze) {
	return 0;
	}

	// Estimation of the blocks that correspond to early reverberations. The
	// estimation is done by analyzing the impulse response. The portions of the
	// impulse response whose energy is not decreasing over its coefficients are
	// considered to be part of the early reverberations. Furthermore, the blocks
	// where the energy is decreasing faster than what it does at the end of the
	// impulse response are also considered to be part of the early
	// reverberations. The estimation is limited to the first
	// kNumSectionsToAnalyze sections.

	RTC_DCHECK_LE(n_sections_, numerators_smooth_.size());
	const float min_numerator_tail =
	*std::min_element(numerators_smooth_.begin() + kNumSectionsToAnalyze,
	numerators_smooth_.begin() + n_sections_);
	int early_reverb_size_minus_1 = 0;
	for (int k = 0; k < kNumSectionsToAnalyze; ++k) {
	if ((numerators_smooth_[k] > numerator_11) \|\|
	(numerators_smooth_[k] < numerator_08 &&
	numerators_smooth_[k] < 0.9f * min_numerator_tail)) {
	early_reverb_size_minus_1 = k;
	}
	}

	return early_reverb_size_minus_1 == 0 ? 0 : early_reverb_size_minus_1 + 1;
	}

	void ReverbDecayEstimator::EarlyReverbLengthEstimator::Dump(
	ApmDataDumper* data_dumper) const {
	data_dumper->DumpRaw("aec3_er_acum_numerator", numerators_smooth_);
	}

	} // namespace webrtc