Add a sparse FIR filter implementation

A Finite Impulse Response filter implementation which takes advantage of sparse coefficients.
The coefficients are assumed to be uniformly distributed and have an initial offset.

BUG=webrtc:3146
R=andrew@webrtc.org, bjornv@webrtc.org, kwiberg@webrtc.org

Review URL: https://webrtc-codereview.appspot.com/49659004

Cr-Commit-Position: refs/heads/master@{#9002}

6 files changed, 347 insertions, 0 deletions

diff --git a/webrtc/common_audio/BUILD.gn b/webrtc/common_audio/BUILD.gn
index 34e1575b24..544a142ee3 100644
--- a/webrtc/common_audio/BUILD.gn
+++ b/webrtc/common_audio/BUILD.gn
@@ -85,6 +85,8 @@ source_set("common_audio") {
     "signal_processing/splitting_filter.c",
     "signal_processing/sqrt_of_one_minus_x_squared.c",
     "signal_processing/vector_scaling_operations.c",
+    "sparse_fir_filter.cc",
+    "sparse_fir_filter.h",
     "vad/include/vad.h",
     "vad/include/webrtc_vad.h",
     "vad/vad.cc",
diff --git a/webrtc/common_audio/common_audio.gyp b/webrtc/common_audio/common_audio.gyp
index c7a4870131..3742f2f3b3 100644
--- a/webrtc/common_audio/common_audio.gyp
+++ b/webrtc/common_audio/common_audio.gyp
@@ -99,6 +99,8 @@
         'signal_processing/splitting_filter.c',
         'signal_processing/sqrt_of_one_minus_x_squared.c',
         'signal_processing/vector_scaling_operations.c',
+        'sparse_fir_filter.cc',
+        'sparse_fir_filter.h',
         'vad/include/vad.h',
         'vad/include/webrtc_vad.h',
         'vad/vad.cc',
@@ -259,6 +261,7 @@
             'ring_buffer_unittest.cc',
             'signal_processing/real_fft_unittest.cc',
             'signal_processing/signal_processing_unittest.cc',
+            'sparse_fir_filter_unittest.cc',
             'vad/vad_core_unittest.cc',
             'vad/vad_filterbank_unittest.cc',
             'vad/vad_gmm_unittest.cc',
diff --git a/webrtc/common_audio/fir_filter_unittest.cc b/webrtc/common_audio/fir_filter_unittest.cc
index 1bd01bea81..13f79d9482 100644
--- a/webrtc/common_audio/fir_filter_unittest.cc
+++ b/webrtc/common_audio/fir_filter_unittest.cc
@@ -16,6 +16,7 @@
 #include "webrtc/base/scoped_ptr.h"
 
 namespace webrtc {
+namespace {
 
 static const float kCoefficients[] = {0.2f, 0.3f, 0.5f, 0.7f, 0.11f};
 static const size_t kCoefficientsLength = sizeof(kCoefficients) /
@@ -34,6 +35,8 @@ void VerifyOutput(const float* expected_output,
                       length * sizeof(expected_output[0])));
 }
 
+}  // namespace
+
 TEST(FIRFilterTest, FilterAsIdentity) {
   const float kCoefficients[] = {1.f, 0.f, 0.f, 0.f, 0.f};
   float output[kInputLength];
diff --git a/webrtc/common_audio/sparse_fir_filter.cc b/webrtc/common_audio/sparse_fir_filter.cc
new file mode 100644
index 0000000000..9b8468d08c
--- /dev/null
+++ b/webrtc/common_audio/sparse_fir_filter.cc
@@ -0,0 +1,60 @@
+/*
+ *  Copyright (c) 2015 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.
+ */
+
+#include "webrtc/common_audio/sparse_fir_filter.h"
+
+#include "webrtc/base/checks.h"
+
+namespace webrtc {
+
+SparseFIRFilter::SparseFIRFilter(const float* nonzero_coeffs,
+                                 size_t num_nonzero_coeffs,
+                                 size_t sparsity,
+                                 size_t offset)
+    : sparsity_(sparsity),
+      offset_(offset),
+      nonzero_coeffs_(nonzero_coeffs, nonzero_coeffs + num_nonzero_coeffs),
+      state_(sparsity_ * (num_nonzero_coeffs - 1) + offset_, 0.f) {
+  CHECK_GE(num_nonzero_coeffs, 1u);
+  CHECK_GE(sparsity, 1u);
+}
+
+void SparseFIRFilter::Filter(const float* in, size_t length, float* out) {
+  // Convolves the input signal |in| with the filter kernel |nonzero_coeffs_|
+  // taking into account the previous state.
+  for (size_t i = 0; i < length; ++i) {
+    out[i] = 0.f;
+    size_t j;
+    for (j = 0; i >= j * sparsity_ + offset_ &&
+                j < nonzero_coeffs_.size(); ++j) {
+      out[i] += in[i - j * sparsity_ - offset_] * nonzero_coeffs_[j];
+    }
+    for (; j < nonzero_coeffs_.size(); ++j) {
+      out[i] += state_[i + (nonzero_coeffs_.size() - j - 1) * sparsity_] *
+                nonzero_coeffs_[j];
+    }
+  }
+
+  // Update current state.
+  if (state_.size() > 0u) {
+    if (length >= state_.size()) {
+      std::memcpy(&state_.front(),
+                  &in[length - state_.size()],
+                  state_.size() * sizeof(*in));
+    } else {
+      std::memmove(&state_.front(),
+                   &state_[length],
+                   (state_.size() - length) * sizeof(state_[0]));
+      std::memcpy(&state_[state_.size() - length], in, length * sizeof(*in));
+    }
+  }
+}
+
+}  // namespace webrtc
diff --git a/webrtc/common_audio/sparse_fir_filter.h b/webrtc/common_audio/sparse_fir_filter.h
new file mode 100644
index 0000000000..f5f3fadcce
--- /dev/null
+++ b/webrtc/common_audio/sparse_fir_filter.h
@@ -0,0 +1,48 @@
+/*
+ *  Copyright (c) 2015 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.
+ */
+
+#ifndef WEBRTC_COMMON_AUDIO_SPARSE_FIR_FILTER_H_
+#define WEBRTC_COMMON_AUDIO_SPARSE_FIR_FILTER_H_
+
+#include <cstring>
+#include <vector>
+
+namespace webrtc {
+
+// A Finite Impulse Response filter implementation which takes advantage of a
+// sparse structure with uniformly distributed non-zero coefficients.
+class SparseFIRFilter {
+ public:
+  // |num_nonzero_coeffs| is the number of non-zero coefficients,
+  // |nonzero_coeffs|. They are assumed to be uniformly distributed every
+  // |sparsity| samples and with an initial |offset|. The rest of the filter
+  // coefficients will be assumed zeros. For example, with sparsity = 3, and
+  // offset = 1 the filter coefficients will be:
+  // B = [0 coeffs[0] 0 0 coeffs[1] 0 0 coeffs[2] ... ]
+  // All initial state values will be zeros.
+  SparseFIRFilter(const float* nonzero_coeffs,
+                  size_t num_nonzero_coeffs,
+                  size_t sparsity,
+                  size_t offset);
+
+  // Filters the |in| data supplied.
+  // |out| must be previously allocated and it must be at least of |length|.
+  void Filter(const float* in, size_t length, float* out);
+
+ private:
+  const size_t sparsity_;
+  const size_t offset_;
+  const std::vector<float> nonzero_coeffs_;
+  std::vector<float> state_;
+};
+
+}  // namespace webrtc
+
+#endif  // WEBRTC_COMMON_AUDIO_SPARSE_FIR_FILTER_H_
diff --git a/webrtc/common_audio/sparse_fir_filter_unittest.cc b/webrtc/common_audio/sparse_fir_filter_unittest.cc
new file mode 100644
index 0000000000..82a53a5287
--- /dev/null
+++ b/webrtc/common_audio/sparse_fir_filter_unittest.cc
@@ -0,0 +1,231 @@
+/*
+ *  Copyright (c) 2015 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.
+ */
+
+#include "webrtc/common_audio/sparse_fir_filter.h"
+
+#include "testing/gtest/include/gtest/gtest.h"
+#include "webrtc/base/arraysize.h"
+#include "webrtc/base/scoped_ptr.h"
+#include "webrtc/common_audio/fir_filter.h"
+
+namespace webrtc {
+namespace {
+
+static const float kCoeffs[] = {0.2f, 0.3f, 0.5f, 0.7f, 0.11f};
+static const float kInput[] =
+    {1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 7.f, 8.f, 9.f, 10.f};
+
+template <size_t N>
+void VerifyOutput(const float (&expected_output)[N], const float (&output)[N]) {
+  EXPECT_EQ(0, memcmp(expected_output, output, sizeof(output)));
+}
+
+}  // namespace
+
+TEST(SparseFIRFilterTest, FilterAsIdentity) {
+  const float kCoeff = 1.f;
+  const size_t kNumCoeff = 1;
+  const size_t kSparsity = 3;
+  const size_t kOffset = 0;
+  float output[arraysize(kInput)];
+  SparseFIRFilter filter(&kCoeff, kNumCoeff, kSparsity, kOffset);
+  filter.Filter(kInput, arraysize(kInput), output);
+  VerifyOutput(kInput, output);
+}
+
+TEST(SparseFIRFilterTest, SameOutputForScalarCoefficientAndDifferentSparsity) {
+  const float kCoeff = 2.f;
+  const size_t kNumCoeff = 1;
+  const size_t kLowSparsity = 1;
+  const size_t kHighSparsity = 7;
+  const size_t kOffset = 0;
+  float low_sparsity_output[arraysize(kInput)];
+  float high_sparsity_output[arraysize(kInput)];
+  SparseFIRFilter low_sparsity_filter(&kCoeff,
+                                      kNumCoeff,
+                                      kLowSparsity,
+                                      kOffset);
+  SparseFIRFilter high_sparsity_filter(&kCoeff,
+                                       kNumCoeff,
+                                       kHighSparsity,
+                                       kOffset);
+  low_sparsity_filter.Filter(kInput, arraysize(kInput), low_sparsity_output);
+  high_sparsity_filter.Filter(kInput, arraysize(kInput), high_sparsity_output);
+  VerifyOutput(low_sparsity_output, high_sparsity_output);
+}
+
+TEST(SparseFIRFilterTest, FilterUsedAsScalarMultiplication) {
+  const float kCoeff = 5.f;
+  const size_t kNumCoeff = 1;
+  const size_t kSparsity = 5;
+  const size_t kOffset = 0;
+  float output[arraysize(kInput)];
+  SparseFIRFilter filter(&kCoeff, kNumCoeff, kSparsity, kOffset);
+  filter.Filter(kInput, arraysize(kInput), output);
+  EXPECT_FLOAT_EQ(5.f, output[0]);
+  EXPECT_FLOAT_EQ(20.f, output[3]);
+  EXPECT_FLOAT_EQ(25.f, output[4]);
+  EXPECT_FLOAT_EQ(50.f, output[arraysize(kInput) - 1]);
+}
+
+TEST(SparseFIRFilterTest, FilterUsedAsInputShifting) {
+  const float kCoeff = 1.f;
+  const size_t kNumCoeff = 1;
+  const size_t kSparsity = 1;
+  const size_t kOffset = 4;
+  float output[arraysize(kInput)];
+  SparseFIRFilter filter(&kCoeff, kNumCoeff, kSparsity, kOffset);
+  filter.Filter(kInput, arraysize(kInput), output);
+  EXPECT_FLOAT_EQ(0.f, output[0]);
+  EXPECT_FLOAT_EQ(0.f, output[3]);
+  EXPECT_FLOAT_EQ(1.f, output[4]);
+  EXPECT_FLOAT_EQ(2.f, output[5]);
+  EXPECT_FLOAT_EQ(6.f, output[arraysize(kInput) - 1]);
+}
+
+TEST(SparseFIRFilterTest, FilterUsedAsArbitraryWeighting) {
+  const size_t kSparsity = 2;
+  const size_t kOffset = 1;
+  float output[arraysize(kInput)];
+  SparseFIRFilter filter(kCoeffs, arraysize(kCoeffs), kSparsity, kOffset);
+  filter.Filter(kInput, arraysize(kInput), output);
+  EXPECT_FLOAT_EQ(0.f, output[0]);
+  EXPECT_FLOAT_EQ(0.9f, output[3]);
+  EXPECT_FLOAT_EQ(1.4f, output[4]);
+  EXPECT_FLOAT_EQ(2.4f, output[5]);
+  EXPECT_FLOAT_EQ(8.61f, output[arraysize(kInput) - 1]);
+}
+
+TEST(SparseFIRFilterTest, FilterInLengthLesserOrEqualToCoefficientsLength) {
+  const size_t kSparsity = 1;
+  const size_t kOffset = 0;
+  float output[arraysize(kInput)];
+  SparseFIRFilter filter(kCoeffs, arraysize(kCoeffs), kSparsity, kOffset);
+  filter.Filter(kInput, 2, output);
+  EXPECT_FLOAT_EQ(0.2f, output[0]);
+  EXPECT_FLOAT_EQ(0.7f, output[1]);
+}
+
+TEST(SparseFIRFilterTest, MultipleFilterCalls) {
+  const size_t kSparsity = 1;
+  const size_t kOffset = 0;
+  float output[arraysize(kInput)];
+  SparseFIRFilter filter(kCoeffs, arraysize(kCoeffs), kSparsity, kOffset);
+  filter.Filter(kInput, 2, output);
+  EXPECT_FLOAT_EQ(0.2f, output[0]);
+  EXPECT_FLOAT_EQ(0.7f, output[1]);
+  filter.Filter(kInput, 2, output);
+  EXPECT_FLOAT_EQ(1.3f, output[0]);
+  EXPECT_FLOAT_EQ(2.4f, output[1]);
+  filter.Filter(kInput, 2, output);
+  EXPECT_FLOAT_EQ(2.81f, output[0]);
+  EXPECT_FLOAT_EQ(2.62f, output[1]);
+  filter.Filter(kInput, 2, output);
+  EXPECT_FLOAT_EQ(2.81f, output[0]);
+  EXPECT_FLOAT_EQ(2.62f, output[1]);
+  filter.Filter(&kInput[3], 3, output);
+  EXPECT_FLOAT_EQ(3.41f, output[0]);
+  EXPECT_FLOAT_EQ(4.12f, output[1]);
+  EXPECT_FLOAT_EQ(6.21f, output[2]);
+  filter.Filter(&kInput[3], 3, output);
+  EXPECT_FLOAT_EQ(8.12f, output[0]);
+  EXPECT_FLOAT_EQ(9.14f, output[1]);
+  EXPECT_FLOAT_EQ(9.45f, output[2]);
+}
+
+TEST(SparseFIRFilterTest, VerifySampleBasedVsBlockBasedFiltering) {
+  const size_t kSparsity = 3;
+  const size_t kOffset = 1;
+  float output_block_based[arraysize(kInput)];
+  SparseFIRFilter filter_block(kCoeffs,
+                               arraysize(kCoeffs),
+                               kSparsity,
+                               kOffset);
+  filter_block.Filter(kInput, arraysize(kInput), output_block_based);
+  float output_sample_based[arraysize(kInput)];
+  SparseFIRFilter filter_sample(kCoeffs,
+                                arraysize(kCoeffs),
+                                kSparsity,
+                                kOffset);
+  for (size_t i = 0; i < arraysize(kInput); ++i)
+    filter_sample.Filter(&kInput[i], 1, &output_sample_based[i]);
+  VerifyOutput(output_block_based, output_sample_based);
+}
+
+TEST(SparseFIRFilterTest, SimpleHighPassFilter) {
+  const size_t kSparsity = 2;
+  const size_t kOffset = 2;
+  const float kHPCoeffs[] = {1.f, -1.f};
+  const float kConstantInput[] =
+      {1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f};
+  float output[arraysize(kConstantInput)];
+  SparseFIRFilter filter(kHPCoeffs, arraysize(kHPCoeffs), kSparsity, kOffset);
+  filter.Filter(kConstantInput, arraysize(kConstantInput), output);
+  EXPECT_FLOAT_EQ(0.f, output[0]);
+  EXPECT_FLOAT_EQ(0.f, output[1]);
+  EXPECT_FLOAT_EQ(1.f, output[2]);
+  EXPECT_FLOAT_EQ(1.f, output[3]);
+  for (size_t i = kSparsity + kOffset; i < arraysize(kConstantInput); ++i)
+    EXPECT_FLOAT_EQ(0.f, output[i]);
+}
+
+TEST(SparseFIRFilterTest, SimpleLowPassFilter) {
+  const size_t kSparsity = 2;
+  const size_t kOffset = 2;
+  const float kLPCoeffs[] = {1.f, 1.f};
+  const float kHighFrequencyInput[] =
+      {1.f, 1.f, -1.f, -1.f, 1.f, 1.f, -1.f, -1.f, 1.f, 1.f};
+  float output[arraysize(kHighFrequencyInput)];
+  SparseFIRFilter filter(kLPCoeffs, arraysize(kLPCoeffs), kSparsity, kOffset);
+  filter.Filter(kHighFrequencyInput, arraysize(kHighFrequencyInput), output);
+  EXPECT_FLOAT_EQ(0.f, output[0]);
+  EXPECT_FLOAT_EQ(0.f, output[1]);
+  EXPECT_FLOAT_EQ(1.f, output[2]);
+  EXPECT_FLOAT_EQ(1.f, output[3]);
+  for (size_t i = kSparsity + kOffset; i < arraysize(kHighFrequencyInput); ++i)
+    EXPECT_FLOAT_EQ(0.f, output[i]);
+}
+
+TEST(SparseFIRFilterTest, SameOutputWhenSwappedCoefficientsAndInput) {
+  const size_t kSparsity = 1;
+  const size_t kOffset = 0;
+  float output[arraysize(kCoeffs)];
+  float output_swapped[arraysize(kCoeffs)];
+  SparseFIRFilter filter(kCoeffs, arraysize(kCoeffs), kSparsity, kOffset);
+  // Use arraysize(kCoeffs) for in_length to get same-length outputs.
+  filter.Filter(kInput, arraysize(kCoeffs), output);
+  SparseFIRFilter filter_swapped(kInput,
+                                 arraysize(kCoeffs),
+                                 kSparsity,
+                                 kOffset);
+  filter_swapped.Filter(kCoeffs, arraysize(kCoeffs), output_swapped);
+  VerifyOutput(output, output_swapped);
+}
+
+TEST(SparseFIRFilterTest, SameOutputAsFIRFilterWhenSparsityOneAndOffsetZero) {
+  const size_t kSparsity = 1;
+  const size_t kOffset = 0;
+  float output[arraysize(kInput)];
+  float sparse_output[arraysize(kInput)];
+  rtc::scoped_ptr<FIRFilter> filter(FIRFilter::Create(kCoeffs,
+                                                      arraysize(kCoeffs),
+                                                      arraysize(kInput)));
+  SparseFIRFilter sparse_filter(kCoeffs,
+                                arraysize(kCoeffs),
+                                kSparsity,
+                                kOffset);
+  filter->Filter(kInput, arraysize(kInput), output);
+  sparse_filter.Filter(kInput, arraysize(kInput), sparse_output);
+  for (size_t i = 0; i < arraysize(kInput); ++i) {
+    EXPECT_FLOAT_EQ(output[i], sparse_output[i]);
+  }
+}
+
+}  // namespace webrtc