Add algorithm for Residual Echo Detector.

webrtc/modules/audio_processing/echo_detector/echo_detector.cc

Index: webrtc/modules/audio_processing/echo_detector/echo_detector.cc
diff --git a/webrtc/modules/audio_processing/echo_detector/echo_detector.cc b/webrtc/modules/audio_processing/echo_detector/echo_detector.cc
index 3885c1fb0716389c9f6d6293761cefaa9b9db2fa..22d31ea7e9b2ff543364ed3eec06f25807a5fa55 100644
--- a/webrtc/modules/audio_processing/echo_detector/echo_detector.cc
+++ b/webrtc/modules/audio_processing/echo_detector/echo_detector.cc
@@ -8,21 +8,141 @@
  *  be found in the AUTHORS file in the root of the source tree.
  */
 
+#include "webrtc/base/checks.h"
 #include "webrtc/modules/audio_processing/echo_detector/echo_detector.h"

[hlundin-webrtc, 2016/10/14 08:00:39]

Aways include the corresponding .h file before any other includes.
https://google.github.io/styleguide/cppguide.html#Names_and_Order_of_Includes

[ivoc, 2016/10/17 16:05:08]

Done.

 
+#include <math.h>

[hlundin-webrtc, 2016/10/14 08:00:39]

System includes go before own project's includes (except echo_detector.h).
https://google.github.io/styleguide/cppguide.html#Names_and_Order_of_Includes

[ivoc, 2016/10/17 16:05:07]

Done.

+
+namespace {
+
+float power(const float* farend, size_t nr_samples) {

[hlundin-webrtc, 2016/10/14 08:00:39]

Power with capital P.

[hlundin-webrtc, 2016/10/14 08:00:39]

With my suggestion in the previous CL to use ArrayView, you can use that here as
well. The code then simplifies to
float Power(rtc::ArrayView<const float> signal) {
  double result = 0.;
  for (auto value : signal) {
    result += value * value;
  }
  return static_cast<float>(result);
}

If you want to go fancy, you can make use of ArrayView's iterators and invoke
std::inner_product. The whole function then collapses to a one-liner:

float Power(rtc::ArrayView<const float> x) {
  return static_cast<float>(rtc::inner_product(x.begin(), x.end(), x.begin(),
0.0));
}

(I am, however, not sure that the intermediate and final sum will be double. You
may want to verify that, if it is important.)

[peah-webrtc, 2016/10/14 08:18:57]

It is not superimportant, but if you instead use the rtc::ArrayView class, this
code could be a bit shorter

float Power(rtc::ArrayView<const float> x) {
  float power = 0.f;
  for (auto x_v : render) {
    power += x_v * x_v;
  }
  return power;
 }
}


....

const float capture_power = Power(rtc::ArrayView<const float>(nearend,
num_samples));

[peah-webrtc, 2016/10/14 08:18:58]

Capital P in the function name.

[peah-webrtc, 2016/10/14 08:18:58]

Since you are using this for both farend and nearend, the name of the input
vector should not be farend.

[kwiberg-webrtc, 2016/10/14 08:58:54]

The docs suggest that, unless you ask it nicely for something else,
std::inner_product will do float * float -> float multiplications and double +
double -> double additions, just like the manual computation you suggest. (Would
it make sense to do the multiplications in double too?)

[hlundin-webrtc, 2016/10/14 09:15:28]

But the current implementation does float * float -> double, then double +
double, and finally cast to float. Does that make sense from a precision
point-of-view, and, can you make std::inner_product do that?

[peah-webrtc, 2016/10/14 13:39:40]

In this case, nr_samples is limited to 160. I really doubt that it is worth
doing the processing in double format.

[kwiberg-webrtc, 2016/10/14 20:33:33]

I was going to say that the current implementation does float * float -> float,
then converts the float to a double and does double += double. But I decided to
look it up, and found out about FLT_EVAL_METHOD:
https://en.wikipedia.org/wiki/C99#IEEE.C2.A0754_floating_point_support.
Basically, the compiler is allowed but not required to use a larger type for
intermediate results, and you can look at the value of FLT_EVAL_METHOD to see
what it claims to be doing.

That Wikipedia page suggests that at least gcc on x86-64 does *not* use a larger
type for intermediate values in float and double expressions, so if you want it,
it's probably best to explicitly cast one of the operands to double.

(Separately from this, the compiler is also allowed to use infinite precision
for intermediates if it wants to. I gather this is mostly used with "fused
multiply and add" instructions, which looks like a good choice for the compiler
to employ here.)

And yes, you should be able to make std::inner_product do that. You may
optionally pass it the addition and multiplication functions you want it to use,
so all you have to do is give it std::plus<double> and std::multiplies<double>
if you want double precision for all intermediate values.

[kwiberg-webrtc, 2016/10/14 20:33:33]

OK, doing it all in float would simplify things. :-) And using
std::inner_product becomes more attractive since you don't have to pass it extra
options.

