Audacity 3.2.0
Biquad.cpp
Go to the documentation of this file.
1/**********************************************************************
2
3Audacity: A Digital Audio Editor
4
5Biquad.cpp
6
7Norm C
8Max 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
20Biquad::Biquad()
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
30void Biquad::Reset()
31{
32 fPrevIn = 0;
33 fPrevPrevIn = 0;
34 fPrevOut = 0;
35 fPrevPrevOut = 0;
36}
37
38void Biquad::Process(const float* pfIn, float* pfOut, int iNumSamples)
39{
40 for (int i = 0; i < iNumSamples; i++)
41 *pfOut++ = ProcessOne(*pfIn++);
42}
43
44const 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
64ArrayOf<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 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
146ArrayOf<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 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
228ArrayOf<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 pBiquad;
306}
307
308void 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
316bool 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
324float Biquad::Calc2D_DistSqr (double fX1, double fY1, double fX2, double fY2)
325{
326 return square (fX1 - fX2) + square (fY1 - fY2);
327}
328
329double 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}
square
#define square(a)
Definition: Biquad.cpp:17
PI
#define PI
Definition: Biquad.cpp:18
DB_TO_LINEAR
#define DB_TO_LINEAR(x)
Definition: MemoryX.h:337
fn
static const auto fn
Definition: WaveformView.cpp:1083
staffpad::audio::simd::sqrt
__finl float_x4 __vecc sqrt(const float_x4 &a)
Definition: SimdTypes_neon.h:144
Biquad::CalcButterworthFilter
static ArrayOf< Biquad > CalcButterworthFilter(int order, double fn, double fc, int type)
Definition: Biquad.cpp:64
Biquad::BilinTransform
static bool BilinTransform(double fSX, double fSY, double *pfZX, double *pfZY)
Definition: Biquad.cpp:316
Biquad::ProcessOne
float ProcessOne(float fIn)
Definition: Biquad.h:36
Biquad::CalcChebyshevType2Filter
static ArrayOf< Biquad > CalcChebyshevType2Filter(int order, double fn, double fc, double ripple, int type)
Definition: Biquad.cpp:228
Biquad::ChebyPoly
static double ChebyPoly(int Order, double NormFreq)
Definition: Biquad.cpp:329
Biquad::s_fChebyCoeffs
static const double s_fChebyCoeffs[MAX_Order][MAX_Order+1]
Definition: Biquad.h:75
Biquad::Biquad
Biquad()
Definition: Biquad.cpp:20
Biquad::Reset
void Reset()
Definition: Biquad.cpp:30
Biquad::fPrevOut
double fPrevOut
Definition: Biquad.h:56
Biquad::fPrevIn
double fPrevIn
Definition: Biquad.h:54
Biquad::A1
@ A1
Denominator coefficient indices.
Definition: Biquad.h:29
Biquad::A2
@ A2
Definition: Biquad.h:29
Biquad::MAX_Order
@ MAX_Order
Definition: Biquad.h:33
Biquad::B1
@ B1
Definition: Biquad.h:27
Biquad::MIN_Order
@ MIN_Order
Possible filter orders for the Calc...Filter(...) functions.
Definition: Biquad.h:32
Biquad::B2
@ B2
Definition: Biquad.h:27
Biquad::B0
@ B0
Numerator coefficient indices.
Definition: Biquad.h:27
Biquad::fPrevPrevIn
double fPrevPrevIn
Definition: Biquad.h:55
Biquad::CalcChebyshevType1Filter
static ArrayOf< Biquad > CalcChebyshevType1Filter(int order, double fn, double fc, double ripple, int type)
Definition: Biquad.cpp:146
Biquad::kLowPass
@ kLowPass
Definition: Biquad.h:61
Biquad::fNumerCoeffs
double fNumerCoeffs[3]
Definition: Biquad.h:52
Biquad::fDenomCoeffs
double fDenomCoeffs[2]
Definition: Biquad.h:53
Biquad::Calc2D_DistSqr
static float Calc2D_DistSqr(double fX1, double fY1, double fX2, double fY2)
Definition: Biquad.cpp:324
Biquad::ComplexDiv
static void ComplexDiv(double fNumerR, double fNumerI, double fDenomR, double fDenomI, double *pfQuotientR, double *pfQuotientI)
Definition: Biquad.cpp:308
Biquad::Process
void Process(const float *pfIn, float *pfOut, int iNumSamples)
Definition: Biquad.cpp:38
Biquad::fPrevPrevOut
double fPrevPrevOut
Definition: Biquad.h:57
Generated by doxygen 1.9.3