Upgrade AudioRendererAlgorithm to use WSOLA,

media/filters/audio_renderer_algorithm_util.cc

+// Copyright (c) 2012 The Chromium Authors. All rights reserved.

[DaleCurtis, 2013/07/16 00:18:40]

2013 and no (c). Here and other files.

[turaj, 2013/07/29 22:09:57]

Done.

+// Use of this source code is governed by a BSD-style license that can be
+// found in the LICENSE file.
+
+#include "media/filters/audio_renderer_algorithm_util.h"
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+
+#include "base/logging.h"
+#include "base/memory/scoped_ptr.h"
+#include "media/base/audio_bus.h"
+
+namespace media {
+
+bool InInterval(int n, interval q) {

[DaleCurtis, 2013/07/16 00:18:40]

Types should have the first letter capitalized; so Interval instead of interval.

[turaj, 2013/07/29 22:09:57]

Done.

+  return n >= q.first && n <= q.second;
+}
+
+float MultiChannelSimilarityMeasure(const float* dot_prod_a_b,
+                                    const float* energy_a,
+                                    const float* energy_b,
+                                    int channels) {
+  float similarity_measure = 0;
+  for (int n = 0; n < channels; ++n) {
+    similarity_measure += dot_prod_a_b[n] / sqrt(energy_a[n] * energy_b[n] +

[DaleCurtis, 2013/07/16 00:18:40]

What is 1e-12 for?

Additionally vector_math.h might be a good home for some of the vector math
primitives used by these functions.

[turaj, 2013/07/29 22:09:57]

It is to prevent dividing by zero. I change it to const. And thank for pointing
to vector_math.h, but so far there are only two functions one implementing "a =
a + b * scale" and "a = b * scale."
I'll use them if I need such functions.
  

On 2013/07/16 00:18:40, DaleCurtis wrote:

+                                                 1e-12);
+  }
+  return similarity_measure;
+}
+
+void MultiChannelDotProduct(const AudioBus* a,
+                            int frame_offset_a,
+                            const AudioBus* b,
+                            int frame_offset_b,
+                            int num_frames,
+                            float* dot_product) {
+  DCHECK(a->channels() == b->channels());

[DaleCurtis, 2013/07/16 00:18:40]

use DCHECK_EQ, DCHECK_GE, DCHECK_LE etc, instead.

[turaj, 2013/07/29 22:09:57]

Done.

+  DCHECK(frame_offset_a >= 0);
+  DCHECK(frame_offset_b >= 0);
+  DCHECK(frame_offset_a + num_frames <= a->frames());
+  DCHECK(frame_offset_b + num_frames <= b->frames());
+
+  memset(dot_product, 0, sizeof(*dot_product) * a->channels());
+  for (int k = 0; k < a->channels(); ++k) {

[ajm, 2013/07/23 18:03:28]

I suppose using i here and j in the inner loop is more natural.

+    const float* ch_a = a->channel(k) + frame_offset_a;
+    const float* ch_b = b->channel(k) + frame_offset_b;
+    for (int n = 0; n < num_frames; ++n) {
+      dot_product[k] += *ch_a++ * *ch_b++;

[ajm, 2013/07/23 18:03:28]

Any reason to prefer *ch_a++ over ch_a[n]?

+    }
+  }
+}
+
+void MultiChannelMovingWindowEnergies(const AudioBus* input,
+                                      int frames_per_window,
+                                      float* energy) {
+  int num_blocks = input->frames() - (frames_per_window - 1);
+  int channels = input->channels();
+
+  for (int k = 0; k < input->channels(); ++k) {

[ajm, 2013/07/23 18:03:28]

Again, I think i, j, k ... is more natural.

+    const float* input_channel = input->channel(k);
+
+    energy[k] = 0;
+    // First window of channel |k|.
+    for(int m = 0; m < frames_per_window; ++m)

[ajm, 2013/07/23 18:03:28]

for (

+      energy[k] += input_channel[m] * input_channel[m];
+
+    const float* slide_out = input_channel;
+    const float* slide_in = &input_channel[frames_per_window];
+    for (int n = 1; n < num_blocks; ++n, ++slide_in, ++slide_out) {

[ajm, 2013/07/23 18:03:28]

This looks fine, but some comments explaining it could be instructive.

+      energy[k + n * channels] = energy[k + (n - 1) * channels] - *slide_out *
+          *slide_out + *slide_in * *slide_in;
+    }
+  }
+}
+
+// Fit the curve f(x) = a * x^2 + b * x + c such that
+// f(-1) = |y[0]|
+// f(0) = |y[1]|
+// f(1) = |y[2]|.
+void CubicInterpol(const float* y_values,
+                   float* extremum,
+                   float* extremum_value) {
+  float a = 0.5f * (y_values[2] + y_values[0]) - y_values[1];
+  float b = 0.5f * (y_values[2] - y_values[0]);
+  float c = y_values[1];
+
+  *extremum = -b / (2.f * a);
+  *extremum_value = a * (*extremum) * (*extremum) + b * (*extremum) + c;
+}
+
+int DecimatedSearch(int decimation,
+                    interval exclude_interval,
+                    const AudioBus* target_block,
+                    const AudioBus* search_segment,
+                    const float* energy_target_block,
+                    const float* energy_candid_blocks) {
+  int channels = search_segment->channels();
+  int block_size = target_block->frames();
+  int num_candid_frames = search_segment->frames() - (block_size - 1);
+  scoped_ptr<float[]> dot_prod(new float[channels]);
+  float similarity[3];  // Three elements for cubic interpolation.
+
+  int n = 0;
+  MultiChannelDotProduct(target_block, 0, search_segment, n, block_size,
+                         dot_prod.get());
+  similarity[0] = MultiChannelSimilarityMeasure(
+      dot_prod.get(), energy_target_block, &energy_candid_blocks[n * channels],
+      channels);
+
+  // Set the starting point as optimal point.
+  float best_similarity = similarity[0];
+  int optimal_index = 0;
+
+  n += decimation;
+  MultiChannelDotProduct(target_block, 0, search_segment, n, block_size,
+                         dot_prod.get());
+  similarity[1] = MultiChannelSimilarityMeasure(
+      dot_prod.get(), energy_target_block, &energy_candid_blocks[n * channels],
+      channels);
+
+  n += decimation;
+  for (; n < num_candid_frames; n += decimation) {
+    MultiChannelDotProduct(target_block, 0, search_segment, n, block_size,
+                           dot_prod.get());
+
+    similarity[2] = MultiChannelSimilarityMeasure(
+        dot_prod.get(), energy_target_block,
+        &energy_candid_blocks[n * channels], channels);
+
+    if (similarity[1] > similarity[0] &&
+        similarity[1] > similarity[2]) {
+      // A local maximum is found. Do a cubic interpolation for a better
+      // estimate of candid maximum.
+      float normalized_candid_optimal_index;
+      float candid_best_similarity;
+      CubicInterpol(similarity, &normalized_candid_optimal_index,
+                    &candid_best_similarity);
+
+      int candid_optimal_index = n - decimation +
+          static_cast<int>(floor(normalized_candid_optimal_index * decimation +
+                                 0.5f));
+      if (candid_best_similarity > best_similarity &&
+          !InInterval(candid_optimal_index, exclude_interval)) {
+        optimal_index = candid_optimal_index;
+        best_similarity = candid_best_similarity;
+      }
+    } else if (n + decimation >= num_candid_frames &&
+        similarity[2] > best_similarity && !InInterval(n, exclude_interval)) {
+      // If this is the end-point and has a better similarity-measure than
+      // optimal, then we accept it as optimal point.
+      optimal_index = n;
+      best_similarity = similarity[2];
+    }
+    memmove(similarity, &similarity[1], 2 * sizeof(*similarity));
+  }
+  return optimal_index;
+}
+
+int PartialSearch(int lim_low,
+                  int lim_high,
+                  interval exclude_interval,
+                  const AudioBus* target_block,
+                  const AudioBus* search_segment,
+                  const float* energy_target_block,
+                  const float* energy_candid_blocks) {
+  int channels = search_segment->channels();
+  int block_size = target_block->frames();
+  scoped_ptr<float[]> dot_prod(new float[channels]);
+
+  float best_similarity = std::numeric_limits<float>::min();
+  int optimal_index = 0;
+
+  for (int n = lim_low; n <= lim_high; ++n) {
+    if (InInterval(n, exclude_interval)) {
+      continue;
+    }
+    MultiChannelDotProduct(target_block, 0, search_segment, n, block_size,
+                           dot_prod.get());
+
+    float similarity = MultiChannelSimilarityMeasure(
+        dot_prod.get(), energy_target_block,
+        &energy_candid_blocks[n * channels], channels);
+
+    if (similarity > best_similarity) {
+      best_similarity = similarity;
+      optimal_index = n;
+    }
+  }
+
+  return optimal_index;
+}
+
+int OptimalIndex(const AudioBus* search_segment,
+                 const AudioBus* target_block,
+                 interval exclude_interval) {
+  int channels = search_segment->channels();
+  DCHECK(channels == target_block->channels());
+  int block_size = target_block->frames();
+  int num_candid_frames = search_segment->frames() - (block_size - 1);
+  const int kSearchDecimation = 5;
+
+  scoped_ptr<float[]> energy_target_frame(new float[channels]);
+  scoped_ptr<float[]> energy_candid_frames(
+      new float[channels * num_candid_frames]);
+
+  // Energy of all candid frames.
+  MultiChannelMovingWindowEnergies(search_segment, block_size,
+                                   energy_candid_frames.get());
+
+  // Energy of target frame.
+  MultiChannelDotProduct(target_block, 0, target_block, 0,
+                         block_size, energy_target_frame.get());
+
+  int optimal_index = DecimatedSearch(kSearchDecimation,
+                                      exclude_interval, target_block,
+                                      search_segment, energy_target_frame.get(),
+                                      energy_candid_frames.get());
+
+  int lim_low = std::max(0, optimal_index - kSearchDecimation);
+  int lim_high = std::min(num_candid_frames - 1,
+                          optimal_index + kSearchDecimation);
+  return PartialSearch(lim_low, lim_high, exclude_interval, target_block,
+                       search_segment, energy_target_frame.get(),
+                       energy_candid_frames.get());
+}
+
+void HannSym(int window_length, float* window) {
+  const float kPi = 3.14159265f;
+  const float scale = 2.f * kPi / static_cast<float>(window_length);
+  for (int n = 0; n < window_length; ++n)
+    window[n] = 0.5 * (1 - cos(n * scale));
+}
+
+}  // namespace media
+
+
+
+