diff --git a/spa/plugins/audioconvert/biquad.c b/spa/plugins/audioconvert/biquad.c index 27c3fee59..68a6d8775 100644 --- a/spa/plugins/audioconvert/biquad.c +++ b/spa/plugins/audioconvert/biquad.c @@ -8,9 +8,6 @@ * found in the LICENSE.WEBKIT file. */ - -#include - #include #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; } } diff --git a/spa/plugins/audioconvert/biquad.h b/spa/plugins/audioconvert/biquad.h index 519a3275b..3344598e6 100644 --- a/spa/plugins/audioconvert/biquad.h +++ b/spa/plugins/audioconvert/biquad.h @@ -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" */ diff --git a/spa/plugins/audioconvert/crossover.c b/spa/plugins/audioconvert/crossover.c index ba3eaf023..537d4937f 100644 --- a/spa/plugins/audioconvert/crossover.c +++ b/spa/plugins/audioconvert/crossover.c @@ -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; diff --git a/spa/plugins/audioconvert/crossover.h b/spa/plugins/audioconvert/crossover.h index 34bf80f80..159d95903 100644 --- a/spa/plugins/audioconvert/crossover.h +++ b/spa/plugins/audioconvert/crossover.h @@ -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_ */