//// AudioReverb.cpp// libraries/audio/src//// Created by Ken Cooke on 10/11/15.// Copyright 2015 High Fidelity, Inc.// Copyright 2026 Overte e.V.// #include "AudioReverb.h" #include <cmath>#include <cstdint>#include <cstring> #ifdef _MSC_VER#include <intrin.h>#define MULHI(a,b) ((int32_t)(__emul(a, b) >> 32))#else#define MULHI(a,b) ((int32_t)(((int64_t)(a) * (int64_t)(b)) >> 32))#endif #ifndef MAX#define MAX(a,b) (((a) > (b)) ? (a) : (b))#endif#ifndef MIN#define MIN(a,b) (((a) < (b)) ? (a) : (b))#endif static const float PHI = 0.6180339887f; // maximum allpass diffusionstatic const float TWOPI = 6.283185307f; static const double PI = 3.14159265358979323846;static const double SQRT2 = 1.41421356237309504880; static const double FIXQ31 = 2147483648.0;static const double FIXQ32 = 4294967296.0; //// Allpass delay modulation//static const int MOD_INTBITS = 4;static const int MOD_FRACBITS = 31 - MOD_INTBITS;static const uint32_t MOD_FRACMASK = (1 << MOD_FRACBITS) - 1;static const float QMOD_TO_FLOAT = 1.0f / (1 << MOD_FRACBITS); //// Reverb delay values, defined for sampleRate=48000 roomSize=100% density=100%//// max path should already be prime, to prevent getPrime(NEXTPOW2(N)) > NEXTPOW2(N)//static const int M_PD0 = 16000; // max predelay = 333ms static const int M_AP0 = 83;static const int M_AP1 = 211;static const int M_AP2 = 311;static const int M_AP3 = 97;static const int M_AP4 = 223;static const int M_AP5 = 293; static const int M_MT0 = 1017;static const int M_MT1 = 4077;static const int M_MT2 = 2039;static const int M_MT3 = 1017;static const int M_MT4 = 2593;static const int M_MT5 = 2039; static const int M_MT1_2 = 600;static const int M_MT4_2 = 600;static const int M_MT1_MAX = MAX(M_MT1, M_MT1_2);static const int M_MT4_MAX = MAX(M_MT4, M_MT4_2); static const int M_AP6 = 817;static const int M_AP7 = 513;static const int M_AP8 = 765;static const int M_AP9 = 465;static const int M_AP10 = 3021;static const int M_AP11 = 2121;static const int M_AP12 = 1705;static const int M_AP13 = 1081;static const int M_AP14 = 3313;static const int M_AP15 = 2205;static const int M_AP16 = 1773;static const int M_AP17 = 981; static const int M_AP7_MAX = M_AP7 + (1 << MOD_INTBITS) + 1; // include max excursionstatic const int M_AP9_MAX = M_AP9 + (1 << MOD_INTBITS) + 1; static const int M_MT6 = 6863;static const int M_MT7 = 7639;static const int M_MT8 = 3019;static const int M_MT9 = 2875; static const int M_LD0 = 8000; // max late delay = 166msstatic const int M_MT6_MAX = MAX(M_MT6, M_LD0);static const int M_MT7_MAX = MAX(M_MT7, M_LD0);static const int M_MT8_MAX = MAX(M_MT8, M_LD0);static const int M_MT9_MAX = MAX(M_MT9, M_LD0); static const int M_AP18 = 131;static const int M_AP19 = 113;static const int M_AP20 = 107;static const int M_AP21 = 127; // Round an integer to the next power-of-two, at compile timeconstexpr uint32_t NEXTPOW2(uint32_t n) { n -= 1; n |= (n >> 1); n |= (n >> 2); n |= (n >> 4); n |= (n >> 8); n |= (n >> 16); n += 1; return n;} //// Filter design tools using analog-matched response.// All filter types approximate the s-plane response, including cutoff > Nyquist.// static float dBToGain(float dB) { return powf(10.0f, dB * (1/20.0f));} static double dBToGain(double dB) { return pow(10.0, dB * (1/20.0));} // Returns the gain of analog (s-plane) peak-filter evaluated at wstatic double analogPeak(double w0, double G, double Q, double w) { double w0sq, wsq, Gsq, Qsq; double num, den; w0sq = w0 * w0; wsq = w * w; Gsq = G * G; Qsq = Q * Q; num = Qsq * (wsq - w0sq) * (wsq - w0sq) + wsq * w0sq * Gsq; den = Qsq * (wsq - w0sq) * (wsq - w0sq) + wsq * w0sq; return sqrt(num / den);} // Returns the gain of analog (s-plane) shelf evaluated at w// non-resonant Q = sqrt(0.5) so 1/Q^2 = 2.0static double analogShelf(double w0, double G, double resonance, int isHigh, double w) { double w0sq, wsq, Qrsq; double num, den; resonance = MIN(MAX(resonance, 0.0), 1.0); Qrsq = 2.0 * pow(G, 1.0 - resonance); // resonant 1/Q^2 w0sq = w0 * w0; wsq = w * w; if (isHigh) { num = (G*wsq - w0sq) * (G*wsq - w0sq) + wsq * w0sq * Qrsq; } else { num = (wsq - G*w0sq) * (wsq - G*w0sq) + wsq * w0sq * Qrsq; } den = wsq * wsq + w0sq * w0sq; return sqrt(num / den);} // Returns the gain of analog (s-plane) highpass/lowpass evaluated at w// Q = sqrt(0.5) = 2nd order Butterworthstatic double analogFilter(double w0, int isHigh, double w) { double w0sq, wsq; double num, den; w0sq = w0 * w0; wsq = w * w; if (isHigh) { num = wsq * wsq; } else { num = w0sq * w0sq; } den = wsq * wsq + w0sq * w0sq; return sqrt(num / den);} //// Biquad Peaking EQ using analog matching.// NOTE: Peak topology becomes ill-conditioned as w0 approaches pi.//static void BQPeakBelowPi(double coef[5], double w0, double dbgain, double Q) { double G, G1, Gsq, G1sq, Qsq; double G11, G00, Gratio; double wpi, wh, w0sq, whsq; double Wsq, B, A; double b0, b1, b2, a0, a1, a2; double temp, scale; int isUnity; // convert cut into boost, invert later G = dBToGain(fabs(dbgain)); // special case near unity gain isUnity = G < 1.001; G = MAX(G, 1.001); // compute the Nyquist gain wpi = w0 + 2.8 * (1.0 - w0/PI); // minimax-like error wpi = MIN(wpi, PI); G1 = analogPeak(w0, G, Q, wpi); G1sq = G1 * G1; Gsq = G * G; Qsq = Q * Q; // compute the analog half-gain frequency temp = G + 2.0 * Qsq - sqrt(Gsq + 4.0 * Qsq * G); wh = sqrt(temp) * w0 / (Q * SQRT2); // prewarp freqs of w0 and wh w0 = tan(0.5 * w0); wh = tan(0.5 * wh); w0sq = w0 * w0; whsq = wh * wh; // compute Wsq, from asymmetry due to G1 G11 = Gsq - G1sq; G00 = Gsq - 1.0; Gratio = G11 / G00; Wsq = w0sq * sqrt(Gratio); // compute B, matching gains at w0 and wh temp = 2.0 * Wsq * (G1 - sqrt(Gratio)); temp += Gsq * whsq * (1.0 - G1sq) / G00; temp += G * Gratio * (w0sq - whsq) * (w0sq - whsq) / whsq; B = sqrt(temp); // compute A, matching gains at w0 and wh temp = 2.0 * Wsq; temp += (2.0 * G11 * w0sq + (G1sq - G) * whsq) / (G - Gsq); temp += Gratio * w0sq * w0sq / (G * whsq); A = sqrt(temp); // design digital filter via bilinear transform b0 = G1 + B + Wsq; b1 = 2.0 * (Wsq - G1); b2 = G1 - B + Wsq; a0 = 1.0 + A + Wsq; a1 = 2.0 * (Wsq - 1.0); a2 = 1.0 - A + Wsq; // for unity gain, ensure poles/zeros in the right place. // needed for smooth interpolation when gain is changed. if (isUnity) { b0 = a0; b1 = a1; b2 = a2; } // invert filter for cut if (dbgain < 0.0) { temp = b0; b0 = a0; a0 = temp; temp = b1; b1 = a1; a1 = temp; temp = b2; b2 = a2; a2 = temp; } // normalize scale = 1.0 / a0; coef[0] = b0 * scale; coef[1] = b1 * scale; coef[2] = b2 * scale; coef[3] = a1 * scale; coef[4] = a2 * scale;} //// Biquad Peaking EQ using a shelf instead of peak.//// This uses a shelf topology, matched to the analog peaking filter located above Nyquist.//// NOTE: the result is close, but not identical to BQPeakBelowPi(), since a pole/zero// pair must jump to the right side of the real axis. Eg. inflection at the peak goes away.// However, interpolation from peak to shelf is well behaved if the switch is made near pi,// as the pole/zero be near (and travel down) the real axis.//static void BQPeakAbovePi(double coef[5], double w0, double dbgain, double Q) { double G, G1, Gsq, Qsq; double wpi, wh, wn, wd; double wna, wda; double gn, gd, gnsq, gdsq; double num, den; double Wnsq, Wdsq, B, A; double b0, b1, b2, a0, a1, a2; double temp, scale; int isUnity; // convert cut into boost, invert later G = dBToGain(fabs(dbgain)); // special case near unity gain isUnity = G < 1.001; G = MAX(G, 1.001); // compute the Nyquist gain wpi = PI; if (w0 < PI) { G1 = G; // use the peak gain } else { G1 = analogPeak(w0, G, Q, wpi); // use Nyquist gain of analog response } Gsq = G * G; Qsq = Q * Q; // compute the analog half-gain frequency temp = G + 2.0 * Qsq - sqrt(Gsq + 4.0 * Qsq * G); wh = sqrt(temp) * w0 / (Q * SQRT2); // approximate wn and wd // use half-gain frequency as mapping wn = 0.5 * wh / sqrt(sqrt((G1))); wd = wn * sqrt(G1); Wnsq = wn * wn; Wdsq = wd * wd; // analog freqs of wn and wd wna = 2.0 * atan(wn); wda = 2.0 * atan(wd); // normalized analog gains at wna and wda temp = 1.0 / G1; gn = temp * analogPeak(w0, G, Q, wna); gd = temp * analogPeak(w0, G, Q, wda); gnsq = gn * gn; gdsq = gd * gd; // compute B, matching gains at wn and wd temp = 1.0 / (wn * wd); den = fabs(gnsq - gdsq); num = gnsq * (Wnsq - Wdsq) * (Wnsq - Wdsq) * (Wnsq + gdsq * Wdsq); B = temp * sqrt(num / den); // compute A, matching gains at wn and wd num = (Wnsq - Wdsq) * (Wnsq - Wdsq) * (Wnsq + gnsq * Wdsq); A = temp * sqrt(num / den); // design digital filter via bilinear transform b0 = G1 * (1.0 + B + Wnsq); b1 = G1 * 2.0 * (Wnsq - 1.0); b2 = G1 * (1.0 - B + Wnsq); a0 = 1.0 + A + Wdsq; a1 = 2.0 * (Wdsq - 1.0); a2 = 1.0 - A + Wdsq; // for unity gain, ensure poles/zeros are in the right place. // allows smooth coefficient interpolation when gain is changed. if (isUnity) { b0 = a0; b1 = a1; b2 = a2; } // invert filter for cut if (dbgain < 0.0) { temp = b0; b0 = a0; a0 = temp; temp = b1; b1 = a1; a1 = temp; temp = b2; b2 = a2; a2 = temp; } // normalize scale = 1.0 / a0; coef[0] = b0 * scale; coef[1] = b1 * scale; coef[2] = b2 * scale; coef[3] = a1 * scale; coef[4] = a2 * scale;} //// Biquad Peaking EQ using analog matching.// Supports full range of w0.//void BQPeak(double coef[5], double w0, double dbgain, double Q) { w0 = MAX(w0, 0.0); // allow w0 > pi Q = MIN(MAX(Q, 1.0e-6), 1.0e+6); // Switch from peak to shelf, just before w0 = pi. // Too early causes a jump in the peak location, which interpolates to lowered gains. // Too late becomes ill-conditioned, and makes no improvement to interpolated response. // 3.14 is a good choice. if (w0 > 3.14) { BQPeakAbovePi(coef, w0, dbgain, Q); } else { BQPeakBelowPi(coef, w0, dbgain, Q); }} //// Biquad Shelf using analog matching.//void BQShelf(double coef[5], double w0, double dbgain, double resonance, int isHigh) { double G, G1; double wpi, wn, wd; double wna, wda; double gn, gd, gnsq, gdsq; double num, den; double Wnsq, Wdsq, B, A; double b0, b1, b2, a0, a1, a2; double temp, scale; int isUnity; w0 = MAX(w0, 0.0); // allow w0 > pi resonance = MIN(MAX(resonance, 0.0), 1.0); // convert cut into boost, invert later G = dBToGain(fabs(dbgain)); // special case near unity gain isUnity = G < 1.001; G = MAX(G, 1.001); // compute the Nyquist gain wpi = w0 + 2.8 * (1.0 - w0/PI); // minimax-like error wpi = MIN(wpi, PI); G1 = analogShelf(w0, G, resonance, isHigh, wpi); // approximate wn and wd if (isHigh) { // use center as mapping wn = 0.5 * w0 / sqrt(sqrt((G * G1))); wd = wn * sqrt(G1); } else { // use wd as mapping wd = 0.5 * w0; wn = wd * sqrt(G/G1); } Wnsq = wn * wn; Wdsq = wd * wd; // analog freqs of wn and wd wna = 2.0 * atan(wn); wda = 2.0 * atan(wd); // normalized analog gains at wna and wda temp = 1.0 / G1; gn = temp * analogShelf(w0, G, resonance, isHigh, wna); gd = temp * analogShelf(w0, G, resonance, isHigh, wda); gnsq = gn * gn; gdsq = gd * gd; // compute B, matching gains at wn and wd temp = 1.0 / (wn * wd); den = fabs(gnsq - gdsq); num = gnsq * (Wnsq - Wdsq) * (Wnsq - Wdsq) * (Wnsq + gdsq * Wdsq); B = temp * sqrt(num / den); // compute A, matching gains at wn and wd num = (Wnsq - Wdsq) * (Wnsq - Wdsq) * (Wnsq + gnsq * Wdsq); A = temp * sqrt(num / den); // design digital filter via bilinear transform b0 = G1 * (1.0 + B + Wnsq); b1 = G1 * 2.0 * (Wnsq - 1.0); b2 = G1 * (1.0 - B + Wnsq); a0 = 1.0 + A + Wdsq; a1 = 2.0 * (Wdsq - 1.0); a2 = 1.0 - A + Wdsq; // for unity gain, ensure poles/zeros in the right place. // needed for smooth interpolation when gain is changed. if (isUnity) { b0 = a0; b1 = a1; b2 = a2; } // invert filter for cut if (dbgain < 0.0) { temp = b0; b0 = a0; a0 = temp; temp = b1; b1 = a1; a1 = temp; temp = b2; b2 = a2; a2 = temp; } // normalize scale = 1.0 / a0; coef[0] = b0 * scale; coef[1] = b1 * scale; coef[2] = b2 * scale; coef[3] = a1 * scale; coef[4] = a2 * scale;} //// Biquad Lowpass/Highpass using analog matching.// Q = sqrt(0.5) = 2nd order Butterworth//void BQFilter(double coef[5], double w0, int isHigh) { double G1;Lines 1–500 of 2055. Load more
// // AudioReverb.cpp // libraries/audio/src // // Created by Ken Cooke on 10/11/15. // Copyright 2015 High Fidelity, Inc. // Copyright 2026 Overte e.V. // #include "AudioReverb.h" #include <cmath> #include <cstdint> #include <cstring> #ifdef _MSC_VER #include <intrin.h> #define MULHI(a,b) ((int32_t)(__emul(a, b) >> 32)) #else #define MULHI(a,b) ((int32_t)(((int64_t)(a) * (int64_t)(b)) >> 32)) #endif #ifndef MAX #define MAX(a,b) (((a) > (b)) ? (a) : (b)) #endif #ifndef MIN #define MIN(a,b) (((a) < (b)) ? (a) : (b)) #endif static const float PHI = 0.6180339887f; // maximum allpass diffusion static const float TWOPI = 6.283185307f; static const double PI = 3.14159265358979323846; static const double SQRT2 = 1.41421356237309504880; static const double FIXQ31 = 2147483648.0; static const double FIXQ32 = 4294967296.0; // // Allpass delay modulation // static const int MOD_INTBITS = 4; static const int MOD_FRACBITS = 31 - MOD_INTBITS; static const uint32_t MOD_FRACMASK = (1 << MOD_FRACBITS) - 1; static const float QMOD_TO_FLOAT = 1.0f / (1 << MOD_FRACBITS); // // Reverb delay values, defined for sampleRate=48000 roomSize=100% density=100% // // max path should already be prime, to prevent getPrime(NEXTPOW2(N)) > NEXTPOW2(N) // static const int M_PD0 = 16000; // max predelay = 333ms static const int M_AP0 = 83; static const int M_AP1 = 211; static const int M_AP2 = 311; static const int M_AP3 = 97; static const int M_AP4 = 223; static const int M_AP5 = 293; static const int M_MT0 = 1017; static const int M_MT1 = 4077; static const int M_MT2 = 2039; static const int M_MT3 = 1017; static const int M_MT4 = 2593; static const int M_MT5 = 2039; static const int M_MT1_2 = 600; static const int M_MT4_2 = 600; static const int M_MT1_MAX = MAX(M_MT1, M_MT1_2); static const int M_MT4_MAX = MAX(M_MT4, M_MT4_2); static const int M_AP6 = 817; static const int M_AP7 = 513; static const int M_AP8 = 765; static const int M_AP9 = 465; static const int M_AP10 = 3021; static const int M_AP11 = 2121; static const int M_AP12 = 1705; static const int M_AP13 = 1081; static const int M_AP14 = 3313; static const int M_AP15 = 2205; static const int M_AP16 = 1773; static const int M_AP17 = 981; static const int M_AP7_MAX = M_AP7 + (1 << MOD_INTBITS) + 1; // include max excursion static const int M_AP9_MAX = M_AP9 + (1 << MOD_INTBITS) + 1; static const int M_MT6 = 6863; static const int M_MT7 = 7639; static const int M_MT8 = 3019; static const int M_MT9 = 2875; static const int M_LD0 = 8000; // max late delay = 166ms static const int M_MT6_MAX = MAX(M_MT6, M_LD0); static const int M_MT7_MAX = MAX(M_MT7, M_LD0); static const int M_MT8_MAX = MAX(M_MT8, M_LD0); static const int M_MT9_MAX = MAX(M_MT9, M_LD0); static const int M_AP18 = 131; static const int M_AP19 = 113; static const int M_AP20 = 107; static const int M_AP21 = 127; // Round an integer to the next power-of-two, at compile time constexpr uint32_t NEXTPOW2(uint32_t n) { n -= 1; n |= (n >> 1); n |= (n >> 2); n |= (n >> 4); n |= (n >> 8); n |= (n >> 16); n += 1; return n; } // // Filter design tools using analog-matched response. // All filter types approximate the s-plane response, including cutoff > Nyquist. // static float dBToGain(float dB) { return powf(10.0f, dB * (1/20.0f)); } static double dBToGain(double dB) { return pow(10.0, dB * (1/20.0)); } // Returns the gain of analog (s-plane) peak-filter evaluated at w static double analogPeak(double w0, double G, double Q, double w) { double w0sq, wsq, Gsq, Qsq; double num, den; w0sq = w0 * w0; wsq = w * w; Gsq = G * G; Qsq = Q * Q; num = Qsq * (wsq - w0sq) * (wsq - w0sq) + wsq * w0sq * Gsq; den = Qsq * (wsq - w0sq) * (wsq - w0sq) + wsq * w0sq; return sqrt(num / den); } // Returns the gain of analog (s-plane) shelf evaluated at w // non-resonant Q = sqrt(0.5) so 1/Q^2 = 2.0 static double analogShelf(double w0, double G, double resonance, int isHigh, double w) { double w0sq, wsq, Qrsq; double num, den; resonance = MIN(MAX(resonance, 0.0), 1.0); Qrsq = 2.0 * pow(G, 1.0 - resonance); // resonant 1/Q^2 w0sq = w0 * w0; wsq = w * w; if (isHigh) { num = (G*wsq - w0sq) * (G*wsq - w0sq) + wsq * w0sq * Qrsq; } else { num = (wsq - G*w0sq) * (wsq - G*w0sq) + wsq * w0sq * Qrsq; } den = wsq * wsq + w0sq * w0sq; return sqrt(num / den); } // Returns the gain of analog (s-plane) highpass/lowpass evaluated at w // Q = sqrt(0.5) = 2nd order Butterworth static double analogFilter(double w0, int isHigh, double w) { double w0sq, wsq; double num, den; w0sq = w0 * w0; wsq = w * w; if (isHigh) { num = wsq * wsq; } else { num = w0sq * w0sq; } den = wsq * wsq + w0sq * w0sq; return sqrt(num / den); } // // Biquad Peaking EQ using analog matching. // NOTE: Peak topology becomes ill-conditioned as w0 approaches pi. // static void BQPeakBelowPi(double coef[5], double w0, double dbgain, double Q) { double G, G1, Gsq, G1sq, Qsq; double G11, G00, Gratio; double wpi, wh, w0sq, whsq; double Wsq, B, A; double b0, b1, b2, a0, a1, a2; double temp, scale; int isUnity; // convert cut into boost, invert later G = dBToGain(fabs(dbgain)); // special case near unity gain isUnity = G < 1.001; G = MAX(G, 1.001); // compute the Nyquist gain wpi = w0 + 2.8 * (1.0 - w0/PI); // minimax-like error wpi = MIN(wpi, PI); G1 = analogPeak(w0, G, Q, wpi); G1sq = G1 * G1; Gsq = G * G; Qsq = Q * Q; // compute the analog half-gain frequency temp = G + 2.0 * Qsq - sqrt(Gsq + 4.0 * Qsq * G); wh = sqrt(temp) * w0 / (Q * SQRT2); // prewarp freqs of w0 and wh w0 = tan(0.5 * w0); wh = tan(0.5 * wh); w0sq = w0 * w0; whsq = wh * wh; // compute Wsq, from asymmetry due to G1 G11 = Gsq - G1sq; G00 = Gsq - 1.0; Gratio = G11 / G00; Wsq = w0sq * sqrt(Gratio); // compute B, matching gains at w0 and wh temp = 2.0 * Wsq * (G1 - sqrt(Gratio)); temp += Gsq * whsq * (1.0 - G1sq) / G00; temp += G * Gratio * (w0sq - whsq) * (w0sq - whsq) / whsq; B = sqrt(temp); // compute A, matching gains at w0 and wh temp = 2.0 * Wsq; temp += (2.0 * G11 * w0sq + (G1sq - G) * whsq) / (G - Gsq); temp += Gratio * w0sq * w0sq / (G * whsq); A = sqrt(temp); // design digital filter via bilinear transform b0 = G1 + B + Wsq; b1 = 2.0 * (Wsq - G1); b2 = G1 - B + Wsq; a0 = 1.0 + A + Wsq; a1 = 2.0 * (Wsq - 1.0); a2 = 1.0 - A + Wsq; // for unity gain, ensure poles/zeros in the right place. // needed for smooth interpolation when gain is changed. if (isUnity) { b0 = a0; b1 = a1; b2 = a2; } // invert filter for cut if (dbgain < 0.0) { temp = b0; b0 = a0; a0 = temp; temp = b1; b1 = a1; a1 = temp; temp = b2; b2 = a2; a2 = temp; } // normalize scale = 1.0 / a0; coef[0] = b0 * scale; coef[1] = b1 * scale; coef[2] = b2 * scale; coef[3] = a1 * scale; coef[4] = a2 * scale; } // // Biquad Peaking EQ using a shelf instead of peak. // // This uses a shelf topology, matched to the analog peaking filter located above Nyquist. // // NOTE: the result is close, but not identical to BQPeakBelowPi(), since a pole/zero // pair must jump to the right side of the real axis. Eg. inflection at the peak goes away. // However, interpolation from peak to shelf is well behaved if the switch is made near pi, // as the pole/zero be near (and travel down) the real axis. // static void BQPeakAbovePi(double coef[5], double w0, double dbgain, double Q) { double G, G1, Gsq, Qsq; double wpi, wh, wn, wd; double wna, wda; double gn, gd, gnsq, gdsq; double num, den; double Wnsq, Wdsq, B, A; double b0, b1, b2, a0, a1, a2; double temp, scale; int isUnity; // convert cut into boost, invert later G = dBToGain(fabs(dbgain)); // special case near unity gain isUnity = G < 1.001; G = MAX(G, 1.001); // compute the Nyquist gain wpi = PI; if (w0 < PI) { G1 = G; // use the peak gain } else { G1 = analogPeak(w0, G, Q, wpi); // use Nyquist gain of analog response } Gsq = G * G; Qsq = Q * Q; // compute the analog half-gain frequency temp = G + 2.0 * Qsq - sqrt(Gsq + 4.0 * Qsq * G); wh = sqrt(temp) * w0 / (Q * SQRT2); // approximate wn and wd // use half-gain frequency as mapping wn = 0.5 * wh / sqrt(sqrt((G1))); wd = wn * sqrt(G1); Wnsq = wn * wn; Wdsq = wd * wd; // analog freqs of wn and wd wna = 2.0 * atan(wn); wda = 2.0 * atan(wd); // normalized analog gains at wna and wda temp = 1.0 / G1; gn = temp * analogPeak(w0, G, Q, wna); gd = temp * analogPeak(w0, G, Q, wda); gnsq = gn * gn; gdsq = gd * gd; // compute B, matching gains at wn and wd temp = 1.0 / (wn * wd); den = fabs(gnsq - gdsq); num = gnsq * (Wnsq - Wdsq) * (Wnsq - Wdsq) * (Wnsq + gdsq * Wdsq); B = temp * sqrt(num / den); // compute A, matching gains at wn and wd num = (Wnsq - Wdsq) * (Wnsq - Wdsq) * (Wnsq + gnsq * Wdsq); A = temp * sqrt(num / den); // design digital filter via bilinear transform b0 = G1 * (1.0 + B + Wnsq); b1 = G1 * 2.0 * (Wnsq - 1.0); b2 = G1 * (1.0 - B + Wnsq); a0 = 1.0 + A + Wdsq; a1 = 2.0 * (Wdsq - 1.0); a2 = 1.0 - A + Wdsq; // for unity gain, ensure poles/zeros are in the right place. // allows smooth coefficient interpolation when gain is changed. if (isUnity) { b0 = a0; b1 = a1; b2 = a2; } // invert filter for cut if (dbgain < 0.0) { temp = b0; b0 = a0; a0 = temp; temp = b1; b1 = a1; a1 = temp; temp = b2; b2 = a2; a2 = temp; } // normalize scale = 1.0 / a0; coef[0] = b0 * scale; coef[1] = b1 * scale; coef[2] = b2 * scale; coef[3] = a1 * scale; coef[4] = a2 * scale; } // // Biquad Peaking EQ using analog matching. // Supports full range of w0. // void BQPeak(double coef[5], double w0, double dbgain, double Q) { w0 = MAX(w0, 0.0); // allow w0 > pi Q = MIN(MAX(Q, 1.0e-6), 1.0e+6); // Switch from peak to shelf, just before w0 = pi. // Too early causes a jump in the peak location, which interpolates to lowered gains. // Too late becomes ill-conditioned, and makes no improvement to interpolated response. // 3.14 is a good choice. if (w0 > 3.14) { BQPeakAbovePi(coef, w0, dbgain, Q); } else { BQPeakBelowPi(coef, w0, dbgain, Q); } } // // Biquad Shelf using analog matching. // void BQShelf(double coef[5], double w0, double dbgain, double resonance, int isHigh) { double G, G1; double wpi, wn, wd; double wna, wda; double gn, gd, gnsq, gdsq; double num, den; double Wnsq, Wdsq, B, A; double b0, b1, b2, a0, a1, a2; double temp, scale; int isUnity; w0 = MAX(w0, 0.0); // allow w0 > pi resonance = MIN(MAX(resonance, 0.0), 1.0); // convert cut into boost, invert later G = dBToGain(fabs(dbgain)); // special case near unity gain isUnity = G < 1.001; G = MAX(G, 1.001); // compute the Nyquist gain wpi = w0 + 2.8 * (1.0 - w0/PI); // minimax-like error wpi = MIN(wpi, PI); G1 = analogShelf(w0, G, resonance, isHigh, wpi); // approximate wn and wd if (isHigh) { // use center as mapping wn = 0.5 * w0 / sqrt(sqrt((G * G1))); wd = wn * sqrt(G1); } else { // use wd as mapping wd = 0.5 * w0; wn = wd * sqrt(G/G1); } Wnsq = wn * wn; Wdsq = wd * wd; // analog freqs of wn and wd wna = 2.0 * atan(wn); wda = 2.0 * atan(wd); // normalized analog gains at wna and wda temp = 1.0 / G1; gn = temp * analogShelf(w0, G, resonance, isHigh, wna); gd = temp * analogShelf(w0, G, resonance, isHigh, wda); gnsq = gn * gn; gdsq = gd * gd; // compute B, matching gains at wn and wd temp = 1.0 / (wn * wd); den = fabs(gnsq - gdsq); num = gnsq * (Wnsq - Wdsq) * (Wnsq - Wdsq) * (Wnsq + gdsq * Wdsq); B = temp * sqrt(num / den); // compute A, matching gains at wn and wd num = (Wnsq - Wdsq) * (Wnsq - Wdsq) * (Wnsq + gnsq * Wdsq); A = temp * sqrt(num / den); // design digital filter via bilinear transform b0 = G1 * (1.0 + B + Wnsq); b1 = G1 * 2.0 * (Wnsq - 1.0); b2 = G1 * (1.0 - B + Wnsq); a0 = 1.0 + A + Wdsq; a1 = 2.0 * (Wdsq - 1.0); a2 = 1.0 - A + Wdsq; // for unity gain, ensure poles/zeros in the right place. // needed for smooth interpolation when gain is changed. if (isUnity) { b0 = a0; b1 = a1; b2 = a2; } // invert filter for cut if (dbgain < 0.0) { temp = b0; b0 = a0; a0 = temp; temp = b1; b1 = a1; a1 = temp; temp = b2; b2 = a2; a2 = temp; } // normalize scale = 1.0 / a0; coef[0] = b0 * scale; coef[1] = b1 * scale; coef[2] = b2 * scale; coef[3] = a1 * scale; coef[4] = a2 * scale; } // // Biquad Lowpass/Highpass using analog matching. // Q = sqrt(0.5) = 2nd order Butterworth // void BQFilter(double coef[5], double w0, int isHigh) { double G1;