Audacity  3.0.3
Biquad.cpp
Go to the documentation of this file.
1 /**********************************************************************
2 
3 Audacity: A Digital Audio Editor
4 
5 Biquad.cpp
6 
7 Norm C
8 Max Maisel
9 
10 ***********************************************************************/
11 
12 #include "Biquad.h"
13 
14 #include <cmath>
15 #include <wx/utils.h>
16 
17 #define square(a) ((a)*(a))
18 #define PI M_PI
19 
21 {
22  fNumerCoeffs[B0] = 1;
23  fNumerCoeffs[B1] = 0;
24  fNumerCoeffs[B2] = 0;
25  fDenomCoeffs[A1] = 0;
26  fDenomCoeffs[A2] = 0;
27  Reset();
28 }
29 
31 {
32  fPrevIn = 0;
33  fPrevPrevIn = 0;
34  fPrevOut = 0;
35  fPrevPrevOut = 0;
36 }
37 
38 void Biquad::Process(float* pfIn, float* pfOut, int iNumSamples)
39 {
40  for (int i = 0; i < iNumSamples; i++)
41  *pfOut++ = ProcessOne(*pfIn++);
42 }
43 
44 const double Biquad::s_fChebyCoeffs[MAX_Order][MAX_Order + 1] =
45 {
46  // For Chebyshev polynomials of the first kind (see http://en.wikipedia.org/wiki/Chebyshev_polynomial)
47  // Coeffs are in the order 0, 1, 2...9
48  { 0, 1}, // order 1
49  {-1, 0, 2}, // order 2 etc.
50  { 0, -3, 0, 4},
51  { 1, 0, -8, 0, 8},
52  { 0, 5, 0, -20, 0, 16},
53  {-1, 0, 18, 0, -48, 0, 32},
54  { 0, -7, 0, 56, 0, -112, 0, 64},
55  { 1, 0, -32, 0, 160, 0, -256, 0, 128},
56  { 0, 9, 0, -120, 0, 432, 0, -576, 0, 256},
57  {-1, 0, 50, 0, -400, 0, 1120, 0, -1280, 0, 512}
58 };
59 
60 // order: filter order
61 // fn: nyquist frequency, i.e. half sample rate
62 // fc: cutoff frequency
63 // subtype: highpass or lowpass
64 ArrayOf<Biquad> Biquad::CalcButterworthFilter(int order, double fn, double fc, int subtype)
65 {
66  ArrayOf<Biquad> pBiquad(size_t((order+1) / 2), true);
67  // Set up the coefficients in all the biquads
68  double fNorm = fc / fn;
69  if (fNorm >= 0.9999)
70  fNorm = 0.9999F;
71  double fC = tan (PI * fNorm / 2);
72  double fDCPoleDistSqr = 1.0F;
73  double fZPoleX, fZPoleY;
74 
75  if ((order & 1) == 0)
76  {
77  // Even order
78  for (int iPair = 0; iPair < order/2; iPair++)
79  {
80  double fSPoleX = fC * cos (PI - (iPair + 0.5) * PI / order);
81  double fSPoleY = fC * sin (PI - (iPair + 0.5) * PI / order);
82  BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
83  pBiquad[iPair].fNumerCoeffs [B0] = 1;
84  if (subtype == kLowPass) // LOWPASS
85  pBiquad[iPair].fNumerCoeffs [B1] = 2;
86  else
87  pBiquad[iPair].fNumerCoeffs [B1] = -2;
88  pBiquad[iPair].fNumerCoeffs [B2] = 1;
89  pBiquad[iPair].fDenomCoeffs [A1] = -2 * fZPoleX;
90  pBiquad[iPair].fDenomCoeffs [A2] = square(fZPoleX) + square(fZPoleY);
91  if (subtype == kLowPass) // LOWPASS
92  fDCPoleDistSqr *= Calc2D_DistSqr (1, 0, fZPoleX, fZPoleY);
93  else
94  fDCPoleDistSqr *= Calc2D_DistSqr (-1, 0, fZPoleX, fZPoleY); // distance from Nyquist
95  }
96  }
97  else
98  {
99  // Odd order - first do the 1st-order section
100  double fSPoleX = -fC;
101  double fSPoleY = 0;
102  BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
103  pBiquad[0].fNumerCoeffs [B0] = 1;
104  if (subtype == kLowPass) // LOWPASS
105  pBiquad[0].fNumerCoeffs [B1] = 1;
106  else
107  pBiquad[0].fNumerCoeffs [B1] = -1;
108  pBiquad[0].fNumerCoeffs [B2] = 0;
109  pBiquad[0].fDenomCoeffs [A1] = -fZPoleX;
110  pBiquad[0].fDenomCoeffs [A2] = 0;
111  if (subtype == kLowPass) // LOWPASS
112  fDCPoleDistSqr = 1 - fZPoleX;
113  else
114  fDCPoleDistSqr = fZPoleX + 1; // dist from Nyquist
115  for (int iPair = 1; iPair <= order/2; iPair++)
116  {
117  double fSPoleX = fC * cos (PI - iPair * PI / order);
118  double fSPoleY = fC * sin (PI - iPair * PI / order);
119  BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
120  pBiquad[iPair].fNumerCoeffs [B0] = 1;
121  if (subtype == kLowPass) // LOWPASS
122  pBiquad[iPair].fNumerCoeffs [B1] = 2;
123  else
124  pBiquad[iPair].fNumerCoeffs [B1] = -2;
125  pBiquad[iPair].fNumerCoeffs [B2] = 1;
126  pBiquad[iPair].fDenomCoeffs [A1] = -2 * fZPoleX;
127  pBiquad[iPair].fDenomCoeffs [A2] = square(fZPoleX) + square(fZPoleY);
128  if (subtype == kLowPass) // LOWPASS
129  fDCPoleDistSqr *= Calc2D_DistSqr (1, 0, fZPoleX, fZPoleY);
130  else
131  fDCPoleDistSqr *= Calc2D_DistSqr (-1, 0, fZPoleX, fZPoleY); // distance from Nyquist
132  }
133  }
134  pBiquad[0].fNumerCoeffs [B0] *= fDCPoleDistSqr / (1 << order); // mult by DC dist from poles, divide by dist from zeroes
135  pBiquad[0].fNumerCoeffs [B1] *= fDCPoleDistSqr / (1 << order);
136  pBiquad[0].fNumerCoeffs [B2] *= fDCPoleDistSqr / (1 << order);
137 
138  return std::move(pBiquad);
139 }
140 
141 // order: filter order
142 // fn: nyquist frequency, i.e. half sample rate
143 // fc: cutoff frequency
144 // ripple: passband ripple in dB
145 // subtype: highpass or lowpass
146 ArrayOf<Biquad> Biquad::CalcChebyshevType1Filter(int order, double fn, double fc, double ripple, int subtype)
147 {
148  ArrayOf<Biquad> pBiquad(size_t((order+1) / 2), true);
149  // Set up the coefficients in all the biquads
150  double fNorm = fc / fn;
151  if (fNorm >= 0.9999)
152  fNorm = 0.9999F;
153  double fC = tan (PI * fNorm / 2);
154  double fDCPoleDistSqr = 1.0F;
155  double fZPoleX, fZPoleY;
156  double fZZeroX;
157  double beta = cos (fNorm*PI);
158 
159  double eps; eps = sqrt (pow (10.0, wxMax(0.001, ripple) / 10.0) - 1);
160  double a; a = log (1 / eps + sqrt(1 / square(eps) + 1)) / order;
161  // Assume even order to start
162  for (int iPair = 0; iPair < order/2; iPair++)
163  {
164  double fSPoleX = -fC * sinh (a) * sin ((2*iPair + 1) * PI / (2 * order));
165  double fSPoleY = fC * cosh (a) * cos ((2*iPair + 1) * PI / (2 * order));
166  BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
167  if (subtype == kLowPass) // LOWPASS
168  {
169  fZZeroX = -1;
170  fDCPoleDistSqr = Calc2D_DistSqr (1, 0, fZPoleX, fZPoleY);
171  fDCPoleDistSqr /= 2*2; // dist from zero at Nyquist
172  }
173  else
174  {
175  // Highpass - do the digital LP->HP transform on the poles and zeroes
176  ComplexDiv (beta - fZPoleX, -fZPoleY, 1 - beta * fZPoleX, -beta * fZPoleY, &fZPoleX, &fZPoleY);
177  fZZeroX = 1;
178  fDCPoleDistSqr = Calc2D_DistSqr (-1, 0, fZPoleX, fZPoleY); // distance from Nyquist
179  fDCPoleDistSqr /= 2*2; // dist from zero at Nyquist
180  }
181  pBiquad[iPair].fNumerCoeffs [B0] = fDCPoleDistSqr;
182  pBiquad[iPair].fNumerCoeffs [B1] = -2 * fZZeroX * fDCPoleDistSqr;
183  pBiquad[iPair].fNumerCoeffs [B2] = fDCPoleDistSqr;
184  pBiquad[iPair].fDenomCoeffs [A1] = -2 * fZPoleX;
185  pBiquad[iPair].fDenomCoeffs [A2] = square(fZPoleX) + square(fZPoleY);
186  }
187  if ((order & 1) == 0)
188  {
189  double fTemp = DB_TO_LINEAR(-wxMax(0.001, ripple)); // at DC the response is down R dB (for even-order)
190  pBiquad[0].fNumerCoeffs [B0] *= fTemp;
191  pBiquad[0].fNumerCoeffs [B1] *= fTemp;
192  pBiquad[0].fNumerCoeffs [B2] *= fTemp;
193  }
194  else
195  {
196  // Odd order - now do the 1st-order section
197  double fSPoleX = -fC * sinh (a);
198  double fSPoleY = 0;
199  BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
200  if (subtype == kLowPass) // LOWPASS
201  {
202  fZZeroX = -1;
203  fDCPoleDistSqr = sqrt(Calc2D_DistSqr (1, 0, fZPoleX, fZPoleY));
204  fDCPoleDistSqr /= 2; // dist from zero at Nyquist
205  }
206  else
207  {
208  // Highpass - do the digital LP->HP transform on the poles and zeroes
209  ComplexDiv (beta - fZPoleX, -fZPoleY, 1 - beta * fZPoleX, -beta * fZPoleY, &fZPoleX, &fZPoleY);
210  fZZeroX = 1;
211  fDCPoleDistSqr = sqrt(Calc2D_DistSqr (-1, 0, fZPoleX, fZPoleY)); // distance from Nyquist
212  fDCPoleDistSqr /= 2; // dist from zero at Nyquist
213  }
214  pBiquad[(order-1)/2].fNumerCoeffs [B0] = fDCPoleDistSqr;
215  pBiquad[(order-1)/2].fNumerCoeffs [B1] = -fZZeroX * fDCPoleDistSqr;
216  pBiquad[(order-1)/2].fNumerCoeffs [B2] = 0;
217  pBiquad[(order-1)/2].fDenomCoeffs [A1] = -fZPoleX;
218  pBiquad[(order-1)/2].fDenomCoeffs [A2] = 0;
219  }
220  return std::move(pBiquad);
221 }
222 
223 // order: filter order
224 // fn: nyquist frequency, i.e. half sample rate
225 // fc: cutoff frequency
226 // ripple: stopband ripple in dB
227 // subtype: highpass or lowpass
228 ArrayOf<Biquad> Biquad::CalcChebyshevType2Filter(int order, double fn, double fc, double ripple, int subtype)
229 {
230  ArrayOf<Biquad> pBiquad(size_t((order+1) / 2), true);
231  // Set up the coefficients in all the biquads
232  double fNorm = fc / fn;
233  if (fNorm >= 0.9999)
234  fNorm = 0.9999F;
235  double fC = tan (PI * fNorm / 2);
236  double fDCPoleDistSqr = 1.0F;
237  double fZPoleX, fZPoleY;
238  double fZZeroX, fZZeroY;
239  double beta = cos (fNorm*PI);
240 
241  double fSZeroX, fSZeroY;
242  double fSPoleX, fSPoleY;
243  double eps = DB_TO_LINEAR(-wxMax(0.001, ripple));
244  double a = log (1 / eps + sqrt(1 / square(eps) + 1)) / order;
245 
246  // Assume even order
247  for (int iPair = 0; iPair < order/2; iPair++)
248  {
249  ComplexDiv (fC, 0, -sinh (a) * sin ((2*iPair + 1) * PI / (2 * order)),
250  cosh (a) * cos ((2*iPair + 1) * PI / (2 * order)),
251  &fSPoleX, &fSPoleY);
252  BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
253  fSZeroX = 0;
254  fSZeroY = fC / cos (((2 * iPair) + 1) * PI / (2 * order));
255  BilinTransform (fSZeroX, fSZeroY, &fZZeroX, &fZZeroY);
256 
257  if (subtype == kLowPass) // LOWPASS
258  {
259  fDCPoleDistSqr = Calc2D_DistSqr (1, 0, fZPoleX, fZPoleY);
260  fDCPoleDistSqr /= Calc2D_DistSqr (1, 0, fZZeroX, fZZeroY);
261  }
262  else
263  {
264  // Highpass - do the digital LP->HP transform on the poles and zeroes
265  ComplexDiv (beta - fZPoleX, -fZPoleY, 1 - beta * fZPoleX, -beta * fZPoleY, &fZPoleX, &fZPoleY);
266  ComplexDiv (beta - fZZeroX, -fZZeroY, 1 - beta * fZZeroX, -beta * fZZeroY, &fZZeroX, &fZZeroY);
267  fDCPoleDistSqr = Calc2D_DistSqr (-1, 0, fZPoleX, fZPoleY); // distance from Nyquist
268  fDCPoleDistSqr /= Calc2D_DistSqr (-1, 0, fZZeroX, fZZeroY);
269  }
270  pBiquad[iPair].fNumerCoeffs [B0] = fDCPoleDistSqr;
271  pBiquad[iPair].fNumerCoeffs [B1] = -2 * fZZeroX * fDCPoleDistSqr;
272  pBiquad[iPair].fNumerCoeffs [B2] = (square(fZZeroX) + square(fZZeroY)) * fDCPoleDistSqr;
273  pBiquad[iPair].fDenomCoeffs [A1] = -2 * fZPoleX;
274  pBiquad[iPair].fDenomCoeffs [A2] = square(fZPoleX) + square(fZPoleY);
275  }
276  // Now, if it's odd order, we have one more to do
277  if (order & 1)
278  {
279  int iPair = (order-1)/2; // we'll do it as a biquad, but it's just first-order
280  ComplexDiv (fC, 0, -sinh (a) * sin ((2*iPair + 1) * PI / (2 * order)),
281  cosh (a) * cos ((2*iPair + 1) * PI / (2 * order)),
282  &fSPoleX, &fSPoleY);
283  BilinTransform (fSPoleX, fSPoleY, &fZPoleX, &fZPoleY);
284  fZZeroX = -1; // in the s-plane, the zero is at infinity
285  fZZeroY = 0;
286  if (subtype == kLowPass) // LOWPASS
287  {
288  fDCPoleDistSqr = sqrt(Calc2D_DistSqr (1, 0, fZPoleX, fZPoleY));
289  fDCPoleDistSqr /= 2;
290  }
291  else
292  {
293  // Highpass - do the digital LP->HP transform on the poles and zeroes
294  ComplexDiv (beta - fZPoleX, -fZPoleY, 1 - beta * fZPoleX, -fZPoleY, &fZPoleX, &fZPoleY);
295  fZZeroX = 1;
296  fDCPoleDistSqr = sqrt(Calc2D_DistSqr (-1, 0, fZPoleX, fZPoleY)); // distance from Nyquist
297  fDCPoleDistSqr /= 2;
298  }
299  pBiquad[iPair].fNumerCoeffs [B0] = fDCPoleDistSqr;
300  pBiquad[iPair].fNumerCoeffs [B1] = -fZZeroX * fDCPoleDistSqr;
301  pBiquad[iPair].fNumerCoeffs [B2] = 0;
302  pBiquad[iPair].fDenomCoeffs [A1] = -fZPoleX;
303  pBiquad[iPair].fDenomCoeffs [A2] = 0;
304  }
305  return std::move(pBiquad);
306 }
307 
308 void Biquad::ComplexDiv (double fNumerR, double fNumerI, double fDenomR, double fDenomI,
309  double* pfQuotientR, double* pfQuotientI)
310 {
311  double fDenom = square(fDenomR) + square(fDenomI);
312  *pfQuotientR = (fNumerR * fDenomR + fNumerI * fDenomI) / fDenom;
313  *pfQuotientI = (fNumerI * fDenomR - fNumerR * fDenomI) / fDenom;
314 }
315 
316 bool Biquad::BilinTransform (double fSX, double fSY, double* pfZX, double* pfZY)
317 {
318  double fDenom = square (1 - fSX) + square (fSY);
319  *pfZX = (1 - square (fSX) - square (fSY)) / fDenom;
320  *pfZY = 2 * fSY / fDenom;
321  return true;
322 }
323 
324 float Biquad::Calc2D_DistSqr (double fX1, double fY1, double fX2, double fY2)
325 {
326  return square (fX1 - fX2) + square (fY1 - fY2);
327 }
328 
329 double Biquad::ChebyPoly(int Order, double NormFreq) // NormFreq = 1 at the f0 point (where response is R dB down)
330 {
331  // Calc cosh (Order * acosh (NormFreq));
332  double x = 1;
333  double fSum = 0;
334  wxASSERT (Order >= MIN_Order && Order <= MAX_Order);
335  for (int i = 0; i <= Order; i++)
336  {
337  fSum += s_fChebyCoeffs [Order-1][i] * x;
338  x *= NormFreq;
339  }
340  return fSum;
341 }
Biquad::fPrevPrevOut
double fPrevPrevOut
Definition: Biquad.h:57
DB_TO_LINEAR
const double MIN_Threshold_Linear DB_TO_LINEAR(MIN_Threshold_dB)
fn
static const auto fn
Definition: WaveformView.cpp:1108
Biquad.h
Biquad::Reset
void Reset()
Definition: Biquad.cpp:30
Biquad::ComplexDiv
static void ComplexDiv(double fNumerR, double fNumerI, double fDenomR, double fDenomI, double *pfQuotientR, double *pfQuotientI)
Definition: Biquad.cpp:308
Biquad::fNumerCoeffs
double fNumerCoeffs[3]
Definition: Biquad.h:52
Biquad::MAX_Order
@ MAX_Order
Definition: Biquad.h:33
Biquad::fPrevPrevIn
double fPrevPrevIn
Definition: Biquad.h:55
Biquad::Calc2D_DistSqr
static float Calc2D_DistSqr(double fX1, double fY1, double fX2, double fY2)
Definition: Biquad.cpp:324
Biquad::fPrevIn
double fPrevIn
Definition: Biquad.h:54
PI
#define PI
Definition: Biquad.cpp:18
Biquad::MIN_Order
@ MIN_Order
Possible filter orders for the Calc...Filter(...) functions.
Definition: Biquad.h:32
Biquad::B1
@ B1
Definition: Biquad.h:27
Biquad::A2
@ A2
Definition: Biquad.h:29
Biquad::s_fChebyCoeffs
static const double s_fChebyCoeffs[MAX_Order][MAX_Order+1]
Definition: Biquad.h:75
Biquad::A1
@ A1
Denominator coefficient indices.
Definition: Biquad.h:29
Biquad::CalcChebyshevType1Filter
static ArrayOf< Biquad > CalcChebyshevType1Filter(int order, double fn, double fc, double ripple, int type)
Definition: Biquad.cpp:146
Biquad::fPrevOut
double fPrevOut
Definition: Biquad.h:56
Biquad::B0
@ B0
Numerator coefficient indices.
Definition: Biquad.h:27
Biquad::B2
@ B2
Definition: Biquad.h:27
Biquad::Process
void Process(float *pfIn, float *pfOut, int iNumSamples)
Definition: Biquad.cpp:38
Biquad::Biquad
Biquad()
Definition: Biquad.cpp:20
Biquad::kLowPass
@ kLowPass
Definition: Biquad.h:61
Biquad::BilinTransform
static bool BilinTransform(double fSX, double fSY, double *pfZX, double *pfZY)
Definition: Biquad.cpp:316
square
#define square(a)
Definition: Biquad.cpp:17
Biquad::ProcessOne
float ProcessOne(float fIn)
Definition: Biquad.h:36
Biquad::ChebyPoly
static double ChebyPoly(int Order, double NormFreq)
Definition: Biquad.cpp:329
Biquad::CalcButterworthFilter
static ArrayOf< Biquad > CalcButterworthFilter(int order, double fn, double fc, int type)
Definition: Biquad.cpp:64
Biquad::CalcChebyshevType2Filter
static ArrayOf< Biquad > CalcChebyshevType2Filter(int order, double fn, double fc, double ripple, int type)
Definition: Biquad.cpp:228
ArrayOf< Biquad >
Biquad::fDenomCoeffs
double fDenomCoeffs[2]
Definition: Biquad.h:53