Autres/audacity
1
0
Fork 0
mirror of https://github.com/cookiengineer/audacity synced 2025-06-16 16:10:06 +02:00
Code Issues Releases Wiki Activity

Consistently use double within Biquad filter coefficient calcultations.

Browse Source 
This avoids catastrophic rounding errors which can even lead to filter
instability in some cases. This is a real problem in the coming IIR envelope
detector.
This commit is contained in:
Max Maisel 2018-08-06 17:44:40 +02:00
parent 86ae0460c9
commit e439e8f4c7
2 changed files with 48 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

69
src/effects/Biquad.cpp
View File

@ -64,20 +64,20 @@ ArrayOf<Biquad> Biquad::CalcButterworthFilter(int order, double fn, double fc, i
{ {
ArrayOf<Biquad> pBiquad(size_t((order+1) / 2), true); ArrayOf<Biquad> pBiquad(size_t((order+1) / 2), true);
// Set up the coefficients in all the biquads // Set up the coefficients in all the biquads
float fNorm = fc / fn; double fNorm = fc / fn;
if (fNorm >= 0.9999) if (fNorm >= 0.9999)
fNorm = 0.9999F; fNorm = 0.9999F;
float fC = tan (PI * fNorm / 2); double fC = tan (PI * fNorm / 2);
float fDCPoleDistSqr = 1.0F; double fDCPoleDistSqr = 1.0F;
float fZPoleX, fZPoleY; double fZPoleX, fZPoleY;
if ((order & 1) == 0) if ((order & 1) == 0)
{ {
// Even order // Even order
for (int iPair = 0; iPair < order/2; iPair++) for (int iPair = 0; iPair < order/2; iPair++)
{ {
float fSPoleX = fC * cos (PI - (iPair + 0.5) * PI / order); double fSPoleX = fC * cos (PI - (iPair + 0.5) * PI / order);
float fSPoleY = fC * sin (PI - (iPair + 0.5) * PI / order); double fSPoleY = fC * sin (PI - (iPair + 0.5) * PI / order);
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY); BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
pBiquad[iPair].fNumerCoeffs [B0] = 1; pBiquad[iPair].fNumerCoeffs [B0] = 1;
if (subtype == kLowPass) // LOWPASS if (subtype == kLowPass) // LOWPASS
@ -96,8 +96,8 @@ ArrayOf<Biquad> Biquad::CalcButterworthFilter(int order, double fn, double fc, i
else else
{ {
// Odd order - first do the 1st-order section // Odd order - first do the 1st-order section
float fSPoleX = -fC; double fSPoleX = -fC;
float fSPoleY = 0; double fSPoleY = 0;
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY); BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
pBiquad[0].fNumerCoeffs [B0] = 1; pBiquad[0].fNumerCoeffs [B0] = 1;
if (subtype == kLowPass) // LOWPASS if (subtype == kLowPass) // LOWPASS
@ -113,8 +113,8 @@ ArrayOf<Biquad> Biquad::CalcButterworthFilter(int order, double fn, double fc, i
fDCPoleDistSqr = fZPoleX + 1; // dist from Nyquist fDCPoleDistSqr = fZPoleX + 1; // dist from Nyquist
for (int iPair = 1; iPair <= order/2; iPair++) for (int iPair = 1; iPair <= order/2; iPair++)
{ {
float fSPoleX = fC * cos (PI - iPair * PI / order); double fSPoleX = fC * cos (PI - iPair * PI / order);
float fSPoleY = fC * sin (PI - iPair * PI / order); double fSPoleY = fC * sin (PI - iPair * PI / order);
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY); BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
pBiquad[iPair].fNumerCoeffs [B0] = 1; pBiquad[iPair].fNumerCoeffs [B0] = 1;
if (subtype == kLowPass) // LOWPASS if (subtype == kLowPass) // LOWPASS
@ -146,22 +146,22 @@ ArrayOf<Biquad> Biquad::CalcChebyshevType1Filter(int order, double fn, double fc
{ {
ArrayOf<Biquad> pBiquad(size_t((order+1) / 2), true); ArrayOf<Biquad> pBiquad(size_t((order+1) / 2), true);
// Set up the coefficients in all the biquads // Set up the coefficients in all the biquads
float fNorm = fc / fn; double fNorm = fc / fn;
if (fNorm >= 0.9999) if (fNorm >= 0.9999)
fNorm = 0.9999F; fNorm = 0.9999F;
float fC = tan (PI * fNorm / 2); double fC = tan (PI * fNorm / 2);
float fDCPoleDistSqr = 1.0F; double fDCPoleDistSqr = 1.0F;
float fZPoleX, fZPoleY; double fZPoleX, fZPoleY;
float fZZeroX; double fZZeroX;
float beta = cos (fNorm*PI); double beta = cos (fNorm*PI);
double eps; eps = sqrt (pow (10.0, wxMax(0.001, ripple) / 10.0) - 1); double eps; eps = sqrt (pow (10.0, wxMax(0.001, ripple) / 10.0) - 1);
double a; a = log (1 / eps + sqrt(1 / square(eps) + 1)) / order; double a; a = log (1 / eps + sqrt(1 / square(eps) + 1)) / order;
// Assume even order to start // Assume even order to start
for (int iPair = 0; iPair < order/2; iPair++) for (int iPair = 0; iPair < order/2; iPair++)
{ {
float fSPoleX = -fC * sinh (a) * sin ((2*iPair + 1) * PI / (2 * order)); double fSPoleX = -fC * sinh (a) * sin ((2*iPair + 1) * PI / (2 * order));
float fSPoleY = fC * cosh (a) * cos ((2*iPair + 1) * PI / (2 * order)); double fSPoleY = fC * cosh (a) * cos ((2*iPair + 1) * PI / (2 * order));
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY); BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
if (subtype == kLowPass) // LOWPASS if (subtype == kLowPass) // LOWPASS
{ {
@ -185,7 +185,7 @@ ArrayOf<Biquad> Biquad::CalcChebyshevType1Filter(int order, double fn, double fc
} }
if ((order & 1) == 0) if ((order & 1) == 0)
{ {
float fTemp = DB_TO_LINEAR(-wxMax(0.001, ripple)); // at DC the response is down R dB (for even-order) double fTemp = DB_TO_LINEAR(-wxMax(0.001, ripple)); // at DC the response is down R dB (for even-order)
pBiquad[0].fNumerCoeffs [B0] *= fTemp; pBiquad[0].fNumerCoeffs [B0] *= fTemp;
pBiquad[0].fNumerCoeffs [B1] *= fTemp; pBiquad[0].fNumerCoeffs [B1] *= fTemp;
pBiquad[0].fNumerCoeffs [B2] *= fTemp; pBiquad[0].fNumerCoeffs [B2] *= fTemp;
@ -193,8 +193,8 @@ ArrayOf<Biquad> Biquad::CalcChebyshevType1Filter(int order, double fn, double fc
else else
{ {
// Odd order - now do the 1st-order section // Odd order - now do the 1st-order section
float fSPoleX = -fC * sinh (a); double fSPoleX = -fC * sinh (a);
float fSPoleY = 0; double fSPoleY = 0;
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY); BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
if (subtype == kLowPass) // LOWPASS if (subtype == kLowPass) // LOWPASS
{ {
@ -228,17 +228,17 @@ ArrayOf<Biquad> Biquad::CalcChebyshevType2Filter(int order, double fn, double fc
{ {
ArrayOf<Biquad> pBiquad(size_t((order+1) / 2), true); ArrayOf<Biquad> pBiquad(size_t((order+1) / 2), true);
// Set up the coefficients in all the biquads // Set up the coefficients in all the biquads
float fNorm = fc / fn; double fNorm = fc / fn;
if (fNorm >= 0.9999) if (fNorm >= 0.9999)
fNorm = 0.9999F; fNorm = 0.9999F;
float fC = tan (PI * fNorm / 2); double fC = tan (PI * fNorm / 2);
float fDCPoleDistSqr = 1.0F; double fDCPoleDistSqr = 1.0F;
float fZPoleX, fZPoleY; double fZPoleX, fZPoleY;
float fZZeroX, fZZeroY; double fZZeroX, fZZeroY;
float beta = cos (fNorm*PI); double beta = cos (fNorm*PI);
float fSZeroX, fSZeroY; double fSZeroX, fSZeroY;
float fSPoleX, fSPoleY; double fSPoleX, fSPoleY;
double eps = DB_TO_LINEAR(-wxMax(0.001, ripple)); double eps = DB_TO_LINEAR(-wxMax(0.001, ripple));
double a = log (1 / eps + sqrt(1 / square(eps) + 1)) / order; double a = log (1 / eps + sqrt(1 / square(eps) + 1)) / order;
@ -304,22 +304,23 @@ ArrayOf<Biquad> Biquad::CalcChebyshevType2Filter(int order, double fn, double fc
return std::move(pBiquad); return std::move(pBiquad);
} }
void Biquad::ComplexDiv (float fNumerR, float fNumerI, float fDenomR, float fDenomI, float* pfQuotientR, float* pfQuotientI) void Biquad::ComplexDiv (double fNumerR, double fNumerI, double fDenomR, double fDenomI,
double* pfQuotientR, double* pfQuotientI)
{ {
float fDenom = square(fDenomR) + square(fDenomI); double fDenom = square(fDenomR) + square(fDenomI);
*pfQuotientR = (fNumerR * fDenomR + fNumerI * fDenomI) / fDenom; *pfQuotientR = (fNumerR * fDenomR + fNumerI * fDenomI) / fDenom;
*pfQuotientI = (fNumerI * fDenomR - fNumerR * fDenomI) / fDenom; *pfQuotientI = (fNumerI * fDenomR - fNumerR * fDenomI) / fDenom;
} }
bool Biquad::BilinTransform (float fSX, float fSY, float* pfZX, float* pfZY) bool Biquad::BilinTransform (double fSX, double fSY, double* pfZX, double* pfZY)
{ {
float fDenom = square (1 - fSX) + square (fSY); double fDenom = square (1 - fSX) + square (fSY);
*pfZX = (1 - square (fSX) - square (fSY)) / fDenom; *pfZX = (1 - square (fSX) - square (fSY)) / fDenom;
*pfZY = 2 * fSY / fDenom; *pfZY = 2 * fSY / fDenom;
return true; return true;
} }
float Biquad::Calc2D_DistSqr (float fX1, float fY1, float fX2, float fY2) float Biquad::Calc2D_DistSqr (double fX1, double fY1, double fX2, double fY2)
{ {
return square (fX1 - fX2) + square (fY1 - fY2); return square (fX1 - fX2) + square (fY1 - fY2);
} }

23
src/effects/Biquad.h
View File

@ -35,7 +35,9 @@ struct Biquad
inline float ProcessOne(float fIn) inline float ProcessOne(float fIn)
{ {
float fOut = fIn * fNumerCoeffs[B0] + // Biquad must use double for all calculations. Otherwise some
// filters may have catastrophic rounding errors!
double fOut = double(fIn) * fNumerCoeffs[B0] +
fPrevIn * fNumerCoeffs[B1] + fPrevIn * fNumerCoeffs[B1] +
fPrevPrevIn * fNumerCoeffs[B2] - fPrevPrevIn * fNumerCoeffs[B2] -
fPrevOut * fDenomCoeffs[A1] - fPrevOut * fDenomCoeffs[A1] -
@ -47,12 +49,12 @@ struct Biquad
return fOut; return fOut;
} }
float fNumerCoeffs[3]; // B0 B1 B2 double fNumerCoeffs[3]; // B0 B1 B2
float fDenomCoeffs[2]; // A1 A2, A0 == 1.0 double fDenomCoeffs[2]; // A1 A2, A0 == 1.0
float fPrevIn; double fPrevIn;
float fPrevPrevIn; double fPrevPrevIn;
float fPrevOut; double fPrevOut;
float fPrevPrevOut; double fPrevPrevOut;
enum kSubTypes enum kSubTypes
{ {
@ -65,9 +67,10 @@ struct Biquad
static ArrayOf<Biquad> CalcChebyshevType1Filter(int order, double fn, double fc, double ripple, int type); static ArrayOf<Biquad> CalcChebyshevType1Filter(int order, double fn, double fc, double ripple, int type);
static ArrayOf<Biquad> CalcChebyshevType2Filter(int order, double fn, double fc, double ripple, int type); static ArrayOf<Biquad> CalcChebyshevType2Filter(int order, double fn, double fc, double ripple, int type);
static void ComplexDiv (float fNumerR, float fNumerI, float fDenomR, float fDenomI, float* pfQuotientR, float* pfQuotientI); static void ComplexDiv (double fNumerR, double fNumerI, double fDenomR, double fDenomI,
static bool BilinTransform (float fSX, float fSY, float* pfZX, float* pfZY); double* pfQuotientR, double* pfQuotientI);
static float Calc2D_DistSqr (float fX1, float fY1, float fX2, float fY2); static bool BilinTransform (double fSX, double fSY, double* pfZX, double* pfZY);
static float Calc2D_DistSqr (double fX1, double fY1, double fX2, double fY2);
static const double s_fChebyCoeffs[MAX_Order][MAX_Order + 1]; static const double s_fChebyCoeffs[MAX_Order][MAX_Order + 1];
static double ChebyPoly(int Order, double NormFreq); static double ChebyPoly(int Order, double NormFreq);