Implement IIRFilter node for WebAudio.

third_party/WebKit/Source/platform/audio/IIRFilter.cpp

Index: third_party/WebKit/Source/platform/audio/IIRFilter.cpp
diff --git a/third_party/WebKit/Source/platform/audio/IIRFilter.cpp b/third_party/WebKit/Source/platform/audio/IIRFilter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4ab4a19444c66f169369f48e9eba8d6dfa767a18
--- /dev/null
+++ b/third_party/WebKit/Source/platform/audio/IIRFilter.cpp
@@ -0,0 +1,134 @@
+// Copyright 2015 The Chromium Authors. All rights reserved.

[tkent, 2016/01/13 03:28:28]

2016

[Raymond Toy, 2016/01/13 18:30:54]

Done.

+// Use of this source code is governed by a BSD-style license that can be
+// found in the LICENSE file.
+
+#include "platform/audio/IIRFilter.h"
+
+#include "wtf/MathExtras.h"
+

[tkent, 2016/01/13 03:28:28]

A blink line is unnecessary.

[Raymond Toy, 2016/01/13 18:30:54]

Done.

+#include <complex>
+
+namespace blink {
+
+// The length of the memory buffers for the IIR filter.  This MUST be a power of two and must be
+// greater than the possible length of the filter coefficients.
+const int kBufferLength = 32;
+static_assert(kBufferLength >= IIRFilter::kMaxOrder + 1,
+    "Internal IIR buffer length must be greater than maximum IIR Filter order.");
+
+IIRFilter::IIRFilter(const AudioDoubleArray* feedforward, const AudioDoubleArray* feedback)
+    : m_bufferIndex(0)
+    , m_feedback(feedback)
+    , m_feedforward(feedforward)
+{
+    // These are guaranteed to be zero-initialized.
+    m_xBuffer.allocate(kBufferLength);
+    m_yBuffer.allocate(kBufferLength);
+}
+
+IIRFilter::~IIRFilter()
+{
+}
+
+void IIRFilter::reset()
+{
+    m_xBuffer.zero();
+    m_yBuffer.zero();
+}
+
+static std::complex<double> evaluatePolynomial(const double* coef, std::complex<double> z, int order)
+{
+    // Use Horner's method to evaluate the polynomial P(z) = sum(coef[k]*z^k, k, 0, order);
+    std::complex<double> result = 0;
+
+    for (int k = order; k >= 0; --k)
+        result = result * z + static_cast<std::complex<double>>(coef[k]);

[tkent, 2016/01/13 03:28:28]

Do we need static_cast<> here?
I think |result * z + std::complex<double>(ceof[k])| works.

[Raymond Toy, 2016/01/13 18:30:54]

Done.

+
+    return result;
+}
+
+void IIRFilter::process(const float* sourceP, float* destP, size_t framesToProcess)
+{
+    // Compute
+    //
+    //   y[n] = sum(b[k] * x[n - k], k = 0, M) - sum(a[k] * y[n - k], k = 1, N)
+    //
+    // where b[k] are the feedforward coefficients and a[k] are the feedback coefficients of the
+    // filter.
+
+    // This is a Direct Form I implementation of an IIR Filter.  Should we consider doing a
+    // different implementation such as Transposed Direct Form II?
+    const double* feedback = m_feedback->data();
+    const double* feedforward = m_feedforward->data();
+
+    ASSERT(feedback);
+    ASSERT(feedforward);
+
+    // Sanity check to see if the feedback coefficients have been scaled appropriately. It must
+    // be EXACTLY 1!
+    ASSERT(feedback[0] == 1);
+
+    int feedbackLength = m_feedback->size();
+    int feedforwardLength = m_feedforward->size();
+    int minLength = std::min(feedbackLength, feedforwardLength);
+
+    double* xBuffer = m_xBuffer.data();
+    double* yBuffer = m_yBuffer.data();
+
+    for (size_t n = 0; n < framesToProcess; ++n) {
+        // To help minimize roundoff, we compute using double's, even though the filter coefficients
+        // only have single precision values.
+        double yn = feedforward[0] * sourceP[n];
+
+        // Run both the feedforward and feedback terms together, when possible.
+        for (int k = 1; k < minLength; ++k) {
+            int n = (m_bufferIndex - k) & (kBufferLength - 1);
+            yn += feedforward[k] * xBuffer[n];
+            yn -= feedback[k] * yBuffer[n];
+        }
+
+        // Handle any remaining feedforward or feedback terms.
+        for (int k = minLength; k < feedforwardLength; ++k)
+            yn += feedforward[k] * xBuffer[(m_bufferIndex - k) & (kBufferLength - 1)];
+
+        for (int k = minLength; k < feedbackLength; ++k)
+            yn -= feedback[k] * yBuffer[(m_bufferIndex - k) & (kBufferLength - 1)];
+
+        // Save the current input and output values in the memory buffers for the next output.
+        m_xBuffer[m_bufferIndex] = sourceP[n];
+        m_yBuffer[m_bufferIndex] = yn;
+
+        m_bufferIndex = (m_bufferIndex + 1) & (kBufferLength - 1);
+
+        destP[n] = yn;
+    }
+}
+
+void IIRFilter::getFrequencyResponse(int nFrequencies, const float* frequency, float* magResponse, float* phaseResponse)
+{
+    // Evaluate the z-transform of the filter at the given normalized frequencies from 0 to 1. (One
+    // corresponds to the Nyquist frequency.)
+    //
+    // The z-tranform of the filter is
+    //
+    // H(z) = sum(b[k]*z^(-k), k, 0, M) / sum(a[k]*z^(-k), k, 0, N);
+    //
+    // The desired frequency response is H(exp(j*omega)), where omega is in [0, 1).
+    //
+    // Let P(x) = sum(c[k]*x^k, k, 0, P) be a polynomial of order P.  Then each of the sums in H(z)
+    // is equivalent to evaluating a polynomial at the point 1/z.
+
+    for (int k = 0; k < nFrequencies; ++k) {
+        // zRecip = 1/z = exp(-j*frequency)
+        double omega = -piDouble * frequency[k];
+        std::complex<double> zRecip = std::complex<double>(cos(omega), sin(omega));
+
+        std::complex<double> numerator = evaluatePolynomial(m_feedforward->data(), zRecip, m_feedforward->size() - 1);
+        std::complex<double> denominator = evaluatePolynomial(m_feedback->data(), zRecip, m_feedback->size() - 1);
+        std::complex<double> response = numerator / denominator;
+        magResponse[k] = static_cast<float>(abs(response));
+        phaseResponse[k] = static_cast<float>(atan2(imag(response), real(response)));
+    }
+}
+
+} // blink