OpenShot Audio Library | OpenShotAudio  0.3.1
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)
static ReferenceCountedArray< IIRCoefficients > designIIRLowpassHighOrderEllipticMethod(FloatType frequency, double sampleRate, FloatType normalisedTransitionWidth, FloatType passbandAmplitudedB, FloatType stopbandAmplitudedB)
static ReferenceCountedArray< IIRCoefficients > designIIRLowpassHighOrderChebyshev2Method(FloatType frequency, double sampleRate, FloatType normalisedTransitionWidth, FloatType passbandAmplitudedB, FloatType stopbandAmplitudedB)
static IIRPolyphaseAllpassStructure designIIRLowpassHalfBandPolyphaseAllpassMethod(FloatType normalisedTransitionWidth, FloatType stopbandAmplitudedB)
static Type decibelsToGain(Type decibels, Type minusInfinityDb=Type(defaultMinusInfinitydB))
Definition: juce_Decibels.h:42
void addArray(const Type *elementsToAdd, int numElementsToAdd)
Definition: juce_Array.h:583
void add(const ElementType &newElement)
Definition: juce_Array.h:418
ObjectClass * add(ObjectClass *newObject)
ReferenceCountedObjectPtr< Coefficients > Ptr
static Matrix identity(size_t size)
Definition: juce_Matrix.cpp:33
static FIRCoefficientsPtr designFIRLowpassLeastSquaresMethod(FloatType frequency, double sampleRate, size_t order, FloatType normalisedTransitionWidth, FloatType stopBandWeight)
void resize(int targetNumItems)
Definition: juce_Array.h:670
static Matrix hankel(const Matrix &vector, size_t size, size_t offset=0)
Definition: juce_Matrix.cpp:62
NumericType * getRawCoefficients() noexcept
static FIRCoefficientsPtr designFIRLowpassTransitionMethod(FloatType frequency, double sampleRate, size_t order, FloatType normalisedTransitionWidth, FloatType spline)
static void ellipticIntegralK(double k, double &K, double &Kp) noexcept
void setUnchecked(int indexToChange, ParameterType newValue)
Definition: juce_Array.h:568
void multiplyWithWindowingTable(FloatType *samples, size_t size) noexcept
static FIRCoefficientsPtr designFIRLowpassWindowMethod(FloatType frequency, double sampleRate, size_t order, WindowingMethod type, FloatType beta=static_cast< FloatType >(2))
int size() const noexcept
Definition: juce_Array.h:215
static Complex< double > asne(Complex< double > w, double k) noexcept
static FIRCoefficientsPtr designFIRLowpassHalfBandEquirippleMethod(FloatType normalisedTransitionWidth, FloatType amplitudedB)
static Complex< double > sne(Complex< double > u, double k) noexcept
static ReferenceCountedArray< IIRCoefficients > designIIRLowpassHighOrderButterworthMethod(FloatType frequency, double sampleRate, FloatType normalisedTransitionWidth, FloatType passbandAmplitudedB, FloatType stopbandAmplitudedB)
static Matrix toeplitz(const Matrix &vector, size_t size)
Definition: juce_Matrix.cpp:44
static Complex< double > cde(Complex< double > u, double k) noexcept
static FIRCoefficientsPtr designFIRLowpassKaiserMethod(FloatType frequency, double sampleRate, FloatType normalisedTransitionWidth, FloatType amplitudedB)
static ReferenceCountedArray< IIRCoefficients > designIIRHighpassHighOrderButterworthMethod(FloatType frequency, double sampleRate, int order)