1 /*
2  *  Copyright (c) 2018 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/stationarity_estimator.h"
12 
13 #include <algorithm>
14 #include <array>
15 
16 #include "api/array_view.h"
17 #include "modules/audio_processing/aec3/aec3_common.h"
18 #include "modules/audio_processing/aec3/spectrum_buffer.h"
19 #include "modules/audio_processing/logging/apm_data_dumper.h"
20 #include "rtc_base/atomic_ops.h"
21 
22 namespace webrtc {
23 
24 namespace {
25 constexpr float kMinNoisePower = 10.f;
26 constexpr int kHangoverBlocks = kNumBlocksPerSecond / 20;
27 constexpr int kNBlocksAverageInitPhase = 20;
28 constexpr int kNBlocksInitialPhase = kNumBlocksPerSecond * 2.;
29 }  // namespace
30 
StationarityEstimator()31 StationarityEstimator::StationarityEstimator()
32     : data_dumper_(
33           new ApmDataDumper(rtc::AtomicOps::Increment(&instance_count_))) {
34   Reset();
35 }
36 
37 StationarityEstimator::~StationarityEstimator() = default;
38 
Reset()39 void StationarityEstimator::Reset() {
40   noise_.Reset();
41   hangovers_.fill(0);
42   stationarity_flags_.fill(false);
43 }
44 
45 // Update just the noise estimator. Usefull until the delay is known
UpdateNoiseEstimator(rtc::ArrayView<const std::array<float,kFftLengthBy2Plus1>> spectrum)46 void StationarityEstimator::UpdateNoiseEstimator(
47     rtc::ArrayView<const std::array<float, kFftLengthBy2Plus1>> spectrum) {
48   noise_.Update(spectrum);
49   data_dumper_->DumpRaw("aec3_stationarity_noise_spectrum", noise_.Spectrum());
50   data_dumper_->DumpRaw("aec3_stationarity_is_block_stationary",
51                         IsBlockStationary());
52 }
53 
UpdateStationarityFlags(const SpectrumBuffer & spectrum_buffer,rtc::ArrayView<const float> render_reverb_contribution_spectrum,int idx_current,int num_lookahead)54 void StationarityEstimator::UpdateStationarityFlags(
55     const SpectrumBuffer& spectrum_buffer,
56     rtc::ArrayView<const float> render_reverb_contribution_spectrum,
57     int idx_current,
58     int num_lookahead) {
59   std::array<int, kWindowLength> indexes;
60   int num_lookahead_bounded = std::min(num_lookahead, kWindowLength - 1);
61   int idx = idx_current;
62 
63   if (num_lookahead_bounded < kWindowLength - 1) {
64     int num_lookback = (kWindowLength - 1) - num_lookahead_bounded;
65     idx = spectrum_buffer.OffsetIndex(idx_current, num_lookback);
66   }
67   // For estimating the stationarity properties of the current frame, the
68   // power for each band is accumulated for several consecutive spectra in the
69   // method EstimateBandStationarity.
70   // In order to avoid getting the indexes of the spectra for every band with
71   // its associated overhead, those indexes are stored in an array and then use
72   // when the estimation is done.
73   indexes[0] = idx;
74   for (size_t k = 1; k < indexes.size(); ++k) {
75     indexes[k] = spectrum_buffer.DecIndex(indexes[k - 1]);
76   }
77   RTC_DCHECK_EQ(
78       spectrum_buffer.DecIndex(indexes[kWindowLength - 1]),
79       spectrum_buffer.OffsetIndex(idx_current, -(num_lookahead_bounded + 1)));
80 
81   for (size_t k = 0; k < stationarity_flags_.size(); ++k) {
82     stationarity_flags_[k] = EstimateBandStationarity(
83         spectrum_buffer, render_reverb_contribution_spectrum, indexes, k);
84   }
85   UpdateHangover();
86   SmoothStationaryPerFreq();
87 }
88 
IsBlockStationary() const89 bool StationarityEstimator::IsBlockStationary() const {
90   float acum_stationarity = 0.f;
91   RTC_DCHECK_EQ(stationarity_flags_.size(), kFftLengthBy2Plus1);
92   for (size_t band = 0; band < stationarity_flags_.size(); ++band) {
93     bool st = IsBandStationary(band);
94     acum_stationarity += static_cast<float>(st);
95   }
96   return ((acum_stationarity * (1.f / kFftLengthBy2Plus1)) > 0.75f);
97 }
98 
EstimateBandStationarity(const SpectrumBuffer & spectrum_buffer,rtc::ArrayView<const float> average_reverb,const std::array<int,kWindowLength> & indexes,size_t band) const99 bool StationarityEstimator::EstimateBandStationarity(
100     const SpectrumBuffer& spectrum_buffer,
101     rtc::ArrayView<const float> average_reverb,
102     const std::array<int, kWindowLength>& indexes,
103     size_t band) const {
104   constexpr float kThrStationarity = 10.f;
105   float acum_power = 0.f;
106   const int num_render_channels =
107       static_cast<int>(spectrum_buffer.buffer[0].size());
108   const float one_by_num_channels = 1.f / num_render_channels;
109   for (auto idx : indexes) {
110     for (int ch = 0; ch < num_render_channels; ++ch) {
111       acum_power += spectrum_buffer.buffer[idx][ch][band] * one_by_num_channels;
112     }
113   }
114   acum_power += average_reverb[band];
115   float noise = kWindowLength * GetStationarityPowerBand(band);
116   RTC_CHECK_LT(0.f, noise);
117   bool stationary = acum_power < kThrStationarity * noise;
118   data_dumper_->DumpRaw("aec3_stationarity_long_ratio", acum_power / noise);
119   return stationary;
120 }
121 
AreAllBandsStationary()122 bool StationarityEstimator::AreAllBandsStationary() {
123   for (auto b : stationarity_flags_) {
124     if (!b)
125       return false;
126   }
127   return true;
128 }
129 
UpdateHangover()130 void StationarityEstimator::UpdateHangover() {
131   bool reduce_hangover = AreAllBandsStationary();
132   for (size_t k = 0; k < stationarity_flags_.size(); ++k) {
133     if (!stationarity_flags_[k]) {
134       hangovers_[k] = kHangoverBlocks;
135     } else if (reduce_hangover) {
136       hangovers_[k] = std::max(hangovers_[k] - 1, 0);
137     }
138   }
139 }
140 
SmoothStationaryPerFreq()141 void StationarityEstimator::SmoothStationaryPerFreq() {
142   std::array<bool, kFftLengthBy2Plus1> all_ahead_stationary_smooth;
143   for (size_t k = 1; k < kFftLengthBy2Plus1 - 1; ++k) {
144     all_ahead_stationary_smooth[k] = stationarity_flags_[k - 1] &&
145                                      stationarity_flags_[k] &&
146                                      stationarity_flags_[k + 1];
147   }
148 
149   all_ahead_stationary_smooth[0] = all_ahead_stationary_smooth[1];
150   all_ahead_stationary_smooth[kFftLengthBy2Plus1 - 1] =
151       all_ahead_stationary_smooth[kFftLengthBy2Plus1 - 2];
152 
153   stationarity_flags_ = all_ahead_stationary_smooth;
154 }
155 
156 int StationarityEstimator::instance_count_ = 0;
157 
NoiseSpectrum()158 StationarityEstimator::NoiseSpectrum::NoiseSpectrum() {
159   Reset();
160 }
161 
162 StationarityEstimator::NoiseSpectrum::~NoiseSpectrum() = default;
163 
Reset()164 void StationarityEstimator::NoiseSpectrum::Reset() {
165   block_counter_ = 0;
166   noise_spectrum_.fill(kMinNoisePower);
167 }
168 
Update(rtc::ArrayView<const std::array<float,kFftLengthBy2Plus1>> spectrum)169 void StationarityEstimator::NoiseSpectrum::Update(
170     rtc::ArrayView<const std::array<float, kFftLengthBy2Plus1>> spectrum) {
171   RTC_DCHECK_LE(1, spectrum[0].size());
172   const int num_render_channels = static_cast<int>(spectrum.size());
173 
174   std::array<float, kFftLengthBy2Plus1> avg_spectrum_data;
175   rtc::ArrayView<const float> avg_spectrum;
176   if (num_render_channels == 1) {
177     avg_spectrum = spectrum[0];
178   } else {
179     // For multiple channels, average the channel spectra before passing to the
180     // noise spectrum estimator.
181     avg_spectrum = avg_spectrum_data;
182     std::copy(spectrum[0].begin(), spectrum[0].end(),
183               avg_spectrum_data.begin());
184     for (int ch = 1; ch < num_render_channels; ++ch) {
185       for (size_t k = 1; k < kFftLengthBy2Plus1; ++k) {
186         avg_spectrum_data[k] += spectrum[ch][k];
187       }
188     }
189 
190     const float one_by_num_channels = 1.f / num_render_channels;
191     for (size_t k = 1; k < kFftLengthBy2Plus1; ++k) {
192       avg_spectrum_data[k] *= one_by_num_channels;
193     }
194   }
195 
196   ++block_counter_;
197   float alpha = GetAlpha();
198   for (size_t k = 0; k < kFftLengthBy2Plus1; ++k) {
199     if (block_counter_ <= kNBlocksAverageInitPhase) {
200       noise_spectrum_[k] += (1.f / kNBlocksAverageInitPhase) * avg_spectrum[k];
201     } else {
202       noise_spectrum_[k] =
203           UpdateBandBySmoothing(avg_spectrum[k], noise_spectrum_[k], alpha);
204     }
205   }
206 }
207 
GetAlpha() const208 float StationarityEstimator::NoiseSpectrum::GetAlpha() const {
209   constexpr float kAlpha = 0.004f;
210   constexpr float kAlphaInit = 0.04f;
211   constexpr float kTiltAlpha = (kAlphaInit - kAlpha) / kNBlocksInitialPhase;
212 
213   if (block_counter_ > (kNBlocksInitialPhase + kNBlocksAverageInitPhase)) {
214     return kAlpha;
215   } else {
216     return kAlphaInit -
217            kTiltAlpha * (block_counter_ - kNBlocksAverageInitPhase);
218   }
219 }
220 
UpdateBandBySmoothing(float power_band,float power_band_noise,float alpha) const221 float StationarityEstimator::NoiseSpectrum::UpdateBandBySmoothing(
222     float power_band,
223     float power_band_noise,
224     float alpha) const {
225   float power_band_noise_updated = power_band_noise;
226   if (power_band_noise < power_band) {
227     RTC_DCHECK_GT(power_band, 0.f);
228     float alpha_inc = alpha * (power_band_noise / power_band);
229     if (block_counter_ > kNBlocksInitialPhase) {
230       if (10.f * power_band_noise < power_band) {
231         alpha_inc *= 0.1f;
232       }
233     }
234     power_band_noise_updated += alpha_inc * (power_band - power_band_noise);
235   } else {
236     power_band_noise_updated += alpha * (power_band - power_band_noise);
237     power_band_noise_updated =
238         std::max(power_band_noise_updated, kMinNoisePower);
239   }
240   return power_band_noise_updated;
241 }
242 
243 }  // namespace webrtc
244