diff options
Diffstat (limited to 'webrtc/modules/audio_processing/aec3/reverb_decay_estimator.cc')
| -rw-r--r-- | webrtc/modules/audio_processing/aec3/reverb_decay_estimator.cc | 409 |
1 files changed, 409 insertions, 0 deletions
diff --git a/webrtc/modules/audio_processing/aec3/reverb_decay_estimator.cc b/webrtc/modules/audio_processing/aec3/reverb_decay_estimator.cc new file mode 100644 index 0000000..f160b83 --- /dev/null +++ b/webrtc/modules/audio_processing/aec3/reverb_decay_estimator.cc @@ -0,0 +1,409 @@ +/* + * 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)) { + 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 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; + } + + 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_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 |