Add SincResampler ported from WebKit.

media/base/sinc_resampler.cc

Index: media/base/sinc_resampler.cc
diff --git a/media/base/sinc_resampler.cc b/media/base/sinc_resampler.cc
new file mode 100644
index 0000000000000000000000000000000000000000..a8337a62effef73c423163b2f6280052950bbbef
--- /dev/null
+++ b/media/base/sinc_resampler.cc
@@ -0,0 +1,193 @@
+// Copyright (c) 2012 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.
+//
+// Input buffer layout, dividing the total buffer into regions (r0 - r5):
+//
+// |----------------|-----------------------------------------|----------------|
+//
+//                                   kBlockSize + kKernelSize / 2
+//                   <--------------------------------------------------------->
+//                                               r0
+//
+//  kKernelSize / 2   kKernelSize / 2         kKernelSize / 2   kKernelSize / 2
+// <---------------> <--------------->       <---------------> <--------------->
+//         r1                r2                      r3                r4
+//
+//                                                     kBlockSize
+//                                     <--------------------------------------->

[Ami GONE FROM CHROMIUM, 2012/07/10 03:23:00]

I hope you have tests covering the boundary cases of these various spans
overlapping in different ways (i.e. testing that the impl is not vulnerable to
</<= type bugs checking endpoints, etc).

[DaleCurtis, 2012/07/10 05:24:48]

Did you have something particular in mind? These are fixed quantities (113-118)
at construction time.  I'll see about moving them into the constructor and if I
can come up with some DCHECKS which enforce their sizing.

I'll also see what kind of effect an off by one error has here.

[Ami GONE FROM CHROMIUM, 2012/07/10 16:48:50]

True, but they depend on the magic constants.  My point is that it's possible
for bugs to hide in the impl that will only be revealed when the magic constants
change, and for bugs to hide in the impl that aren't visible b/c the unittest
doesn't happen to hit all boundary conditions with the particular inputs
involved (i/o ratio, chunksize, requested frame counts).

[DaleCurtis, 2012/07/10 21:38:55]

I think the current unit test covers this adequately, an off by one error in any
of these results in several orders of magnitude of RMS error in the results:
i.e. off by one in r5 results in RMS error of 1.99998 vs 0.00022.

I've moved the initialization of these regions to the constructor, made them
const, and adding some basic DCHECKS for layout.