[hlundin-webrtc, 2016/10/14 21:15:46]

It seems to me that storing the intermediate as float or double doesn't matter.
I don't have strong opinions whether to use std::inner_product or simply stick
to what you have, but passing special operators to std::inner_product is not
attractive.

[ivoc, 2016/10/17 16:05:07]

Thanks for all the comments, this seems to be a popular topic for discussion :-)
I decided to use innter_product and do everything in floats.

[hlundin-webrtc, 2016/10/17 18:49:39]

Good call!

+  double result = 0.;

[peah-webrtc, 2016/10/14 08:18:58]

I think you can use float for the power computation.

[ivoc, 2016/10/17 16:05:07]

Done.

+  for (size_t i = 0; i < nr_samples; i++) {

[peah-webrtc, 2016/10/14 08:18:58]

++i is faster

[hlundin-webrtc, 2016/10/14 21:15:46]

Not for a simple scalar.
https://google.github.io/styleguide/cppguide.html#Preincrement_and_Predecrement
But it doesn't hurt to be consistent and always use pre-increment.

+    result += farend[i] * farend[i];
+  }
+  return static_cast<float>(result);

[peah-webrtc, 2016/10/14 08:18:58]

If you use float, you don't need the cast.

+}
+
+}  // namespace
+
 namespace webrtc {
 
-void EchoDetector::BufferFarend(const float* /*farend*/,
-                                size_t /*num_samples*/) {
-  // TODO(ivoc): Add implementation.
+void MeanVarianceEstimator::UpdateEstimate(float value) {
+  mean_ = (1.f - alpha_) * mean_ + alpha_ * value;
+  variance_ =
+      (1.f - alpha_) * variance_ + alpha_ * (value - mean_) * (value - mean_);
+}
+
+float MeanVarianceEstimator::SigmaEstimate() {
+  // Variance is >= 0, so this is okay.

[hlundin-webrtc, 2016/10/14 08:00:39]

DCHECK this assumption.

[ivoc, 2016/10/17 16:05:08]

Done.

+  return sqrtf(variance_);
+}
+
+float MeanVarianceEstimator::MeanEstimate() {
+  return mean_;
+}
+
+void MeanVarianceEstimator::Clear() {
+  mean_ = 0.f;
+  variance_ = 0.f;
+}
+
+void CovarianceEstimator::UpdateCovarianceEstimate(float x,

[peah-webrtc, 2016/10/14 08:18:58]

I don't think that covariance is fully correct here. What is computed is the
normalized covariance/cross-correlation. Would it be possibly to change the
naming to reflect what it is being computed?

[ivoc, 2016/10/17 16:05:07]

Done.

+                                                   float x_mean,
+                                                   float x_sigma,
+                                                   float y,
+                                                   float y_mean,
+                                                   float y_sigma) {
+  covariance_ =
+      (1.f - alpha_) * covariance_ + alpha_ * (x - x_mean) * (y - y_mean);

[peah-webrtc, 2016/10/14 08:18:58]

The style guide mentions that abbreviations should be avoided in naming. Would
it be possible to rename PCC? What about CrossCorrelation or normalize
cross-correlation?

[ivoc, 2016/10/17 16:05:07]

I chose normalized_cross_correlation.

+  PCC_ = covariance_ / (x_sigma * y_sigma + .0001f);
+}
+
+float CovarianceEstimator::PCCEstimate() {

[peah-webrtc, 2016/10/14 08:18:58]

Please change this method name to avoid the abbreviation.

[ivoc, 2016/10/17 16:05:07]

Done.

+  return PCC_;
+}
+
+void CovarianceEstimator::Clear() {
+  covariance_ = 0.f;
+  PCC_ = 0.f;
+}
+
+CircularBuffer::CircularBuffer(int size) : buffer_(size), size_(size) {}
+CircularBuffer::~CircularBuffer() {}

[hlundin-webrtc, 2016/10/14 08:00:39]

Replace {} with "= default".

[peah-webrtc, 2016/10/14 08:18:57]

Alternatively you can use = default; 

[ivoc, 2016/10/17 16:05:07]

Done.

+
+void CircularBuffer::Insert(float value) {
+  buffer_[next_insertion_index_] = value;

[peah-webrtc, 2016/10/14 08:18:57]

Please add a DCHECK that next_insertion_index_ < buffer_.size()

[ivoc, 2016/10/17 16:05:07]

Done.

+  next_insertion_index_++;

[peah-webrtc, 2016/10/14 08:18:57]

++next_insertion_index_ is faster.

[ivoc, 2016/10/17 16:05:07]

For scalars I don't think it is, but you're right that it's probably better to
use it everywhere.

+  next_insertion_index_ %= size_;

[peah-webrtc, 2016/10/14 08:18:58]

buffer_.size() here and elsewhere.

[ivoc, 2016/10/17 16:05:07]

Done.

+}
+
+float CircularBuffer::GetValue(int delay) {
+  RTC_DCHECK(delay >= 0 && delay < size_);

[hlundin-webrtc, 2016/10/14 08:00:39]

Separate into two DCHECKs to get better printouts. Also, use DCHECK_GE/LT.

[ivoc, 2016/10/17 16:05:07]

Since I changed delay to size_t, the >=0 check doesn't seem useful anymore. I
did switch to LT for the < size() check.

+  const int index = (size_ + next_insertion_index_ - delay - 1) % size_;

[hlundin-webrtc, 2016/10/14 08:00:39]

size_t

[ivoc, 2016/10/17 16:05:07]

Done.

+  return buffer_[index];

[peah-webrtc, 2016/10/14 08:18:58]

Please add a DCHECK that index < buffer_.size() and that index >= 0

[ivoc, 2016/10/17 16:05:07]

Done.

+}
+
+void CircularBuffer::Clear() {
+  std::fill(buffer_.begin(), buffer_.end(), 0.f);
+  next_insertion_index_ = 0;
+}
+
+EchoDetector::EchoDetector()
+    : render_buffer_(render_buffer_size_),

[peah-webrtc, 2016/10/14 08:18:58]

Do you need to have these in the initializer list? Can't you initialize them
inside at the declaration instead? That would shorten the code and also place
the initialization closer to the constant they are initialized to (which makes
it easier to see what they are initialized to).

[ivoc, 2016/10/17 16:05:07]

Good idea!

+      render_power_(lookback_periods_),
+      render_power_mean_(lookback_periods_),
+      render_power_sigma_(lookback_periods_) {}
+
+EchoDetector::~EchoDetector() {}

[hlundin-webrtc, 2016/10/14 08:00:39]

= default

[peah-webrtc, 2016/10/14 08:18:58]

Alternatively you can use = default; 

[ivoc, 2016/10/17 16:05:07]

Done.

+
+void EchoDetector::BufferFarend(const float* farend, size_t num_samples) {

[peah-webrtc, 2016/10/14 08:18:57]

After some discussions I think we will be changing farend to render and nearend
to capture for all other submodules. So probably it would make sense to do that
there as well.

[ivoc, 2016/10/17 16:05:07]

Done.

+  float pow = power(farend, num_samples);

[peah-webrtc, 2016/10/14 08:18:57]

If you change the name of the power method to Power, you don't need to
abbreviate the power with pow.

[ivoc, 2016/10/17 16:05:07]

Done.

+  render_buffer_.Insert(pow);
+  render_buffer_delay_ =
+      std::min(render_buffer_delay_ + 1, render_buffer_size_ - 1);

[peah-webrtc, 2016/10/14 08:18:57]

Why is this minimum needed? 

+}
+
+void EchoDetector::Process(const float* nearend, size_t num_samples) {
+  const float cap_pow = power(nearend, num_samples);

[hlundin-webrtc, 2016/10/14 08:00:39]

cap_* sounds like a value capped by something.

[peah-webrtc, 2016/10/14 08:18:57]

I think that along the style guide abbreviation rule it would be good not to
abbreviate capture as cap. Same thing with render and ren

[ivoc, 2016/10/17 16:05:08]

Done.

+  capture_variance_.UpdateEstimate(cap_pow);
+  const float cap_mean = capture_variance_.MeanEstimate();
+  const float cap_sigma = capture_variance_.SigmaEstimate();
+  float echo_likelihood = 0.f;
+
+  const float latest_ren_pow = render_buffer_.GetValue(render_buffer_delay_);

[peah-webrtc, 2016/10/14 08:18:58]

This delay is only used to access the most recently added value in the render
buffer, right? 
Since the render buffer knows that the delay and the latest value that was
added, it is a bit of double bookkeeping keeping track of that here as well. 
What about adding a getter function to the render_buffer that does this? Or am I
missing something?

[ivoc, 2016/10/17 16:05:07]

As discussed offline, I moved the buffer size into the circular buffer class,
and converted the other uses of the circular buffer to regular vectors.

+  render_buffer_delay_ = std::max(0, render_buffer_delay_ - 1);
+  render_variance_.UpdateEstimate(latest_ren_pow);
+  render_power_.Insert(latest_ren_pow);
+  render_power_mean_.Insert(render_variance_.MeanEstimate());
+  render_power_sigma_.Insert(render_variance_.SigmaEstimate());
+
+  for (int delay = 0; delay < lookback_periods_; delay++) {

[hlundin-webrtc, 2016/10/14 08:00:39]

size_t

[peah-webrtc, 2016/10/14 08:18:57]

size_t ?

[peah-webrtc, 2016/10/14 08:18:57]

size_t

[peah-webrtc, 2016/10/14 08:18:58]

++delay

[peah-webrtc, 2016/10/14 08:18:58]

This name is not fully clear to me. What is meant by periods? Do you mean
frames?

[peah-webrtc, 2016/10/14 08:18:58]

Use covariance_.size()

[ivoc, 2016/10/17 16:05:07]

Done.

+    const float ren_pow = render_power_.GetValue(delay);
+    const float ren_mean = render_power_mean_.GetValue(delay);
+    const float ren_sigma = render_power_sigma_.GetValue(delay);
+    covariance_[delay].UpdateCovarianceEstimate(cap_pow, cap_mean, cap_sigma,
+                                                ren_pow, ren_mean, ren_sigma);
+    echo_likelihood =

[hlundin-webrtc, 2016/10/14 08:00:39]

No need for the local variable echo_likelihood. You could just as well update
echo_likelihood_ here.

[ivoc, 2016/10/17 16:05:08]

Done.

+        std::max(echo_likelihood, covariance_[delay].PCCEstimate());
+  }
+  echo_likelihood_ = echo_likelihood;
 }
 
-void EchoDetector::Process(const float* /*nearend*/, size_t /*num_samples*/) {
-  // TODO(ivoc): Add implementation.
+void EchoDetector::Initialize(int sample_rate_hz) {

[hlundin-webrtc, 2016/10/14 08:00:39]

sample_rate_hz not used? Delete it.

[ivoc, 2016/10/17 16:05:08]

Done.

+  render_buffer_.Clear();
+  render_power_.Clear();
+  render_power_mean_.Clear();
+  render_power_sigma_.Clear();
+  render_variance_.Clear();
+  capture_variance_.Clear();
+  for (auto& cov : covariance_) {
+    cov.Clear();
+  }
+  echo_likelihood_ = 0.f;
 }
 
-void EchoDetector::Initialize(int /*sample_rate_hz*/) {
-  // TODO(ivoc): Add implementation.
+float EchoDetector::EchoLikelihood() {

[peah-webrtc, 2016/10/14 08:18:57]

const

[ivoc, 2016/10/17 16:05:08]

Done.

+  return echo_likelihood_;
 }
 
 }  // namespace webrtc