OpenShot Library | OpenShotAudio  0.2.2
juce_FilterDesign.cpp
1 /*
2  ==============================================================================
3 
4  This file is part of the JUCE library.
5  Copyright (c) 2017 - ROLI Ltd.
6 
7  JUCE is an open source library subject to commercial or open-source
8  licensing.
9 
10  By using JUCE, you agree to the terms of both the JUCE 5 End-User License
11  Agreement and JUCE 5 Privacy Policy (both updated and effective as of the
12  27th April 2017).
13 
14  End User License Agreement: www.juce.com/juce-5-licence
15  Privacy Policy: www.juce.com/juce-5-privacy-policy
16 
17  Or: You may also use this code under the terms of the GPL v3 (see
18  www.gnu.org/licenses).
19 
20  JUCE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY, AND ALL WARRANTIES, WHETHER
21  EXPRESSED OR IMPLIED, INCLUDING MERCHANTABILITY AND FITNESS FOR PURPOSE, ARE
22  DISCLAIMED.
23 
24  ==============================================================================
25 */
26 
27 namespace juce
28 {
29 namespace dsp
30 {
31 
32 template <typename FloatType>
35  double sampleRate, size_t order,
36  WindowingMethod type,
37  FloatType beta)
38 {
39  jassert (sampleRate > 0);
40  jassert (frequency > 0 && frequency <= sampleRate * 0.5);
41 
42  auto* result = new typename FIR::Coefficients<FloatType> (order + 1u);
43 
44  auto* c = result->getRawCoefficients();
45  auto normalisedFrequency = frequency / sampleRate;
46 
47  for (size_t i = 0; i <= order; ++i)
48  {
49  if (i == order * 0.5)
50  {
51  c[i] = static_cast<FloatType> (normalisedFrequency * 2);
52  }
53  else
54  {
55  auto indice = MathConstants<double>::pi * (static_cast<double> (i) - 0.5 * static_cast<double> (order));
56  c[i] = static_cast<FloatType> (std::sin (2.0 * indice * normalisedFrequency) / indice);
57  }
58  }
59 
60  WindowingFunction<FloatType> theWindow (order + 1, type, false, beta);
61  theWindow.multiplyWithWindowingTable (c, order + 1);
62 
63  return *result;
64 }
65 
66 template <typename FloatType>
68  FilterDesign<FloatType>::designFIRLowpassKaiserMethod (FloatType frequency, double sampleRate,
69  FloatType normalisedTransitionWidth,
70  FloatType amplitudedB)
71 {
72  jassert (sampleRate > 0);
73  jassert (frequency > 0 && frequency <= sampleRate * 0.5);
74  jassert (normalisedTransitionWidth > 0 && normalisedTransitionWidth <= 0.5);
75  jassert (amplitudedB >= -100 && amplitudedB <= 0);
76 
77  FloatType beta = 0;
78 
79  if (amplitudedB < -50)
80  beta = static_cast<FloatType> (0.1102 * (-amplitudedB - 8.7));
81  else if (amplitudedB <= -21)
82  beta = static_cast<FloatType> (0.5842 * std::pow (-amplitudedB - 21, 0.4) + 0.07886 * (-amplitudedB - 21));
83 
84  int order = amplitudedB < -21 ? roundToInt (std::ceil ((-amplitudedB - 7.95) / (2.285 * normalisedTransitionWidth * MathConstants<double>::twoPi)))
85  : roundToInt (std::ceil (5.79 / (normalisedTransitionWidth * MathConstants<double>::twoPi)));
86 
87  jassert (order >= 0);
88 
89  return designFIRLowpassWindowMethod (frequency, sampleRate, static_cast<size_t> (order),
91 }
92 
93 
94 template <typename FloatType>
96  FilterDesign<FloatType>::designFIRLowpassTransitionMethod (FloatType frequency, double sampleRate, size_t order,
97  FloatType normalisedTransitionWidth, FloatType spline)
98 {
99  jassert (sampleRate > 0);
100  jassert (frequency > 0 && frequency <= sampleRate * 0.5);
101  jassert (normalisedTransitionWidth > 0 && normalisedTransitionWidth <= 0.5);
102  jassert (spline >= 1.0 && spline <= 4.0);
103 
104  auto normalisedFrequency = frequency / static_cast<FloatType> (sampleRate);
105 
106  auto* result = new typename FIR::Coefficients<FloatType> (order + 1u);
107  auto* c = result->getRawCoefficients();
108 
109  for (size_t i = 0; i <= order; ++i)
110  {
111  if (i == order / 2 && order % 2 == 0)
112  {
113  c[i] = static_cast<FloatType> (2 * normalisedFrequency);
114  }
115  else
116  {
117  auto indice = MathConstants<double>::pi * (i - 0.5 * order);
118  auto indice2 = MathConstants<double>::pi * normalisedTransitionWidth * (i - 0.5 * order) / spline;
119  c[i] = static_cast<FloatType> (std::sin (2 * indice * normalisedFrequency)
120  / indice * std::pow (std::sin (indice2) / indice2, spline));
121  }
122  }
123 
124  return *result;
125 }
126 
127 template <typename FloatType>
130  double sampleRate, size_t order,
131  FloatType normalisedTransitionWidth,
132  FloatType stopBandWeight)
133 {
134  jassert (sampleRate > 0);
135  jassert (frequency > 0 && frequency <= sampleRate * 0.5);
136  jassert (normalisedTransitionWidth > 0 && normalisedTransitionWidth <= 0.5);
137  jassert (stopBandWeight >= 1.0 && stopBandWeight <= 100.0);
138 
139  auto normalisedFrequency = static_cast<double> (frequency) / sampleRate;
140 
141  auto wp = MathConstants<double>::twoPi * (static_cast<double> (normalisedFrequency - normalisedTransitionWidth / 2.0));
142  auto ws = MathConstants<double>::twoPi * (static_cast<double> (normalisedFrequency + normalisedTransitionWidth / 2.0));
143 
144  auto N = order + 1;
145 
146  auto* result = new typename FIR::Coefficients<FloatType> (static_cast<size_t> (N));
147  auto* c = result->getRawCoefficients();
148 
149  if (N % 2 == 1)
150  {
151  // Type I
152  auto M = (N - 1) / 2;
153 
154  Matrix<double> b (M + 1, 1),
155  q (2 * M + 1, 1);
156 
157  auto sinc = [](double x) { return x == 0 ? 1 : std::sin (x * MathConstants<double>::pi)
158  / (MathConstants<double>::pi * x); };
159 
160  auto factorp = wp / MathConstants<double>::pi;
161  auto factors = ws / MathConstants<double>::pi;
162 
163  for (size_t i = 0; i <= M; ++i)
164  b (i, 0) = factorp * sinc (factorp * i);
165 
166  q (0, 0) = factorp + stopBandWeight * (1.0 - factors);
167 
168  for (size_t i = 1; i <= 2 * M; ++i)
169  q (i, 0) = factorp * sinc (factorp * i) - stopBandWeight * factors * sinc (factors * i);
170 
171  auto Q1 = Matrix<double>::toeplitz (q, M + 1);
172  auto Q2 = Matrix<double>::hankel (q, M + 1, 0);
173 
174  Q1 += Q2; Q1 *= 0.5;
175 
176  Q1.solve (b);
177 
178  c[M] = static_cast<FloatType> (b (0, 0));
179 
180  for (size_t i = 1; i <= M; ++i)
181  {
182  c[M - i] = static_cast<FloatType> (b (i, 0) * 0.5);
183  c[M + i] = static_cast<FloatType> (b (i, 0) * 0.5);
184  }
185  }
186  else
187  {
188  // Type II
189  auto M = N / 2;
190 
191  Matrix<double> b (M, 1);
192  Matrix<double> qp (2 * M, 1);
193  Matrix<double> qs (2 * M, 1);
194 
195  auto sinc = [](double x) { return x == 0 ? 1 : std::sin (x * MathConstants<double>::pi)
196  / (MathConstants<double>::pi * x); };
197 
198  auto factorp = wp / MathConstants<double>::pi;
199  auto factors = ws / MathConstants<double>::pi;
200 
201  for (size_t i = 0; i < M; ++i)
202  b (i, 0) = factorp * sinc (factorp * (i + 0.5));
203 
204  for (size_t i = 0; i < 2 * M; ++i)
205  {
206  qp (i, 0) = 0.25 * factorp * sinc (factorp * i);
207  qs (i, 0) = -0.25 * stopBandWeight * factors * sinc (factors * i);
208  }
209 
210  auto Q1p = Matrix<double>::toeplitz (qp, M);
211  auto Q2p = Matrix<double>::hankel (qp, M, 1);
212  auto Q1s = Matrix<double>::toeplitz (qs, M);
213  auto Q2s = Matrix<double>::hankel (qs, M, 1);
214 
215  auto Id = Matrix<double>::identity (M);
216  Id *= (0.25 * stopBandWeight);
217 
218  Q1p += Q2p;
219  Q1s += Q2s;
220  Q1s += Id;
221 
222  auto& Q = Q1s;
223  Q += Q1p;
224 
225  Q.solve (b);
226 
227  for (size_t i = 0; i < M; ++i)
228  {
229  c[M - i - 1] = static_cast<FloatType> (b (i, 0) * 0.25);
230  c[M + i] = static_cast<FloatType> (b (i, 0) * 0.25);
231  }
232  }
233 
234  return *result;
235 }
236 
237 template <typename FloatType>
240  FloatType amplitudedB)
241 {
242  jassert (normalisedTransitionWidth > 0 && normalisedTransitionWidth <= 0.5);
243  jassert (amplitudedB >= -300 && amplitudedB <= -10);
244 
245  auto wpT = (0.5 - normalisedTransitionWidth) * MathConstants<double>::pi;
246 
247  auto n = roundToInt (std::ceil ((amplitudedB - 18.18840664 * wpT + 33.64775300) / (18.54155181 * wpT - 29.13196871)));
248  auto kp = (n * wpT - 1.57111377 * n + 0.00665857) / (-1.01927560 * n + 0.37221484);
249  auto A = (0.01525753 * n + 0.03682344 + 9.24760314 / (double) n) * kp + 1.01701407 + 0.73512298 / (double) n;
250  auto B = (0.00233667 * n - 1.35418408 + 5.75145813 / (double) n) * kp + 1.02999650 - 0.72759508 / (double) n;
251 
254 
255  auto diff = (hn.size() - hnm.size()) / 2;
256 
257  for (int i = 0; i < diff; ++i)
258  {
259  hnm.add (0.0);
260  hnm.insert (0, 0.0);
261  }
262 
263  auto hh = hn;
264 
265  for (int i = 0; i < hn.size(); ++i)
266  hh.setUnchecked (i, A * hh[i] + B * hnm[i]);
267 
268  auto* result = new typename FIR::Coefficients<FloatType> (static_cast<size_t> (hh.size()));
269  auto* c = result->getRawCoefficients();
270 
271  for (int i = 0; i < hh.size(); ++i)
272  c[i] = (float) hh[i];
273 
274  double NN;
275 
276  if (n % 2 == 0)
277  {
278  NN = 2.0 * result->getMagnitudeForFrequency (0.5, 1.0);
279  }
280  else
281  {
282  auto w01 = std::sqrt (kp * kp + (1 - kp * kp) * std::pow (std::cos (MathConstants<double>::pi / (2.0 * n + 1.0)), 2.0));
283  auto om01 = std::acos (-w01);
284 
285  NN = -2.0 * result->getMagnitudeForFrequency (om01 / MathConstants<double>::twoPi, 1.0);
286  }
287 
288  for (int i = 0; i < hh.size(); ++i)
289  c[i] = static_cast<FloatType> ((A * hn[i] + B * hnm[i]) / NN);
290 
291  c[2 * n + 1] = static_cast<FloatType> (0.5);
292 
293  return *result;
294 }
295 
296 template <typename FloatType>
298 {
299  Array<double> alpha;
300  alpha.resize (2 * n + 1);
301 
302  alpha.setUnchecked (2 * n, 1.0 / std::pow (1.0 - kp * kp, n));
303 
304  if (n > 0)
305  alpha.setUnchecked (2 * n - 2, -(2 * n * kp * kp + 1) * alpha[2 * n]);
306 
307  if (n > 1)
308  alpha.setUnchecked (2 * n - 4, -(4 * n + 1 + (n - 1) * (2 * n - 1) * kp * kp) / (2.0 * n) * alpha[2 * n - 2]
309  - (2 * n + 1) * ((n + 1) * kp * kp + 1) / (2.0 * n) * alpha[2 * n]);
310 
311  for (int k = n; k >= 3; --k)
312  {
313  auto c1 = (3 * (n*(n + 2) - k * (k - 2)) + 2 * k - 3 + 2 * (k - 2)*(2 * k - 3) * kp * kp) * alpha[2 * k - 4];
314  auto c2 = (3 * (n*(n + 2) - (k - 1) * (k + 1)) + 2 * (2 * k - 1) + 2 * k*(2 * k - 1) * kp * kp) * alpha[2 * k - 2];
315  auto c3 = (n * (n + 2) - (k - 1) * (k + 1)) * alpha[2 * k];
316  auto c4 = (n * (n + 2) - (k - 3) * (k - 1));
317 
318  alpha.setUnchecked (2 * k - 6, -(c1 + c2 + c3) / c4);
319  }
320 
321  Array<double> ai;
322  ai.resize (2 * n + 1 + 1);
323 
324  for (int k = 0; k <= n; ++k)
325  ai.setUnchecked (2 * k + 1, alpha[2 * k] / (2.0 * k + 1.0));
326 
327  Array<double> hn;
328  hn.resize (2 * n + 1 + 2 * n + 1 + 1);
329 
330  for (int k = 0; k <= n; ++k)
331  {
332  hn.setUnchecked (2 * n + 1 + (2 * k + 1), 0.5 * ai[2 * k + 1]);
333  hn.setUnchecked (2 * n + 1 - (2 * k + 1), 0.5 * ai[2 * k + 1]);
334  }
335 
336  return hn;
337 }
338 
339 template <typename FloatType>
342  FloatType normalisedTransitionWidth,
343  FloatType passbandAmplitudedB,
344  FloatType stopbandAmplitudedB)
345 {
346  return designIIRLowpassHighOrderGeneralMethod (0, frequency, sampleRate, normalisedTransitionWidth,
347  passbandAmplitudedB, stopbandAmplitudedB);
348 }
349 
350 template <typename FloatType>
353  FloatType normalisedTransitionWidth,
354  FloatType passbandAmplitudedB,
355  FloatType stopbandAmplitudedB)
356 {
357  return designIIRLowpassHighOrderGeneralMethod (1, frequency, sampleRate, normalisedTransitionWidth,
358  passbandAmplitudedB, stopbandAmplitudedB);
359 }
360 
361 template <typename FloatType>
364  FloatType normalisedTransitionWidth,
365  FloatType passbandAmplitudedB,
366  FloatType stopbandAmplitudedB)
367 {
368  return designIIRLowpassHighOrderGeneralMethod (2, frequency, sampleRate, normalisedTransitionWidth,
369  passbandAmplitudedB, stopbandAmplitudedB);
370 }
371 
372 template <typename FloatType>
375  FloatType normalisedTransitionWidth,
376  FloatType passbandAmplitudedB,
377  FloatType stopbandAmplitudedB)
378 {
379  return designIIRLowpassHighOrderGeneralMethod (3, frequency, sampleRate, normalisedTransitionWidth,
380  passbandAmplitudedB, stopbandAmplitudedB);
381 }
382 
383 template <typename FloatType>
385  FilterDesign<FloatType>::designIIRLowpassHighOrderGeneralMethod (int type, FloatType frequency, double sampleRate,
386  FloatType normalisedTransitionWidth,
387  FloatType passbandAmplitudedB,
388  FloatType stopbandAmplitudedB)
389 {
390  jassert (sampleRate > 0);
391  jassert (frequency > 0 && frequency <= sampleRate * 0.5);
392  jassert (normalisedTransitionWidth > 0 && normalisedTransitionWidth <= 0.5);
393  jassert (passbandAmplitudedB > -20 && passbandAmplitudedB < 0);
394  jassert (stopbandAmplitudedB > -300 && stopbandAmplitudedB < -20);
395 
396  auto normalisedFrequency = frequency / sampleRate;
397 
398  auto fp = normalisedFrequency - normalisedTransitionWidth / 2;
399  auto fs = normalisedFrequency + normalisedTransitionWidth / 2;
400 
401  double Ap = passbandAmplitudedB;
402  double As = stopbandAmplitudedB;
403  auto Gp = Decibels::decibelsToGain (Ap, -300.0);
404  auto Gs = Decibels::decibelsToGain (As, -300.0);
405  auto epsp = std::sqrt (1.0 / (Gp * Gp) - 1.0);
406  auto epss = std::sqrt (1.0 / (Gs * Gs) - 1.0);
407 
408  auto omegap = std::tan (MathConstants<double>::pi * fp);
409  auto omegas = std::tan (MathConstants<double>::pi * fs);
410  constexpr auto halfPi = MathConstants<double>::halfPi;
411 
412  auto k = omegap / omegas;
413  auto k1 = epsp / epss;
414 
415  int N;
416 
417  if (type == 0)
418  {
419  N = roundToInt (std::ceil (std::log (1.0 / k1) / std::log (1.0 / k)));
420  }
421  else if (type == 1 || type == 2)
422  {
423  N = roundToInt (std::ceil (std::acosh (1.0 / k1) / std::acosh (1.0 / k)));
424  }
425  else
426  {
427  double K, Kp, K1, K1p;
428 
431 
432  N = roundToInt (std::ceil ((K1p * K) / (K1 * Kp)));
433  }
434 
435  const int r = N % 2;
436  const int L = (N - r) / 2;
437  const double H0 = (type == 1 || type == 3) ? std::pow (Gp, 1.0 - r) : 1.0;
438 
439  Array<Complex<double>> pa, za;
440  Complex<double> j (0, 1);
441 
442  if (type == 0)
443  {
444  if (r == 1)
445  pa.add (-omegap * std::pow (epsp, -1.0 / (double) N));
446 
447  for (int i = 1; i <= L; ++i)
448  {
449  auto ui = (2 * i - 1.0) / (double) N;
450  pa.add (omegap * std::pow (epsp, -1.0 / (double) N) * j * exp (ui * halfPi * j));
451  }
452  }
453  else if (type == 1)
454  {
455  auto v0 = std::asinh (1.0 / epsp) / (N * halfPi);
456 
457  if (r == 1)
458  pa.add (-omegap * std::sinh (v0 * halfPi));
459 
460  for (int i = 1; i <= L; ++i)
461  {
462  auto ui = (2 * i - 1.0) / (double) N;
463  pa.add (omegap * j * std::cos ((ui - j * v0) * halfPi));
464  }
465  }
466  else if (type == 2)
467  {
468  auto v0 = std::asinh (epss) / (N * halfPi);
469 
470  if (r == 1)
471  pa.add(-1.0 / (k / omegap * std::sinh (v0 * halfPi)));
472 
473  for (int i = 1; i <= L; ++i)
474  {
475  auto ui = (2 * i - 1.0) / (double) N;
476 
477  pa.add (1.0 / (k / omegap * j * std::cos ((ui - j * v0) * halfPi)));
478  za.add (1.0 / (k / omegap * j * std::cos (ui * halfPi)));
479  }
480  }
481  else
482  {
483  auto v0 = -j * (SpecialFunctions::asne (j / epsp, k1) / (double) N);
484 
485  if (r == 1)
486  pa.add (omegap * j * SpecialFunctions::sne (j * v0, k));
487 
488  for (int i = 1; i <= L; ++i)
489  {
490  auto ui = (2 * i - 1.0) / (double) N;
491  auto zetai = SpecialFunctions::cde (ui, k);
492 
493  pa.add (omegap * j * SpecialFunctions::cde (ui - j * v0, k));
494  za.add (omegap * j / (k * zetai));
495  }
496  }
497 
498  Array<Complex<double>> p, z, g;
499 
500  if (r == 1)
501  {
502  p.add ((1.0 + pa[0]) / (1.0 - pa[0]));
503  g.add (0.5 * (1.0 - p[0]));
504  }
505 
506  for (int i = 0; i < L; ++i)
507  {
508  p.add ((1.0 + pa[i + r]) / (1.0 - pa[i + r]));
509  z.add (za.size() == 0 ? -1.0 : (1.0 + za[i]) / (1.0 - za[i]));
510  g.add ((1.0 - p[i + r]) / (1.0 - z[i]));
511  }
512 
514 
515  if (r == 1)
516  {
517  auto b0 = static_cast<FloatType> (H0 * std::real (g[0]));
518  auto b1 = b0;
519  auto a1 = static_cast<FloatType> (-std::real (p[0]));
520 
521  cascadedCoefficients.add (new IIR::Coefficients<FloatType> (b0, b1, 1.0f, a1));
522  }
523 
524  for (int i = 0; i < L; ++i)
525  {
526  auto gain = std::pow (std::abs (g[i + r]), 2.0);
527 
528  auto b0 = static_cast<FloatType> (gain);
529  auto b1 = static_cast<FloatType> (std::real (-z[i] - std::conj (z[i])) * gain);
530  auto b2 = static_cast<FloatType> (std::real ( z[i] * std::conj (z[i])) * gain);
531 
532  auto a1 = static_cast<FloatType> (std::real (-p[i+r] - std::conj (p[i + r])));
533  auto a2 = static_cast<FloatType> (std::real ( p[i+r] * std::conj (p[i + r])));
534 
535  cascadedCoefficients.add (new IIR::Coefficients<FloatType> (b0, b1, b2, 1, a1, a2));
536  }
537 
538  return cascadedCoefficients;
539 }
540 
541 template <typename FloatType>
544  double sampleRate, int order)
545 {
546  jassert (sampleRate > 0);
547  jassert (frequency > 0 && frequency <= sampleRate * 0.5);
548  jassert (order > 0);
549 
551 
552  if (order % 2 == 1)
553  {
554  arrayFilters.add (*IIR::Coefficients<FloatType>::makeFirstOrderLowPass (sampleRate, frequency));
555 
556  for (auto i = 0; i < order / 2; ++i)
557  {
558  auto Q = 1.0 / (2.0 * std::cos ((i + 1.0) * MathConstants<double>::pi / order));
559  arrayFilters.add (*IIR::Coefficients<FloatType>::makeLowPass (sampleRate, frequency,
560  static_cast<FloatType> (Q)));
561  }
562  }
563  else
564  {
565  for (auto i = 0; i < order / 2; ++i)
566  {
567  auto Q = 1.0 / (2.0 * std::cos ((2.0 * i + 1.0) * MathConstants<double>::pi / (order * 2.0)));
568  arrayFilters.add (*IIR::Coefficients<FloatType>::makeLowPass (sampleRate, frequency,
569  static_cast<FloatType> (Q)));
570  }
571  }
572 
573  return arrayFilters;
574 }
575 
576 template <typename FloatType>
579  double sampleRate, int order)
580 {
581  jassert (sampleRate > 0);
582  jassert (frequency > 0 && frequency <= sampleRate * 0.5);
583  jassert (order > 0);
584 
586 
587  if (order % 2 == 1)
588  {
589  arrayFilters.add (*IIR::Coefficients<FloatType>::makeFirstOrderHighPass (sampleRate, frequency));
590 
591  for (auto i = 0; i < order / 2; ++i)
592  {
593  auto Q = 1.0 / (2.0 * std::cos ((i + 1.0) * MathConstants<double>::pi / order));
594  arrayFilters.add (*IIR::Coefficients<FloatType>::makeHighPass (sampleRate, frequency,
595  static_cast<FloatType> (Q)));
596  }
597  }
598  else
599  {
600  for (auto i = 0; i < order / 2; ++i)
601  {
602  auto Q = 1.0 / (2.0 * std::cos ((2.0 * i + 1.0) * MathConstants<double>::pi / (order * 2.0)));
603  arrayFilters.add (*IIR::Coefficients<FloatType>::makeHighPass (sampleRate, frequency,
604  static_cast<FloatType> (Q)));
605  }
606  }
607 
608  return arrayFilters;
609 }
610 
611 template <typename FloatType>
614  FloatType stopbandAmplitudedB)
615 {
616  jassert (normalisedTransitionWidth > 0 && normalisedTransitionWidth <= 0.5);
617  jassert (stopbandAmplitudedB > -300 && stopbandAmplitudedB < -10);
618 
619  const double wt = MathConstants<double>::twoPi * normalisedTransitionWidth;
620  const double ds = Decibels::decibelsToGain (stopbandAmplitudedB, static_cast<FloatType> (-300.0));
621 
622  auto k = std::pow (std::tan ((MathConstants<double>::pi - wt) / 4), 2.0);
623  auto kp = std::sqrt (1.0 - k * k);
624  auto e = (1 - std::sqrt (kp)) / (1 + std::sqrt (kp)) * 0.5;
625  auto q = e + 2 * std::pow (e, 5.0) + 15 * std::pow (e, 9.0) + 150 * std::pow (e, 13.0);
626 
627  auto k1 = ds * ds / (1 - ds * ds);
628  int n = roundToInt (std::ceil (std::log (k1 * k1 / 16) / std::log (q)));
629 
630  if (n % 2 == 0)
631  ++n;
632 
633  if (n == 1)
634  n = 3;
635 
636  auto q1 = std::pow (q, (double) n);
637  k1 = 4 * std::sqrt (q1);
638 
639  const int N = (n - 1) / 2;
640  Array<double> ai;
641 
642  for (int i = 1; i <= N; ++i)
643  {
644  double num = 0.0;
645  double delta = 1.0;
646  int m = 0;
647 
648  while (std::abs (delta) > 1e-100)
649  {
650  delta = std::pow (-1, m) * std::pow (q, m * (m + 1))
651  * std::sin ((2 * m + 1) * MathConstants<double>::pi * i / (double) n);
652  num += delta;
653  m++;
654  }
655 
656  num *= 2 * std::pow (q, 0.25);
657 
658  double den = 0.0;
659  delta = 1.0;
660  m = 1;
661 
662  while (std::abs (delta) > 1e-100)
663  {
664  delta = std::pow (-1, m) * std::pow (q, m * m)
665  * std::cos (m * MathConstants<double>::twoPi * i / (double) n);
666  den += delta;
667  ++m;
668  }
669 
670  den = 1 + 2 * den;
671 
672  auto wi = num / den;
673  auto api = std::sqrt ((1 - wi * wi * k) * (1 - wi * wi / k)) / (1 + wi * wi);
674 
675  ai.add ((1 - api) / (1 + api));
676  }
677 
679 
680  for (int i = 0; i < N; i += 2)
681  structure.directPath.add (new IIR::Coefficients<FloatType> (static_cast<FloatType> (ai[i]),
682  0, 1, 1, 0, static_cast<FloatType> (ai[i])));
683 
684  structure.delayedPath.add (new IIR::Coefficients<FloatType> (0, 1, 1, 0));
685 
686  for (int i = 1; i < N; i += 2)
687  structure.delayedPath.add (new IIR::Coefficients<FloatType> (static_cast<FloatType> (ai[i]),
688  0, 1, 1, 0, static_cast<FloatType> (ai[i])));
689 
690  structure.alpha.addArray (ai);
691 
692  return structure;
693 }
694 
695 
696 template struct FilterDesign<float>;
697 template struct FilterDesign<double>;
698 
699 } // namespace dsp
700 } // namespace juce
static ReferenceCountedArray< IIRCoefficients > designIIRLowpassHighOrderChebyshev1Method(FloatType frequency, double sampleRate, FloatType normalisedTransitionWidth, FloatType passbandAmplitudedB, FloatType stopbandAmplitudedB)
This method returns an array of IIR::Coefficients, made to be used in cascaded IIRFilters, providing a minimum phase low-pass filter without any ripple in the stop band only.
A set of coefficients for use in an FIRFilter object.
static ReferenceCountedArray< IIRCoefficients > designIIRLowpassHighOrderEllipticMethod(FloatType frequency, double sampleRate, FloatType normalisedTransitionWidth, FloatType passbandAmplitudedB, FloatType stopbandAmplitudedB)
This method returns an array of IIR::Coefficients, made to be used in cascaded IIR::Filters, providing a minimum phase low-pass filter with ripples in both the pass band and in the stop band.
static ReferenceCountedArray< IIRCoefficients > designIIRLowpassHighOrderChebyshev2Method(FloatType frequency, double sampleRate, FloatType normalisedTransitionWidth, FloatType passbandAmplitudedB, FloatType stopbandAmplitudedB)
This method returns an array of IIR::Coefficients, made to be used in cascaded IIRFilters, providing a minimum phase low-pass filter without any ripple in the pass band only.
static IIRPolyphaseAllpassStructure designIIRLowpassHalfBandPolyphaseAllpassMethod(FloatType normalisedTransitionWidth, FloatType stopbandAmplitudedB)
This method generates arrays of IIR::Coefficients for a low-pass filter, with a cutoff frequency at h...
static Type decibelsToGain(Type decibels, Type minusInfinityDb=Type(defaultMinusInfinitydB))
Converts a dBFS value to its equivalent gain level.
Definition: juce_Decibels.h:46
void addArray(const Type *elementsToAdd, int numElementsToAdd)
Adds elements from an array to the end of this array.
Definition: juce_Array.h:587
void add(const ElementType &newElement)
Appends a new element at the end of the array.
Definition: juce_Array.h:422
ObjectClass * add(ObjectClass *newObject)
Appends a new object to the end of the array.
ReferenceCountedObjectPtr< Coefficients > Ptr
The Coefficients structure is ref-counted, so this is a handy type that can be used as a pointer to o...
static Matrix identity(size_t size)
Creates the identity matrix.
Definition: juce_Matrix.cpp:33
static FIRCoefficientsPtr designFIRLowpassLeastSquaresMethod(FloatType frequency, double sampleRate, size_t order, FloatType normalisedTransitionWidth, FloatType stopBandWeight)
This method generates a FIR::Coefficients for a low-pass filter, by minimizing the average error betw...
void resize(int targetNumItems)
This will enlarge or shrink the array to the given number of elements, by adding or removing items fr...
Definition: juce_Array.h:674
static Matrix hankel(const Matrix &vector, size_t size, size_t offset=0)
Creates a squared size x size Hankel Matrix from a vector with an optional offset.
Definition: juce_Matrix.cpp:62
NumericType * getRawCoefficients() noexcept
Returns a raw data pointer to the coefficients.
static FIRCoefficientsPtr designFIRLowpassTransitionMethod(FloatType frequency, double sampleRate, size_t order, FloatType normalisedTransitionWidth, FloatType spline)
This method is also a variant of the function designFIRLowpassWindowMethod, using a rectangular windo...
static void ellipticIntegralK(double k, double &K, double &Kp) noexcept
Computes the complete elliptic integral of the first kind K for a given double value k...
void setUnchecked(int indexToChange, ParameterType newValue)
Replaces an element with a new value without doing any bounds-checking.
Definition: juce_Array.h:572
A set of coefficients for use in an Filter object.
void multiplyWithWindowingTable(FloatType *samples, size_t size) noexcept
Multiplies the content of a buffer with the given window.
A class which provides multiple windowing functions useful for filter design and spectrum analyzers...
static FIRCoefficientsPtr designFIRLowpassWindowMethod(FloatType frequency, double sampleRate, size_t order, WindowingMethod type, FloatType beta=static_cast< FloatType >(2))
This method generates a FIR::Coefficients for a low-pass filter, using the windowing design method...
int size() const noexcept
Returns the current number of elements in the array.
Definition: juce_Array.h:219
static Complex< double > asne(Complex< double > w, double k) noexcept
Computes the inverse of the Jacobian elliptic function sn for the elliptic modulus k and the quarter-...
static FIRCoefficientsPtr designFIRLowpassHalfBandEquirippleMethod(FloatType normalisedTransitionWidth, FloatType amplitudedB)
This method generates a FIR::Coefficients for a low-pass filter, with a cutoff frequency at half band...
A smart-pointer class which points to a reference-counted object.
static Complex< double > sne(Complex< double > u, double k) noexcept
Computes the Jacobian elliptic function sn for the elliptic modulus k and the quarter-period units u...
static ReferenceCountedArray< IIRCoefficients > designIIRLowpassHighOrderButterworthMethod(FloatType frequency, double sampleRate, FloatType normalisedTransitionWidth, FloatType passbandAmplitudedB, FloatType stopbandAmplitudedB)
This method returns an array of IIR::Coefficients, made to be used in cascaded IIRFilters, providing a minimum phase low-pass filter without any ripple in the pass band and in the stop band.
The structure returned by the function designIIRLowpassHalfBandPolyphaseAllpassMethod.
Commonly used mathematical constants.
static Matrix toeplitz(const Matrix &vector, size_t size)
Creates a Toeplitz Matrix from a vector with a given squared size.
Definition: juce_Matrix.cpp:44
static Complex< double > cde(Complex< double > u, double k) noexcept
Computes the Jacobian elliptic function cd for the elliptic modulus k and the quarter-period units u...
General matrix and vectors class, meant for classic math manipulation such as additions, multiplications, and linear systems of equations solving.
Definition: juce_Matrix.h:45
This class provides a set of functions which generates FIR::Coefficients and IIR::Coefficients, of high-order low-pass filters.
static FIRCoefficientsPtr designFIRLowpassKaiserMethod(FloatType frequency, double sampleRate, FloatType normalisedTransitionWidth, FloatType amplitudedB)
This a variant of the function designFIRLowpassWindowMethod, which allows the user to specify a trans...
Holds a list of objects derived from ReferenceCountedObject, or which implement basic reference-count...
static ReferenceCountedArray< IIRCoefficients > designIIRHighpassHighOrderButterworthMethod(FloatType frequency, double sampleRate, int order)
This method returns an array of IIR::Coefficients, made to be used in cascaded IIRFilters, providing a minimum phase high-pass filter without any ripple in the pass band and in the stop band.