#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)
#define __VFORCE_H
#include <Accelerate/Accelerate.h>
#endif
namespace WebCore {
#if USE(ACCELERATE)
constexpr int kBufferSize = 1024;
#endif
Biquad::Biquad()
{
#if USE(ACCELERATE)
m_inputBuffer.resize(kBufferSize + 2);
m_outputBuffer.resize(kBufferSize + 2);
#endif
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);
setNormalizedCoefficients(0, 1, 0, 0, 1, 0, 0);
reset(); }
Biquad::~Biquad() = default;
void Biquad::process(const float* sourceP, float* destP, size_t framesToProcess)
{
if (hasSampleAccurateValues()) {
size_t n = framesToProcess;
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) {
float x = *sourceP++;
float y = b0[k] * x + b1[k] * x1 + b2[k] * x2 - a1[k] * y1 - a2[k] * y2;
*destP++ = y;
x2 = x1;
x1 = x;
y2 = y1;
y1 = y;
}
m_x1 = DenormalDisabler::flushDenormalFloatToZero(x1);
m_x2 = DenormalDisabler::flushDenormalFloatToZero(x2);
m_y1 = DenormalDisabler::flushDenormalFloatToZero(y1);
m_y2 = DenormalDisabler::flushDenormalFloatToZero(y2);
return;
}
#if USE(ACCELERATE)
auto* inputP = m_inputBuffer.data();
auto* outputP = m_outputBuffer.data();
inputP[0] = m_x2;
inputP[1] = m_x1;
outputP[0] = m_y2;
outputP[1] = m_y1;
processFast(sourceP, destP, framesToProcess);
ASSERT(framesToProcess >= 2);
m_x1 = inputP[1];
m_x2 = inputP[0];
m_y1 = destP[framesToProcess - 1];
m_y2 = destP[framesToProcess - 2];
#else
size_t n = framesToProcess;
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--) {
float x = *sourceP++;
float y = b0 * x + b1 * x1 + b2 * x2 - a1 * y1 - a2 * y2;
*destP++ = y;
x2 = x1;
x1 = x;
y2 = y1;
y1 = y;
}
m_x1 = DenormalDisabler::flushDenormalFloatToZero(x1);
m_x2 = DenormalDisabler::flushDenormalFloatToZero(x2);
m_y1 = DenormalDisabler::flushDenormalFloatToZero(y1);
m_y2 = DenormalDisabler::flushDenormalFloatToZero(y2);
#endif
}
#if USE(ACCELERATE)
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;
size_t n = framesToProcess;
while (n > 0) {
int framesThisTime = n < kBufferSize ? n : kBufferSize;
for (int i = 0; i < framesThisTime; ++i)
input2P[i] = *sourceP++;
processSliceFast(inputP, outputP, filterCoefficients, framesThisTime);
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)
{
vDSP_deq22D(sourceP, 1, coefficientsP, destP, 1, 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)
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)
{
cutoff = std::max(0.0, std::min(cutoff, 1.0));
if (cutoff == 1) {
setNormalizedCoefficients(index, 1, 0, 0, 1, 0, 0);
} else if (cutoff > 0) {
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 {
setNormalizedCoefficients(index, 0, 0, 0, 1, 0, 0);
}
}
void Biquad::setHighpassParams(size_t index, double cutoff, double resonance)
{
cutoff = std::max(0.0, std::min(cutoff, 1.0));
if (cutoff == 1) {
setNormalizedCoefficients(index, 0, 0, 0, 1, 0, 0);
} else if (cutoff > 0) {
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 {
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)
{
frequency = std::max(0.0, std::min(frequency, 1.0));
double A = pow(10.0, dbGain / 40);
if (frequency == 1) {
setNormalizedCoefficients(index, A * A, 0, 0, 1, 0, 0);
} else if (frequency > 0) {
double w0 = piDouble * frequency;
double S = 1; 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 {
setNormalizedCoefficients(index, 1, 0, 0, 1, 0, 0);
}
}
void Biquad::setHighShelfParams(size_t index, double frequency, double dbGain)
{
frequency = std::max(0.0, std::min(frequency, 1.0));
double A = pow(10.0, dbGain / 40);
if (frequency == 1) {
setNormalizedCoefficients(index, 1, 0, 0, 1, 0, 0);
} else if (frequency > 0) {
double w0 = piDouble * frequency;
double S = 1; 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 {
setNormalizedCoefficients(index, A * A, 0, 0, 1, 0, 0);
}
}
void Biquad::setPeakingParams(size_t index, double frequency, double Q, double dbGain)
{
frequency = std::max(0.0, std::min(frequency, 1.0));
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 {
setNormalizedCoefficients(index, A * A, 0, 0, 1, 0, 0);
}
} else {
setNormalizedCoefficients(index, 1, 0, 0, 1, 0, 0);
}
}
void Biquad::setAllpassParams(size_t index, double frequency, double Q)
{
frequency = std::max(0.0, std::min(frequency, 1.0));
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 {
setNormalizedCoefficients(index, -1, 0, 0, 1, 0, 0);
}
} else {
setNormalizedCoefficients(index, 1, 0, 0, 1, 0, 0);
}
}
void Biquad::setNotchParams(size_t index, double frequency, double Q)
{
frequency = std::max(0.0, std::min(frequency, 1.0));
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 {
setNormalizedCoefficients(index, 0, 0, 0, 1, 0, 0);
}
} else {
setNormalizedCoefficients(index, 1, 0, 0, 1, 0, 0);
}
}
void Biquad::setBandpassParams(size_t index, double frequency, double Q)
{
frequency = std::max(0.0, frequency);
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 {
setNormalizedCoefficients(index, 1, 0, 0, 1, 0, 0);
}
} else {
setNormalizedCoefficients(index, 0, 0, 0, 1, 0, 0);
}
}
void Biquad::getFrequencyResponse(unsigned nFrequencies, const float* frequency, float* magResponse, float* phaseResponse)
{
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) {
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)
{
return (n - 2) * std::log(r) + std::log(std::fabs(c1 * (n + 1) * r * r + c2)) - logEPS;
}
static double rootFinder(double low, double high, double logEPS, double c1, double c2, double r)
{
constexpr double accuracyThreshold = 0.5;
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);
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) {
high = root;
fHigh = fr;
side = -1;
} else if (fLow * fr > 0) {
low = root;
fLow = fr;
if (side == 1)
fHigh /= 2;
side = 1;
} else {
break;
}
}
ASSERT(iteration < maxIterations);
return root;
}
double Biquad::tailFrame(size_t coefIndex, double maxFrame)
{
constexpr double maxTailAmplitude = 1 / 32768.0;
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) {
double rplus = (-a1 + std::sqrt(discrim)) / 2;
double rminus = (-a1 - std::sqrt(discrim)) / 2;
double r1 = std::fabs(rplus) >= std::fabs(rminus) ? rplus : rminus;
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));
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) {
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;
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) {
tailFrame = 0;
} else {
ASSERT(std::isfinite(tailFrame));
}
}
return tailFrame;
}
double r = -a1 / 2;
double tailFrame = 0;
if (!r) {
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));
if (!c1 && !c2)
tailFrame = 0;
else {
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;
}
}
#endif // ENABLE(WEB_AUDIO)