29template <
typename FloatType>
32 double sampleRate,
size_t order,
36 jassert (sampleRate > 0);
37 jassert (frequency > 0 && frequency <= sampleRate * 0.5);
41 auto* c = result->getRawCoefficients();
42 auto normalisedFrequency = frequency / sampleRate;
44 for (
size_t i = 0; i <= order; ++i)
48 c[i] =
static_cast<FloatType
> (normalisedFrequency * 2);
53 c[i] =
static_cast<FloatType
> (std::sin (2.0 * indice * normalisedFrequency) / indice);
63template <
typename FloatType>
66 FloatType normalisedTransitionWidth,
67 FloatType amplitudedB)
69 jassert (sampleRate > 0);
70 jassert (frequency > 0 && frequency <= sampleRate * 0.5);
71 jassert (normalisedTransitionWidth > 0 && normalisedTransitionWidth <= 0.5);
72 jassert (amplitudedB >= -100 && amplitudedB <= 0);
76 if (amplitudedB < -50)
77 beta =
static_cast<FloatType
> (0.1102 * (-amplitudedB - 8.7));
78 else if (amplitudedB <= -21)
79 beta =
static_cast<FloatType
> (0.5842 * std::pow (-amplitudedB - 21, 0.4) + 0.07886 * (-amplitudedB - 21));
81 int order = amplitudedB < -21 ? roundToInt (std::ceil ((-amplitudedB - 7.95) / (2.285 * normalisedTransitionWidth *
MathConstants<double>::twoPi)))
86 return designFIRLowpassWindowMethod (frequency, sampleRate,
static_cast<size_t> (order),
91template <
typename FloatType>
94 FloatType normalisedTransitionWidth, FloatType spline)
96 jassert (sampleRate > 0);
97 jassert (frequency > 0 && frequency <= sampleRate * 0.5);
98 jassert (normalisedTransitionWidth > 0 && normalisedTransitionWidth <= 0.5);
99 jassert (spline >= 1.0 && spline <= 4.0);
101 auto normalisedFrequency = frequency /
static_cast<FloatType
> (sampleRate);
104 auto* c = result->getRawCoefficients();
106 for (
size_t i = 0; i <= order; ++i)
108 if (i == order / 2 && order % 2 == 0)
110 c[i] =
static_cast<FloatType
> (2 * normalisedFrequency);
116 c[i] =
static_cast<FloatType
> (std::sin (2 * indice * normalisedFrequency)
117 / indice * std::pow (std::sin (indice2) / indice2, spline));
124template <
typename FloatType>
127 double sampleRate,
size_t order,
128 FloatType normalisedTransitionWidth,
129 FloatType stopBandWeight)
131 jassert (sampleRate > 0);
132 jassert (frequency > 0 && frequency <= sampleRate * 0.5);
133 jassert (normalisedTransitionWidth > 0 && normalisedTransitionWidth <= 0.5);
134 jassert (stopBandWeight >= 1.0 && stopBandWeight <= 100.0);
136 auto normalisedFrequency =
static_cast<double> (frequency) / sampleRate;
144 auto* c = result->getRawCoefficients();
146 auto sinc = [] (
double x)
154 auto M = (N - 1) / 2;
162 for (
size_t i = 0; i <= M; ++i)
163 b (i, 0) = factorp * sinc (factorp * (
double) i);
165 q (0, 0) = factorp + stopBandWeight * (1.0 - factors);
167 for (
size_t i = 1; i <= 2 * M; ++i)
168 q (i, 0) = factorp * sinc (factorp * (
double) i) - stopBandWeight * factors * sinc (factors * (
double) i);
177 c[M] =
static_cast<FloatType
> (b (0, 0));
179 for (
size_t i = 1; i <= M; ++i)
181 c[M - i] =
static_cast<FloatType
> (b (i, 0) * 0.5);
182 c[M + i] =
static_cast<FloatType
> (b (i, 0) * 0.5);
197 for (
size_t i = 0; i < M; ++i)
198 b (i, 0) = factorp * sinc (factorp * ((
double) i + 0.5));
200 for (
size_t i = 0; i < 2 * M; ++i)
202 qp (i, 0) = 0.25 * factorp * sinc (factorp * (
double) i);
203 qs (i, 0) = -0.25 * stopBandWeight * factors * sinc (factors * (
double) i);
212 Id *= (0.25 * stopBandWeight);
223 for (
size_t i = 0; i < M; ++i)
225 c[M - i - 1] =
static_cast<FloatType
> (b (i, 0) * 0.25);
226 c[M + i] =
static_cast<FloatType
> (b (i, 0) * 0.25);
233template <
typename FloatType>
236 FloatType amplitudedB)
238 jassert (normalisedTransitionWidth > 0 && normalisedTransitionWidth <= 0.5);
239 jassert (amplitudedB >= -300 && amplitudedB <= -10);
243 auto n = roundToInt (std::ceil ((amplitudedB - 18.18840664 * wpT + 33.64775300) / (18.54155181 * wpT - 29.13196871)));
244 auto kp = (n * wpT - 1.57111377 * n + 0.00665857) / (-1.01927560 * n + 0.37221484);
245 auto A = (0.01525753 * n + 0.03682344 + 9.24760314 / (double) n) * kp + 1.01701407 + 0.73512298 / (double) n;
246 auto B = (0.00233667 * n - 1.35418408 + 5.75145813 / (double) n) * kp + 1.02999650 - 0.72759508 / (double) n;
251 auto diff = (hn.size() - hnm.size()) / 2;
253 for (
int i = 0; i < diff; ++i)
261 for (
int i = 0; i < hn.size(); ++i)
262 hh.setUnchecked (i, A * hh[i] + B * hnm[i]);
265 auto* c = result->getRawCoefficients();
267 for (
int i = 0; i < hh.size(); ++i)
268 c[i] = (
float) hh[i];
273 return 2.0 * result->getMagnitudeForFrequency (0.5, 1.0);
277 if (std::abs (w01) > 1.0)
278 return 2.0 * result->getMagnitudeForFrequency (0.5, 1.0);
280 auto om01 = std::acos (-w01);
284 for (
int i = 0; i < hh.size(); ++i)
285 c[i] =
static_cast<FloatType
> ((A * hn[i] + B * hnm[i]) / NN);
287 c[2 * n + 1] =
static_cast<FloatType
> (0.5);
292template <
typename FloatType>
298 alpha.
setUnchecked (2 * n, 1.0 / std::pow (1.0 - kp * kp, n));
301 alpha.
setUnchecked (2 * n - 2, -(2 * n * kp * kp + 1) * alpha[2 * n]);
304 alpha.
setUnchecked (2 * n - 4, -(4 * n + 1 + (n - 1) * (2 * n - 1) * kp * kp) / (2.0 * n) * alpha[2 * n - 2]
305 - (2 * n + 1) * ((n + 1) * kp * kp + 1) / (2.0 * n) * alpha[2 * n]);
307 for (
int k = n; k >= 3; --k)
309 auto c1 = (3 * (n*(n + 2) - k * (k - 2)) + 2 * k - 3 + 2 * (k - 2)*(2 * k - 3) * kp * kp) * alpha[2 * k - 4];
310 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];
311 auto c3 = (n * (n + 2) - (k - 1) * (k + 1)) * alpha[2 * k];
312 auto c4 = (n * (n + 2) - (k - 3) * (k - 1));
318 ai.
resize (2 * n + 1 + 1);
320 for (
int k = 0; k <= n; ++k)
321 ai.
setUnchecked (2 * k + 1, alpha[2 * k] / (2.0 * k + 1.0));
324 hn.
resize (2 * n + 1 + 2 * n + 1 + 1);
326 for (
int k = 0; k <= n; ++k)
328 hn.
setUnchecked (2 * n + 1 + (2 * k + 1), 0.5 * ai[2 * k + 1]);
329 hn.
setUnchecked (2 * n + 1 - (2 * k + 1), 0.5 * ai[2 * k + 1]);
335template <
typename FloatType>
336ReferenceCountedArray<IIR::Coefficients<FloatType>>
338 FloatType normalisedTransitionWidth,
339 FloatType passbandAmplitudedB,
340 FloatType stopbandAmplitudedB)
342 return designIIRLowpassHighOrderGeneralMethod (0, frequency, sampleRate, normalisedTransitionWidth,
343 passbandAmplitudedB, stopbandAmplitudedB);
346template <
typename FloatType>
349 FloatType normalisedTransitionWidth,
350 FloatType passbandAmplitudedB,
351 FloatType stopbandAmplitudedB)
353 return designIIRLowpassHighOrderGeneralMethod (1, frequency, sampleRate, normalisedTransitionWidth,
354 passbandAmplitudedB, stopbandAmplitudedB);
357template <
typename FloatType>
360 FloatType normalisedTransitionWidth,
361 FloatType passbandAmplitudedB,
362 FloatType stopbandAmplitudedB)
364 return designIIRLowpassHighOrderGeneralMethod (2, frequency, sampleRate, normalisedTransitionWidth,
365 passbandAmplitudedB, stopbandAmplitudedB);
368template <
typename FloatType>
371 FloatType normalisedTransitionWidth,
372 FloatType passbandAmplitudedB,
373 FloatType stopbandAmplitudedB)
375 return designIIRLowpassHighOrderGeneralMethod (3, frequency, sampleRate, normalisedTransitionWidth,
376 passbandAmplitudedB, stopbandAmplitudedB);
379template <
typename FloatType>
382 FloatType normalisedTransitionWidth,
383 FloatType passbandAmplitudedB,
384 FloatType stopbandAmplitudedB)
386 jassert (0 < sampleRate);
387 jassert (0 < frequency && frequency <= sampleRate * 0.5);
388 jassert (0 < normalisedTransitionWidth && normalisedTransitionWidth <= 0.5);
389 jassert (-20 < passbandAmplitudedB && passbandAmplitudedB < 0);
390 jassert (-300 < stopbandAmplitudedB && stopbandAmplitudedB < -20);
392 auto normalisedFrequency = frequency / sampleRate;
394 auto fp = normalisedFrequency - normalisedTransitionWidth / 2;
395 jassert (0.0 < fp && fp < 0.5);
397 auto fs = normalisedFrequency + normalisedTransitionWidth / 2;
398 jassert (0.0 < fs && fs < 0.5);
400 double Ap = passbandAmplitudedB;
401 double As = stopbandAmplitudedB;
404 auto epsp = std::sqrt (1.0 / (Gp * Gp) - 1.0);
405 auto epss = std::sqrt (1.0 / (Gs * Gs) - 1.0);
411 auto k = omegap / omegas;
412 auto k1 = epsp / epss;
418 N = roundToInt (std::ceil (std::log (1.0 / k1) / std::log (1.0 / k)));
420 else if (type == 1 || type == 2)
422 N = roundToInt (std::ceil (std::acosh (1.0 / k1) / std::acosh (1.0 / k)));
426 double K, Kp, K1, K1p;
431 N = roundToInt (std::ceil ((K1p * K) / (K1 * Kp)));
435 const int L = (N - r) / 2;
436 const double H0 = (type == 1 || type == 3) ? std::pow (Gp, 1.0 - r) : 1.0;
438 Array<Complex<double>> pa, za;
439 Complex<double> j (0, 1);
444 pa.add (-omegap * std::pow (epsp, -1.0 / (
double) N));
446 for (
int i = 1; i <= L; ++i)
448 auto ui = (2 * i - 1.0) / (
double) N;
449 pa.add (omegap * std::pow (epsp, -1.0 / (
double) N) * j * exp (ui * halfPi * j));
454 auto v0 = std::asinh (1.0 / epsp) / (N * halfPi);
457 pa.add (-omegap * std::sinh (v0 * halfPi));
459 for (
int i = 1; i <= L; ++i)
461 auto ui = (2 * i - 1.0) / (
double) N;
462 pa.add (omegap * j * std::cos ((ui - j * v0) * halfPi));
467 auto v0 = std::asinh (epss) / (N * halfPi);
470 pa.add (-1.0 / (k / omegap * std::sinh (v0 * halfPi)));
472 for (
int i = 1; i <= L; ++i)
474 auto ui = (2 * i - 1.0) / (
double) N;
476 pa.add (1.0 / (k / omegap * j * std::cos ((ui - j * v0) * halfPi)));
477 za.add (1.0 / (k / omegap * j * std::cos (ui * halfPi)));
487 for (
int i = 1; i <= L; ++i)
489 auto ui = (2 * i - 1.0) / (
double) N;
493 za.add (omegap * j / (k * zetai));
497 Array<Complex<double>> p, z, g;
501 p.add ((1.0 + pa[0]) / (1.0 - pa[0]));
502 g.add (0.5 * (1.0 - p[0]));
505 for (
int i = 0; i < L; ++i)
507 p.add ((1.0 + pa[i + r]) / (1.0 - pa[i + r]));
508 z.add (za.size() == 0 ? -1.0 : (1.0 + za[i]) / (1.0 - za[i]));
509 g.add ((1.0 - p[i + r]) / (1.0 - z[i]));
512 ReferenceCountedArray<IIR::Coefficients<FloatType>> cascadedCoefficients;
516 auto b0 =
static_cast<FloatType
> (H0 * std::real (g[0]));
518 auto a1 =
static_cast<FloatType
> (-std::real (p[0]));
520 cascadedCoefficients.add (
new IIR::Coefficients<FloatType> (b0, b1, 1.0f, a1));
523 for (
int i = 0; i < L; ++i)
525 auto gain = std::pow (std::abs (g[i + r]), 2.0);
527 auto b0 =
static_cast<FloatType
> (gain);
528 auto b1 =
static_cast<FloatType
> (std::real (-z[i] - std::conj (z[i])) * gain);
529 auto b2 =
static_cast<FloatType
> (std::real ( z[i] * std::conj (z[i])) * gain);
531 auto a1 =
static_cast<FloatType
> (std::real (-p[i+r] - std::conj (p[i + r])));
532 auto a2 =
static_cast<FloatType
> (std::real ( p[i+r] * std::conj (p[i + r])));
534 cascadedCoefficients.add (
new IIR::Coefficients<FloatType> (b0, b1, b2, 1, a1, a2));
537 return cascadedCoefficients;
540template <
typename FloatType>
541ReferenceCountedArray<IIR::Coefficients<FloatType>>
543 double sampleRate,
int order)
545 jassert (sampleRate > 0);
546 jassert (frequency > 0 && frequency <= sampleRate * 0.5);
555 for (
int i = 0; i < order / 2; ++i)
559 static_cast<FloatType
> (Q)));
564 for (
int i = 0; i < order / 2; ++i)
568 static_cast<FloatType
> (Q)));
575template <
typename FloatType>
578 double sampleRate,
int order)
580 jassert (sampleRate > 0);
581 jassert (frequency > 0 && frequency <= sampleRate * 0.5);
590 for (
int i = 0; i < order / 2; ++i)
594 static_cast<FloatType
> (Q)));
599 for (
int i = 0; i < order / 2; ++i)
603 static_cast<FloatType
> (Q)));
610template <
typename FloatType>
613 FloatType stopbandAmplitudedB)
615 jassert (normalisedTransitionWidth > 0 && normalisedTransitionWidth <= 0.5);
616 jassert (stopbandAmplitudedB > -300 && stopbandAmplitudedB < -10);
622 auto kp = std::sqrt (1.0 - k * k);
623 auto e = (1 - std::sqrt (kp)) / (1 + std::sqrt (kp)) * 0.5;
624 auto q = e + 2 * std::pow (e, 5.0) + 15 * std::pow (e, 9.0) + 150 * std::pow (e, 13.0);
626 auto k1 = ds * ds / (1 - ds * ds);
627 int n = roundToInt (std::ceil (std::log (k1 * k1 / 16) / std::log (q)));
635 auto q1 = std::pow (q, (
double) n);
636 k1 = 4 * std::sqrt (q1);
638 const int N = (n - 1) / 2;
641 for (
int i = 1; i <= N; ++i)
647 while (std::abs (delta) > 1e-100)
649 delta = std::pow (-1, m) * std::pow (q, m * (m + 1))
655 num *= 2 * std::pow (q, 0.25);
661 while (std::abs (delta) > 1e-100)
663 delta = std::pow (-1, m) * std::pow (q, m * m)
672 auto api = std::sqrt ((1 - wi * wi * k) * (1 - wi * wi / k)) / (1 + wi * wi);
674 ai.
add ((1 - api) / (1 + api));
679 for (
int i = 0; i < N; i += 2)
681 0, 1, 1, 0,
static_cast<FloatType
> (ai[i])));
685 for (
int i = 1; i < N; i += 2)
687 0, 1, 1, 0,
static_cast<FloatType
> (ai[i])));
void setUnchecked(int indexToChange, ParameterType newValue)
void addArray(const Type *elementsToAdd, int numElementsToAdd)
void add(const ElementType &newElement)
void resize(int targetNumItems)
static Type decibelsToGain(Type decibels, Type minusInfinityDb=Type(defaultMinusInfinitydB))
ObjectClass * add(ObjectClass *newObject)
static Matrix hankel(const Matrix &vector, size_t size, size_t offset=0)
static Matrix identity(size_t size)
static Matrix toeplitz(const Matrix &vector, size_t size)
void multiplyWithWindowingTable(FloatType *samples, size_t size) const noexcept
ReferenceCountedObjectPtr< Coefficients > Ptr
static ReferenceCountedArray< IIRCoefficients > designIIRLowpassHighOrderButterworthMethod(FloatType frequency, double sampleRate, FloatType normalisedTransitionWidth, FloatType passbandAmplitudedB, FloatType stopbandAmplitudedB)
static FIRCoefficientsPtr designFIRLowpassLeastSquaresMethod(FloatType frequency, double sampleRate, size_t order, FloatType normalisedTransitionWidth, FloatType stopBandWeight)
static FIRCoefficientsPtr designFIRLowpassKaiserMethod(FloatType frequency, double sampleRate, FloatType normalisedTransitionWidth, FloatType amplitudedB)
static IIRPolyphaseAllpassStructure designIIRLowpassHalfBandPolyphaseAllpassMethod(FloatType normalisedTransitionWidth, FloatType stopbandAmplitudedB)
static ReferenceCountedArray< IIRCoefficients > designIIRLowpassHighOrderChebyshev1Method(FloatType frequency, double sampleRate, FloatType normalisedTransitionWidth, FloatType passbandAmplitudedB, FloatType stopbandAmplitudedB)
static ReferenceCountedArray< IIRCoefficients > designIIRHighpassHighOrderButterworthMethod(FloatType frequency, double sampleRate, int order)
static FIRCoefficientsPtr designFIRLowpassTransitionMethod(FloatType frequency, double sampleRate, size_t order, FloatType normalisedTransitionWidth, FloatType spline)
static ReferenceCountedArray< IIRCoefficients > designIIRLowpassHighOrderChebyshev2Method(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 FIRCoefficientsPtr designFIRLowpassHalfBandEquirippleMethod(FloatType normalisedTransitionWidth, FloatType amplitudedB)
static FIRCoefficientsPtr designFIRLowpassWindowMethod(FloatType frequency, double sampleRate, size_t order, WindowingMethod type, FloatType beta=static_cast< FloatType >(2))
static Complex< double > sne(Complex< double > u, double k) noexcept
static Complex< double > cde(Complex< double > u, double k) noexcept
static Complex< double > asne(Complex< double > w, double k) noexcept
static void ellipticIntegralK(double k, double &K, double &Kp) noexcept