/* * 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 #include #include #include #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 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& 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 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 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)) { RTC_DCHECK_GT(config.filter.refined.length_blocks, static_cast(kEarlyReverbMinSizeBlocks)); } ReverbDecayEstimator::~ReverbDecayEstimator() = default; void ReverbDecayEstimator::Update(rtc::ArrayView filter, const absl::optional& filter_quality, int filter_delay_blocks, bool usable_linear_filter, bool stationary_signal) { const int filter_size = static_cast(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 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 filter) { auto h = rtc::ArrayView( filter.begin() + block_to_analyze_ * kFftLengthBy2, kFftLengthBy2); // Compute squared filter coeffiecients for the block to analyze_; std::array 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_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(numerators_.size() - 1)); float x_value = static_cast(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