Add new variance update option and unittests for intelligibility

webrtc/modules/audio_processing/intelligibility/intelligibility_utils.cc

 /*
  *  Copyright (c) 2014 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.
  */
 
@@ skipping 78 matching lines @@
       running_mean_sq_(new complex<float>[freqs]()),
       sub_running_mean_(new complex<float>[freqs]()),
       sub_running_mean_sq_(new complex<float>[freqs]()),
       variance_(new float[freqs]()),
       conj_sum_(new float[freqs]()),
       freqs_(freqs),
       window_size_(window_size),
       decay_(decay),
       history_cursor_(0),
       count_(0),
-      array_mean_(0.0f) {
+      array_mean_(0.0f),
+      buffer_full_(false) {
   history_.reset(new rtc::scoped_ptr<complex<float>[]>[freqs_]());
   for (int i = 0; i < freqs_; ++i) {
     history_[i].reset(new complex<float>[window_size_]());
   }
   subhistory_.reset(new rtc::scoped_ptr<complex<float>[]>[freqs_]());
   for (int i = 0; i < freqs_; ++i) {
     subhistory_[i].reset(new complex<float>[window_size_]());
   }
   subhistory_sq_.reset(new rtc::scoped_ptr<complex<float>[]>[freqs_]());
   for (int i = 0; i < freqs_; ++i) {
     subhistory_sq_[i].reset(new complex<float>[window_size_]());
   }
   switch (type) {
     case kStepInfinite:
       step_func_ = &VarianceArray::InfiniteStep;
       break;
     case kStepDecaying:
       step_func_ = &VarianceArray::DecayStep;
       break;
     case kStepWindowed:
       step_func_ = &VarianceArray::WindowedStep;
       break;
     case kStepBlocked:
       step_func_ = &VarianceArray::BlockedStep;
       break;
+    case kStepBlockBasedMovingAverage:
+      step_func_ = &VarianceArray::BlockBasedMovingAverage;
+      break;
   }
 }
 
 // Compute the variance with Welford's algorithm, adding some fudge to
 // the input in case of all-zeroes.
 void VarianceArray::InfiniteStep(const complex<float>* data, bool skip_fudge) {
   array_mean_ = 0.0f;
   ++count_;
   for (int i = 0; i < freqs_; ++i) {
     complex<float> sample = data[i];
@@ skipping 81 matching lines @@
   history_cursor_ = (history_cursor_ + 1) % window_size_;
   ++count_;
 }
 
 // Variance with a window of blocks. Within each block, the variances are
 // recomputed from scratch at every stp, using |Var(X) = E(X^2) - E^2(X)|.
 // Once a block is filled with kWindowBlockSize samples, it is added to the
 // history window and a new block is started. The variances for the window
 // are recomputed from scratch at each of these transitions.
 void VarianceArray::BlockedStep(const complex<float>* data, bool /*dummy*/) {
-  int blocks = min(window_size_, history_cursor_);
+  int blocks = min(window_size_, history_cursor_ + 1);
   for (int i = 0; i < freqs_; ++i) {
     AddToMean(data[i], count_ + 1, &sub_running_mean_[i]);
     AddToMean(data[i] * std::conj(data[i]), count_ + 1,
               &sub_running_mean_sq_[i]);
     subhistory_[i][history_cursor_ % window_size_] = sub_running_mean_[i];
     subhistory_sq_[i][history_cursor_ % window_size_] = sub_running_mean_sq_[i];
 
     variance_[i] =
         (NewMean(running_mean_sq_[i], sub_running_mean_sq_[i], blocks) -
          NewMean(running_mean_[i], sub_running_mean_[i], blocks) *
              std::conj(NewMean(running_mean_[i], sub_running_mean_[i], blocks)))
             .real();
     if (count_ == kWindowBlockSize - 1) {
       sub_running_mean_[i] = complex<float>(0.0f, 0.0f);
       sub_running_mean_sq_[i] = complex<float>(0.0f, 0.0f);
       running_mean_[i] = complex<float>(0.0f, 0.0f);
       running_mean_sq_[i] = complex<float>(0.0f, 0.0f);
       for (int j = 0; j < min(window_size_, history_cursor_); ++j) {
-        AddToMean(subhistory_[i][j], j, &running_mean_[i]);
-        AddToMean(subhistory_sq_[i][j], j, &running_mean_sq_[i]);
+        AddToMean(subhistory_[i][j], j + 1, &running_mean_[i]);
+        AddToMean(subhistory_sq_[i][j], j + 1, &running_mean_sq_[i]);
       }
       ++history_cursor_;
     }
   }
   ++count_;
   if (count_ == kWindowBlockSize) {
     count_ = 0;
   }
 }
 
+// Recomputes variances from scratch each window based on previous window.
+void VarianceArray::BlockBasedMovingAverage(

[turaj, 2015/06/26 00:32:58]

There is a concern, which is not proven, that keeping track of sum of the
elements in a circular buffer in the following fashion, i.e. subtract the oldest
and add the newest, might diverge from the true sum, over time. Although it
might not be the case for this application as a call typically is not so long. 

Could you please add a TODO that we add a counter so that after a while
|running_mean_| and |running_mean_sq_| should be replaces with their "true"
value by summing over all the element?

[ekm, 2015/06/26 19:07:09]

Interesting. This is just do to floating point errors? and we can't expect the
total error to converge to 0?

[turaj, 2015/06/29 17:33:35]

Something like that, but I'm not sure.

+    const std::complex<float>* data, bool /*dummy*/) {
+  for (int i = 0; i < freqs_; ++i) {
+    sub_running_mean_[i] += data[i];
+    sub_running_mean_sq_[i] += data[i] * std::conj(data[i]);
+  }
+  ++count_;
+
+  // TODO(ekmeyerson) make kWindowBlockSize nonconstant to allow
+  // experimentation with different block size,window size pairs.
+  if (count_ >= kWindowBlockSize) {
+    count_ = 0;
+
+    for (int i = 0; i < freqs_; ++i) {
+      running_mean_[i] -= subhistory_[i][history_cursor_];
+      running_mean_sq_[i] -= subhistory_sq_[i][history_cursor_];
+
+      float scale = 1.f / kWindowBlockSize;
+      subhistory_[i][history_cursor_] = sub_running_mean_[i] * scale;
+      subhistory_sq_[i][history_cursor_] = sub_running_mean_sq_[i] * scale;
+
+      sub_running_mean_[i] = std::complex<float>(0.0f, 0.0f);
+      sub_running_mean_sq_[i] = std::complex<float>(0.0f, 0.0f);
+
+      running_mean_[i] += subhistory_[i][history_cursor_];
+      running_mean_sq_[i] += subhistory_sq_[i][history_cursor_];
+
+      scale = 1.f / (buffer_full_ ? window_size_ : history_cursor_ + 1);
+      variance_[i] = std::real(running_mean_sq_[i] * scale - running_mean_[i] *
+                              scale * std::conj(running_mean_[i]) * scale);
+    }
+
+    ++history_cursor_;
+    if (history_cursor_ >= window_size_) {
+      buffer_full_ = true;
+      history_cursor_ = 0;
+    }
+  }
+}
+
 void VarianceArray::Clear() {
   memset(running_mean_.get(), 0, sizeof(*running_mean_.get()) * freqs_);
   memset(running_mean_sq_.get(), 0, sizeof(*running_mean_sq_.get()) * freqs_);
   memset(variance_.get(), 0, sizeof(*variance_.get()) * freqs_);
   memset(conj_sum_.get(), 0, sizeof(*conj_sum_.get()) * freqs_);
   history_cursor_ = 0;
   count_ = 0;
   array_mean_ = 0.0f;
 }
 
@@ skipping 24 matching lines @@
       factor = 1.0f;
     }
     out_block[i] = factor * in_block[i];
     current_[i] = UpdateFactor(target_[i], current_[i], change_limit_);
   }
 }
 
 }  // namespace intelligibility
 
 }  // namespace webrtc