mirror of
https://github.com/cookiengineer/audacity
synced 2025-05-01 08:09:41 +02:00
Refactored filter calculation out of ScienFilter into Biquad.
This allows easy code reuse in other effects.
This commit is contained in:
Max Maisel
parent
6da48db127
commit
064ddb5a54
@ -10,8 +10,11 @@ Max Maisel
|
||||
***********************************************************************/
|
||||
|
||||
#include "Biquad.h"
|
||||
#include "Audacity.h"
|
||||
#include <cmath>
|
||||
|
||||
#define square(a) ((a)*(a))
|
||||
#define PI M_PI
|
||||
|
||||
Biquad::Biquad()
|
||||
{
|
||||
@ -39,6 +42,284 @@ void Biquad::Process(int iNumSamples)
|
||||
*pfOut++ = ProcessOne(*pfIn++);
|
||||
}
|
||||
|
||||
const double Biquad::s_fChebyCoeffs[MAX_Order][MAX_Order + 1] =
|
||||
{
|
||||
// For Chebyshev polynomials of the first kind (see http://en.wikipedia.org/wiki/Chebyshev_polynomial)
|
||||
// Coeffs are in the order 0, 1, 2...9
|
||||
{ 0, 1}, // order 1
|
||||
{-1, 0, 2}, // order 2 etc.
|
||||
{ 0, -3, 0, 4},
|
||||
{ 1, 0, -8, 0, 8},
|
||||
{ 0, 5, 0, -20, 0, 16},
|
||||
{-1, 0, 18, 0, -48, 0, 32},
|
||||
{ 0, -7, 0, 56, 0, -112, 0, 64},
|
||||
{ 1, 0, -32, 0, 160, 0, -256, 0, 128},
|
||||
{ 0, 9, 0, -120, 0, 432, 0, -576, 0, 256},
|
||||
{-1, 0, 50, 0, -400, 0, 1120, 0, -1280, 0, 512}
|
||||
};
|
||||
|
||||
// order: filter order
|
||||
// fn: nyquist frequency, i.e. half sample rate
|
||||
// fc: cutoff frequency
|
||||
// subtype: highpass or lowpass
|
||||
ArrayOf<Biquad> Biquad::CalcButterworthFilter(int order, double fn, double fc, int subtype)
|
||||
{
|
||||
ArrayOf<Biquad> pBiquad(size_t((order+1) / 2), true);
|
||||
// Set up the coefficients in all the biquads
|
||||
float fNorm = fc / fn;
|
||||
if (fNorm >= 0.9999)
|
||||
fNorm = 0.9999F;
|
||||
float fC = tan (PI * fNorm / 2);
|
||||
float fDCPoleDistSqr = 1.0F;
|
||||
float fZPoleX, fZPoleY;
|
||||
|
||||
if ((order & 1) == 0)
|
||||
{
|
||||
// Even order
|
||||
for (int iPair = 0; iPair < order/2; iPair++)
|
||||
{
|
||||
float fSPoleX = fC * cos (PI - (iPair + 0.5) * PI / order);
|
||||
float fSPoleY = fC * sin (PI - (iPair + 0.5) * PI / order);
|
||||
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
|
||||
pBiquad[iPair].fNumerCoeffs [B0] = 1;
|
||||
if (subtype == kLowPass) // LOWPASS
|
||||
pBiquad[iPair].fNumerCoeffs [B1] = 2;
|
||||
else
|
||||
pBiquad[iPair].fNumerCoeffs [B1] = -2;
|
||||
pBiquad[iPair].fNumerCoeffs [B2] = 1;
|
||||
pBiquad[iPair].fDenomCoeffs [A1] = -2 * fZPoleX;
|
||||
pBiquad[iPair].fDenomCoeffs [A2] = square(fZPoleX) + square(fZPoleY);
|
||||
if (subtype == kLowPass) // LOWPASS
|
||||
fDCPoleDistSqr *= Calc2D_DistSqr (1, 0, fZPoleX, fZPoleY);
|
||||
else
|
||||
fDCPoleDistSqr *= Calc2D_DistSqr (-1, 0, fZPoleX, fZPoleY); // distance from Nyquist
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Odd order - first do the 1st-order section
|
||||
float fSPoleX = -fC;
|
||||
float fSPoleY = 0;
|
||||
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
|
||||
pBiquad[0].fNumerCoeffs [B0] = 1;
|
||||
if (subtype == kLowPass) // LOWPASS
|
||||
pBiquad[0].fNumerCoeffs [B1] = 1;
|
||||
else
|
||||
pBiquad[0].fNumerCoeffs [B1] = -1;
|
||||
pBiquad[0].fNumerCoeffs [B2] = 0;
|
||||
pBiquad[0].fDenomCoeffs [A1] = -fZPoleX;
|
||||
pBiquad[0].fDenomCoeffs [A2] = 0;
|
||||
if (subtype == kLowPass) // LOWPASS
|
||||
fDCPoleDistSqr = 1 - fZPoleX;
|
||||
else
|
||||
fDCPoleDistSqr = fZPoleX + 1; // dist from Nyquist
|
||||
for (int iPair = 1; iPair <= order/2; iPair++)
|
||||
{
|
||||
float fSPoleX = fC * cos (PI - iPair * PI / order);
|
||||
float fSPoleY = fC * sin (PI - iPair * PI / order);
|
||||
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
|
||||
pBiquad[iPair].fNumerCoeffs [B0] = 1;
|
||||
if (subtype == kLowPass) // LOWPASS
|
||||
pBiquad[iPair].fNumerCoeffs [B1] = 2;
|
||||
else
|
||||
pBiquad[iPair].fNumerCoeffs [B1] = -2;
|
||||
pBiquad[iPair].fNumerCoeffs [B2] = 1;
|
||||
pBiquad[iPair].fDenomCoeffs [A1] = -2 * fZPoleX;
|
||||
pBiquad[iPair].fDenomCoeffs [A2] = square(fZPoleX) + square(fZPoleY);
|
||||
if (subtype == kLowPass) // LOWPASS
|
||||
fDCPoleDistSqr *= Calc2D_DistSqr (1, 0, fZPoleX, fZPoleY);
|
||||
else
|
||||
fDCPoleDistSqr *= Calc2D_DistSqr (-1, 0, fZPoleX, fZPoleY); // distance from Nyquist
|
||||
}
|
||||
}
|
||||
pBiquad[0].fNumerCoeffs [B0] *= fDCPoleDistSqr / (1 << order); // mult by DC dist from poles, divide by dist from zeroes
|
||||
pBiquad[0].fNumerCoeffs [B1] *= fDCPoleDistSqr / (1 << order);
|
||||
pBiquad[0].fNumerCoeffs [B2] *= fDCPoleDistSqr / (1 << order);
|
||||
|
||||
return std::move(pBiquad);
|
||||
}
|
||||
|
||||
// order: filter order
|
||||
// fn: nyquist frequency, i.e. half sample rate
|
||||
// fc: cutoff frequency
|
||||
// ripple: passband ripple in dB
|
||||
// subtype: highpass or lowpass
|
||||
ArrayOf<Biquad> Biquad::CalcChebyshevType1Filter(int order, double fn, double fc, double ripple, int subtype)
|
||||
{
|
||||
ArrayOf<Biquad> pBiquad(size_t((order+1) / 2), true);
|
||||
// Set up the coefficients in all the biquads
|
||||
float fNorm = fc / fn;
|
||||
if (fNorm >= 0.9999)
|
||||
fNorm = 0.9999F;
|
||||
float fC = tan (PI * fNorm / 2);
|
||||
float fDCPoleDistSqr = 1.0F;
|
||||
float fZPoleX, fZPoleY;
|
||||
float fZZeroX;
|
||||
float beta = cos (fNorm*PI);
|
||||
|
||||
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;
|
||||
// Assume even order to start
|
||||
for (int iPair = 0; iPair < order/2; iPair++)
|
||||
{
|
||||
float fSPoleX = -fC * sinh (a) * sin ((2*iPair + 1) * PI / (2 * order));
|
||||
float fSPoleY = fC * cosh (a) * cos ((2*iPair + 1) * PI / (2 * order));
|
||||
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
|
||||
if (subtype == kLowPass) // LOWPASS
|
||||
{
|
||||
fZZeroX = -1;
|
||||
fDCPoleDistSqr = Calc2D_DistSqr (1, 0, fZPoleX, fZPoleY);
|
||||
fDCPoleDistSqr /= 2*2; // dist from zero at Nyquist
|
||||
}
|
||||
else
|
||||
{
|
||||
// Highpass - do the digital LP->HP transform on the poles and zeroes
|
||||
ComplexDiv (beta - fZPoleX, -fZPoleY, 1 - beta * fZPoleX, -beta * fZPoleY, &fZPoleX, &fZPoleY);
|
||||
fZZeroX = 1;
|
||||
fDCPoleDistSqr = Calc2D_DistSqr (-1, 0, fZPoleX, fZPoleY); // distance from Nyquist
|
||||
fDCPoleDistSqr /= 2*2; // dist from zero at Nyquist
|
||||
}
|
||||
pBiquad[iPair].fNumerCoeffs [B0] = fDCPoleDistSqr;
|
||||
pBiquad[iPair].fNumerCoeffs [B1] = -2 * fZZeroX * fDCPoleDistSqr;
|
||||
pBiquad[iPair].fNumerCoeffs [B2] = fDCPoleDistSqr;
|
||||
pBiquad[iPair].fDenomCoeffs [A1] = -2 * fZPoleX;
|
||||
pBiquad[iPair].fDenomCoeffs [A2] = square(fZPoleX) + square(fZPoleY);
|
||||
}
|
||||
if ((order & 1) == 0)
|
||||
{
|
||||
float 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 [B1] *= fTemp;
|
||||
pBiquad[0].fNumerCoeffs [B2] *= fTemp;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Odd order - now do the 1st-order section
|
||||
float fSPoleX = -fC * sinh (a);
|
||||
float fSPoleY = 0;
|
||||
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
|
||||
if (subtype == kLowPass) // LOWPASS
|
||||
{
|
||||
fZZeroX = -1;
|
||||
fDCPoleDistSqr = sqrt(Calc2D_DistSqr (1, 0, fZPoleX, fZPoleY));
|
||||
fDCPoleDistSqr /= 2; // dist from zero at Nyquist
|
||||
}
|
||||
else
|
||||
{
|
||||
// Highpass - do the digital LP->HP transform on the poles and zeroes
|
||||
ComplexDiv (beta - fZPoleX, -fZPoleY, 1 - beta * fZPoleX, -beta * fZPoleY, &fZPoleX, &fZPoleY);
|
||||
fZZeroX = 1;
|
||||
fDCPoleDistSqr = sqrt(Calc2D_DistSqr (-1, 0, fZPoleX, fZPoleY)); // distance from Nyquist
|
||||
fDCPoleDistSqr /= 2; // dist from zero at Nyquist
|
||||
}
|
||||
pBiquad[(order-1)/2].fNumerCoeffs [B0] = fDCPoleDistSqr;
|
||||
pBiquad[(order-1)/2].fNumerCoeffs [B1] = -fZZeroX * fDCPoleDistSqr;
|
||||
pBiquad[(order-1)/2].fNumerCoeffs [B2] = 0;
|
||||
pBiquad[(order-1)/2].fDenomCoeffs [A1] = -fZPoleX;
|
||||
pBiquad[(order-1)/2].fDenomCoeffs [A2] = 0;
|
||||
}
|
||||
return std::move(pBiquad);
|
||||
}
|
||||
|
||||
// order: filter order
|
||||
// fn: nyquist frequency, i.e. half sample rate
|
||||
// fc: cutoff frequency
|
||||
// ripple: stopband ripple in dB
|
||||
// subtype: highpass or lowpass
|
||||
ArrayOf<Biquad> Biquad::CalcChebyshevType2Filter(int order, double fn, double fc, double ripple, int subtype)
|
||||
{
|
||||
ArrayOf<Biquad> pBiquad(size_t((order+1) / 2), true);
|
||||
// Set up the coefficients in all the biquads
|
||||
float fNorm = fc / fn;
|
||||
if (fNorm >= 0.9999)
|
||||
fNorm = 0.9999F;
|
||||
float fC = tan (PI * fNorm / 2);
|
||||
float fDCPoleDistSqr = 1.0F;
|
||||
float fZPoleX, fZPoleY;
|
||||
float fZZeroX, fZZeroY;
|
||||
float beta = cos (fNorm*PI);
|
||||
|
||||
float fSZeroX, fSZeroY;
|
||||
float fSPoleX, fSPoleY;
|
||||
double eps = DB_TO_LINEAR(-wxMax(0.001, ripple));
|
||||
double a = log (1 / eps + sqrt(1 / square(eps) + 1)) / order;
|
||||
|
||||
// Assume even order
|
||||
for (int iPair = 0; iPair < order/2; iPair++)
|
||||
{
|
||||
ComplexDiv (fC, 0, -sinh (a) * sin ((2*iPair + 1) * PI / (2 * order)),
|
||||
cosh (a) * cos ((2*iPair + 1) * PI / (2 * order)),
|
||||
&fSPoleX, &fSPoleY);
|
||||
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
|
||||
fSZeroX = 0;
|
||||
fSZeroY = fC / cos (((2 * iPair) + 1) * PI / (2 * order));
|
||||
BilinTransform (fSZeroX, fSZeroY, &fZZeroX, &fZZeroY);
|
||||
|
||||
if (subtype == kLowPass) // LOWPASS
|
||||
{
|
||||
fDCPoleDistSqr = Calc2D_DistSqr (1, 0, fZPoleX, fZPoleY);
|
||||
fDCPoleDistSqr /= Calc2D_DistSqr (1, 0, fZZeroX, fZZeroY);
|
||||
}
|
||||
else
|
||||
{
|
||||
// Highpass - do the digital LP->HP transform on the poles and zeroes
|
||||
ComplexDiv (beta - fZPoleX, -fZPoleY, 1 - beta * fZPoleX, -beta * fZPoleY, &fZPoleX, &fZPoleY);
|
||||
ComplexDiv (beta - fZZeroX, -fZZeroY, 1 - beta * fZZeroX, -beta * fZZeroY, &fZZeroX, &fZZeroY);
|
||||
fDCPoleDistSqr = Calc2D_DistSqr (-1, 0, fZPoleX, fZPoleY); // distance from Nyquist
|
||||
fDCPoleDistSqr /= Calc2D_DistSqr (-1, 0, fZZeroX, fZZeroY);
|
||||
}
|
||||
pBiquad[iPair].fNumerCoeffs [B0] = fDCPoleDistSqr;
|
||||
pBiquad[iPair].fNumerCoeffs [B1] = -2 * fZZeroX * fDCPoleDistSqr;
|
||||
pBiquad[iPair].fNumerCoeffs [B2] = (square(fZZeroX) + square(fZZeroY)) * fDCPoleDistSqr;
|
||||
pBiquad[iPair].fDenomCoeffs [A1] = -2 * fZPoleX;
|
||||
pBiquad[iPair].fDenomCoeffs [A2] = square(fZPoleX) + square(fZPoleY);
|
||||
}
|
||||
// Now, if it's odd order, we have one more to do
|
||||
if (order & 1)
|
||||
{
|
||||
int iPair = (order-1)/2; // we'll do it as a biquad, but it's just first-order
|
||||
ComplexDiv (fC, 0, -sinh (a) * sin ((2*iPair + 1) * PI / (2 * order)),
|
||||
cosh (a) * cos ((2*iPair + 1) * PI / (2 * order)),
|
||||
&fSPoleX, &fSPoleY);
|
||||
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
|
||||
fZZeroX = -1; // in the s-plane, the zero is at infinity
|
||||
fZZeroY = 0;
|
||||
if (subtype == kLowPass) // LOWPASS
|
||||
{
|
||||
fDCPoleDistSqr = sqrt(Calc2D_DistSqr (1, 0, fZPoleX, fZPoleY));
|
||||
fDCPoleDistSqr /= 2;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Highpass - do the digital LP->HP transform on the poles and zeroes
|
||||
ComplexDiv (beta - fZPoleX, -fZPoleY, 1 - beta * fZPoleX, -fZPoleY, &fZPoleX, &fZPoleY);
|
||||
fZZeroX = 1;
|
||||
fDCPoleDistSqr = sqrt(Calc2D_DistSqr (-1, 0, fZPoleX, fZPoleY)); // distance from Nyquist
|
||||
fDCPoleDistSqr /= 2;
|
||||
}
|
||||
pBiquad[iPair].fNumerCoeffs [B0] = fDCPoleDistSqr;
|
||||
pBiquad[iPair].fNumerCoeffs [B1] = -fZZeroX * fDCPoleDistSqr;
|
||||
pBiquad[iPair].fNumerCoeffs [B2] = 0;
|
||||
pBiquad[iPair].fDenomCoeffs [A1] = -fZPoleX;
|
||||
pBiquad[iPair].fDenomCoeffs [A2] = 0;
|
||||
}
|
||||
return std::move(pBiquad);
|
||||
}
|
||||
|
||||
double Biquad::ChebyPoly(int Order, double NormFreq) // NormFreq = 1 at the f0 point (where response is R dB down)
|
||||
{
|
||||
// Calc cosh (Order * acosh (NormFreq));
|
||||
double x = 1;
|
||||
double fSum = 0;
|
||||
wxASSERT (Order >= MIN_Order && Order <= MAX_Order);
|
||||
for (int i = 0; i <= Order; i++)
|
||||
{
|
||||
fSum += s_fChebyCoeffs [Order-1][i] * x;
|
||||
x *= NormFreq;
|
||||
}
|
||||
return fSum;
|
||||
}
|
||||
|
||||
void ComplexDiv (float fNumerR, float fNumerI, float fDenomR, float fDenomI, float* pfQuotientR, float* pfQuotientI)
|
||||
{
|
||||
float fDenom = square(fDenomR) + square(fDenomI);
|
||||
|
||||
@ -12,6 +12,7 @@ Max Maisel
|
||||
#ifndef __BIQUAD_H__
|
||||
#define __BIQUAD_H__
|
||||
|
||||
#include "MemoryX.h"
|
||||
|
||||
/// \brief Represents a biquad digital filter.
|
||||
struct Biquad
|
||||
@ -25,7 +26,11 @@ struct Biquad
|
||||
/// Numerator coefficient indices
|
||||
B0=0, B1, B2,
|
||||
/// Denominator coefficient indices
|
||||
A1=0, A2
|
||||
A1=0, A2,
|
||||
|
||||
/// Possible filter orders for the Calc...Filter(...) functions
|
||||
MIN_Order = 1,
|
||||
MAX_Order = 10
|
||||
};
|
||||
|
||||
inline float ProcessOne(float fIn)
|
||||
@ -50,6 +55,20 @@ struct Biquad
|
||||
float fPrevPrevIn;
|
||||
float fPrevOut;
|
||||
float fPrevPrevOut;
|
||||
|
||||
enum kSubTypes
|
||||
{
|
||||
kLowPass,
|
||||
kHighPass,
|
||||
nSubTypes
|
||||
};
|
||||
|
||||
static ArrayOf<Biquad> CalcButterworthFilter(int order, double fn, double fc, 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 const double s_fChebyCoeffs[MAX_Order][MAX_Order + 1];
|
||||
static double ChebyPoly(int Order, double NormFreq);
|
||||
};
|
||||
|
||||
void ComplexDiv (float fNumerR, float fNumerI, float fDenomR, float fDenomI, float* pfQuotientR, float* pfQuotientI);
|
||||
|
||||
@ -105,9 +105,9 @@ static const EnumValueSymbol kTypeStrings[nTypes] =
|
||||
|
||||
enum kSubTypes
|
||||
{
|
||||
kLowPass,
|
||||
kHighPass,
|
||||
nSubTypes
|
||||
kLowPass = Biquad::kLowPass,
|
||||
kHighPass = Biquad::kHighPass,
|
||||
nSubTypes = Biquad::nSubTypes
|
||||
};
|
||||
|
||||
static const EnumValueSymbol kSubTypeStrings[nSubTypes] =
|
||||
@ -127,22 +127,6 @@ Param( Cutoff, float, wxT("Cutoff"), 1000.0, 1.0, FLT_MAX
|
||||
Param( Passband, float, wxT("PassbandRipple"), 1.0, 0.0, 100.0, 1 );
|
||||
Param( Stopband, float, wxT("StopbandRipple"), 30.0, 0.0, 100.0, 1 );
|
||||
|
||||
static const double s_fChebyCoeffs[MAX_Order][MAX_Order + 1] =
|
||||
{
|
||||
// For Chebyshev polynomials of the first kind (see http://en.wikipedia.org/wiki/Chebyshev_polynomial)
|
||||
// Coeffs are in the order 0, 1, 2...9
|
||||
{ 0, 1}, // order 1
|
||||
{-1, 0, 2}, // order 2 etc.
|
||||
{ 0, -3, 0, 4},
|
||||
{ 1, 0, -8, 0, 8},
|
||||
{ 0, 5, 0, -20, 0, 16},
|
||||
{-1, 0, 18, 0, -48, 0, 32},
|
||||
{ 0, -7, 0, 56, 0, -112, 0, 64},
|
||||
{ 1, 0, -32, 0, 160, 0, -256, 0, 128},
|
||||
{ 0, 9, 0, -120, 0, 432, 0, -576, 0, 256},
|
||||
{-1, 0, 50, 0, -400, 0, 1120, 0, -1280, 0, 512}
|
||||
};
|
||||
|
||||
//----------------------------------------------------------------------------
|
||||
// EffectScienFilter
|
||||
//----------------------------------------------------------------------------
|
||||
@ -161,7 +145,6 @@ BEGIN_EVENT_TABLE(EffectScienFilter, wxEvtHandler)
|
||||
END_EVENT_TABLE()
|
||||
|
||||
EffectScienFilter::EffectScienFilter()
|
||||
: mpBiquad{ size_t( MAX_Order / 2 ), true }
|
||||
{
|
||||
mOrder = DEF_Order;
|
||||
mFilterType = DEF_Type;
|
||||
@ -621,230 +604,20 @@ bool EffectScienFilter::TransferGraphLimitsFromWindow()
|
||||
return true;
|
||||
}
|
||||
|
||||
bool EffectScienFilter::CalcFilter()
|
||||
void EffectScienFilter::CalcFilter()
|
||||
{
|
||||
// Set up the coefficients in all the biquads
|
||||
float fNorm = mCutoff / mNyquist;
|
||||
if (fNorm >= 0.9999)
|
||||
fNorm = 0.9999F;
|
||||
float fC = tan (PI * fNorm / 2);
|
||||
float fDCPoleDistSqr = 1.0F;
|
||||
float fZPoleX, fZPoleY;
|
||||
float fZZeroX, fZZeroY;
|
||||
float beta = cos (fNorm*PI);
|
||||
switch (mFilterType)
|
||||
{
|
||||
case kButterworth: // Butterworth
|
||||
if ((mOrder & 1) == 0)
|
||||
{
|
||||
// Even order
|
||||
for (int iPair = 0; iPair < mOrder/2; iPair++)
|
||||
{
|
||||
float fSPoleX = fC * cos (PI - (iPair + 0.5) * PI / mOrder);
|
||||
float fSPoleY = fC * sin (PI - (iPair + 0.5) * PI / mOrder);
|
||||
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
|
||||
mpBiquad[iPair].fNumerCoeffs [0] = 1;
|
||||
if (mFilterSubtype == kLowPass) // LOWPASS
|
||||
mpBiquad[iPair].fNumerCoeffs [1] = 2;
|
||||
else
|
||||
mpBiquad[iPair].fNumerCoeffs [1] = -2;
|
||||
mpBiquad[iPair].fNumerCoeffs [2] = 1;
|
||||
mpBiquad[iPair].fDenomCoeffs [0] = -2 * fZPoleX;
|
||||
mpBiquad[iPair].fDenomCoeffs [1] = square(fZPoleX) + square(fZPoleY);
|
||||
if (mFilterSubtype == kLowPass) // LOWPASS
|
||||
fDCPoleDistSqr *= Calc2D_DistSqr (1, 0, fZPoleX, fZPoleY);
|
||||
else
|
||||
fDCPoleDistSqr *= Calc2D_DistSqr (-1, 0, fZPoleX, fZPoleY); // distance from Nyquist
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
// Odd order - first do the 1st-order section
|
||||
float fSPoleX = -fC;
|
||||
float fSPoleY = 0;
|
||||
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
|
||||
mpBiquad[0].fNumerCoeffs [0] = 1;
|
||||
if (mFilterSubtype == kLowPass) // LOWPASS
|
||||
mpBiquad[0].fNumerCoeffs [1] = 1;
|
||||
else
|
||||
mpBiquad[0].fNumerCoeffs [1] = -1;
|
||||
mpBiquad[0].fNumerCoeffs [2] = 0;
|
||||
mpBiquad[0].fDenomCoeffs [0] = -fZPoleX;
|
||||
mpBiquad[0].fDenomCoeffs [1] = 0;
|
||||
if (mFilterSubtype == kLowPass) // LOWPASS
|
||||
fDCPoleDistSqr = 1 - fZPoleX;
|
||||
else
|
||||
fDCPoleDistSqr = fZPoleX + 1; // dist from Nyquist
|
||||
for (int iPair = 1; iPair <= mOrder/2; iPair++)
|
||||
{
|
||||
fSPoleX = fC * cos (PI - iPair * PI / mOrder);
|
||||
fSPoleY = fC * sin (PI - iPair * PI / mOrder);
|
||||
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
|
||||
mpBiquad[iPair].fNumerCoeffs [0] = 1;
|
||||
if (mFilterSubtype == kLowPass) // LOWPASS
|
||||
mpBiquad[iPair].fNumerCoeffs [1] = 2;
|
||||
else
|
||||
mpBiquad[iPair].fNumerCoeffs [1] = -2;
|
||||
mpBiquad[iPair].fNumerCoeffs [2] = 1;
|
||||
mpBiquad[iPair].fDenomCoeffs [0] = -2 * fZPoleX;
|
||||
mpBiquad[iPair].fDenomCoeffs [1] = square(fZPoleX) + square(fZPoleY);
|
||||
if (mFilterSubtype == kLowPass) // LOWPASS
|
||||
fDCPoleDistSqr *= Calc2D_DistSqr (1, 0, fZPoleX, fZPoleY);
|
||||
else
|
||||
fDCPoleDistSqr *= Calc2D_DistSqr (-1, 0, fZPoleX, fZPoleY); // distance from Nyquist
|
||||
}
|
||||
}
|
||||
mpBiquad[0].fNumerCoeffs [0] *= fDCPoleDistSqr / (1 << mOrder); // mult by DC dist from poles, divide by dist from zeroes
|
||||
mpBiquad[0].fNumerCoeffs [1] *= fDCPoleDistSqr / (1 << mOrder);
|
||||
mpBiquad[0].fNumerCoeffs [2] *= fDCPoleDistSqr / (1 << mOrder);
|
||||
case kButterworth:
|
||||
mpBiquad = Biquad::CalcButterworthFilter(mOrder, mNyquist, mCutoff, mFilterSubtype);
|
||||
break;
|
||||
|
||||
case kChebyshevTypeI: // Chebyshev Type 1
|
||||
double eps; eps = sqrt (pow (10.0, wxMax(0.001, mRipple) / 10.0) - 1);
|
||||
double a; a = log (1 / eps + sqrt(1 / square(eps) + 1)) / mOrder;
|
||||
// Assume even order to start
|
||||
for (int iPair = 0; iPair < mOrder/2; iPair++)
|
||||
{
|
||||
float fSPoleX = -fC * sinh (a) * sin ((2*iPair + 1) * PI / (2 * mOrder));
|
||||
float fSPoleY = fC * cosh (a) * cos ((2*iPair + 1) * PI / (2 * mOrder));
|
||||
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
|
||||
if (mFilterSubtype == kLowPass) // LOWPASS
|
||||
{
|
||||
fZZeroX = -1;
|
||||
fDCPoleDistSqr = Calc2D_DistSqr (1, 0, fZPoleX, fZPoleY);
|
||||
fDCPoleDistSqr /= 2*2; // dist from zero at Nyquist
|
||||
}
|
||||
else
|
||||
{
|
||||
// Highpass - do the digital LP->HP transform on the poles and zeroes
|
||||
ComplexDiv (beta - fZPoleX, -fZPoleY, 1 - beta * fZPoleX, -beta * fZPoleY, &fZPoleX, &fZPoleY);
|
||||
fZZeroX = 1;
|
||||
fDCPoleDistSqr = Calc2D_DistSqr (-1, 0, fZPoleX, fZPoleY); // distance from Nyquist
|
||||
fDCPoleDistSqr /= 2*2; // dist from zero at Nyquist
|
||||
}
|
||||
mpBiquad[iPair].fNumerCoeffs [0] = fDCPoleDistSqr;
|
||||
mpBiquad[iPair].fNumerCoeffs [1] = -2 * fZZeroX * fDCPoleDistSqr;
|
||||
mpBiquad[iPair].fNumerCoeffs [2] = fDCPoleDistSqr;
|
||||
mpBiquad[iPair].fDenomCoeffs [0] = -2 * fZPoleX;
|
||||
mpBiquad[iPair].fDenomCoeffs [1] = square(fZPoleX) + square(fZPoleY);
|
||||
}
|
||||
if ((mOrder & 1) == 0)
|
||||
{
|
||||
float fTemp = DB_TO_LINEAR(-wxMax(0.001, mRipple)); // at DC the response is down R dB (for even-order)
|
||||
mpBiquad[0].fNumerCoeffs [0] *= fTemp;
|
||||
mpBiquad[0].fNumerCoeffs [1] *= fTemp;
|
||||
mpBiquad[0].fNumerCoeffs [2] *= fTemp;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Odd order - now do the 1st-order section
|
||||
float fSPoleX = -fC * sinh (a);
|
||||
float fSPoleY = 0;
|
||||
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
|
||||
if (mFilterSubtype == kLowPass) // LOWPASS
|
||||
{
|
||||
fZZeroX = -1;
|
||||
fDCPoleDistSqr = sqrt(Calc2D_DistSqr (1, 0, fZPoleX, fZPoleY));
|
||||
fDCPoleDistSqr /= 2; // dist from zero at Nyquist
|
||||
}
|
||||
else
|
||||
{
|
||||
// Highpass - do the digital LP->HP transform on the poles and zeroes
|
||||
ComplexDiv (beta - fZPoleX, -fZPoleY, 1 - beta * fZPoleX, -beta * fZPoleY, &fZPoleX, &fZPoleY);
|
||||
fZZeroX = 1;
|
||||
fDCPoleDistSqr = sqrt(Calc2D_DistSqr (-1, 0, fZPoleX, fZPoleY)); // distance from Nyquist
|
||||
fDCPoleDistSqr /= 2; // dist from zero at Nyquist
|
||||
}
|
||||
mpBiquad[(mOrder-1)/2].fNumerCoeffs [0] = fDCPoleDistSqr;
|
||||
mpBiquad[(mOrder-1)/2].fNumerCoeffs [1] = -fZZeroX * fDCPoleDistSqr;
|
||||
mpBiquad[(mOrder-1)/2].fNumerCoeffs [2] = 0;
|
||||
mpBiquad[(mOrder-1)/2].fDenomCoeffs [0] = -fZPoleX;
|
||||
mpBiquad[(mOrder-1)/2].fDenomCoeffs [1] = 0;
|
||||
}
|
||||
case kChebyshevTypeI:
|
||||
mpBiquad = Biquad::CalcChebyshevType1Filter(mOrder, mNyquist, mCutoff, mRipple, mFilterSubtype);
|
||||
break;
|
||||
|
||||
case kChebyshevTypeII: // Chebyshev Type 2
|
||||
float fSZeroX, fSZeroY;
|
||||
float fSPoleX, fSPoleY;
|
||||
eps = DB_TO_LINEAR(-wxMax(0.001, mStopbandRipple));
|
||||
a = log (1 / eps + sqrt(1 / square(eps) + 1)) / mOrder;
|
||||
|
||||
// Assume even order
|
||||
for (int iPair = 0; iPair < mOrder/2; iPair++)
|
||||
{
|
||||
ComplexDiv (fC, 0, -sinh (a) * sin ((2*iPair + 1) * PI / (2 * mOrder)),
|
||||
cosh (a) * cos ((2*iPair + 1) * PI / (2 * mOrder)),
|
||||
&fSPoleX, &fSPoleY);
|
||||
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
|
||||
fSZeroX = 0;
|
||||
fSZeroY = fC / cos (((2 * iPair) + 1) * PI / (2 * mOrder));
|
||||
BilinTransform (fSZeroX, fSZeroY, &fZZeroX, &fZZeroY);
|
||||
|
||||
if (mFilterSubtype == kLowPass) // LOWPASS
|
||||
{
|
||||
fDCPoleDistSqr = Calc2D_DistSqr (1, 0, fZPoleX, fZPoleY);
|
||||
fDCPoleDistSqr /= Calc2D_DistSqr (1, 0, fZZeroX, fZZeroY);
|
||||
}
|
||||
else
|
||||
{
|
||||
// Highpass - do the digital LP->HP transform on the poles and zeroes
|
||||
ComplexDiv (beta - fZPoleX, -fZPoleY, 1 - beta * fZPoleX, -beta * fZPoleY, &fZPoleX, &fZPoleY);
|
||||
ComplexDiv (beta - fZZeroX, -fZZeroY, 1 - beta * fZZeroX, -beta * fZZeroY, &fZZeroX, &fZZeroY);
|
||||
fDCPoleDistSqr = Calc2D_DistSqr (-1, 0, fZPoleX, fZPoleY); // distance from Nyquist
|
||||
fDCPoleDistSqr /= Calc2D_DistSqr (-1, 0, fZZeroX, fZZeroY);
|
||||
}
|
||||
mpBiquad[iPair].fNumerCoeffs [0] = fDCPoleDistSqr;
|
||||
mpBiquad[iPair].fNumerCoeffs [1] = -2 * fZZeroX * fDCPoleDistSqr;
|
||||
mpBiquad[iPair].fNumerCoeffs [2] = (square(fZZeroX) + square(fZZeroY)) * fDCPoleDistSqr;
|
||||
mpBiquad[iPair].fDenomCoeffs [0] = -2 * fZPoleX;
|
||||
mpBiquad[iPair].fDenomCoeffs [1] = square(fZPoleX) + square(fZPoleY);
|
||||
}
|
||||
// Now, if it's odd order, we have one more to do
|
||||
if (mOrder & 1)
|
||||
{
|
||||
int iPair = (mOrder-1)/2; // we'll do it as a biquad, but it's just first-order
|
||||
ComplexDiv (fC, 0, -sinh (a) * sin ((2*iPair + 1) * PI / (2 * mOrder)),
|
||||
cosh (a) * cos ((2*iPair + 1) * PI / (2 * mOrder)),
|
||||
&fSPoleX, &fSPoleY);
|
||||
BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
|
||||
fZZeroX = -1; // in the s-plane, the zero is at infinity
|
||||
fZZeroY = 0;
|
||||
if (mFilterSubtype == kLowPass) // LOWPASS
|
||||
{
|
||||
fDCPoleDistSqr = sqrt(Calc2D_DistSqr (1, 0, fZPoleX, fZPoleY));
|
||||
fDCPoleDistSqr /= 2;
|
||||
}
|
||||
else
|
||||
{
|
||||
// Highpass - do the digital LP->HP transform on the poles and zeroes
|
||||
ComplexDiv (beta - fZPoleX, -fZPoleY, 1 - beta * fZPoleX, -fZPoleY, &fZPoleX, &fZPoleY);
|
||||
fZZeroX = 1;
|
||||
fDCPoleDistSqr = sqrt(Calc2D_DistSqr (-1, 0, fZPoleX, fZPoleY)); // distance from Nyquist
|
||||
fDCPoleDistSqr /= 2;
|
||||
}
|
||||
mpBiquad[iPair].fNumerCoeffs [0] = fDCPoleDistSqr;
|
||||
mpBiquad[iPair].fNumerCoeffs [1] = -fZZeroX * fDCPoleDistSqr;
|
||||
mpBiquad[iPair].fNumerCoeffs [2] = 0;
|
||||
mpBiquad[iPair].fDenomCoeffs [0] = -fZPoleX;
|
||||
mpBiquad[iPair].fDenomCoeffs [1] = 0;
|
||||
}
|
||||
case kChebyshevTypeII:
|
||||
mpBiquad = Biquad::CalcChebyshevType2Filter(mOrder, mNyquist, mCutoff, mStopbandRipple, mFilterSubtype);
|
||||
break;
|
||||
}
|
||||
return true;
|
||||
}
|
||||
|
||||
double EffectScienFilter::ChebyPoly(int Order, double NormFreq) // NormFreq = 1 at the f0 point (where response is R dB down)
|
||||
{
|
||||
// Calc cosh (Order * acosh (NormFreq));
|
||||
double x = 1;
|
||||
double fSum = 0;
|
||||
wxASSERT (Order >= MIN_Order && Order <= MAX_Order);
|
||||
for (int i = 0; i <= Order; i++)
|
||||
{
|
||||
fSum += s_fChebyCoeffs [Order-1][i] * x;
|
||||
x *= NormFreq;
|
||||
}
|
||||
return fSum;
|
||||
}
|
||||
|
||||
float EffectScienFilter::FilterMagnAtFreq(float Freq)
|
||||
@ -882,14 +655,17 @@ float EffectScienFilter::FilterMagnAtFreq(float Freq)
|
||||
|
||||
case kChebyshevTypeI: // Chebyshev Type 1
|
||||
double eps; eps = sqrt(pow (10.0, wxMax(0.001, mRipple)/10.0) - 1);
|
||||
double chebyPolyVal;
|
||||
switch (mFilterSubtype)
|
||||
{
|
||||
case 0: // lowpass
|
||||
default:
|
||||
Magn = sqrt (1 / (1 + square(eps) * square(ChebyPoly(mOrder, FreqWarped/CutoffWarped))));
|
||||
chebyPolyVal = Biquad::ChebyPoly(mOrder, FreqWarped/CutoffWarped);
|
||||
Magn = sqrt (1 / (1 + square(eps) * square(chebyPolyVal)));
|
||||
break;
|
||||
case 1:
|
||||
Magn = sqrt (1 / (1 + square(eps) * square(ChebyPoly(mOrder, CutoffWarped/FreqWarped))));
|
||||
chebyPolyVal = Biquad::ChebyPoly(mOrder, CutoffWarped/FreqWarped);
|
||||
Magn = sqrt (1 / (1 + square(eps) * square(chebyPolyVal)));
|
||||
break;
|
||||
}
|
||||
break;
|
||||
@ -900,10 +676,12 @@ float EffectScienFilter::FilterMagnAtFreq(float Freq)
|
||||
{
|
||||
case kLowPass: // lowpass
|
||||
default:
|
||||
Magn = sqrt (1 / (1 + 1 / (square(eps) * square(ChebyPoly(mOrder, CutoffWarped/FreqWarped)))));
|
||||
chebyPolyVal = Biquad::ChebyPoly(mOrder, CutoffWarped/FreqWarped);
|
||||
Magn = sqrt (1 / (1 + 1 / (square(eps) * square(chebyPolyVal))));
|
||||
break;
|
||||
case kHighPass:
|
||||
Magn = sqrt (1 / (1 + 1 / (square(eps) * square(ChebyPoly(mOrder, FreqWarped/CutoffWarped)))));
|
||||
chebyPolyVal = Biquad::ChebyPoly(mOrder, FreqWarped/CutoffWarped);
|
||||
Magn = sqrt (1 / (1 + 1 / (square(eps) * square(chebyPolyVal))));
|
||||
break;
|
||||
}
|
||||
break;
|
||||
|
||||
@ -69,12 +69,8 @@ private:
|
||||
// EffectScienFilter implementation
|
||||
|
||||
bool TransferGraphLimitsFromWindow();
|
||||
bool CalcFilter();
|
||||
double ChebyPoly (int Order, double NormFreq);
|
||||
void CalcFilter();
|
||||
float FilterMagnAtFreq(float Freq);
|
||||
|
||||
bool CalcFilterCoeffs (void);
|
||||
|
||||
void EnableDisableRippleCtl (int FilterType);
|
||||
|
||||
void OnSize( wxSizeEvent & evt );
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user