Adding full initial version of delay estimation functionality in echo canceller 3

webrtc/modules/audio_processing/aec3/correlator.cc

+/*

[hlundin-webrtc, 2017/02/01 08:30:00]

This file needs:
1. Reference to or comments containing equations to explain the basic algorithm.
2. Shorter variable names, matching equations.
3. More comments on the sub-steps.
4. More DCHECK for documenting the expected outputs of various sub-steps.

[peah-webrtc, 2017/02/02 14:04:45]

Done.

+ *  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/correlator.h"
+
+#include <algorithm>
+
+#include "webrtc/modules/audio_processing/aec3/aec3_constants.h"
+#include "webrtc/modules/audio_processing/include/audio_processing.h"
+
+namespace webrtc {
+
+Correlator::IndexedBuffer::IndexedBuffer(size_t size) : data(size, 0.f) {}
+
+Correlator::IndexedBuffer::~IndexedBuffer() = default;
+
+Correlator::Correlator(ApmDataDumper* data_dumper,
+                       size_t window_size_sub_blocks,
+                       size_t num_alignment_shifts,
+                       size_t alignment_shift_sub_blocks)
+    : data_dumper_(data_dumper),
+      correlator_intra_lag_shift_(alignment_shift_sub_blocks * kSubBlockSize),
+      correlations_(
+          num_alignment_shifts + 1,
+          std::vector<float>(window_size_sub_blocks * kSubBlockSize, 0.f)),
+      lag_estimates_(num_alignment_shifts + 1),
+      render_buffer_(kSubBlockSize *
+                     (alignment_shift_sub_blocks * num_alignment_shifts + 1 +
+                      window_size_sub_blocks +
+                      1)) {
+  RTC_DCHECK(data_dumper);
+  RTC_DCHECK_EQ(0, render_buffer_.data.size() % kSubBlockSize);
+  RTC_DCHECK_LT(0, window_size_sub_blocks);
+}
+
+Correlator::~Correlator() = default;
+
+void Correlator::Update(rtc::ArrayView<const float> render,
+                        rtc::ArrayView<const float> capture) {
+  RTC_DCHECK_EQ(kSubBlockSize, capture.size());
+  RTC_DCHECK_EQ(render.size(), capture.size());
+  const float signal_strength_threshold =
+      correlations_[0].size() * 150.f * 150.f;
+
+  // Insert the new x sub block into x_buffer.

[aleloi, 2017/01/27 15:37:47]

Is x_sub and x_buffer references to other variable names? If so, please update.

[peah-webrtc, 2017/02/02 14:04:45]

x sub block refers to x subblock.

Changed the name to x_buffer_.

Done.

+  render_buffer_.index =
+      (render_buffer_.index + kSubBlockSize) % render_buffer_.data.size();

[aleloi, 2017/01/27 15:37:47]

Will the buffer always have room for kSubBlockSize samples more? Can you add a
DCHECK?

[peah-webrtc, 2017/02/02 14:04:45]

Done.

+  std::copy(render.begin(), render.end(),

[hlundin-webrtc, 2017/02/01 08:30:00]

RTC_DCHECK_LE(render_buffer_.index + render.size(), render_buffer_.data.size());

[peah-webrtc, 2017/02/02 14:04:45]

Done.

+            render_buffer_.data.begin() + render_buffer_.index);
+
+  // Update the correlations for all sub-correlators.
+  size_t alignment_shift = 0;
+  for (size_t sub_correlator = 0; sub_correlator < correlations_.size();
+       ++sub_correlator) {
+    float error = 0.f;
+    bool correlations_updated = false;
+    size_t render_0_index =

[hlundin-webrtc, 2017/02/01 08:30:00]

What is this? Can the name be shorter?

[peah-webrtc, 2017/02/02 14:04:45]

This is the index of the first value in the render buffer to use for the matched
filtering. I changed the name to better reflect that.

PTAL.

+        (render_buffer_.index - alignment_shift + render_buffer_.data.size()) %
+        render_buffer_.data.size();
+
+    // Update the correlation estimates.
+    RTC_DCHECK_GT(render_buffer_.data.size(),
+                  render_0_index + kSubBlockSize - 1);
+    for (size_t i = 0; i < kSubBlockSize; ++i, ++render_0_index) {
+      const size_t render_buffer_part1_size =

[hlundin-webrtc, 2017/02/01 08:30:00]

Please, comment on the fact that this is due to processing across the boundary
of a circular buffer.
Also, the name render_buffer_part1_size is very long and reduces readability.

[peah-webrtc, 2017/02/02 14:04:45]

Done.

+          std::min(correlations_[sub_correlator].size(), render_0_index + 1);
+
+      // Compute the correlation estimate update factor.
+      float s = 0.f;

[hlundin-webrtc, 2017/02/01 08:30:00]

For the sake of my own understanding, I'm going to match this with the notation
in
https://en.wikipedia.org/wiki/Least_mean_squares_filter#Normalised_least_mean....

render <=> x
correlations <=> h
s <=> h^T*x
capture <=> d
update_factor <=> e = d - h^T*x
normalized_update_factor <=> mu * e / x^T*x

[peah-webrtc, 2017/02/02 14:04:45]

Done.

+      float signal_strength = 0.f;
+      for (size_t lag = 0; lag < render_buffer_part1_size; ++lag) {
+        const float render_k = render_buffer_.data[render_0_index - lag];
+        signal_strength += render_k * render_k;
+        s += correlations_[sub_correlator][lag] * render_k;

[hlundin-webrtc, 2017/02/01 08:30:00]

Consider calculating s only when you know you'll need it, i.e., when
signal_strength is large enough.

[peah-webrtc, 2017/02/02 14:04:45]

I don't think that's possible. In order to compute the matched filter accuracy,
s is needed.

WDYT?

+      }
+
+      for (size_t lag = render_buffer_part1_size,
+                  render_buffer_index = render_buffer_.data.size() - 1;
+           lag < correlations_[sub_correlator].size();
+           ++lag, --render_buffer_index) {
+        const float render_k = render_buffer_.data[render_buffer_index];
+        signal_strength += render_k * render_k;
+        s += correlations_[sub_correlator][lag] * render_k;
+      }
+
+      const float update_factor =
+          std::min(32767.f, std::max(-32768.f, capture[i] - s));
+      error += update_factor * update_factor;
+
+      // Update the correlation.
+      if (signal_strength > signal_strength_threshold) {
+        correlations_updated = true;
+        const float normalized_update_factor =
+            0.7f * update_factor / signal_strength;
+
+        for (size_t lag = 0; lag < render_buffer_part1_size; ++lag) {
+          correlations_[sub_correlator][lag] +=
+              normalized_update_factor *
+              render_buffer_.data[render_0_index - lag];
+        }
+
+        for (size_t lag = render_buffer_part1_size,
+                    render_buffer_index = render_buffer_.data.size() - 1;
+             lag < correlations_[sub_correlator].size();
+             ++lag, --render_buffer_index) {
+          correlations_[sub_correlator][lag] +=
+              normalized_update_factor *
+              render_buffer_.data[render_buffer_index];
+        }
+      }
+    }
+
+    const float error_anchor =

[hlundin-webrtc, 2017/02/01 08:30:00]

// Energy of capture.

[peah-webrtc, 2017/02/02 14:04:45]

Yes, but I prefer to name it differently as it is used as an anchor for the
matched filter error. To me that name makes the computation of the LagEstimate
information more readable. WDYT?

+        std::accumulate(capture.begin(), capture.end(), 0.f,

[hlundin-webrtc, 2017/02/01 08:30:00]

This should be the same as: 
std::inner_product(capture.begin(), capture.end(), capture.begin(), 0.f);

[peah-webrtc, 2017/02/02 14:04:45]

Nice!

Done.

+                        [](float a, float b) -> float { return a + b * b; });
+    const size_t lag_estimate = std::distance(

[hlundin-webrtc, 2017/02/01 08:30:00]

// Distance from start of sub_correlator to strongest (squared) correlation
peak.

[peah-webrtc, 2017/02/02 14:04:45]

Done.

+        correlations_[sub_correlator].begin(),
+        std::max_element(
+            correlations_[sub_correlator].begin(),
+            correlations_[sub_correlator].end(),
+            [](float a, float b) -> bool { return a * a < b * b; }));
+

[hlundin-webrtc, 2017/02/01 08:30:00]

What can you say about lag_estimate. DCHECK any lower and upper bounds.

[peah-webrtc, 2017/02/02 14:04:45]

There are no bounds, except from 0 and filters_[n].size(). Since the max_element
method ensures that the lag_estimate is within these bounds, no DCHECKS are
required.

WDYT?

+    lag_estimates_[sub_correlator].Set(
+        error_anchor - error, error < 0.3f * error_anchor,

[hlundin-webrtc, 2017/02/01 08:30:00]

Explain the thresholding with 0.3.

[peah-webrtc, 2017/02/02 14:04:45]

It is a tuning parameter. I made a constant and adjusted the naming to describe
what it does. I'm not sure that we should describe it more than that in the
text.

WDYT?

+        lag_estimate + alignment_shift, correlations_updated);
+
+    // TODO(peah): Remove once development of EchoCanceller3 is fully done.
+    RTC_DCHECK_EQ(4, correlations_.size());
+    switch (sub_correlator) {
+      case 0:
+        data_dumper_->DumpRaw("aec3_correlator_0_h", correlations_[0]);
+        break;
+      case 1:
+        data_dumper_->DumpRaw("aec3_correlator_1_h", correlations_[1]);
+        break;
+      case 2:
+        data_dumper_->DumpRaw("aec3_correlator_2_h", correlations_[2]);
+        break;
+      case 3:
+        data_dumper_->DumpRaw("aec3_correlator_3_h", correlations_[3]);
+        break;
+      default:
+        RTC_DCHECK(false);
+    }
+
+    alignment_shift += correlator_intra_lag_shift_;
+  }
+}
+
+}  // namespace webrtc