Finalization of the first version of EchoCanceller 3

webrtc/modules/audio_processing/aec3/adaptive_fir_filter.cc

Index: webrtc/modules/audio_processing/aec3/adaptive_fir_filter.cc
diff --git a/webrtc/modules/audio_processing/aec3/adaptive_fir_filter.cc b/webrtc/modules/audio_processing/aec3/adaptive_fir_filter.cc
new file mode 100644
index 0000000000000000000000000000000000000000..9d97fefe52e295528fd0e25d22b234f153cb9816
--- /dev/null
+++ b/webrtc/modules/audio_processing/aec3/adaptive_fir_filter.cc
@@ -0,0 +1,159 @@
+/*
+ *  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/adaptive_fir_filter.h"
+
+#include <algorithm>
+#include <functional>
+
+#include "webrtc/base/checks.h"
+#include "webrtc/modules/audio_processing/aec3/fft_data.h"
+
+namespace webrtc {
+
+namespace {
+
+// Constrains an adaptive filter partition.

[hlundin-webrtc, 2017/02/13 21:36:59]

Comment on what this does mathematically.
IFFT of H -> h
Zero out upper part(?) of h
FFT h -> H

[peah-webrtc, 2017/02/20 07:37:14]

Done.

+void Constrain(const Aec3Fft& fft, FftData* H) {
+  std::array<float, kFftLength> h;
+  fft.Ifft(*H, &h);
+  constexpr float kScale = 1.0f / kFftLengthBy2;
+  std::for_each(h.begin(), h.begin() + kFftLengthBy2,
+                [&](float& a) { a *= kScale; });

[hlundin-webrtc, 2017/02/13 21:36:59]

Is the default capture [&] necessary for the lambda to pick up kScale?

[peah-webrtc, 2017/02/20 07:37:13]

No, it was not it turned out.

Done.

+  std::fill(h.begin() + kFftLengthBy2, h.end(), 0.f);
+  fft.Fft(&h, H);
+}
+
+// Adapts the filter partitions.

[hlundin-webrtc, 2017/02/13 21:36:59]

And what is the math?

[hlundin-webrtc, 2017/02/13 21:36:59]

What is the relationship between X, G, and H?

[peah-webrtc, 2017/02/20 07:37:13]

I included the math.

PTAL if you think that is sufficient!

[peah-webrtc, 2017/02/20 07:37:13]

Done.

+void AdaptPartitions(const FftBuffer& X_buffer,
+                     const FftData& G,
+                     rtc::ArrayView<FftData> H) {
+  rtc::ArrayView<const FftData> X_buffer_data = X_buffer.Buffer();
+  size_t index = X_buffer.Position();
+  for (auto& H_j : H) {
+    const FftData& X = X_buffer_data[index];
+    for (size_t k = 0; k < kFftLengthBy2Plus1; ++k) {
+      H_j.re[k] += X.re[k] * G.re[k] + X.im[k] * G.im[k];

[hlundin-webrtc, 2017/02/13 21:36:59]

Please, comment on why the signs are this way. Normally, one would expect
x * y = (xr + i*xi) * (yr + i*yi) = (xr*yr - xi*yi) + i*(xr*yi + xi*yr)

[peah-webrtc, 2017/02/20 07:37:14]

It is the conjugate, I think. I added that to the equation.

PTAL

[hlundin-webrtc, 2017/02/21 09:40:03]

Acknowledged.

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

Acknowledged.

+      H_j.im[k] += X.re[k] * G.im[k] - X.im[k] * G.re[k];
+    }
+
+    index = index < (X_buffer_data.size() - 1) ? index + 1 : 0;
+  }
+}
+
+// Produces the filter output.
+void ApplyFilter(const FftBuffer& X_buffer,
+                 rtc::ArrayView<const FftData> H,
+                 FftData* S) {
+  S->re.fill(0.f);
+  S->im.fill(0.f);
+
+  rtc::ArrayView<const FftData> X_buffer_data = X_buffer.Buffer();
+  size_t index = X_buffer.Position();
+  for (auto& H_j : H) {
+    const FftData& X = X_buffer_data[index];
+    for (size_t k = 0; k < kFftLengthBy2Plus1; ++k) {
+      S->re[k] += X.re[k] * H_j.re[k] - X.im[k] * H_j.im[k];
+      S->im[k] += X.re[k] * H_j.im[k] + X.im[k] * H_j.re[k];
+    }
+    index = index < (X_buffer_data.size() - 1) ? index + 1 : 0;
+  }
+}
+
+// Computes and stores the frequency response of the filter.
+void UpdateFrequencyResponse(
+    rtc::ArrayView<const FftData> H,
+    std::vector<std::array<float, kFftLengthBy2Plus1>>* H2) {
+  for (size_t k = 0; k < H.size(); ++k) {
+    std::transform(H[k].re.begin(), H[k].re.end(), H[k].im.begin(),
+                   (*H2)[k].begin(),
+                   [](float a, float b) { return a * a + b * b; });
+  }
+}
+
+// Computes and stores the echo return loss of the filter.

[hlundin-webrtc, 2017/02/13 21:36:59]

Show me the math.

[peah-webrtc, 2017/02/20 07:37:14]

I added a comment describing what it is. PTAL

+void UpdateErlEstimator(
+    const std::vector<std::array<float, kFftLengthBy2Plus1>>& H2,
+    std::array<float, kFftLengthBy2Plus1>* erl) {
+  erl->fill(0.f);
+  for (auto& H2_j : H2) {
+    std::transform(H2_j.begin(), H2_j.end(), erl->begin(), erl->begin(),
+                   std::plus<float>());
+  }
+}
+
+// Resets the filter.
+void ResetFilter(rtc::ArrayView<FftData> H) {
+  for (auto& H_j : H) {
+    H_j.Clear();
+  }
+}
+
+}  // namespace
+
+int AdaptiveFirFilter::instance_index_ = -1;
+
+AdaptiveFirFilter::AdaptiveFirFilter(size_t size_partitions,
+                                     bool use_filter_statistics,
+                                     ApmDataDumper* data_dumper)
+    : data_dumper_(data_dumper), H_(size_partitions) {
+  RTC_DCHECK(data_dumper_);
+  ResetFilter(H_);
+
+  if (use_filter_statistics) {
+    H2_.reset(new std::vector<std::array<float, kFftLengthBy2Plus1>>(
+        size_partitions, std::array<float, kFftLengthBy2Plus1>()));
+    for (auto H2_k : *H2_) {
+      H2_k.fill(0.f);
+    }
+
+    erl_.reset(new std::array<float, kFftLengthBy2Plus1>());
+    erl_->fill(0.f);
+  }
+  ++instance_index_;

[hlundin-webrtc, 2017/02/13 21:36:59]

Atomic operation?

[peah-webrtc, 2017/02/20 07:37:13]

It could be removed. 

Done.

+}
+
+AdaptiveFirFilter::~AdaptiveFirFilter() = default;
+
+void AdaptiveFirFilter::HandleEchoPathChange() {
+  ResetFilter(H_);
+  if (H2_) {
+    for (auto H2_k : *H2_) {
+      H2_k.fill(0.f);
+    }
+    RTC_DCHECK(erl_);
+    erl_->fill(0.f);
+  }
+}
+
+void AdaptiveFirFilter::Filter(const FftBuffer& X_buffer, FftData* S) const {
+  ApplyFilter(X_buffer, H_, S);

[hlundin-webrtc, 2017/02/13 21:36:59]

Why the extra function? You can implement it here.

[peah-webrtc, 2017/02/20 07:37:14]

The reason was that I was planning for it to be optimized. That optimization has
now been added.

PTAL.

[hlundin-webrtc, 2017/02/21 09:40:03]

Acknowledged.

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

Acknowledged.

+}
+
+void AdaptiveFirFilter::Adapt(const FftBuffer& X_buffer, const FftData& G) {
+  // Adapt the filter.
+  AdaptPartitions(X_buffer, G, H_);
+
+  // Constrain the filter partitions in a cyclic manner.
+  Constrain(fft_, &H_[partition_to_constrain_]);
+  partition_to_constrain_ = partition_to_constrain_ < (H_.size() - 1)
+                                ? partition_to_constrain_ + 1
+                                : 0;
+
+  // Optionally update the frequency response and echo return loss for the
+  // filter.
+  if (H2_) {
+    RTC_DCHECK(erl_);
+    UpdateFrequencyResponse(H_, H2_.get());
+    UpdateErlEstimator(*H2_, erl_.get());
+  }
+}
+
+}  // namespace webrtc