/*
 * Copyright (C) 2010 Google 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.
 * 3.  Neither the name of Apple Inc. ("Apple") nor the names of
 *     its contributors may be used to endorse or promote products derived
 *     from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY APPLE 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 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 "Biquad.h"

#include "AudioUtilities.h"
#include "DenormalDisabler.h"
#include <algorithm>
#include <stdio.h>
#include <wtf/MathExtras.h>

#if USE(ACCELERATE)
// Work around a bug where VForce.h forward declares std::complex in a way that's incompatible with libc++ complex.
#define __VFORCE_H
#include <Accelerate/Accelerate.h>
#endif

namespace WebCore {

#if USE(ACCELERATE)
constexpr int kBufferSize = 1024;
#endif

Biquad::Biquad()
{
#if USE(ACCELERATE)
    // Allocate two samples more for filter history
    m_inputBuffer.resize(kBufferSize + 2);
    m_outputBuffer.resize(kBufferSize + 2);
#endif

    // Allocate enough space for the a-rate filter coefficients to handle a
    // rendering quantum of 128 frames.
    m_b0.resize(AudioUtilities::renderQuantumSize);
    m_b1.resize(AudioUtilities::renderQuantumSize);
    m_b2.resize(AudioUtilities::renderQuantumSize);
    m_a1.resize(AudioUtilities::renderQuantumSize);
    m_a2.resize(AudioUtilities::renderQuantumSize);

    // Initialize as pass-thru (straight-wire, no filter effect)
    setNormalizedCoefficients(0, 1, 0, 0, 1, 0, 0);

    reset(); // clear filter memory
}

Biquad::~Biquad() = default;

void Biquad::process(const float* sourceP, float* destP, size_t framesToProcess)
{
    if (hasSampleAccurateValues()) {
        size_t n = framesToProcess;

        // Create local copies of member variables
        double x1 = m_x1;
        double x2 = m_x2;
        double y1 = m_y1;
        double y2 = m_y2;

        auto* b0 = m_b0.data();
        auto* b1 = m_b1.data();
        auto* b2 = m_b2.data();
        auto* a1 = m_a1.data();
        auto* a2 = m_a2.data();

        for (size_t k = 0; k < n; ++k) {
            // FIXME: this can be optimized by pipelining the multiply adds...
            float x = *sourceP++;
            float y = b0[k] * x + b1[k] * x1 + b2[k] * x2 - a1[k] * y1 - a2[k] * y2;

            *destP++ = y;

            // Update state variables
            x2 = x1;
            x1 = x;
            y2 = y1;
            y1 = y;
        }

        // Local variables back to member. Flush denormals here so we
        // don't slow down the inner loop above.
        m_x1 = DenormalDisabler::flushDenormalFloatToZero(x1);
        m_x2 = DenormalDisabler::flushDenormalFloatToZero(x2);
        m_y1 = DenormalDisabler::flushDenormalFloatToZero(y1);
        m_y2 = DenormalDisabler::flushDenormalFloatToZero(y2);

        // There is an assumption here that once we have sample accurate values we
        // can never go back to not having sample accurate values. This is
        // currently true in the way AudioParamTimline is implemented: once an
        // event is inserted, sample accurate processing is always enabled.
        //
        // If so, then we never have to update the state variables for the MACOSX
        // path. The structure of the state variable in these cases aren't well
        // documented so it's not clear how to update them anyway.
        return;
    }

#if USE(ACCELERATE)
    auto* inputP = m_inputBuffer.data();
    auto* outputP = m_outputBuffer.data();

    // Set up filter state. This is needed in case we're switching from
    // filtering with variable coefficients (i.e., with automations) to
    // fixed coefficients (without automations).
    inputP[0] = m_x2;
    inputP[1] = m_x1;
    outputP[0] = m_y2;
    outputP[1] = m_y1;

    // Use vecLib if available
    processFast(sourceP, destP, framesToProcess);

    ASSERT(framesToProcess >= 2);

    // Copy the last inputs and outputs to the filter memory variables.
    // This is needed because the next rendering quantum might be an
    // automation which needs the history to continue correctly. Because
    // sourceP and destP can be the same block of memory, we can't read from
    // sourceP to get the last inputs. Fortunately, processFast has put the
    // last inputs in input[0] and input[1].
    m_x1 = inputP[1];
    m_x2 = inputP[0];
    m_y1 = destP[framesToProcess - 1];
    m_y2 = destP[framesToProcess - 2];
#else
    size_t n = framesToProcess;

    // Create local copies of member variables
    double x1 = m_x1;
    double x2 = m_x2;
    double y1 = m_y1;
    double y2 = m_y2;

    double b0 = m_b0[0];
    double b1 = m_b1[0];
    double b2 = m_b2[0];
    double a1 = m_a1[0];
    double a2 = m_a2[0];

    while (n--) {
        // FIXME: this can be optimized by pipelining the multiply adds...
        float x = *sourceP++;
        float y = b0 * x + b1 * x1 + b2 * x2 - a1 * y1 - a2 * y2;

        *destP++ = y;

        // Update state variables
        x2 = x1;
        x1 = x;
        y2 = y1;
        y1 = y;
    }

    // Local variables back to member. Flush denormals here so we
    // don't slow down the inner loop above.
    m_x1 = DenormalDisabler::flushDenormalFloatToZero(x1);
    m_x2 = DenormalDisabler::flushDenormalFloatToZero(x2);
    m_y1 = DenormalDisabler::flushDenormalFloatToZero(y1);
    m_y2 = DenormalDisabler::flushDenormalFloatToZero(y2);
#endif
}

#if USE(ACCELERATE)

// Here we have optimized version using Accelerate.framework

void Biquad::processFast(const float* sourceP, float* destP, size_t framesToProcess)
{
    double filterCoefficients[5];
    filterCoefficients[0] = m_b0[0];
    filterCoefficients[1] = m_b1[0];
    filterCoefficients[2] = m_b2[0];
    filterCoefficients[3] = m_a1[0];
    filterCoefficients[4] = m_a2[0];

    double* inputP = m_inputBuffer.data();
    double* outputP = m_outputBuffer.data();

    double* input2P = inputP + 2;
    double* output2P = outputP + 2;

    // Break up processing into smaller slices (kBufferSize) if necessary.

    size_t n = framesToProcess;

    while (n > 0) {
        int framesThisTime = n < kBufferSize ? n : kBufferSize;

        // Copy input to input buffer
        for (int i = 0; i < framesThisTime; ++i)
            input2P[i] = *sourceP++;

        processSliceFast(inputP, outputP, filterCoefficients, framesThisTime);

        // Copy output buffer to output (converts float -> double).
        for (int i = 0; i < framesThisTime; ++i)
            *destP++ = static_cast<float>(output2P[i]);

        n -= framesThisTime;
    }
}

void Biquad::processSliceFast(double* sourceP, double* destP, double* coefficientsP, size_t framesToProcess)
{
    // Use double-precision for filter stability
    vDSP_deq22D(sourceP, 1, coefficientsP, destP, 1, framesToProcess);

    // Save history. Note that sourceP and destP reference m_inputBuffer and m_outputBuffer respectively.
    // These buffers are allocated (in the constructor) with space for two extra samples so it's OK to access
    // array values two beyond framesToProcess.
    sourceP[0] = sourceP[framesToProcess - 2 + 2];
    sourceP[1] = sourceP[framesToProcess - 1 + 2];
    destP[0] = destP[framesToProcess - 2 + 2];
    destP[1] = destP[framesToProcess - 1 + 2];
}

#endif // USE(ACCELERATE)


void Biquad::reset()
{
#if USE(ACCELERATE)
    // Two extra samples for filter history
    double* inputP = m_inputBuffer.data();
    inputP[0] = 0;
    inputP[1] = 0;

    double* outputP = m_outputBuffer.data();
    outputP[0] = 0;
    outputP[1] = 0;
#endif

    m_x1 = m_x2 = m_y1 = m_y2 = 0;

}

void Biquad::setLowpassParams(size_t index, double cutoff, double resonance)
{
    // Limit cutoff to 0 to 1.
    cutoff = std::max(0.0, std::min(cutoff, 1.0));
    
    if (cutoff == 1) {
        // When cutoff is 1, the z-transform is 1.
        setNormalizedCoefficients(index, 1, 0, 0, 1, 0, 0);
    } else if (cutoff > 0) {
        // Compute biquad coefficients for lowpass filter
        resonance = pow(10.0, 0.05 * resonance);

        double theta = piDouble * cutoff;
        double alpha = sin(theta) / (2 * resonance);
        double cosw = cos(theta);
        double beta = (1 - cosw) / 2;

        double b0 = beta;
        double b1 = 2 * beta;
        double b2 = beta;

        double a0 = 1 + alpha;
        double a1 = -2 * cosw;
        double a2 = 1 - alpha;

        setNormalizedCoefficients(index, b0, b1, b2, a0, a1, a2);
    } else {
        // When cutoff is zero, nothing gets through the filter, so set
        // coefficients up correctly.
        setNormalizedCoefficients(index, 0, 0, 0, 1, 0, 0);
    }
}

void Biquad::setHighpassParams(size_t index, double cutoff, double resonance)
{
    // Limit cutoff to 0 to 1.
    cutoff = std::max(0.0, std::min(cutoff, 1.0));

    if (cutoff == 1) {
        // The z-transform is 0.
        setNormalizedCoefficients(index, 0, 0, 0, 1, 0, 0);
    } else if (cutoff > 0) {
        // Compute biquad coefficients for highpass filter
        resonance = pow(10.0, 0.05 * resonance);

        double theta = piDouble * cutoff;
        double alpha = sin(theta) / (2 * resonance);
        double cosw = cos(theta);
        double beta = (1 + cosw) / 2;

        double b0 = beta;
        double b1 = -(2 * beta);
        double b2 = beta;

        double a0 = 1 + alpha;
        double a1 = -2 * cosw;
        double a2 = 1 - alpha;

        setNormalizedCoefficients(index, b0, b1, b2, a0, a1, a2);
    } else {
      // When cutoff is zero, we need to be careful because the above
      // gives a quadratic divided by the same quadratic, with poles
      // and zeros on the unit circle in the same place. When cutoff
      // is zero, the z-transform is 1.
        setNormalizedCoefficients(index, 1, 0, 0, 1, 0, 0);
    }
}

void Biquad::setNormalizedCoefficients(size_t index, double b0, double b1, double b2, double a0, double a1, double a2)
{
    double a0Inverse = 1 / a0;
    
    m_b0[index] = b0 * a0Inverse;
    m_b1[index] = b1 * a0Inverse;
    m_b2[index] = b2 * a0Inverse;
    m_a1[index] = a1 * a0Inverse;
    m_a2[index] = a2 * a0Inverse;
}

void Biquad::setLowShelfParams(size_t index, double frequency, double dbGain)
{
    // Clip frequencies to between 0 and 1, inclusive.
    frequency = std::max(0.0, std::min(frequency, 1.0));
    
    double A = pow(10.0, dbGain / 40);

    if (frequency == 1) {
        // The z-transform is a constant gain.
        setNormalizedCoefficients(index, A * A, 0, 0, 1, 0, 0);
    } else if (frequency > 0) {
        double w0 = piDouble * frequency;
        double S = 1; // filter slope (1 is max value)
        double alpha = 0.5 * sin(w0) * sqrt((A + 1 / A) * (1 / S - 1) + 2);
        double k = cos(w0);
        double k2 = 2 * sqrt(A) * alpha;
        double aPlusOne = A + 1;
        double aMinusOne = A - 1;

        double b0 = A * (aPlusOne - aMinusOne * k + k2);
        double b1 = 2 * A * (aMinusOne - aPlusOne * k);
        double b2 = A * (aPlusOne - aMinusOne * k - k2);
        double a0 = aPlusOne + aMinusOne * k + k2;
        double a1 = -2 * (aMinusOne + aPlusOne * k);
        double a2 = aPlusOne + aMinusOne * k - k2;

        setNormalizedCoefficients(index, b0, b1, b2, a0, a1, a2);
    } else {
        // When frequency is 0, the z-transform is 1.
        setNormalizedCoefficients(index, 1, 0, 0, 1, 0, 0);
    }
}

void Biquad::setHighShelfParams(size_t index, double frequency, double dbGain)
{
    // Clip frequencies to between 0 and 1, inclusive.
    frequency = std::max(0.0, std::min(frequency, 1.0));

    double A = pow(10.0, dbGain / 40);

    if (frequency == 1) {
        // The z-transform is 1.
        setNormalizedCoefficients(index, 1, 0, 0, 1, 0, 0);
    } else if (frequency > 0) {
        double w0 = piDouble * frequency;
        double S = 1; // filter slope (1 is max value)
        double alpha = 0.5 * sin(w0) * sqrt((A + 1 / A) * (1 / S - 1) + 2);
        double k = cos(w0);
        double k2 = 2 * sqrt(A) * alpha;
        double aPlusOne = A + 1;
        double aMinusOne = A - 1;

        double b0 = A * (aPlusOne + aMinusOne * k + k2);
        double b1 = -2 * A * (aMinusOne + aPlusOne * k);
        double b2 = A * (aPlusOne + aMinusOne * k - k2);
        double a0 = aPlusOne - aMinusOne * k + k2;
        double a1 = 2 * (aMinusOne - aPlusOne * k);
        double a2 = aPlusOne - aMinusOne * k - k2;

        setNormalizedCoefficients(index, b0, b1, b2, a0, a1, a2);
    } else {
        // When frequency = 0, the filter is just a gain, A^2.
        setNormalizedCoefficients(index, A * A, 0, 0, 1, 0, 0);
    }
}

void Biquad::setPeakingParams(size_t index, double frequency, double Q, double dbGain)
{
    // Clip frequencies to between 0 and 1, inclusive.
    frequency = std::max(0.0, std::min(frequency, 1.0));

    // Don't let Q go negative, which causes an unstable filter.
    Q = std::max(0.0, Q);

    double A = pow(10.0, dbGain / 40);

    if (frequency > 0 && frequency < 1) {
        if (Q > 0) {
            double w0 = piDouble * frequency;
            double alpha = sin(w0) / (2 * Q);
            double k = cos(w0);

            double b0 = 1 + alpha * A;
            double b1 = -2 * k;
            double b2 = 1 - alpha * A;
            double a0 = 1 + alpha / A;
            double a1 = -2 * k;
            double a2 = 1 - alpha / A;

            setNormalizedCoefficients(index, b0, b1, b2, a0, a1, a2);
        } else {
            // When Q = 0, the above formulas have problems. If we look at
            // the z-transform, we can see that the limit as Q->0 is A^2, so
            // set the filter that way.
            setNormalizedCoefficients(index, A * A, 0, 0, 1, 0, 0);
        }
    } else {
        // When frequency is 0 or 1, the z-transform is 1.
        setNormalizedCoefficients(index, 1, 0, 0, 1, 0, 0);
    }
}

void Biquad::setAllpassParams(size_t index, double frequency, double Q)
{
    // Clip frequencies to between 0 and 1, inclusive.
    frequency = std::max(0.0, std::min(frequency, 1.0));

    // Don't let Q go negative, which causes an unstable filter.
    Q = std::max(0.0, Q);

    if (frequency > 0 && frequency < 1) {
        if (Q > 0) {
            double w0 = piDouble * frequency;
            double alpha = sin(w0) / (2 * Q);
            double k = cos(w0);

            double b0 = 1 - alpha;
            double b1 = -2 * k;
            double b2 = 1 + alpha;
            double a0 = 1 + alpha;
            double a1 = -2 * k;
            double a2 = 1 - alpha;

            setNormalizedCoefficients(index, b0, b1, b2, a0, a1, a2);
        } else {
            // When Q = 0, the above formulas have problems. If we look at
            // the z-transform, we can see that the limit as Q->0 is -1, so
            // set the filter that way.
            setNormalizedCoefficients(index, -1, 0, 0, 1, 0, 0);
        }
    } else {
        // When frequency is 0 or 1, the z-transform is 1.
        setNormalizedCoefficients(index, 1, 0, 0, 1, 0, 0);
    }
}

void Biquad::setNotchParams(size_t index, double frequency, double Q)
{
    // Clip frequencies to between 0 and 1, inclusive.
    frequency = std::max(0.0, std::min(frequency, 1.0));

    // Don't let Q go negative, which causes an unstable filter.
    Q = std::max(0.0, Q);

    if (frequency > 0 && frequency < 1) {
        if (Q > 0) {
            double w0 = piDouble * frequency;
            double alpha = sin(w0) / (2 * Q);
            double k = cos(w0);

            double b0 = 1;
            double b1 = -2 * k;
            double b2 = 1;
            double a0 = 1 + alpha;
            double a1 = -2 * k;
            double a2 = 1 - alpha;

            setNormalizedCoefficients(index, b0, b1, b2, a0, a1, a2);
        } else {
            // When Q = 0, the above formulas have problems. If we look at
            // the z-transform, we can see that the limit as Q->0 is 0, so
            // set the filter that way.
            setNormalizedCoefficients(index, 0, 0, 0, 1, 0, 0);
        }
    } else {
        // When frequency is 0 or 1, the z-transform is 1.
        setNormalizedCoefficients(index, 1, 0, 0, 1, 0, 0);
    }
}

void Biquad::setBandpassParams(size_t index, double frequency, double Q)
{
    // No negative frequencies allowed.
    frequency = std::max(0.0, frequency);

    // Don't let Q go negative, which causes an unstable filter.
    Q = std::max(0.0, Q);

    if (frequency > 0 && frequency < 1) {
        double w0 = piDouble * frequency;
        if (Q > 0) {
            double alpha = sin(w0) / (2 * Q);
            double k = cos(w0);
    
            double b0 = alpha;
            double b1 = 0;
            double b2 = -alpha;
            double a0 = 1 + alpha;
            double a1 = -2 * k;
            double a2 = 1 - alpha;

            setNormalizedCoefficients(index, b0, b1, b2, a0, a1, a2);
        } else {
            // When Q = 0, the above formulas have problems. If we look at
            // the z-transform, we can see that the limit as Q->0 is 1, so
            // set the filter that way.
            setNormalizedCoefficients(index, 1, 0, 0, 1, 0, 0);
        }
    } else {
        // When the cutoff is zero, the z-transform approaches 0, if Q
        // > 0. When both Q and cutoff are zero, the z-transform is
        // pretty much undefined. What should we do in this case?
        // For now, just make the filter 0. When the cutoff is 1, the
        // z-transform also approaches 0.
        setNormalizedCoefficients(index, 0, 0, 0, 1, 0, 0);
    }
}

void Biquad::getFrequencyResponse(unsigned nFrequencies, const float* frequency, float* magResponse, float* phaseResponse)
{
    // Evaluate the Z-transform of the filter at given normalized
    // frequency from 0 to 1. (1 corresponds to the Nyquist
    // frequency.)
    //
    // The z-transform of the filter is
    //
    // H(z) = (b0 + b1*z^(-1) + b2*z^(-2))/(1 + a1*z^(-1) + a2*z^(-2))
    //
    // Evaluate as
    //
    // b0 + (b1 + b2*z1)*z1
    // --------------------
    // 1 + (a1 + a2*z1)*z1
    //
    // with z1 = 1/z and z = exp(j*pi*frequency). Hence z1 = exp(-j*pi*frequency)

    // Make local copies of the coefficients as a micro-optimization.
    double b0 = m_b0[0];
    double b1 = m_b1[0];
    double b2 = m_b2[0];
    double a1 = m_a1[0];
    double a2 = m_a2[0];
    
    for (unsigned k = 0; k < nFrequencies; ++k) {
        if (frequency[k] < 0 || frequency[k] > 1) {
            // Out-of-bounds frequencies should return NaN.
            magResponse[k] = std::nanf("");
            phaseResponse[k] = std::nanf("");
        } else {
            double omega = -piDouble * frequency[k];
            std::complex<double> z = std::complex<double>(cos(omega), sin(omega));
            std::complex<double> numerator = b0 + (b1 + b2 * z) * z;
            std::complex<double> denominator = std::complex<double>(1, 0) + (a1 + a2 * z) * z;
            std::complex<double> response = numerator / denominator;
            magResponse[k] = static_cast<float>(abs(response));
            phaseResponse[k] = static_cast<float>(atan2(imag(response), real(response)));
        }
    }
}

static double repeatedRootResponse(double n, double c1, double c2, double r, double logEPS)
{
    // The response is h(n) = r^(n-2)*[c1*(n+1)*r^2+c2]. We're looking
    // for n such that |h(n)| = eps. Equivalently, we want a root
    // of the equation log(|h(n)|) - log(eps) = 0 or
    //
    //   (n-2)*log(r) + log(|c1*(n+1)*r^2+c2|) - log(eps)
    //
    // This helps with finding a nuemrical solution because this
    // approximately linearizes the response for large n.

    return (n - 2) * std::log(r) + std::log(std::fabs(c1 * (n + 1) * r * r + c2)) - logEPS;
}


// Regula Falsi root finder, Illinois variant
// (https://en.wikipedia.org/wiki/False_position_method#The_Illinois_algorithm).
//
// This finds a root of the repeated root response where the root is
// assumed to lie between |low| and |high|. The response is given by
// |c1|, |c2|, and |r| as determined by |RepeatedRootResponse|.
// |logEPS| is the log the the maximum allowed amplitude in the
// response.
static double rootFinder(double low, double high, double logEPS, double c1, double c2, double r)
{
    // Desired accuray of the root (in frames). This doesn't need to be
    // super-accurate, so half frame is good enough, and should be less
    // than 1 because the algorithm may prematurely terminate.
    constexpr double accuracyThreshold = 0.5;
    // Max number of iterations to do. If we haven't converged by now,
    // just return whatever we've found.
    constexpr int maxIterations = 10;

    int side = 0;
    double root = 0;
    double fLow = repeatedRootResponse(low, c1, c2, r, logEPS);
    double fHigh = repeatedRootResponse(high, c1, c2, r, logEPS);

    // The function values must be finite and have opposite signs!
    ASSERT(std::isfinite(fLow));
    ASSERT(std::isfinite(fHigh));
    ASSERT(fLow * fHigh <= 0);

    int iteration;
    for (iteration = 0; iteration < maxIterations; ++iteration) {
        root = (fLow * high - fHigh * low) / (fLow - fHigh);
        if (fabs(high - low) < accuracyThreshold * std::fabs(high + low))
            break;
        double fr = repeatedRootResponse(root, c1, c2, r, logEPS);

        ASSERT(std::isfinite(fr));

        if (fr * fHigh > 0) {
            // fr and fHigh have same sign. Copy root to fHigh
            high = root;
            fHigh = fr;
            side = -1;
        } else if (fLow * fr > 0) {
            // fr and fLow have same sign. Copy root to fLow
            low = root;
            fLow = fr;
            if (side == 1)
                fHigh /= 2;
            side = 1;
        } else {
            // fLow * fr looks like zero, so assume we've converged.
            break;
        }
    }

    // Want to know if the max number of iterations is ever exceeded so
    // we can understand why that happened.
    ASSERT(iteration < maxIterations);

    return root;
}

double Biquad::tailFrame(size_t coefIndex, double maxFrame)
{
    // The Biquad filter is given by
    //
    //   H(z) = (b0 + b1/z + b2/z^2)/(1 + a1/z + a2/z^2).
    //
    // To compute the tail time, compute the impulse response, h(n), of
    // H(z), which we can do analytically. From this impulse response,
    // find the value n0 where |h(n)| <= eps for n >= n0.
    //
    // Assume first that the two poles of H(z) are not repeated, say r1
    // and r2. Then, we can compute a partial fraction expansion of
    // H(z):
    //
    //   H(z) = (b0+b1/z+b2/z^2)/[(1-r1/z)*(1-r2/z)]
    //        = b0 + C2/(z-r2) - C1/(z-r1)
    //
    //  where
    //    C2 = (b0*r2^2+b1*r2+b2)/(r2-r1)
    //    C1 = (b0*r1^2+b1*r1+b2)/(r2-r1)
    //
    // Expand H(z) then this in powers of 1/z gives:
    //
    //   H(z) = b0 -(C2/r2+C1/r1) + sum(C2*r2^(i-1)/z^i + C1*r1^(i-1)/z^i)
    //
    // Thus, for n > 1 (we don't care about small n),
    //
    //   h(n) = C2*r2^(n-1) + C1*r1^(n-1)
    //
    // We need to find n0 such that |h(n)| < eps for n > n0.
    //
    // Case 1: r1 and r2 are real and distinct, with |r1|>=|r2|.
    //
    // Then
    //
    //   h(n) = C1*r1^(n-1)*(1 + C2/C1*(r2/r1)^(n-1))
    //
    // so
    //
    //   |h(n)| = |C1|*|r1|^(n-1)*|1+C2/C1*(r2/r1)^(n-1)|
    //          <= |C1|*|r1|^(n-1)*[1 + |C2/C1|*|r2/r1|^(n-1)]
    //          <= |C1|*|r1|^(n-1)*[1 + |C2/C1|]
    //
    // so
    //
    //   |h(n)| = |C1|*|r1|^(n-1)*|1+C2/C1*(r2/r1)^(n-1)|
    //          <= |C1|*|r1|^(n-1)*[1 + |C2/C1|*|r2/r1|^(n-1)]
    //          <= |C1|*|r1|^(n-1)*[1 + |C2/C1|]
    //
    // by using the triangle inequality and the fact that |r2|<=|r1|.
    // And we want |h(n)|<=eps which is true if
    //
    //   |C1|*|r1|^(n-1)*[1 + |C2/C1|] <= eps
    //
    // or
    //
    //   n >= 1 + log(eps/C)/log(|r1|)
    //
    // where C = |C1|*[1+|C2/C1|] = |C1| + |C2|.
    //
    // Case 2: r1 and r2 are complex
    //
    // Thne we can write r1=r*exp(i*p) and r2=r*exp(-i*p). So,
    //
    //   |h(n)| = |C2*r^(n-1)*exp(-i*p*(n-1)) + C1*r^(n-1)*exp(i*p*(n-1))|
    //          = |C1|*r^(n-1)*|1 + C2/C1*exp(-i*p*(n-1))/exp(i*n*(n-1))|
    //          <= |C1|*r^(n-1)*[1 + |C2/C1|]
    //
    // Again, this is easily solved to give
    //
    //   n >= 1 + log(eps/C)/log(r)
    //
    // where C = |C1|*[1+|C2/C1|] = |C1| + |C2|.
    //
    // Case 3: Repeated roots, r1=r2=r.
    //
    // In this case,
    //
    //   H(z) = (b0+b1/z+b2/z^2)/[(1-r/z)^2
    //
    // Expanding this in powers of 1/z gives:
    //
    //   H(z) = C1*sum((i+1)*r^i/z^i) - C2 * sum(r^(i-2)/z^i) + b2/r^2
    //        = b2/r^2 + sum([C1*(i+1)*r^i + C2*r^(i-2)]/z^i)
    // where
    //   C1 = (b0*r^2+b1*r+b2)/r^2
    //   C2 = b1*r+2*b2
    //
    // Thus, the impulse response is
    //
    //   h(n) = C1*(n+1)*r^n + C2*r^(n-2)
    //        = r^(n-2)*[C1*(n+1)*r^2+C2]
    //
    // So
    //
    //   |h(n)| = |r|^(n-2)*|C1*(n+1)*r^2+C2|
    //
    // To find n such that |h(n)| < eps, we need a numerical method in
    // general, so there's no real reason to simplify this or use other
    // approximations. Just solve |h(n)|=eps directly.
    //
    // Thus, for an set of filter coefficients, we can compute the tail
    // time.
    //

    // 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 double maxTailAmplitude = 1 / 32768.0;

    // Find the roots of 1+a1/z+a2/z^2 = 0. Or equivalently,
    // z^2+a1*z+a2 = 0. From the quadratic formula the roots are
    // (-a1+/-sqrt(a1^2-4*a2))/2.

    double a1 = m_a1[coefIndex];
    double a2 = m_a2[coefIndex];
    double b0 = m_b0[coefIndex];
    double b1 = m_b1[coefIndex];
    double b2 = m_b2[coefIndex];

    double discrim = a1 * a1 - 4 * a2;
    if (discrim > 0) {
        // Compute the real roots so that r1 has the largest magnitude.
        double rplus = (-a1 + std::sqrt(discrim)) / 2;
        double rminus = (-a1 - std::sqrt(discrim)) / 2;
        double r1 = std::fabs(rplus) >= std::fabs(rminus) ? rplus : rminus;
        // Use the fact that a2 = r1*r2
        double r2 = a2 / r1;

        double c1 = (b0 * r1 * r1 + b1 * r1 + b2) / (r2 - r1);
        double c2 = (b0 * r2 * r2 + b1 * r2 + b2) / (r2 - r1);

        ASSERT(std::isfinite(r1));
        ASSERT(std::isfinite(r2));
        ASSERT(std::isfinite(c1));
        ASSERT(std::isfinite(c2));

        // It's possible for maxTailAmplitude to be greater than c1 + c2.
        // This may produce a negative tail frame. Just clamp the tail
        // frame to 0.
        double tailFrame = clampTo(1 + std::log(maxTailAmplitude / (std::fabs(c1) + std::fabs(c2))) / std::log(std::fabs(r1)), 0);

        ASSERT(std::isfinite(tailFrame));
        return tailFrame;
    }

    if (discrim < 0) {
        // Two complex roots.
        // One root is -a1/2 + i*sqrt(-discrim)/2.
        double x = -a1 / 2;
        double y = std::sqrt(-discrim) / 2;
        std::complex<double> r1(x, y);
        std::complex<double> r2(x, -y);
        double r = std::hypot(x, y);

        ASSERT(std::isfinite(r));

        double tailFrame = 0;
        // It's possible for r to be 1. (LPF with Q very large can cause this.)
        if (r == 1)
            tailFrame = maxFrame;
        else {
            double c1 = std::abs((b0 * r1 * r1 + b1 * r1 + b2) / (r2 - r1));
            double c2 = std::abs((b0 * r2 * r2 + b1 * r2 + b2) / (r2 - r1));

            ASSERT(std::isfinite(c1));
            ASSERT(std::isfinite(c2));
            tailFrame = 1 + std::log(maxTailAmplitude / (c1 + c2)) / std::log(r);
            if (!c1 && !c2) {
                // If c1 = c2 = 0, then H(z) = b0. Hence, there's no tail
                // because this is just a wire from input to output.
                tailFrame = 0;
            } else {
                // Otherwise, check that the tail has finite length. Not
                // strictly necessary, but we want to know if this ever
                // happens.
                ASSERT(std::isfinite(tailFrame));
            }
        }
        return tailFrame;
    }

    // Repeated roots. This should be pretty rare because all the
    // coefficients need to be just the right values to get a
    // discriminant of exactly zero.
    double r = -a1 / 2;

    double tailFrame = 0;
    if (!r) {
        // Double pole at 0. This just delays the signal by 2 frames,
        // so set the tail frame to 2.
        tailFrame = 2;
    } else {
        double c1 = (b0 * r * r + b1 * r + b2) / (r * r);
        double c2 = b1 * r + 2 * b2;

        ASSERT(std::isfinite(c1));
        ASSERT(std::isfinite(c2));

        // It can happen that c1=c2=0. This basically means that H(z) =
        // constant, which is the limiting case for several of the
        // biquad filters.
        if (!c1 && !c2)
            tailFrame = 0;
        else {
            // The function c*(n+1)*r^n is not monotonic, but it's easy to
            // find the max point since the derivative is
            // c*r^n*(1+(n+1)*log(r)). This has a root at
            // -(1+log(r))/log(r). so we can start our search from that
            // point to maxFrames.

            double low = std::clamp(-(1 + std::log(r)) / std::log(r), 1.0, static_cast<double>(maxFrame - 1));
            double high = maxFrame;

            ASSERT(std::isfinite(low));
            ASSERT(std::isfinite(high));

            tailFrame = rootFinder(low, high, std::log(maxTailAmplitude), c1, c2, r);
        }
    }

    return tailFrame;
}

} // namespace WebCore

#endif // ENABLE(WEB_AUDIO)