1 /*
2  *  Copyright (c) 2017 The WebRTC project authors. All Rights Reserved.
3  *
4  *  Use of this source code is governed by a BSD-style license
5  *  that can be found in the LICENSE file in the root of the source
6  *  tree. An additional intellectual property rights grant can be found
7  *  in the file PATENTS.  All contributing project authors may
8  *  be found in the AUTHORS file in the root of the source tree.
9  */
10 
11 #include "modules/audio_processing/aec3/erl_estimator.h"
12 
13 #include <algorithm>
14 #include <numeric>
15 
16 #include "rtc_base/checks.h"
17 
18 namespace webrtc {
19 
20 namespace {
21 
22 constexpr float kMinErl = 0.01f;
23 constexpr float kMaxErl = 1000.f;
24 
25 }  // namespace
26 
ErlEstimator(size_t startup_phase_length_blocks_)27 ErlEstimator::ErlEstimator(size_t startup_phase_length_blocks_)
28     : startup_phase_length_blocks__(startup_phase_length_blocks_) {
29   erl_.fill(kMaxErl);
30   hold_counters_.fill(0);
31   erl_time_domain_ = kMaxErl;
32   hold_counter_time_domain_ = 0;
33 }
34 
35 ErlEstimator::~ErlEstimator() = default;
36 
Reset()37 void ErlEstimator::Reset() {
38   blocks_since_reset_ = 0;
39 }
40 
Update(const std::vector<bool> & converged_filters,rtc::ArrayView<const std::array<float,kFftLengthBy2Plus1>> render_spectra,rtc::ArrayView<const std::array<float,kFftLengthBy2Plus1>> capture_spectra)41 void ErlEstimator::Update(
42     const std::vector<bool>& converged_filters,
43     rtc::ArrayView<const std::array<float, kFftLengthBy2Plus1>> render_spectra,
44     rtc::ArrayView<const std::array<float, kFftLengthBy2Plus1>>
45         capture_spectra) {
46   const size_t num_capture_channels = converged_filters.size();
47   RTC_DCHECK_EQ(capture_spectra.size(), num_capture_channels);
48 
49   // Corresponds to WGN of power -46 dBFS.
50   constexpr float kX2Min = 44015068.0f;
51 
52   const auto first_converged_iter =
53       std::find(converged_filters.begin(), converged_filters.end(), true);
54   const bool any_filter_converged =
55       first_converged_iter != converged_filters.end();
56 
57   if (++blocks_since_reset_ < startup_phase_length_blocks__ ||
58       !any_filter_converged) {
59     return;
60   }
61 
62   // Use the maximum spectrum across capture and the maximum across render.
63   std::array<float, kFftLengthBy2Plus1> max_capture_spectrum_data;
64   std::array<float, kFftLengthBy2Plus1> max_capture_spectrum =
65       capture_spectra[/*channel=*/0];
66   if (num_capture_channels > 1) {
67     // Initialize using the first channel with a converged filter.
68     const size_t first_converged =
69         std::distance(converged_filters.begin(), first_converged_iter);
70     RTC_DCHECK_GE(first_converged, 0);
71     RTC_DCHECK_LT(first_converged, num_capture_channels);
72     max_capture_spectrum_data = capture_spectra[first_converged];
73 
74     for (size_t ch = first_converged + 1; ch < num_capture_channels; ++ch) {
75       if (!converged_filters[ch]) {
76         continue;
77       }
78       for (size_t k = 0; k < kFftLengthBy2Plus1; ++k) {
79         max_capture_spectrum_data[k] =
80             std::max(max_capture_spectrum_data[k], capture_spectra[ch][k]);
81       }
82     }
83     max_capture_spectrum = max_capture_spectrum_data;
84   }
85 
86   const size_t num_render_channels = render_spectra.size();
87   std::array<float, kFftLengthBy2Plus1> max_render_spectrum_data;
88   rtc::ArrayView<const float, kFftLengthBy2Plus1> max_render_spectrum =
89       render_spectra[/*channel=*/0];
90   if (num_render_channels > 1) {
91     std::copy(render_spectra[0].begin(), render_spectra[0].end(),
92               max_render_spectrum_data.begin());
93     for (size_t ch = 1; ch < num_render_channels; ++ch) {
94       for (size_t k = 0; k < kFftLengthBy2Plus1; ++k) {
95         max_render_spectrum_data[k] =
96             std::max(max_render_spectrum_data[k], render_spectra[ch][k]);
97       }
98     }
99     max_render_spectrum = max_render_spectrum_data;
100   }
101 
102   const auto& X2 = max_render_spectrum;
103   const auto& Y2 = max_capture_spectrum;
104 
105   // Update the estimates in a maximum statistics manner.
106   for (size_t k = 1; k < kFftLengthBy2; ++k) {
107     if (X2[k] > kX2Min) {
108       const float new_erl = Y2[k] / X2[k];
109       if (new_erl < erl_[k]) {
110         hold_counters_[k - 1] = 1000;
111         erl_[k] += 0.1f * (new_erl - erl_[k]);
112         erl_[k] = std::max(erl_[k], kMinErl);
113       }
114     }
115   }
116 
117   std::for_each(hold_counters_.begin(), hold_counters_.end(),
118                 [](int& a) { --a; });
119   std::transform(hold_counters_.begin(), hold_counters_.end(), erl_.begin() + 1,
120                  erl_.begin() + 1, [](int a, float b) {
121                    return a > 0 ? b : std::min(kMaxErl, 2.f * b);
122                  });
123 
124   erl_[0] = erl_[1];
125   erl_[kFftLengthBy2] = erl_[kFftLengthBy2 - 1];
126 
127   // Compute ERL over all frequency bins.
128   const float X2_sum = std::accumulate(X2.begin(), X2.end(), 0.0f);
129 
130   if (X2_sum > kX2Min * X2.size()) {
131     const float Y2_sum = std::accumulate(Y2.begin(), Y2.end(), 0.0f);
132     const float new_erl = Y2_sum / X2_sum;
133     if (new_erl < erl_time_domain_) {
134       hold_counter_time_domain_ = 1000;
135       erl_time_domain_ += 0.1f * (new_erl - erl_time_domain_);
136       erl_time_domain_ = std::max(erl_time_domain_, kMinErl);
137     }
138   }
139 
140   --hold_counter_time_domain_;
141   erl_time_domain_ = (hold_counter_time_domain_ > 0)
142                          ? erl_time_domain_
143                          : std::min(kMaxErl, 2.f * erl_time_domain_);
144 }
145 
146 }  // namespace webrtc
147