Upgrade AudioRendererAlgorithm to use WSOLA,

media/filters/wsola_internals.cc

+// Copyright 2013 The Chromium Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style license that can be
+// found in the LICENSE file.
+
+#include "media/filters/wsola_internals.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 {
+
+namespace internal {
+
+bool InInterval(int n, Interval q) {
+  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) {
+  const float kEpsilon = 1e-12;
+  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] +
+                                                 kEpsilon);
+  }
+  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_EQ(a->channels(), b->channels());
+  DCHECK_GE(frame_offset_a, 0);
+  DCHECK_GE(frame_offset_b, 0);
+  DCHECK_LE(frame_offset_a + num_frames, a->frames());
+  DCHECK_LE(frame_offset_b + num_frames, b->frames());
+
+  memset(dot_product, 0, sizeof(*dot_product) * a->channels());
+  for (int k = 0; k < a->channels(); ++k) {
+    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++;
+    }
+  }
+}
+
+void MultiChannelMovingBlockEnergies(const AudioBus* input,
+                                     int frames_per_block,
+                                     float* energy) {
+  int num_blocks = input->frames() - (frames_per_block - 1);
+  int channels = input->channels();
+
+  for (int k = 0; k < input->channels(); ++k) {
+    const float* input_channel = input->channel(k);
+
+    energy[k] = 0;
+
+    // First block of channel |k|.
+    for (int m = 0; m < frames_per_block; ++m) {
+      energy[k] += input_channel[m] * input_channel[m];
+    }
+
+    const float* slide_out = input_channel;
+    const float* slide_in = &input_channel[frames_per_block];

[DaleCurtis, 2013/08/13 21:11:04]

just input_channel + frames_per_block ?

[turaj, 2013/08/16 22:13:56]

Done.

+    for (int n = 1; n < num_blocks; ++n, ++slide_in, ++slide_out) {
+      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 CubicInterpolation(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];
+
+  DCHECK_NE(a, 0);
+  *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_candidate_blocks) {
+  int channels = search_segment->channels();
+  int block_size = target_block->frames();
+  int num_candidate_blocks = 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_candidate_blocks[n * channels], channels);
+
+  // Set the starting point as optimal point.
+  float best_similarity = similarity[0];
+  int optimal_index = 0;
+
+  n += decimation;
+  if (n >= num_candidate_blocks) {
+    return 0;
+  }
+
+  MultiChannelDotProduct(target_block, 0, search_segment, n, block_size,
+                         dot_prod.get());
+  similarity[1] = MultiChannelSimilarityMeasure(
+      dot_prod.get(), energy_target_block,
+      &energy_candidate_blocks[n * channels], channels);
+
+  n += decimation;
+  if (n >= num_candidate_blocks) {
+    // We cannot do any more sampling. Compare these two values and return the
+    // optimal index.
+    return similarity[1] > similarity[0] ? decimation : 0;
+  }
+
+  for (; n < num_candidate_blocks; n += decimation) {
+    MultiChannelDotProduct(target_block, 0, search_segment, n, block_size,
+                           dot_prod.get());
+
+    similarity[2] = MultiChannelSimilarityMeasure(
+        dot_prod.get(), energy_target_block,
+        &energy_candidate_blocks[n * channels], channels);
+
+    if ((similarity[1] > similarity[0] && similarity[1] >= similarity[2]) ||
+        (similarity[1] >= similarity[0] && similarity[1] > similarity[2])) {
+      // A local maximum is found. Do a cubic interpolation for a better
+      // estimate of candidate maximum.
+      float normalized_candidate_index;
+      float candidate_similarity;
+      CubicInterpolation(similarity, &normalized_candidate_index,
+                         &candidate_similarity);
+
+      int candidate_index = n - decimation + static_cast<int>(

[DaleCurtis, 2013/08/13 21:11:04]

no need for cast+floor

[turaj, 2013/08/16 22:13:56]

Done.

[DaleCurtis, 2013/08/19 22:15:23]

Integer value, so just use static_cast instead of floor().

[turaj, 2013/08/21 01:01:19]

Done.

+          floor(normalized_candidate_index * decimation +  0.5f));
+      if (candidate_similarity > best_similarity &&
+          !InInterval(candidate_index, exclude_interval)) {
+        optimal_index = candidate_index;
+        best_similarity = candidate_similarity;
+      }
+    } else if (n + decimation >= num_candidate_blocks &&
+        similarity[2] > best_similarity && !InInterval(n, exclude_interval)) {

[DaleCurtis, 2013/08/13 21:11:04]

should align with (

[turaj, 2013/08/16 22:13:56]

Done.

+      // 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));

[DaleCurtis, 2013/08/13 21:11:04]

memcpy?

[turaj, 2013/08/16 22:13:56]

I thought overlapping in any direction is undefined with memcpy, do I
misinterpret the manual?

On 2013/08/13 21:11:04, DaleCurtis wrote:

+  }
+  return optimal_index;
+}
+
+int FullSearch(int lim_low,
+               int lim_high,
+               Interval exclude_interval,
+               const AudioBus* target_block,
+               const AudioBus* search_block,
+               const float* energy_target_block,
+               const float* energy_candidate_blocks) {
+  int channels = search_block->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_block, n, block_size,
+                           dot_prod.get());
+
+    float similarity = MultiChannelSimilarityMeasure(
+        dot_prod.get(), energy_target_block,
+        &energy_candidate_blocks[n * channels], channels);
+
+    if (similarity > best_similarity) {
+      best_similarity = similarity;
+      optimal_index = n;
+    }
+  }
+
+  return optimal_index;
+}
+
+int OptimalIndex(const AudioBus* search_block,
+                 const AudioBus* target_block,
+                 Interval exclude_interval) {
+  int channels = search_block->channels();
+  DCHECK(channels == target_block->channels());

[DaleCurtis, 2013/08/13 21:11:04]

DCHECK_EQ

[turaj, 2013/08/16 22:13:56]

Done.

+  int target_size = target_block->frames();
+  int num_candidate_blocks = search_block->frames() - (target_size - 1);
+  const int kSearchDecimation = 5;

[DaleCurtis, 2013/08/13 21:11:04]

Comment.

[turaj, 2013/08/16 22:13:56]

Done.

+
+  scoped_ptr<float[]> energy_target_block(new float[channels]);

[DaleCurtis, 2013/08/13 21:11:04]

Bummer about these allocations, but don't worry about them until we have some
performance numbers on this feature.

[turaj, 2013/08/16 22:13:56]

We can define local variables of a maximum expected size. It that is
significantly cheaper.

On 2013/08/13 21:11:04, DaleCurtis wrote:

[DaleCurtis, 2013/08/19 22:15:23]

We should save that for a follow up optimization pass if necessary.

+  scoped_ptr<float[]> energy_candidate_blocks(
+      new float[channels * num_candidate_blocks]);
+
+  // Energy of all candid frames.
+  MultiChannelMovingBlockEnergies(search_block, target_size,
+                                  energy_candidate_blocks.get());
+
+  // Energy of target frame.
+  MultiChannelDotProduct(target_block, 0, target_block, 0,
+                         target_size, energy_target_block.get());
+
+  int optimal_index = DecimatedSearch(kSearchDecimation,

[DaleCurtis, 2013/08/13 21:11:04]

Have you run the performance numbers w/ and w/o decimated search?  I.e. if you
benchmark the unittests, what kind of speedup does it actually provide?

[turaj, 2013/08/16 22:13:56]

No I haven't, but I can do that. Computationally, beside the allocation, this
should be where most CPU is spend.

On 2013/08/13 21:11:04, DaleCurtis wrote:

+                                      exclude_interval, target_block,
+                                      search_block, energy_target_block.get(),
+                                      energy_candidate_blocks.get());
+
+  int lim_low = std::max(0, optimal_index - kSearchDecimation);
+  int lim_high = std::min(num_candidate_blocks - 1,
+                          optimal_index + kSearchDecimation);
+  return FullSearch(lim_low, lim_high, exclude_interval, target_block,
+                    search_block, energy_target_block.get(),
+                    energy_candidate_blocks.get());
+}
+
+void GetSymmetricHanningWindow(int window_length, float* window) {
+  const float kPi = 3.14159265f;

[DaleCurtis, 2013/08/13 21:11:04]

use M_PI, see top of sinc_resampler_unittest.cc for important #include notes,
essentially this needs to be done:

// MSVC++ requires this to be set before any other includes to get M_PI.
#define _USE_MATH_DEFINES

#include <cmath>

[turaj, 2013/08/16 22:13:56]

Done.

+  const float scale = 2.f * kPi / static_cast<float>(window_length);

[DaleCurtis, 2013/08/13 21:11:04]

double

[turaj, 2013/08/16 22:13:56]

Done.

+  for (int n = 0; n < window_length; ++n)
+    window[n] = 0.5 * (1 - cos(n * scale));

[DaleCurtis, 2013/08/13 21:11:04]

1.0  since cos takes double.

[turaj, 2013/08/16 22:13:56]

Can I use cosf()? just because |window| is float.

On 2013/08/13 21:11:04, DaleCurtis wrote:

[DaleCurtis, 2013/08/19 22:15:23]

Yeah that's probably better than double everywhere.

[turaj, 2013/08/21 01:01:19]

Then does it make sense to have

const float kPi = 3.141592f;

On 2013/08/19 22:15:23, DaleCurtis wrote:

+}
+
+}  // namespace internal
+
+}  // namespace media
+