+//                                                         r5
+//
+// The algorithm:
+//
+// 1) Consume input frames into r0 (r1 is zero-initialized).
+// 2) Position kernel centered at start of r0 (r2) and generate output frames
+//    until kernel is centered at start of r4 or we've finished generating all
+//    the output frames.
+// 3) Copy r3 to r1 and r4 to r2.
+// 4) Consume input frames into r5 (zero-pad if we run out of input).
+// 5) Goto (2) until all of input is consumed.
+//
+// Note: we're glossing over how the sub-sample handling works with
+// |virtual_source_idx_|, etc.
+
+// MSVC++ requires this to be set before any other includes to get M_PI.
+#define _USE_MATH_DEFINES
+
+#include "media/base/sinc_resampler.h"
+
+#include <cmath>
+
+#include "base/logging.h"
+
+namespace media {
+
+SincResampler::SincResampler(double io_sample_rate_ratio, const ReadCB& read_cb)
+    : io_sample_rate_ratio_(io_sample_rate_ratio),
+      virtual_source_idx_(0),
+      buffer_primed_(false),
+      read_cb_(read_cb),
+      // TODO(dalecurtis): When we switch to AVX/SSE optimization, we'll need to
+      // allocate with 32-byte alignment and ensure they're sized % 32 bytes.
+      kernel_storage_(new float[kKernelStorageSize]),
+      input_buffer_(new float[kBufferSize]) {
+  DCHECK_EQ(kKernelSize % 2, 0) << "kKernelSize must be even!";
+  DCHECK_GT(kBlockSize, kKernelSize)
+      << "kBlockSize must be greater than kKernelSize!";
+
+  memset(kernel_storage_.get(), 0,
+         sizeof(*kernel_storage_.get()) * kKernelStorageSize);
+  memset(input_buffer_.get(), 0, sizeof(*input_buffer_.get()) * kBufferSize);
+
+  InitializeKernel();
+}
+
+SincResampler::~SincResampler() {}
+
+void SincResampler::InitializeKernel() {
+  // Blackman window parameters.
+  double alpha = 0.16;

[Ami GONE FROM CHROMIUM, 2012/07/10 03:23:00]

static const this and the next three?

[DaleCurtis, 2012/07/10 21:38:55]

Done.

+  double a0 = 0.5 * (1.0 - alpha);
+  double a1 = 0.5;
+  double a2 = 0.5 * alpha;
+
+  // |sinc_scale_factor| is basically the normalized cutoff frequency of the
+  // low-pass filter.
+  double sinc_scale_factor =
+      io_sample_rate_ratio_ > 1.0 ? 1.0 / io_sample_rate_ratio_ : 1.0;
+
+  // The sinc function is an idealized brick-wall filter, but since we're
+  // windowing it the transition from pass to stop does not happen right away.
+  // So we should adjust the lowpass filter cutoff slightly downward to avoid
+  // some aliasing at the very high-end.
+  // TODO(crogers): this value is empirical and to be more exact should vary
+  // depending on kKernelSize.
+  sinc_scale_factor *= 0.9;
+
+  // Generates a set of windowed sinc() kernels.
+  // We generate a range of sub-sample offsets from 0.0 to 1.0.
+  for (int offset_idx = 0; offset_idx <= kKernelOffsetCount; ++offset_idx) {
+    double subsample_offset =
+        static_cast<double>(offset_idx) / kKernelOffsetCount;
+
+    for (int i = 0; i < kKernelSize; ++i) {
+      // Compute the sinc with offset.
+      double s =
+          sinc_scale_factor * M_PI * (i - kKernelSize / 2 - subsample_offset);
+      double sinc = (!s ? 1.0 : sin(s) / s) * sinc_scale_factor;
+
+      // Compute Blackman window, matching the offset of the sinc().
+      double x = (i - subsample_offset) / kKernelSize;
+      double window = a0 - a1 * cos(2.0 * M_PI * x) + a2 * cos(4.0 * M_PI * x);
+
+      // Window the sinc() function and store at the correct offset.
+      kernel_storage_[i + offset_idx * kKernelSize] = sinc * window;
+    }
+  }
+}
+
+void SincResampler::Resample(float* destination, int frames) {
+  int remaining_frames = frames;
+
+  // Setup various region pointers in the buffer (see diagram above).
+  float* r0 = input_buffer_.get() + kKernelSize / 2;
+  float* r1 = input_buffer_.get();
+  float* r2 = r0;
+  float* r3 = r0 + kBlockSize - kKernelSize / 2;
+  float* r4 = r0 + kBlockSize;
+  float* r5 = r0 + kKernelSize / 2;
+
+  // Step (1) -- Prime the input buffer at the start of the input stream.
+  if (!buffer_primed_) {
+    read_cb_.Run(r0, kBlockSize + kKernelSize / 2);
+    buffer_primed_ = true;
+  }
+
+  // Step (2) -- Resample!
+  while (remaining_frames) {
+    while (virtual_source_idx_ < kBlockSize) {
+      // |virtual_source_idx_| lies in between two kernel offsets so figure out
+      // what they are.
+      int source_idx = static_cast<int>(virtual_source_idx_);
+      double subsample_remainder = virtual_source_idx_ - source_idx;
+
+      double virtual_offset_idx = subsample_remainder * kKernelOffsetCount;
+      int offset_idx = static_cast<int>(virtual_offset_idx);
+
+      float* k1 = kernel_storage_.get() + offset_idx * kKernelSize;
+      float* k2 = k1 + kKernelSize;
+
+      // Initialize input pointer based on quantized |virtual_source_idx_|.
+      float* input_ptr = r1 + source_idx;
+
+      // We'll compute "convolutions" for the two kernels which straddle
+      // |virtual_source_idx_|.
+      float sum1 = 0;
+      float sum2 = 0;
+
+      // Figure out how much to weight each kernel's "convolution".
+      double kernel_interpolation_factor = virtual_offset_idx - offset_idx;
+
+      // Generate a single output sample.
+      int n = kKernelSize;
+      float input;
+      // TODO(dalecurtis): For initial commit, I've ripped out all the SSE
+      // optimizations, these definitely need to go back in before release.
+      while (n--) {
+        input = *input_ptr++;
+        sum1 += input * *k1++;
+        sum2 += input * *k2++;
+      }
+
+      // Linearly interpolate the two "convolutions".
+      double result = (1.0 - kernel_interpolation_factor) * sum1
+          + kernel_interpolation_factor * sum2;
+
+      *destination++ = result;
+
+      // Advance the virtual index.
+      virtual_source_idx_ += io_sample_rate_ratio_;
+
+      if (!--remaining_frames)
+        return;

[Ami GONE FROM CHROMIUM, 2012/07/10 16:48:50]

What happens if virtual_source_idx_ just passed kBlockSize 3 lines up, and we
return here?  Does the next invocation of Resample() read out of bounds?
(this is the sort of thing I was talking about in terms of testing boundary
conditions)

[Chris Rogers, 2012/07/10 17:37:11]

I think this should be OK because the next time Resample() is called, the
while() on line 128 will fall through and we'll wrap it

[DaleCurtis, 2012/07/10 21:38:55]

Yes, as Chris noted the state will fall through to the next call and everything
works as normal.

+    }
+
+    // Wrap back around to the start.
+    virtual_source_idx_ -= kBlockSize;
+
+    // Step (3) Copy r3 to r1 and r4 to r2.
+    // This wraps the last input frames back to the start of the buffer.
+    memcpy(r1, r3, sizeof(*input_buffer_.get()) * (kKernelSize / 2));
+    memcpy(r2, r4, sizeof(*input_buffer_.get()) * (kKernelSize / 2));
+
+    // Step (4)
+    // Refresh the buffer with more input.
+    read_cb_.Run(r5, kBlockSize);
+  }
+}
+
+int SincResampler::ChunkSize() {
+  return kBlockSize / io_sample_rate_ratio_;
+}
+
+}  // namespace media