Source/WebCore/platform/audio/IIRFilter.cpp

/*
 * Copyright 2016 The Chromium Authors. All rights reserved.
 * Copyright (C) 2020, Apple Inc. All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1.  Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2.  Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY APPLE INC. AND ITS CONTRIBUTORS ``AS IS'' AND ANY
 * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL APPLE INC. OR ITS CONTRIBUTORS BE LIABLE FOR ANY
 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 * ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

#include "config.h"

#if ENABLE(WEB_AUDIO)
#include "IIRFilter.h"

#include "AudioArray.h"
#include "AudioUtilities.h"
#include "VectorMath.h"
#include <algorithm>
#include <complex>
#include <wtf/MathExtras.h>

namespace WebCore {

// 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.
constexpr int bufferLength = 32;
static_assert(bufferLength >= IIRFilter::maxOrder + 1, "Internal IIR buffer length must be greater than maximum IIR Filter order.");

static std::complex<double> evaluatePolynomial(const double* coefficients, 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 + std::complex<double>(coefficients[k]);
    return result;
}

IIRFilter::IIRFilter(const Vector<double>& feedforward, const Vector<double>& feedback)
    : m_feedforward(feedforward)
    , m_feedback(feedback)
{
    m_xBuffer.fill(0, bufferLength);
    m_yBuffer.fill(0, bufferLength);
}

void IIRFilter::reset()
{
    m_xBuffer.fill(0);
    m_yBuffer.fill(0);
    m_bufferIndex = 0;
}

void IIRFilter::process(const float* source, float* destination, 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] * source[n];

        // Run both the feedforward and feedback terms together, when possible.
        for (int k = 1; k < minLength; ++k) {
            int m = (m_bufferIndex - k) & (bufferLength - 1);
            yn += feedforward[k] * xBuffer[m];
            yn -= feedback[k] * yBuffer[m];
        }

        // Handle any remaining feedforward or feedback terms.
        for (int k = minLength; k < feedforwardLength; ++k)
            yn += feedforward[k] * xBuffer[(m_bufferIndex - k) & (bufferLength - 1)];
        for (int k = minLength; k < feedbackLength; ++k)
            yn -= feedback[k] * yBuffer[(m_bufferIndex - k) & (bufferLength - 1)];

        // Save the current input and output values in the memory buffers for the
        // next output.
        m_xBuffer[m_bufferIndex] = source[n];
        m_yBuffer[m_bufferIndex] = yn;

        m_bufferIndex = (m_bufferIndex + 1) & (bufferLength - 1);

        destination[n] = yn;
    }
}

void IIRFilter::getFrequencyResponse(unsigned length, 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 (unsigned k = 0; k < length; ++k) {
        if (frequency[k] < 0 || frequency[k] > 1) {
            // Out-of-bounds frequencies should return NaN.
            magResponse[k] = std::nanf("");
            phaseResponse[k] = std::nanf("");
        } else {
            // zRecip = 1/z = exp(-j*frequency)
            double omega = -piDouble * frequency[k];
            auto zRecip = std::complex<double>(cos(omega), sin(omega));

            auto numerator = evaluatePolynomial(m_feedforward.data(), zRecip, m_feedforward.size() - 1);
            auto denominator = evaluatePolynomial(m_feedback.data(), zRecip, m_feedback.size() - 1);
            auto response = numerator / denominator;
            magResponse[k] = static_cast<float>(abs(response));
            phaseResponse[k] = static_cast<float>(atan2(imag(response), real(response)));
        }
    }
}

double IIRFilter::tailTime(double sampleRate, bool isFilterStable)
{
    // The maximum tail time. This is somewhat arbitrary, but we're assuming that
    // no one is going to expect the IIRFilter to produce an output after this
    // much time after the inputs have stopped.
    constexpr double maxTailTime = 10;

    // If the maximum amplitude of the impulse response is less than this, we
    // assume that we've reached the tail of the response. Currently, this means
    // that the impulse is less than 1 bit of a 16-bit PCM value.
    constexpr float maxTailAmplitude = 1 / 32768.0;

    // If filter is not stable, just return max tail. Since the filter is not
    // stable, the impulse response won't converge to zero, so we don't need to
    // find the impulse response to find the actual tail time.
    if (!isFilterStable || !sampleRate)
        return maxTailTime;

    // How to compute the tail time? We're going to filter an impulse
    // for |maxTailTime| seconds, in blocks of kRenderQuantumFrames at
    // a time. The maximum magnitude of this block is saved. After all
    // of the samples have been computed, find the last block with a
    // maximum magnitude greater than |kMaxTaileAmplitude|. That block
    // index + 1 will be the tail time. We don't need to be
    // super-accurate in computing the tail time since we process on
    // blocks, block accuracy is good enough, and the value just needs
    // to be larger than the "real" tail time, so we don't prematurely
    // zero out the output of the node.

    // Number of render quanta needed to reach the max tail time.
    int numberOfBlocks = std::ceil(sampleRate * maxTailTime / AudioUtilities::renderQuantumSize);
    RELEASE_ASSERT(numberOfBlocks);

    // Input and output buffers for filtering.
    AudioFloatArray input(AudioUtilities::renderQuantumSize);
    AudioFloatArray output(AudioUtilities::renderQuantumSize);

    // Array to hold the max magnitudes
    AudioFloatArray magnitudes(numberOfBlocks);

    // Create the impulse input signal.
    input[0] = 1;

    // Process the first block and get the max magnitude of the output.
    process(input.data(), output.data(), AudioUtilities::renderQuantumSize);
    magnitudes[0] = VectorMath::maximumMagnitude(output.data(), AudioUtilities::renderQuantumSize);

    // Process the rest of the signal, getting the max magnitude of the
    // output for each block.
    input[0] = 0;

    for (int k = 1; k < numberOfBlocks; ++k) {
        process(input.data(), output.data(), AudioUtilities::renderQuantumSize);
        magnitudes[k] = VectorMath::maximumMagnitude(output.data(), AudioUtilities::renderQuantumSize);
    }

    // Done computing the impulse response; reset the state so the actual node
    // starts in the expected initial state.
    reset();

    // Find the last block with amplitude greater than the threshold.
    int index = numberOfBlocks - 1;
    for (int k = index; k >= 0; --k) {
        if (magnitudes[k] > maxTailAmplitude) {
            index = k;
            break;
        }
    }

    // The magnitude first become lower than the threshold at the next block.
    // Compute the corresponding time value value; that's the tail time.
    return (index + 1) * AudioUtilities::renderQuantumSize / sampleRate;
}

} // namespace WebCore

#endif // ENABLE(WEB_AUDIO)