Transparency increasing tuning for AEC3

webrtc/modules/audio_processing/aec3/suppression_gain.cc

Index: webrtc/modules/audio_processing/aec3/suppression_gain.cc
diff --git a/webrtc/modules/audio_processing/aec3/suppression_gain.cc b/webrtc/modules/audio_processing/aec3/suppression_gain.cc
index 86af60f316fa867e6f2276c8aac7311fd867fa44..8918385127ab0c03a3c6ec6a86df6d33c79bfbc8 100644
--- a/webrtc/modules/audio_processing/aec3/suppression_gain.cc
+++ b/webrtc/modules/audio_processing/aec3/suppression_gain.cc
@@ -25,183 +25,244 @@
 namespace webrtc {
 namespace {
 
-void GainPostProcessing(std::array<float, kFftLengthBy2Plus1>* gain_squared) {
+// Adjust the gains according to the presense of known external filters.

[ivoc, 2017/05/22 14:58:00]

I think it should be "presence".

[peah-webrtc, 2017/05/22 21:52:47]

Done.

+void AdjustForExternalFilters(std::array<float, kFftLengthBy2Plus1>* gain) {
   // Limit the low frequency gains to avoid the impact of the high-pass filter
   // on the lower-frequency gain influencing the overall achieved gain.
-  (*gain_squared)[1] = std::min((*gain_squared)[1], (*gain_squared)[2]);
-  (*gain_squared)[0] = (*gain_squared)[1];
+  (*gain)[0] = (*gain)[1] = std::min((*gain)[1], (*gain)[2]);
 
   // Limit the high frequency gains to avoid the impact of the anti-aliasing
   // filter on the upper-frequency gains influencing the overall achieved
   // gain. TODO(peah): Update this when new anti-aliasing filters are
   // implemented.
   constexpr size_t kAntiAliasingImpactLimit = (64 * 2000) / 8000;
-  std::for_each(gain_squared->begin() + kAntiAliasingImpactLimit,
-                gain_squared->end() - 1,
-                [gain_squared, kAntiAliasingImpactLimit](float& a) {
-                  a = std::min(a, (*gain_squared)[kAntiAliasingImpactLimit]);
-                });
-  (*gain_squared)[kFftLengthBy2] = (*gain_squared)[kFftLengthBy2Minus1];
+  const float min_upper_gain = (*gain)[kAntiAliasingImpactLimit];
+  std::for_each(
+      gain->begin() + kAntiAliasingImpactLimit, gain->end() - 1,
+      [min_upper_gain](float& a) { a = std::min(a, min_upper_gain); });
+  (*gain)[kFftLengthBy2] = (*gain)[kFftLengthBy2Minus1];
 }
 
-constexpr int kNumIterations = 2;
-constexpr float kEchoMaskingMargin = 1.f / 20.f;
-constexpr float kBandMaskingFactor = 1.f / 10.f;
-constexpr float kTimeMaskingFactor = 1.f / 10.f;
+// Computes the gain to apply for the bands beyond the first band.
+float UpperBandsGain(
+    bool saturated_echo,
+    const std::vector<std::vector<float>>& render,
+    const std::array<float, kFftLengthBy2Plus1>& low_band_gain) {
+  if (render.size() == 1) {

[ivoc, 2017/05/22 14:58:00]

Should there be a DHECK for render.size() > 0?

[peah-webrtc, 2017/05/22 21:52:47]

Good point!
Done.

+    return 1.f;
+  }
 
-// TODO(peah): Add further optimizations, in particular for the divisions.
-void ComputeGains(
-    Aec3Optimization optimization,
-    const std::array<float, kFftLengthBy2Plus1>& nearend_power,
-    const std::array<float, kFftLengthBy2Plus1>& residual_echo_power,
-    const std::array<float, kFftLengthBy2Plus1>& comfort_noise_power,
-    float strong_nearend_margin,
-    std::array<float, kFftLengthBy2Minus1>* previous_gain_squared,
-    std::array<float, kFftLengthBy2Minus1>* previous_masker,
+  const float gain_below_8_khz =
+      *std::min_element(low_band_gain.begin() + 32, low_band_gain.end());

[ivoc, 2017/05/22 14:58:00]

I think the constant (32) should be derived from kFftLengthBy2, either by
introducing kFftLengthBy4, or by using kFftLengthBy2 / 2.

[peah-webrtc, 2017/05/22 21:52:47]

Done.

+
+  // Always attenuate the upper bands when there is saturated echo.
+  if (saturated_echo) {
+    return std::min(0.001f, gain_below_8_khz);
+  }
+
+  // Compute the upper and lower band energies.
+  auto sum_of_squares = [](float a, float b) { return a + b * b; };

[aleloi, 2017/05/22 20:41:45]

const auto

[peah-webrtc, 2017/05/22 21:52:47]

Done.

+  const float low_band_energy =
+      std::accumulate(render[0].begin(), render[0].end(), 0.f, sum_of_squares);
+  float high_band_energy = 0.f;
+  for (size_t k = 1; k < render.size(); ++k) {
+    const float energy = std::accumulate(render[k].begin(), render[k].end(),
+                                         0.f, sum_of_squares);
+    high_band_energy = std::max(high_band_energy, energy);
+  }
+
+  // If there is more power in the lower frequencies than the upper frequencies,
+  // or if the power in upper frequencies is low, do not bound the gain in the
+  // upper bands.
+  float anti_howling_gain;
+  constexpr float kThreshold = kSubBlockSize * 10.f * 10.f;
+  if (high_band_energy < std::max(low_band_energy, kThreshold)) {
+    anti_howling_gain = 1.f;
+  } else {
+    // In all other cases, bound the gain for upper frequencies.
+    RTC_DCHECK_LE(low_band_energy, high_band_energy);
+    RTC_DCHECK_NE(0.f, high_band_energy);
+    anti_howling_gain = 0.01f * sqrtf(low_band_energy / high_band_energy);
+  }
+
+  // Choose the gain as the minimum of the lower and upper gains.
+  return std::min(gain_below_8_khz, anti_howling_gain);
+}
+
+// Limits the gain increase.
+void UpdateMaxGainIncrease(
+    size_t no_saturation_counter,
+    bool low_noise_render,
+    const std::array<float, kFftLengthBy2Plus1>& last_echo,
+    const std::array<float, kFftLengthBy2Plus1>& echo,
+    const std::array<float, kFftLengthBy2Plus1>& last_gain,
+    const std::array<float, kFftLengthBy2Plus1>& new_gain,
+    std::array<float, kFftLengthBy2Plus1>* gain_increase) {
+  float max_increasing;
+  float max_decreasing;
+  float rate_increasing;
+  float rate_decreasing;
+  float min_increasing;
+  float min_decreasing;
+
+  if (low_noise_render) {
+    max_increasing = 8.f;
+    max_decreasing = 8.f;
+    rate_increasing = 2.f;
+    rate_decreasing = 2.f;
+    min_increasing = 4.f;
+    min_decreasing = 4.f;
+  } else if (no_saturation_counter > 10) {
+    max_increasing = 4.f;
+    max_decreasing = 4.f;
+    rate_increasing = 2.f;
+    rate_decreasing = 2.f;
+    min_increasing = 1.2f;
+    min_decreasing = 2.f;

[ivoc, 2017/05/22 14:58:00]

Just checking, is this 2.f correct or a typo? (it's the only value that breaks
the pattern)

[peah-webrtc, 2017/05/22 21:52:47]

No, it should be correct. The idea when there is saturation is that the minimum
should be lower if the echo is increasing.

+  } else {
+    max_increasing = 1.2f;
+    max_decreasing = 1.2f;
+    rate_increasing = 1.5f;
+    rate_decreasing = 1.5f;
+    min_increasing = 1.f;
+    min_decreasing = 1.f;
+  }
+
+  for (size_t k = 0; k < new_gain.size(); ++k) {
+    if (echo[k] > last_echo[k]) {
+      (*gain_increase)[k] =
+          new_gain[k] > last_gain[k]
+              ? std::min(max_increasing, (*gain_increase)[k] * rate_increasing)
+              : min_increasing;
+    } else {
+      (*gain_increase)[k] =
+          new_gain[k] > last_gain[k]
+              ? std::min(max_decreasing, (*gain_increase)[k] * rate_decreasing)
+              : min_decreasing;
+    }
+  }
+}
+
+// Computes the gain to reduce the echo to a non audible level.
+void GainToNoAudibleEcho(
+    bool low_noise_render,
+    bool saturated_echo,
+    const std::array<float, kFftLengthBy2Plus1>& nearend,
+    const std::array<float, kFftLengthBy2Plus1>& echo,
+    const std::array<float, kFftLengthBy2Plus1>& masker,
+    const std::array<float, kFftLengthBy2Plus1>& min_gain,
+    const std::array<float, kFftLengthBy2Plus1>& max_gain,
+    const std::array<float, kFftLengthBy2Plus1>& one_by_echo,
     std::array<float, kFftLengthBy2Plus1>* gain) {
-  std::array<float, kFftLengthBy2Minus1> masker;
-  std::array<float, kFftLengthBy2Minus1> same_band_masker;
-  std::array<float, kFftLengthBy2Minus1> one_by_residual_echo_power;
-  std::array<bool, kFftLengthBy2Minus1> strong_nearend;
-  std::array<float, kFftLengthBy2Plus1> neighboring_bands_masker;
-  std::array<float, kFftLengthBy2Plus1>* gain_squared = gain;
-  aec3::VectorMath math(optimization);
-
-  // Precompute 1/residual_echo_power.
-  std::transform(residual_echo_power.begin() + 1, residual_echo_power.end() - 1,
-                 one_by_residual_echo_power.begin(),
-                 [](float a) { return a > 0.f ? 1.f / a : -1.f; });
-
-  // Precompute indicators for bands with strong nearend.
-  std::transform(
-      residual_echo_power.begin() + 1, residual_echo_power.end() - 1,
-      nearend_power.begin() + 1, strong_nearend.begin(),
-      [&](float a, float b) { return a <= strong_nearend_margin * b; });
-
-  //  Precompute masker for the same band.
-  std::transform(comfort_noise_power.begin() + 1, comfort_noise_power.end() - 1,
-                 previous_masker->begin(), same_band_masker.begin(),
-                 [&](float a, float b) { return a + kTimeMaskingFactor * b; });
-
-  for (int k = 0; k < kNumIterations; ++k) {
-    if (k == 0) {
-      // Add masker from the same band.
-      std::copy(same_band_masker.begin(), same_band_masker.end(),
-                masker.begin());
+  constexpr float kEchoMaskingMargin = 1.f / 100.f;
+  const float echo_masking_margin =
+      low_noise_render ? 0.01f : kEchoMaskingMargin;

[ivoc, 2017/05/22 14:58:00]

This seems to have no effect, kEchoMaskingMargin is 0.01f, so it gets the same
value regardless of the value of low_noise_render.

[peah-webrtc, 2017/05/22 21:52:47]

Oups! Good find! Removed that.

Done.

+  const float nearend_masking_margin =
+      low_noise_render ? 2.f : (saturated_echo ? 0.001f : 0.01f);
+
+  for (size_t k = 0; k < gain->size(); ++k) {
+    RTC_DCHECK_LE(0.f, nearend_masking_margin * nearend[k]);
+    if (echo[k] <= nearend_masking_margin * nearend[k]) {
+      (*gain)[k] = 1.f;
     } else {
-      // Add masker for neighboring bands.
-      math.Multiply(nearend_power, *gain_squared, neighboring_bands_masker);
-      math.Accumulate(comfort_noise_power, neighboring_bands_masker);
-      std::transform(
-          neighboring_bands_masker.begin(), neighboring_bands_masker.end() - 2,
-          neighboring_bands_masker.begin() + 2, masker.begin(),
-          [&](float a, float b) { return kBandMaskingFactor * (a + b); });
-
-      // Add masker from the same band.
-      math.Accumulate(same_band_masker, masker);
+      (*gain)[k] = echo_masking_margin * masker[k] * one_by_echo[k];
     }
 
-    // Compute new gain as:
-    // G2(t,f) = (comfort_noise_power(t,f) + G2(t-1)*nearend_power(t-1)) *
-    //           kTimeMaskingFactor
-    //           * kEchoMaskingMargin / residual_echo_power(t,f).
-    // or
-    // G2(t,f) = ((comfort_noise_power(t,f) + G2(t-1) *
-    //             nearend_power(t-1)) * kTimeMaskingFactor +
-    //            (comfort_noise_power(t, f-1) + comfort_noise_power(t, f+1) +
-    //             (G2(t,f-1)*nearend_power(t, f-1) +
-    //              G2(t,f+1)*nearend_power(t, f+1)) *
-    //             kTimeMaskingFactor) * kBandMaskingFactor)
-    //           * kEchoMaskingMargin / residual_echo_power(t,f).
-    std::transform(
-        masker.begin(), masker.end(), one_by_residual_echo_power.begin(),
-        gain_squared->begin() + 1, [&](float a, float b) {
-          return b >= 0 ? std::min(kEchoMaskingMargin * a * b, 1.f) : 1.f;
-        });
-
-    // Limit gain for bands with strong nearend.
-    std::transform(gain_squared->begin() + 1, gain_squared->end() - 1,
-                   strong_nearend.begin(), gain_squared->begin() + 1,
-                   [](float a, bool b) { return b ? 1.f : a; });
-
-    // Limit the allowed gain update over time.
-    std::transform(gain_squared->begin() + 1, gain_squared->end() - 1,
-                   previous_gain_squared->begin(), gain_squared->begin() + 1,
-                   [](float a, float b) {
-                     return b < 0.001f ? std::min(a, 0.001f)
-                                       : std::min(a, b * 2.f);
-                   });
-
-    // Process the gains to avoid artefacts caused by gain realization in the
-    // filterbank and impact of external pre-processing of the signal.
-    GainPostProcessing(gain_squared);
+    (*gain)[k] = std::min(std::max((*gain)[k], min_gain[k]), max_gain[k]);
   }
+}
 
-  std::copy(gain_squared->begin() + 1, gain_squared->end() - 1,
-            previous_gain_squared->begin());
-
-  math.Multiply(
-      rtc::ArrayView<const float>(&(*gain_squared)[1], previous_masker->size()),
-      rtc::ArrayView<const float>(&nearend_power[1], previous_masker->size()),
-      *previous_masker);
-  math.Accumulate(rtc::ArrayView<const float>(&comfort_noise_power[1],
-                                              previous_masker->size()),
-                  *previous_masker);
-  math.Sqrt(*gain);
+// Computes the signal output power that masks the echo signal.
+void MaskingPower(const std::array<float, kFftLengthBy2Plus1>& nearend,
+                  const std::array<float, kFftLengthBy2Plus1>& comfort_noise,
+                  const std::array<float, kFftLengthBy2Plus1>& last_masker,
+                  const std::array<float, kFftLengthBy2Plus1>& gain,
+                  std::array<float, kFftLengthBy2Plus1>* masker) {
+  std::array<float, kFftLengthBy2Plus1> side_band_masker;
+  for (size_t k = 0; k < gain.size(); ++k) {
+    side_band_masker[k] = nearend[k] * gain[k] + comfort_noise[k];
+    (*masker)[k] = comfort_noise[k] + 0.1f * last_masker[k];
+  }
+  for (size_t k = 1; k < gain.size() - 1; ++k) {

[ivoc, 2017/05/22 14:58:00]

It's possible to merge this loop with the previous one, which I expect to be a
bit faster, but I'm not sure if this code is called often enough for it to be
worth it, so I'll leave that up to you. You could do it like this (using a
lambda would be useful if you decide to do this): 

Handle k == 0
for (k = 1; k<gain.size() - 1; ++k) {
  ...
}
Handle k == gain.size() - 1

[peah-webrtc, 2017/05/22 21:52:47]

I considered that, but it makes the code somewhat longer so I decided to put the
loops as this. I think they should all fit into the cache line as they are
fairly short.

+    (*masker)[k] += 0.1f * (side_band_masker[k - 1] + side_band_masker[k + 1]);
+  }
 }
 
 }  // namespace
 
-// Computes an upper bound on the gain to apply for high frequencies.
-float HighFrequencyGainBound(bool saturated_echo,
-                             const std::vector<std::vector<float>>& render) {
-  if (render.size() == 1) {
-    return 1.f;
+// TODO(peah): Add further optimizations, in particular for the divisions.
+void SuppressionGain::LowerBandGain(
+    bool low_noise_render,
+    bool saturated_echo,
+    const std::array<float, kFftLengthBy2Plus1>& nearend,
+    const std::array<float, kFftLengthBy2Plus1>& echo,
+    const std::array<float, kFftLengthBy2Plus1>& comfort_noise,
+    std::array<float, kFftLengthBy2Plus1>* gain) {
+  // Count the number of blocks since saturation.
+  no_saturation_counter_ = saturated_echo ? 0 : no_saturation_counter_ + 1;
+
+  // Precompute 1/echo.
+  std::array<float, kFftLengthBy2Plus1> one_by_echo;
+  std::transform(echo.begin(), echo.end(), one_by_echo.begin(),
+                 [](float a) { return a > 0.f ? 1.f / a : 1.f; });

[ivoc, 2017/05/22 14:58:00]

Why is the default value 1.f? Doesn't that introduce a discontinuity? i.e. if a
== 0.f: one_by_echo will be 1.f, but if a == 0.001f, one_by_echo will be 1000.f.

[aleloi, 2017/05/22 20:41:45]

Another related issue: perhaps we should avoid dividing by echo when the echo is
< a_small_threshold_number, and not when it's == 0. To avoid very large numbers
and discontinuities as Ivo pointed out.

[peah-webrtc, 2017/05/22 21:52:47]

This is not as big a problem as it may seem. If 1/echo is large, the resulting
gain will exceed 1 which means it then will be corrected to 1.
It is a bit unfortunate to have to perform the division far from where it is
used, as it then would be clear that this is a no-issue. 

[peah-webrtc, 2017/05/22 21:52:47]

To avoid that we should have one_by_echo=inf if a == 0. However, in case the
denominator is zero, the value is anyway never used. The construct is only there
to ensure that something is returned instead of dividing by zero.

I thought of DCHECKing for ensuring that it is really never used, but that is
not straightforward to do. I'll add a comment instead

+
+  // Compute the minimum gain as the attenuating gain to put the signal just
+  // above the zero sample values.
+  std::array<float, kFftLengthBy2Plus1> min_gain;
+  const float min_echo_power = low_noise_render ? 192.f : 64.f;
+  if (no_saturation_counter_ > 10) {
+    for (size_t k = 0; k < nearend.size(); ++k) {
+      const float denom = std::min(nearend[k], echo[k]);
+      min_gain[k] = denom > 0.f ? min_echo_power / denom : 1.f;
+      min_gain[k] = std::min(min_gain[k], 1.f);
+    }
+  } else {
+    min_gain.fill(0.f);
   }
 
-  // Always attenuate the upper bands when there is saturated echo.
-  if (saturated_echo) {
-    return 0.001f;
+  // Compute the maximum gain by limiting the gain increase from the previous
+  // gain.
+  std::array<float, kFftLengthBy2Plus1> max_gain;
+  for (size_t k = 0; k < gain->size(); ++k) {
+    max_gain[k] =
+        std::min(std::max(last_gain_[k] * gain_increase_[k], 0.001f), 1.f);
   }
 
-  // Compute the upper and lower band energies.
-  float low_band_energy =
-      std::accumulate(render[0].begin(), render[0].end(), 0.f,
-                      [](float a, float b) -> float { return a + b * b; });
-  float high_band_energies = 0.f;
-  for (size_t k = 1; k < render.size(); ++k) {
-    high_band_energies = std::max(
-        high_band_energies,
-        std::accumulate(render[k].begin(), render[k].end(), 0.f,
-                        [](float a, float b) -> float { return a + b * b; }));
+  // Iteratively compute the gain required to attenuate the echo to a non
+  // noticeable level.
+  gain->fill(0.f);
+  for (int k = 0; k < 2; ++k) {
+    std::array<float, kFftLengthBy2Plus1> masker;
+    MaskingPower(nearend, comfort_noise, last_masker_, *gain, &masker);
+    GainToNoAudibleEcho(low_noise_render, saturated_echo, nearend, echo, masker,
+                        min_gain, max_gain, one_by_echo, gain);
+    AdjustForExternalFilters(gain);
   }
 
-  // If there is more power in the lower frequencies than the upper frequencies,
-  // or if the power in upper frequencies  is low, do not bound the gain in the
-  // upper bands.
-  if (high_band_energies < low_band_energy ||
-      high_band_energies < kSubBlockSize * 10.f * 10.f) {
-    return 1.f;
-  }
+  // Update the allowed maximum gain increase.
+  UpdateMaxGainIncrease(no_saturation_counter_, low_noise_render, last_echo_,
+                        echo, last_gain_, *gain, &gain_increase_);
 
-  // In all other cases, bound the gain for upper frequencies.
-  RTC_DCHECK_LE(low_band_energy, high_band_energies);
-  return 0.01f * sqrtf(low_band_energy / high_band_energies);
+  // Store data required for the gain computation of the next block.
+  std::copy(echo.begin(), echo.end(), last_echo_.begin());
+  std::copy(gain->begin(), gain->end(), last_gain_.begin());
+  MaskingPower(nearend, comfort_noise, last_masker_, *gain, &last_masker_);
+  aec3::VectorMath(optimization_).Sqrt(*gain);
 }
 
 SuppressionGain::SuppressionGain(Aec3Optimization optimization)
     : optimization_(optimization) {
-  previous_gain_squared_.fill(1.f);
-  previous_masker_.fill(0.f);
+  last_gain_.fill(1.f);
+  last_masker_.fill(0.f);
+  gain_increase_.fill(1.f);
 }
 
 void SuppressionGain::GetGain(
-    const std::array<float, kFftLengthBy2Plus1>& nearend_power,
-    const std::array<float, kFftLengthBy2Plus1>& residual_echo_power,
-    const std::array<float, kFftLengthBy2Plus1>& comfort_noise_power,
+    const std::array<float, kFftLengthBy2Plus1>& nearend,
+    const std::array<float, kFftLengthBy2Plus1>& echo,
+    const std::array<float, kFftLengthBy2Plus1>& comfort_noise,
     bool saturated_echo,
     const std::vector<std::vector<float>>& render,
-    size_t num_capture_bands,
     bool force_zero_gain,
     float* high_bands_gain,
     std::array<float, kFftLengthBy2Plus1>* low_band_gain) {
@@ -209,32 +270,42 @@ void SuppressionGain::GetGain(
   RTC_DCHECK(low_band_gain);
 
   if (force_zero_gain) {
-    previous_gain_squared_.fill(0.f);
-    std::copy(comfort_noise_power.begin() + 1, comfort_noise_power.end() - 1,
-              previous_masker_.begin());
+    last_gain_.fill(0.f);
+    std::copy(comfort_noise.begin(), comfort_noise.end(),
+              last_masker_.begin() + 1);
     low_band_gain->fill(0.f);
+    gain_increase_.fill(1.f);
     *high_bands_gain = 0.f;
     return;
   }
 
-  // Choose margin to use.
-  const float margin = saturated_echo ? 0.001f : 0.01f;
-  ComputeGains(optimization_, nearend_power, residual_echo_power,
-               comfort_noise_power, margin, &previous_gain_squared_,
-               &previous_masker_, low_band_gain);
+  bool low_noise_render = low_render_detector_.Detect(render);
 
-  if (num_capture_bands > 1) {
-    // Compute the gain for upper frequencies.
-    const float min_high_band_gain =
-        HighFrequencyGainBound(saturated_echo, render);
-    *high_bands_gain =
-        *std::min_element(low_band_gain->begin() + 32, low_band_gain->end());
+  // Compute gain for the lower band.
+  LowerBandGain(low_noise_render, saturated_echo, nearend, echo, comfort_noise,
+                low_band_gain);
 
-    *high_bands_gain = std::min(*high_bands_gain, min_high_band_gain);
+  // Compute the gain for the upper bands.
+  *high_bands_gain = UpperBandsGain(saturated_echo, render, *low_band_gain);
+}
 
-  } else {
-    *high_bands_gain = 1.f;
+// Detects when the render signal can be considered to have low power and
+// consist of stationary noise.
+bool SuppressionGain::LowNoiseRenderDetector::Detect(
+    const std::vector<std::vector<float>>& render) {
+  float x2_sum = 0.f;
+  float x2_max = 0.f;
+  for (auto x_k : render[0]) {
+    const float x2 = x_k * x_k;
+    x2_sum += x2;
+    x2_max = std::max(x2_max, x2);
   }
+
+  constexpr float kThreshold = 50.f * 50.f * 64.f;
+  const bool low_noise_render =
+      average_power_ < kThreshold && x2_max < 3 * average_power_;
+  average_power_ = average_power_ * 0.9f + x2_sum * 0.1f;
+  return low_noise_render;
 }
 
 }  // namespace webrtc