Finalization of the first version of EchoCanceller 3

webrtc/modules/audio_processing/aec3/suppression_gain.cc

+/*
+ *  Copyright (c) 2017 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 "webrtc/modules/audio_processing/aec3/suppression_gain.h"
+
+#include "webrtc/typedefs.h"
+#if defined(WEBRTC_ARCH_X86_FAMILY)
+#include <emmintrin.h>
+#endif
+#include <math.h>
+#include <algorithm>
+#include <functional>
+
+namespace webrtc {
+namespace {
+
+constexpr int kNumIterations = 2;
+constexpr float kEchoMaskingMargin = 1.f / 10.f;
+constexpr float kBandMaskingFactor = 1.f / 3.f;
+constexpr float kTimeMaskingFactor = 1.f / 10.f;
+
+}  // namespace
+
+#if defined(WEBRTC_ARCH_X86_FAMILY)
+
+// Optimized SSE2 code for the gain computation.
+// TODO(peah): Add further optimizations, in particular for the divisions.
+void ComputeGains_SSE2(
+    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, kFftLengthBy2 - 1>* previous_gain_squared,

[hlundin-webrtc, 2017/02/21 16:34:03]

You already have kFftLengthBy2Plus1 defined. Why not add kFftLengthBy2Minus1?
You would use it a lot in this file at least.

[peah-webrtc, 2017/02/21 23:00:41]

Nice!
Done.

+    std::array<float, kFftLengthBy2 - 1>* previous_masker,
+    std::array<float, kFftLengthBy2Plus1>* gain) {
+  std::array<float, kFftLengthBy2 - 1> masker;
+  std::array<float, kFftLengthBy2 - 1> same_band_masker;
+  std::array<float, kFftLengthBy2 - 1> one_by_residual_echo_power;
+  std::array<bool, kFftLengthBy2 - 1> strong_nearend;
+  std::array<float, kFftLengthBy2Plus1> neighboring_bands_masker;
+
+  // 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());
+    } else {
+      // Add masker for neightboring bands.
+      std::transform(nearend_power.begin(), nearend_power.end(), gain->begin(),
+                     neighboring_bands_masker.begin(),
+                     std::multiplies<float>());
+      std::transform(neighboring_bands_masker.begin(),
+                     neighboring_bands_masker.end(),
+                     comfort_noise_power.begin(),
+                     neighboring_bands_masker.begin(), std::plus<float>());
+      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.
+      std::transform(same_band_masker.begin(), same_band_masker.end(),
+                     masker.begin(), masker.begin(), std::plus<float>());
+    }
+
+    // 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->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->begin() + 1, gain->end() - 1, strong_nearend.begin(),
+                   gain->begin() + 1,
+                   [](float a, bool b) { return b ? 1.f : a; });
+
+    // Limit the allowed gain update over time.
+    std::transform(
+        gain->begin() + 1, gain->end() - 1, previous_gain_squared->begin(),
+        gain->begin() + 1, [](float a, float b) {
+          return b < 0.0001f ? std::min(a, 0.0001f) : std::min(a, b * 2.f);
+        });
+
+    (*gain)[0] = (*gain)[1];
+    (*gain)[kFftLengthBy2] = (*gain)[kFftLengthBy2 - 1];
+  }
+
+  std::copy(gain->begin() + 1, gain->end() - 1, previous_gain_squared->begin());
+
+  std::transform(gain->begin() + 1, gain->end() - 1, nearend_power.begin() + 1,
+                 previous_masker->begin(), std::multiplies<float>());
+  std::transform(previous_masker->begin(), previous_masker->end(),
+                 comfort_noise_power.begin() + 1, previous_masker->begin(),
+                 std::plus<float>());
+
+  for (size_t k = 0; k < kFftLengthBy2; k += 4) {
+    __m128 g = _mm_loadu_ps(&(*gain)[k]);
+    g = _mm_sqrt_ps(g);
+    _mm_storeu_ps(&(*gain)[k], g);
+  }
+
+  (*gain)[kFftLengthBy2] = sqrtf((*gain)[kFftLengthBy2]);
+}
+
+#endif
+
+void ComputeGains(
+    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, kFftLengthBy2 - 1>* previous_gain_squared,
+    std::array<float, kFftLengthBy2 - 1>* previous_masker,
+    std::array<float, kFftLengthBy2Plus1>* gain) {
+  std::array<float, kFftLengthBy2 - 1> masker;
+  std::array<float, kFftLengthBy2 - 1> same_band_masker;
+  std::array<float, kFftLengthBy2 - 1> one_by_residual_echo_power;
+  std::array<bool, kFftLengthBy2 - 1> strong_nearend;
+  std::array<float, kFftLengthBy2Plus1> neighboring_bands_masker;
+
+  // 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; });

[hlundin-webrtc, 2017/02/21 16:34:03]

Is there any risk that we get _really_ large results here, when a is very close
to 0? Should you compare with some threshold larger than 0? Just a thought, it
is up to you.
EDIT: I see on line 211 that the subsequent processing with this result is
capped by 1. So I guess this is not an issue.

[peah-webrtc, 2017/02/21 23:00:41]

Yes, the division should actually be done as it is. But I think the separation
of doing this in separate steps is making this more difficult that it actually
is. This code should be optimized in upcoming CLs. Let's rewrite this part then.

+
+  // 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());
+    } else {
+      // Add masker for neightboring bands.

[hlundin-webrtc, 2017/02/21 16:34:03]

neighboring

[peah-webrtc, 2017/02/21 23:00:41]

Done.

+      std::transform(nearend_power.begin(), nearend_power.end(), gain->begin(),
+                     neighboring_bands_masker.begin(),
+                     std::multiplies<float>());
+      std::transform(neighboring_bands_masker.begin(),
+                     neighboring_bands_masker.end(),
+                     comfort_noise_power.begin(),
+                     neighboring_bands_masker.begin(), std::plus<float>());
+      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.
+      std::transform(same_band_masker.begin(), same_band_masker.end(),
+                     masker.begin(), masker.begin(), std::plus<float>());
+    }
+
+    // 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->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->begin() + 1, gain->end() - 1, strong_nearend.begin(),
+                   gain->begin() + 1,
+                   [](float a, bool b) { return b ? 1.f : a; });
+
+    // Limit the allowed gain update over time.
+    std::transform(
+        gain->begin() + 1, gain->end() - 1, previous_gain_squared->begin(),
+        gain->begin() + 1, [](float a, float b) {
+          return b < 0.0001f ? std::min(a, 0.0001f) : std::min(a, b * 2.f);
+        });
+
+    (*gain)[0] = (*gain)[1];
+    (*gain)[kFftLengthBy2] = (*gain)[kFftLengthBy2 - 1];
+  }
+
+  std::copy(gain->begin() + 1, gain->end() - 1, previous_gain_squared->begin());

[hlundin-webrtc, 2017/02/21 16:34:03]

It's a bit confusing that you copy to "previous_gain_squared" from "gain"
(presumably not squared). Can you align the naming somewhat?

[peah-webrtc, 2017/02/21 23:00:41]

I simplified this with an alias.
PTAL.

[hlundin-webrtc, 2017/02/22 15:24:05]

Better.

[peah-webrtc, 2017/02/22 23:51:36]

Acknowledged.

+
+  std::transform(gain->begin() + 1, gain->end() - 1, nearend_power.begin() + 1,
+                 previous_masker->begin(), std::multiplies<float>());
+  std::transform(previous_masker->begin(), previous_masker->end(),
+                 comfort_noise_power.begin() + 1, previous_masker->begin(),
+                 std::plus<float>());
+
+  std::transform(gain->begin(), gain->end(), gain->begin(),
+                 [](float a) { return sqrtf(a); });
+}
+
+SuppressionGain::SuppressionGain(Aec3Optimization optimization)
+    : optimization_(optimization) {
+  previous_gain_squared_.fill(1.f);
+  previous_masker_.fill(0.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,
+    float strong_nearend_margin,
+    std::array<float, kFftLengthBy2Plus1>* gain) {
+  RTC_DCHECK(gain);
+  if (optimization_ == Aec3Optimization::kNone) {

[hlundin-webrtc, 2017/02/21 16:34:03]

I think you should use switch instead:
switch (optimization_) {
#if defined(WEBRTC_ARCH_X86_FAMILY)
  case Aec3Optimization::kSse2: {
    ComputeGains_SSE2(...);
  }
#endif
  default: {
    ComputeGains(...);
  }
}

This has the benefit that if you add another type of optimization, there is no
risk of exiting this method without having computed anything at all. As an
alternative, you can omit the default case, and explicitly name both
alternatives. That way, the compiler will help you find all places that needs
attention when you add another optimization type.

[peah-webrtc, 2017/02/21 23:00:41]

Great one!
I also changed to this in the other places using optimizations.

[hlundin-webrtc, 2017/02/22 15:24:05]

That was the way I hoped it would pan out. :)

[peah-webrtc, 2017/02/23 11:18:32]

:-)
Acknowledged.

+    ComputeGains(nearend_power, residual_echo_power, comfort_noise_power,
+                 strong_nearend_margin, &previous_gain_squared_,
+                 &previous_masker_, gain);
+  }
+#if defined(WEBRTC_ARCH_X86_FAMILY)
+  if (optimization_ == Aec3Optimization::kSse2) {
+    ComputeGains_SSE2(nearend_power, residual_echo_power, comfort_noise_power,
+                      strong_nearend_margin, &previous_gain_squared_,
+                      &previous_masker_, gain);
+  }
+#endif
+}
+
+}  // namespace webrtc