audioconvert: use the same biquad as filter-chain

This commit is contained in:
Wim Taymans 2024-11-27 16:47:53 +01:00
parent 98365470c7
commit 0a0b1d594b
4 changed files with 338 additions and 68 deletions

View file

@ -8,9 +8,6 @@
* found in the LICENSE.WEBKIT file.
*/
#include <spa/utils/defs.h>
#include <math.h>
#include "biquad.h"
@ -18,9 +15,8 @@
#define M_PI 3.14159265358979323846
#endif
#ifndef M_SQRT2
#define M_SQRT2 1.41421356237309504880
#endif
/* S = 1 in Q */
#define BIQUAD_SHELVING_DEFAULT_Q 0.707106781186548
static void set_coefficient(struct biquad *bq, double b0, double b1, double b2,
double a0, double a1, double a2)
@ -33,82 +29,346 @@ static void set_coefficient(struct biquad *bq, double b0, double b1, double b2,
bq->a2 = (float)(a2 * a0_inv);
}
static void biquad_lowpass(struct biquad *bq, double cutoff)
static void biquad_lowpass(struct biquad *bq, double cutoff, double resonance)
{
/* Limit cutoff to 0 to 1. */
cutoff = SPA_CLAMP(cutoff, 0.0, 1.0);
cutoff = fmax(0.0, fmin(cutoff, 1.0));
if (cutoff >= 1.0) {
/* When cutoff is 1, the z-transform is 1. */
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
} else if (cutoff > 0) {
/* Compute biquad coefficients for lowpass filter */
double theta = M_PI * cutoff;
double sn = 0.5 * M_SQRT2 * sin(theta);
double beta = 0.5 * (1 - sn) / (1 + sn);
double gamma_coeff = (0.5 + beta) * cos(theta);
double alpha = 0.25 * (0.5 + beta - gamma_coeff);
double b0 = 2 * alpha;
double b1 = 2 * 2 * alpha;
double b2 = 2 * alpha;
double a1 = 2 * -gamma_coeff;
double a2 = 2 * beta;
set_coefficient(bq, b0, b1, b2, 1, a1, a2);
} else {
/* When cutoff is zero, nothing gets through the filter, so set
if (cutoff == 1 || cutoff == 0) {
/* When cutoff is 1, the z-transform is 1.
* When cutoff is zero, nothing gets through the filter, so set
* coefficients up correctly.
*/
set_coefficient(bq, 0, 0, 0, 1, 0, 0);
set_coefficient(bq, cutoff, 0, 0, 1, 0, 0);
return;
}
/* Compute biquad coefficients for lowpass filter */
resonance = fmax(0.0, resonance); /* can't go negative */
double g = pow(10.0, 0.05 * resonance);
double d = sqrt((4 - sqrt(16 - 16 / (g * g))) / 2);
double theta = M_PI * cutoff;
double sn = 0.5 * d * sin(theta);
double beta = 0.5 * (1 - sn) / (1 + sn);
double gamma = (0.5 + beta) * cos(theta);
double alpha = 0.25 * (0.5 + beta - gamma);
double b0 = 2 * alpha;
double b1 = 2 * 2 * alpha;
double b2 = 2 * alpha;
double a1 = 2 * -gamma;
double a2 = 2 * beta;
set_coefficient(bq, b0, b1, b2, 1, a1, a2);
}
static void biquad_highpass(struct biquad *bq, double cutoff)
static void biquad_highpass(struct biquad *bq, double cutoff, double resonance)
{
/* Limit cutoff to 0 to 1. */
cutoff = SPA_CLAMP(cutoff, 0.0, 1.0);
cutoff = fmax(0.0, fmin(cutoff, 1.0));
if (cutoff >= 1.0) {
/* The z-transform is 0. */
set_coefficient(bq, 0, 0, 0, 1, 0, 0);
} else if (cutoff > 0) {
/* Compute biquad coefficients for highpass filter */
double theta = M_PI * cutoff;
double sn = 0.5 * M_SQRT2 * sin(theta);
double beta = 0.5 * (1 - sn) / (1 + sn);
double gamma_coeff = (0.5 + beta) * cos(theta);
double alpha = 0.25 * (0.5 + beta + gamma_coeff);
double b0 = 2 * alpha;
double b1 = 2 * -2 * alpha;
double b2 = 2 * alpha;
double a1 = 2 * -gamma_coeff;
double a2 = 2 * beta;
set_coefficient(bq, b0, b1, b2, 1, a1, a2);
} else {
if (cutoff == 1 || cutoff == 0) {
/* When cutoff is one, the z-transform is 0. */
/* 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.
*/
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
set_coefficient(bq, 1 - cutoff, 0, 0, 1, 0, 0);
return;
}
/* Compute biquad coefficients for highpass filter */
resonance = fmax(0.0, resonance); /* can't go negative */
double g = pow(10.0, 0.05 * resonance);
double d = sqrt((4 - sqrt(16 - 16 / (g * g))) / 2);
double theta = M_PI * cutoff;
double sn = 0.5 * d * sin(theta);
double beta = 0.5 * (1 - sn) / (1 + sn);
double gamma = (0.5 + beta) * cos(theta);
double alpha = 0.25 * (0.5 + beta + gamma);
double b0 = 2 * alpha;
double b1 = 2 * -2 * alpha;
double b2 = 2 * alpha;
double a1 = 2 * -gamma;
double a2 = 2 * beta;
set_coefficient(bq, b0, b1, b2, 1, a1, a2);
}
void biquad_set(struct biquad *bq, enum biquad_type type, double freq)
static void biquad_bandpass(struct biquad *bq, double frequency, double Q)
{
/* No negative frequencies allowed. */
frequency = fmax(0.0, frequency);
/* Don't let Q go negative, which causes an unstable filter. */
Q = fmax(0.0, Q);
if (frequency <= 0 || frequency >= 1) {
/* 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.
*/
set_coefficient(bq, 0, 0, 0, 1, 0, 0);
return;
}
if (Q <= 0) {
/* 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.
*/
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
return;
}
double w0 = M_PI * frequency;
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;
set_coefficient(bq, b0, b1, b2, a0, a1, a2);
}
static void biquad_lowshelf(struct biquad *bq, double frequency, double Q,
double db_gain)
{
/* Clip frequencies to between 0 and 1, inclusive. */
frequency = fmax(0.0, fmin(frequency, 1.0));
double A = pow(10.0, db_gain / 40);
if (frequency == 1) {
/* The z-transform is a constant gain. */
set_coefficient(bq, A * A, 0, 0, 1, 0, 0);
return;
}
if (frequency <= 0) {
/* When frequency is 0, the z-transform is 1. */
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
return;
}
/* Set Q to an equivalent value to S = 1 if not specified */
if (Q <= 0)
Q = BIQUAD_SHELVING_DEFAULT_Q;
double w0 = M_PI * frequency;
double alpha = sin(w0) / (2 * Q);
double k = cos(w0);
double k2 = 2 * sqrt(A) * alpha;
double a_plus_one = A + 1;
double a_minus_one = A - 1;
double b0 = A * (a_plus_one - a_minus_one * k + k2);
double b1 = 2 * A * (a_minus_one - a_plus_one * k);
double b2 = A * (a_plus_one - a_minus_one * k - k2);
double a0 = a_plus_one + a_minus_one * k + k2;
double a1 = -2 * (a_minus_one + a_plus_one * k);
double a2 = a_plus_one + a_minus_one * k - k2;
set_coefficient(bq, b0, b1, b2, a0, a1, a2);
}
static void biquad_highshelf(struct biquad *bq, double frequency, double Q,
double db_gain)
{
/* Clip frequencies to between 0 and 1, inclusive. */
frequency = fmax(0.0, fmin(frequency, 1.0));
double A = pow(10.0, db_gain / 40);
if (frequency == 1) {
/* The z-transform is 1. */
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
return;
}
if (frequency <= 0) {
/* When frequency = 0, the filter is just a gain, A^2. */
set_coefficient(bq, A * A, 0, 0, 1, 0, 0);
return;
}
/* Set Q to an equivalent value to S = 1 if not specified */
if (Q <= 0)
Q = BIQUAD_SHELVING_DEFAULT_Q;
double w0 = M_PI * frequency;
double alpha = sin(w0) / (2 * Q);
double k = cos(w0);
double k2 = 2 * sqrt(A) * alpha;
double a_plus_one = A + 1;
double a_minus_one = A - 1;
double b0 = A * (a_plus_one + a_minus_one * k + k2);
double b1 = -2 * A * (a_minus_one + a_plus_one * k);
double b2 = A * (a_plus_one + a_minus_one * k - k2);
double a0 = a_plus_one - a_minus_one * k + k2;
double a1 = 2 * (a_minus_one - a_plus_one * k);
double a2 = a_plus_one - a_minus_one * k - k2;
set_coefficient(bq, b0, b1, b2, a0, a1, a2);
}
static void biquad_peaking(struct biquad *bq, double frequency, double Q,
double db_gain)
{
/* Clip frequencies to between 0 and 1, inclusive. */
frequency = fmax(0.0, fmin(frequency, 1.0));
/* Don't let Q go negative, which causes an unstable filter. */
Q = fmax(0.0, Q);
double A = pow(10.0, db_gain / 40);
if (frequency <= 0 || frequency >= 1) {
/* When frequency is 0 or 1, the z-transform is 1. */
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
return;
}
if (Q <= 0) {
/* 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.
*/
set_coefficient(bq, A * A, 0, 0, 1, 0, 0);
return;
}
double w0 = M_PI * 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;
set_coefficient(bq, b0, b1, b2, a0, a1, a2);
}
static void biquad_notch(struct biquad *bq, double frequency, double Q)
{
/* Clip frequencies to between 0 and 1, inclusive. */
frequency = fmax(0.0, fmin(frequency, 1.0));
/* Don't let Q go negative, which causes an unstable filter. */
Q = fmax(0.0, Q);
if (frequency <= 0 || frequency >= 1) {
/* When frequency is 0 or 1, the z-transform is 1. */
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
return;
}
if (Q <= 0) {
/* 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.
*/
set_coefficient(bq, 0, 0, 0, 1, 0, 0);
return;
}
double w0 = M_PI * 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;
set_coefficient(bq, b0, b1, b2, a0, a1, a2);
}
static void biquad_allpass(struct biquad *bq, double frequency, double Q)
{
/* Clip frequencies to between 0 and 1, inclusive. */
frequency = fmax(0.0, fmin(frequency, 1.0));
/* Don't let Q go negative, which causes an unstable filter. */
Q = fmax(0.0, Q);
if (frequency <= 0 || frequency >= 1) {
/* When frequency is 0 or 1, the z-transform is 1. */
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
return;
}
if (Q <= 0) {
/* 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.
*/
set_coefficient(bq, -1, 0, 0, 1, 0, 0);
return;
}
double w0 = M_PI * 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;
set_coefficient(bq, b0, b1, b2, a0, a1, a2);
}
void biquad_set(struct biquad *bq, enum biquad_type type, double freq, double Q,
double gain)
{
/* Clear history values. */
bq->type = type;
bq->x1 = 0;
bq->x2 = 0;
switch (type) {
case BQ_NONE:
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
break;
case BQ_LOWPASS:
biquad_lowpass(bq, freq);
biquad_lowpass(bq, freq, Q);
break;
case BQ_HIGHPASS:
biquad_highpass(bq, freq);
biquad_highpass(bq, freq, Q);
break;
case BQ_BANDPASS:
biquad_bandpass(bq, freq, Q);
break;
case BQ_LOWSHELF:
biquad_lowshelf(bq, freq, Q, gain);
break;
case BQ_HIGHSHELF:
biquad_highshelf(bq, freq, Q, gain);
break;
case BQ_PEAKING:
biquad_peaking(bq, freq, Q, gain);
break;
case BQ_NOTCH:
biquad_notch(bq, freq, Q);
break;
case BQ_ALLPASS:
biquad_allpass(bq, freq, Q);
break;
case BQ_NONE:
case BQ_RAW:
/* Default is an identity filter. */
set_coefficient(bq, 1, 0, 0, 1, 0, 0);
break;
}
}

View file

@ -10,6 +10,20 @@
extern "C" {
#endif
/* The type of the biquad filters */
enum biquad_type {
BQ_NONE,
BQ_LOWPASS,
BQ_HIGHPASS,
BQ_BANDPASS,
BQ_LOWSHELF,
BQ_HIGHSHELF,
BQ_PEAKING,
BQ_NOTCH,
BQ_ALLPASS,
BQ_RAW,
};
/* The biquad filter parameters. The transfer function H(z) is (b0 + b1 * z^(-1)
* + b2 * z^(-2)) / (1 + a1 * z^(-1) + a2 * z^(-2)). The previous two inputs
* are stored in x1 and x2, and the previous two outputs are stored in y1 and
@ -19,15 +33,10 @@ extern "C" {
* float is used during the actual filtering for faster computation.
*/
struct biquad {
enum biquad_type type;
float b0, b1, b2;
float a1, a2;
};
/* The type of the biquad filters */
enum biquad_type {
BQ_NONE,
BQ_LOWPASS,
BQ_HIGHPASS,
float x1, x2;
};
/* Initialize a biquad filter parameters from its type and parameters.
@ -36,8 +45,11 @@ enum biquad_type {
* type - The type of the biquad filter.
* frequency - The value should be in the range [0, 1]. It is relative to
* half of the sampling rate.
* Q - Quality factor. See Web Audio API for details.
* gain - The value is in dB. See Web Audio API for details.
*/
void biquad_set(struct biquad *bq, enum biquad_type type, double freq);
void biquad_set(struct biquad *bq, enum biquad_type type, double freq, double Q,
double gain);
#ifdef __cplusplus
} /* extern "C" */

View file

@ -10,7 +10,7 @@
void lr4_set(struct lr4 *lr4, enum biquad_type type, float freq)
{
biquad_set(&lr4->bq, type, freq);
biquad_set(&lr4->bq, type, freq, 0, 0);
lr4->x1 = 0;
lr4->x2 = 0;
lr4->y1 = 0;

View file

@ -25,6 +25,4 @@ struct lr4 {
void lr4_set(struct lr4 *lr4, enum biquad_type type, float freq);
void lr4_process(struct lr4 *lr4, float *dst, const float *src, const float vol, int samples);
#endif /* CROSSOVER_H_ */