Loading...
Searching...
No Matches
juce_FilterDesign.cpp
/*
==============================================================================
This file is part of the JUCE library.
Copyright (c) 2022 - Raw Material Software Limited
JUCE is an open source library subject to commercial or open-source
licensing.
By using JUCE, you agree to the terms of both the JUCE 7 End-User License
Agreement and JUCE Privacy Policy.
End User License Agreement: www.juce.com/juce-7-licence
Privacy Policy: www.juce.com/juce-privacy-policy
Or: You may also use this code under the terms of the GPL v3 (see
www.gnu.org/licenses).
JUCE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY, AND ALL WARRANTIES, WHETHER
EXPRESSED OR IMPLIED, INCLUDING MERCHANTABILITY AND FITNESS FOR PURPOSE, ARE
DISCLAIMED.
==============================================================================
*/
namespace juce::dsp
{
template <typename FloatType>
typename FIR::Coefficients<FloatType>::Ptr
FilterDesign<FloatType>::designFIRLowpassWindowMethod (FloatType frequency,
double sampleRate, size_t order,
WindowingMethod type,
FloatType beta)
{
jassert (sampleRate > 0);
jassert (frequency > 0 && frequency <= sampleRate * 0.5);
auto* result = new typename FIR::Coefficients<FloatType> (order + 1u);
auto* c = result->getRawCoefficients();
auto normalisedFrequency = frequency / sampleRate;
for (size_t i = 0; i <= order; ++i)
{
if (i == order / 2)
{
c[i] = static_cast<FloatType> (normalisedFrequency * 2);
}
else
{
auto indice = MathConstants<double>::pi * (static_cast<double> (i) - 0.5 * static_cast<double> (order));
c[i] = static_cast<FloatType> (std::sin (2.0 * indice * normalisedFrequency) / indice);
}
}
WindowingFunction<FloatType> theWindow (order + 1, type, false, beta);
theWindow.multiplyWithWindowingTable (c, order + 1);
return *result;
}
template <typename FloatType>
typename FIR::Coefficients<FloatType>::Ptr
FilterDesign<FloatType>::designFIRLowpassKaiserMethod (FloatType frequency, double sampleRate,
FloatType normalisedTransitionWidth,
FloatType amplitudedB)
{
jassert (sampleRate > 0);
jassert (frequency > 0 && frequency <= sampleRate * 0.5);
jassert (normalisedTransitionWidth > 0 && normalisedTransitionWidth <= 0.5);
jassert (amplitudedB >= -100 && amplitudedB <= 0);
FloatType beta = 0;
if (amplitudedB < -50)
beta = static_cast<FloatType> (0.1102 * (-amplitudedB - 8.7));
else if (amplitudedB <= -21)
beta = static_cast<FloatType> (0.5842 * std::pow (-amplitudedB - 21, 0.4) + 0.07886 * (-amplitudedB - 21));
int order = amplitudedB < -21 ? roundToInt (std::ceil ((-amplitudedB - 7.95) / (2.285 * normalisedTransitionWidth * MathConstants<double>::twoPi)))
: roundToInt (std::ceil (5.79 / (normalisedTransitionWidth * MathConstants<double>::twoPi)));
jassert (order >= 0);
return designFIRLowpassWindowMethod (frequency, sampleRate, static_cast<size_t> (order),
WindowingFunction<FloatType>::kaiser, beta);
}
template <typename FloatType>
typename FIR::Coefficients<FloatType>::Ptr
FilterDesign<FloatType>::designFIRLowpassTransitionMethod (FloatType frequency, double sampleRate, size_t order,
FloatType normalisedTransitionWidth, FloatType spline)
{
jassert (sampleRate > 0);
jassert (frequency > 0 && frequency <= sampleRate * 0.5);
jassert (normalisedTransitionWidth > 0 && normalisedTransitionWidth <= 0.5);
jassert (spline >= 1.0 && spline <= 4.0);
auto normalisedFrequency = frequency / static_cast<FloatType> (sampleRate);
auto* result = new typename FIR::Coefficients<FloatType> (order + 1u);
auto* c = result->getRawCoefficients();
for (size_t i = 0; i <= order; ++i)
{
if (i == order / 2 && order % 2 == 0)
{
c[i] = static_cast<FloatType> (2 * normalisedFrequency);
}
else
{
auto indice = MathConstants<double>::pi * ((double) i - 0.5 * (double) order);
auto indice2 = MathConstants<double>::pi * normalisedTransitionWidth * ((double) i - 0.5 * (double) order) / spline;
c[i] = static_cast<FloatType> (std::sin (2 * indice * normalisedFrequency)
/ indice * std::pow (std::sin (indice2) / indice2, spline));
}
}
return *result;
}
template <typename FloatType>
typename FIR::Coefficients<FloatType>::Ptr
FilterDesign<FloatType>::designFIRLowpassLeastSquaresMethod (FloatType frequency,
double sampleRate, size_t order,
FloatType normalisedTransitionWidth,
FloatType stopBandWeight)
{
jassert (sampleRate > 0);
jassert (frequency > 0 && frequency <= sampleRate * 0.5);
jassert (normalisedTransitionWidth > 0 && normalisedTransitionWidth <= 0.5);
jassert (stopBandWeight >= 1.0 && stopBandWeight <= 100.0);
auto normalisedFrequency = static_cast<double> (frequency) / sampleRate;
auto wp = MathConstants<double>::twoPi * (static_cast<double> (normalisedFrequency - normalisedTransitionWidth / 2.0));
auto ws = MathConstants<double>::twoPi * (static_cast<double> (normalisedFrequency + normalisedTransitionWidth / 2.0));
auto N = order + 1;
auto* result = new typename FIR::Coefficients<FloatType> (static_cast<size_t> (N));
auto* c = result->getRawCoefficients();
auto sinc = [] (double x)
{
return approximatelyEqual (x, 0.0) ? 1 : std::sin (x * MathConstants<double>::pi) / (MathConstants<double>::pi * x);
};
if (N % 2 == 1)
{
// Type I
auto M = (N - 1) / 2;
Matrix<double> b (M + 1, 1),
q (2 * M + 1, 1);
auto factorp = wp / MathConstants<double>::pi;
auto factors = ws / MathConstants<double>::pi;
for (size_t i = 0; i <= M; ++i)
b (i, 0) = factorp * sinc (factorp * (double) i);
q (0, 0) = factorp + stopBandWeight * (1.0 - factors);
for (size_t i = 1; i <= 2 * M; ++i)
q (i, 0) = factorp * sinc (factorp * (double) i) - stopBandWeight * factors * sinc (factors * (double) i);
auto Q1 = Matrix<double>::toeplitz (q, M + 1);
auto Q2 = Matrix<double>::hankel (q, M + 1, 0);
Q1 += Q2; Q1 *= 0.5;
Q1.solve (b);
c[M] = static_cast<FloatType> (b (0, 0));
for (size_t i = 1; i <= M; ++i)
{
c[M - i] = static_cast<FloatType> (b (i, 0) * 0.5);
c[M + i] = static_cast<FloatType> (b (i, 0) * 0.5);
}
}
else
{
// Type II
auto M = N / 2;
Matrix<double> b (M, 1);
Matrix<double> qp (2 * M, 1);
Matrix<double> qs (2 * M, 1);
auto factorp = wp / MathConstants<double>::pi;
auto factors = ws / MathConstants<double>::pi;
for (size_t i = 0; i < M; ++i)
b (i, 0) = factorp * sinc (factorp * ((double) i + 0.5));
for (size_t i = 0; i < 2 * M; ++i)
{
qp (i, 0) = 0.25 * factorp * sinc (factorp * (double) i);
qs (i, 0) = -0.25 * stopBandWeight * factors * sinc (factors * (double) i);
}
auto Q1p = Matrix<double>::toeplitz (qp, M);
auto Q2p = Matrix<double>::hankel (qp, M, 1);
auto Q1s = Matrix<double>::toeplitz (qs, M);
auto Q2s = Matrix<double>::hankel (qs, M, 1);
auto Id = Matrix<double>::identity (M);
Id *= (0.25 * stopBandWeight);
Q1p += Q2p;
Q1s += Q2s;
Q1s += Id;
auto& Q = Q1s;
Q += Q1p;
Q.solve (b);
for (size_t i = 0; i < M; ++i)
{
c[M - i - 1] = static_cast<FloatType> (b (i, 0) * 0.25);
c[M + i] = static_cast<FloatType> (b (i, 0) * 0.25);
}
}
return *result;
}
template <typename FloatType>
typename FIR::Coefficients<FloatType>::Ptr
FilterDesign<FloatType>::designFIRLowpassHalfBandEquirippleMethod (FloatType normalisedTransitionWidth,
FloatType amplitudedB)
{
jassert (normalisedTransitionWidth > 0 && normalisedTransitionWidth <= 0.5);
jassert (amplitudedB >= -300 && amplitudedB <= -10);
auto wpT = (0.5 - normalisedTransitionWidth) * MathConstants<double>::pi;
auto n = roundToInt (std::ceil ((amplitudedB - 18.18840664 * wpT + 33.64775300) / (18.54155181 * wpT - 29.13196871)));
auto kp = (n * wpT - 1.57111377 * n + 0.00665857) / (-1.01927560 * n + 0.37221484);
auto A = (0.01525753 * n + 0.03682344 + 9.24760314 / (double) n) * kp + 1.01701407 + 0.73512298 / (double) n;
auto B = (0.00233667 * n - 1.35418408 + 5.75145813 / (double) n) * kp + 1.02999650 - 0.72759508 / (double) n;
auto hn = FilterDesign<FloatType>::getPartialImpulseResponseHn (n, kp);
auto hnm = FilterDesign<FloatType>::getPartialImpulseResponseHn (n - 1, kp);
auto diff = (hn.size() - hnm.size()) / 2;
for (int i = 0; i < diff; ++i)
{
hnm.add (0.0);
hnm.insert (0, 0.0);
}
auto hh = hn;
for (int i = 0; i < hn.size(); ++i)
hh.setUnchecked (i, A * hh[i] + B * hnm[i]);
auto* result = new typename FIR::Coefficients<FloatType> (static_cast<size_t> (hh.size()));
auto* c = result->getRawCoefficients();
for (int i = 0; i < hh.size(); ++i)
c[i] = (float) hh[i];
auto NN = [&]
{
if (n % 2 == 0)
return 2.0 * result->getMagnitudeForFrequency (0.5, 1.0);
auto w01 = std::sqrt (kp * kp + (1 - kp * kp) * std::pow (std::cos (MathConstants<double>::pi / (2.0 * n + 1.0)), 2.0));
if (std::abs (w01) > 1.0)
return 2.0 * result->getMagnitudeForFrequency (0.5, 1.0);
auto om01 = std::acos (-w01);
return -2.0 * result->getMagnitudeForFrequency (om01 / MathConstants<double>::twoPi, 1.0);
}();
for (int i = 0; i < hh.size(); ++i)
c[i] = static_cast<FloatType> ((A * hn[i] + B * hnm[i]) / NN);
c[2 * n + 1] = static_cast<FloatType> (0.5);
return *result;
}
template <typename FloatType>
Array<double> FilterDesign<FloatType>::getPartialImpulseResponseHn (int n, double kp)
{
Array<double> alpha;
alpha.resize (2 * n + 1);
alpha.setUnchecked (2 * n, 1.0 / std::pow (1.0 - kp * kp, n));
if (n > 0)
alpha.setUnchecked (2 * n - 2, -(2 * n * kp * kp + 1) * alpha[2 * n]);
if (n > 1)
alpha.setUnchecked (2 * n - 4, -(4 * n + 1 + (n - 1) * (2 * n - 1) * kp * kp) / (2.0 * n) * alpha[2 * n - 2]
- (2 * n + 1) * ((n + 1) * kp * kp + 1) / (2.0 * n) * alpha[2 * n]);
for (int k = n; k >= 3; --k)
{
auto c1 = (3 * (n*(n + 2) - k * (k - 2)) + 2 * k - 3 + 2 * (k - 2)*(2 * k - 3) * kp * kp) * alpha[2 * k - 4];
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];
auto c3 = (n * (n + 2) - (k - 1) * (k + 1)) * alpha[2 * k];
auto c4 = (n * (n + 2) - (k - 3) * (k - 1));
alpha.setUnchecked (2 * k - 6, -(c1 + c2 + c3) / c4);
}
Array<double> ai;
ai.resize (2 * n + 1 + 1);
for (int k = 0; k <= n; ++k)
ai.setUnchecked (2 * k + 1, alpha[2 * k] / (2.0 * k + 1.0));
Array<double> hn;
hn.resize (2 * n + 1 + 2 * n + 1 + 1);
for (int k = 0; k <= n; ++k)
{
hn.setUnchecked (2 * n + 1 + (2 * k + 1), 0.5 * ai[2 * k + 1]);
hn.setUnchecked (2 * n + 1 - (2 * k + 1), 0.5 * ai[2 * k + 1]);
}
return hn;
}
template <typename FloatType>
ReferenceCountedArray<IIR::Coefficients<FloatType>>
FilterDesign<FloatType>::designIIRLowpassHighOrderButterworthMethod (FloatType frequency, double sampleRate,
FloatType normalisedTransitionWidth,
FloatType passbandAmplitudedB,
FloatType stopbandAmplitudedB)
{
return designIIRLowpassHighOrderGeneralMethod (0, frequency, sampleRate, normalisedTransitionWidth,
passbandAmplitudedB, stopbandAmplitudedB);
}
template <typename FloatType>
ReferenceCountedArray<IIR::Coefficients<FloatType>>
FilterDesign<FloatType>::designIIRLowpassHighOrderChebyshev1Method (FloatType frequency, double sampleRate,
FloatType normalisedTransitionWidth,
FloatType passbandAmplitudedB,
FloatType stopbandAmplitudedB)
{
return designIIRLowpassHighOrderGeneralMethod (1, frequency, sampleRate, normalisedTransitionWidth,
passbandAmplitudedB, stopbandAmplitudedB);
}
template <typename FloatType>
ReferenceCountedArray<IIR::Coefficients<FloatType>>
FilterDesign<FloatType>::designIIRLowpassHighOrderChebyshev2Method (FloatType frequency, double sampleRate,
FloatType normalisedTransitionWidth,
FloatType passbandAmplitudedB,
FloatType stopbandAmplitudedB)
{
return designIIRLowpassHighOrderGeneralMethod (2, frequency, sampleRate, normalisedTransitionWidth,
passbandAmplitudedB, stopbandAmplitudedB);
}
template <typename FloatType>
ReferenceCountedArray<IIR::Coefficients<FloatType>>
FilterDesign<FloatType>::designIIRLowpassHighOrderEllipticMethod (FloatType frequency, double sampleRate,
FloatType normalisedTransitionWidth,
FloatType passbandAmplitudedB,
FloatType stopbandAmplitudedB)
{
return designIIRLowpassHighOrderGeneralMethod (3, frequency, sampleRate, normalisedTransitionWidth,
passbandAmplitudedB, stopbandAmplitudedB);
}
template <typename FloatType>
ReferenceCountedArray<IIR::Coefficients<FloatType>>
FilterDesign<FloatType>::designIIRLowpassHighOrderGeneralMethod (int type, FloatType frequency, double sampleRate,
FloatType normalisedTransitionWidth,
FloatType passbandAmplitudedB,
FloatType stopbandAmplitudedB)
{
jassert (0 < sampleRate);
jassert (0 < frequency && frequency <= sampleRate * 0.5);
jassert (0 < normalisedTransitionWidth && normalisedTransitionWidth <= 0.5);
jassert (-20 < passbandAmplitudedB && passbandAmplitudedB < 0);
jassert (-300 < stopbandAmplitudedB && stopbandAmplitudedB < -20);
auto normalisedFrequency = frequency / sampleRate;
auto fp = normalisedFrequency - normalisedTransitionWidth / 2;
jassert (0.0 < fp && fp < 0.5);
auto fs = normalisedFrequency + normalisedTransitionWidth / 2;
jassert (0.0 < fs && fs < 0.5);
double Ap = passbandAmplitudedB;
double As = stopbandAmplitudedB;
auto Gp = Decibels::decibelsToGain (Ap, -300.0);
auto Gs = Decibels::decibelsToGain (As, -300.0);
auto epsp = std::sqrt (1.0 / (Gp * Gp) - 1.0);
auto epss = std::sqrt (1.0 / (Gs * Gs) - 1.0);
auto omegap = std::tan (MathConstants<double>::pi * fp);
auto omegas = std::tan (MathConstants<double>::pi * fs);
constexpr auto halfPi = MathConstants<double>::halfPi;
auto k = omegap / omegas;
auto k1 = epsp / epss;
int N;
if (type == 0)
{
N = roundToInt (std::ceil (std::log (1.0 / k1) / std::log (1.0 / k)));
}
else if (type == 1 || type == 2)
{
N = roundToInt (std::ceil (std::acosh (1.0 / k1) / std::acosh (1.0 / k)));
}
else
{
double K, Kp, K1, K1p;
SpecialFunctions::ellipticIntegralK (k, K, Kp);
SpecialFunctions::ellipticIntegralK (k1, K1, K1p);
N = roundToInt (std::ceil ((K1p * K) / (K1 * Kp)));
}
const int r = N % 2;
const int L = (N - r) / 2;
const double H0 = (type == 1 || type == 3) ? std::pow (Gp, 1.0 - r) : 1.0;
Array<Complex<double>> pa, za;
Complex<double> j (0, 1);
if (type == 0)
{
if (r == 1)
pa.add (-omegap * std::pow (epsp, -1.0 / (double) N));
for (int i = 1; i <= L; ++i)
{
auto ui = (2 * i - 1.0) / (double) N;
pa.add (omegap * std::pow (epsp, -1.0 / (double) N) * j * exp (ui * halfPi * j));
}
}
else if (type == 1)
{
auto v0 = std::asinh (1.0 / epsp) / (N * halfPi);
if (r == 1)
pa.add (-omegap * std::sinh (v0 * halfPi));
for (int i = 1; i <= L; ++i)
{
auto ui = (2 * i - 1.0) / (double) N;
pa.add (omegap * j * std::cos ((ui - j * v0) * halfPi));
}
}
else if (type == 2)
{
auto v0 = std::asinh (epss) / (N * halfPi);
if (r == 1)
pa.add (-1.0 / (k / omegap * std::sinh (v0 * halfPi)));
for (int i = 1; i <= L; ++i)
{
auto ui = (2 * i - 1.0) / (double) N;
pa.add (1.0 / (k / omegap * j * std::cos ((ui - j * v0) * halfPi)));
za.add (1.0 / (k / omegap * j * std::cos (ui * halfPi)));
}
}
else
{
auto v0 = -j * (SpecialFunctions::asne (j / epsp, k1) / (double) N);
if (r == 1)
pa.add (omegap * j * SpecialFunctions::sne (j * v0, k));
for (int i = 1; i <= L; ++i)
{
auto ui = (2 * i - 1.0) / (double) N;
auto zetai = SpecialFunctions::cde (ui, k);
pa.add (omegap * j * SpecialFunctions::cde (ui - j * v0, k));
za.add (omegap * j / (k * zetai));
}
}
Array<Complex<double>> p, z, g;
if (r == 1)
{
p.add ((1.0 + pa[0]) / (1.0 - pa[0]));
g.add (0.5 * (1.0 - p[0]));
}
for (int i = 0; i < L; ++i)
{
p.add ((1.0 + pa[i + r]) / (1.0 - pa[i + r]));
z.add (za.size() == 0 ? -1.0 : (1.0 + za[i]) / (1.0 - za[i]));
g.add ((1.0 - p[i + r]) / (1.0 - z[i]));
}
ReferenceCountedArray<IIR::Coefficients<FloatType>> cascadedCoefficients;
if (r == 1)
{
auto b0 = static_cast<FloatType> (H0 * std::real (g[0]));
auto b1 = b0;
auto a1 = static_cast<FloatType> (-std::real (p[0]));
cascadedCoefficients.add (new IIR::Coefficients<FloatType> (b0, b1, 1.0f, a1));
}
for (int i = 0; i < L; ++i)
{
auto gain = std::pow (std::abs (g[i + r]), 2.0);
auto b0 = static_cast<FloatType> (gain);
auto b1 = static_cast<FloatType> (std::real (-z[i] - std::conj (z[i])) * gain);
auto b2 = static_cast<FloatType> (std::real ( z[i] * std::conj (z[i])) * gain);
auto a1 = static_cast<FloatType> (std::real (-p[i+r] - std::conj (p[i + r])));
auto a2 = static_cast<FloatType> (std::real ( p[i+r] * std::conj (p[i + r])));
cascadedCoefficients.add (new IIR::Coefficients<FloatType> (b0, b1, b2, 1, a1, a2));
}
return cascadedCoefficients;
}
template <typename FloatType>
ReferenceCountedArray<IIR::Coefficients<FloatType>>
FilterDesign<FloatType>::designIIRLowpassHighOrderButterworthMethod (FloatType frequency,
double sampleRate, int order)
{
jassert (sampleRate > 0);
jassert (frequency > 0 && frequency <= sampleRate * 0.5);
jassert (order > 0);
ReferenceCountedArray<IIR::Coefficients<FloatType>> arrayFilters;
if (order % 2 == 1)
{
arrayFilters.add (*IIR::Coefficients<FloatType>::makeFirstOrderLowPass (sampleRate, frequency));
for (int i = 0; i < order / 2; ++i)
{
auto Q = 1.0 / (2.0 * std::cos ((i + 1.0) * MathConstants<double>::pi / order));
arrayFilters.add (*IIR::Coefficients<FloatType>::makeLowPass (sampleRate, frequency,
static_cast<FloatType> (Q)));
}
}
else
{
for (int i = 0; i < order / 2; ++i)
{
auto Q = 1.0 / (2.0 * std::cos ((2.0 * i + 1.0) * MathConstants<double>::pi / (order * 2.0)));
arrayFilters.add (*IIR::Coefficients<FloatType>::makeLowPass (sampleRate, frequency,
static_cast<FloatType> (Q)));
}
}
return arrayFilters;
}
template <typename FloatType>
ReferenceCountedArray<IIR::Coefficients<FloatType>>
FilterDesign<FloatType>::designIIRHighpassHighOrderButterworthMethod (FloatType frequency,
double sampleRate, int order)
{
jassert (sampleRate > 0);
jassert (frequency > 0 && frequency <= sampleRate * 0.5);
jassert (order > 0);
ReferenceCountedArray<IIR::Coefficients<FloatType>> arrayFilters;
if (order % 2 == 1)
{
arrayFilters.add (*IIR::Coefficients<FloatType>::makeFirstOrderHighPass (sampleRate, frequency));
for (int i = 0; i < order / 2; ++i)
{
auto Q = 1.0 / (2.0 * std::cos ((i + 1.0) * MathConstants<double>::pi / order));
arrayFilters.add (*IIR::Coefficients<FloatType>::makeHighPass (sampleRate, frequency,
static_cast<FloatType> (Q)));
}
}
else
{
for (int i = 0; i < order / 2; ++i)
{
auto Q = 1.0 / (2.0 * std::cos ((2.0 * i + 1.0) * MathConstants<double>::pi / (order * 2.0)));
arrayFilters.add (*IIR::Coefficients<FloatType>::makeHighPass (sampleRate, frequency,
static_cast<FloatType> (Q)));
}
}
return arrayFilters;
}
template <typename FloatType>
typename FilterDesign<FloatType>::IIRPolyphaseAllpassStructure
FilterDesign<FloatType>::designIIRLowpassHalfBandPolyphaseAllpassMethod (FloatType normalisedTransitionWidth,
FloatType stopbandAmplitudedB)
{
jassert (normalisedTransitionWidth > 0 && normalisedTransitionWidth <= 0.5);
jassert (stopbandAmplitudedB > -300 && stopbandAmplitudedB < -10);
const double wt = MathConstants<double>::twoPi * normalisedTransitionWidth;
const double ds = Decibels::decibelsToGain (stopbandAmplitudedB, static_cast<FloatType> (-300.0));
auto k = std::pow (std::tan ((MathConstants<double>::pi - wt) / 4), 2.0);
auto kp = std::sqrt (1.0 - k * k);
auto e = (1 - std::sqrt (kp)) / (1 + std::sqrt (kp)) * 0.5;
auto q = e + 2 * std::pow (e, 5.0) + 15 * std::pow (e, 9.0) + 150 * std::pow (e, 13.0);
auto k1 = ds * ds / (1 - ds * ds);
int n = roundToInt (std::ceil (std::log (k1 * k1 / 16) / std::log (q)));
if (n % 2 == 0)
++n;
if (n == 1)
n = 3;
auto q1 = std::pow (q, (double) n);
k1 = 4 * std::sqrt (q1);
const int N = (n - 1) / 2;
Array<double> ai;
for (int i = 1; i <= N; ++i)
{
double num = 0.0;
double delta = 1.0;
int m = 0;
while (std::abs (delta) > 1e-100)
{
delta = std::pow (-1, m) * std::pow (q, m * (m + 1))
* std::sin ((2 * m + 1) * MathConstants<double>::pi * i / (double) n);
num += delta;
m++;
}
num *= 2 * std::pow (q, 0.25);
double den = 0.0;
delta = 1.0;
m = 1;
while (std::abs (delta) > 1e-100)
{
delta = std::pow (-1, m) * std::pow (q, m * m)
* std::cos (m * MathConstants<double>::twoPi * i / (double) n);
den += delta;
++m;
}
den = 1 + 2 * den;
auto wi = num / den;
auto api = std::sqrt ((1 - wi * wi * k) * (1 - wi * wi / k)) / (1 + wi * wi);
ai.add ((1 - api) / (1 + api));
}
IIRPolyphaseAllpassStructure structure;
for (int i = 0; i < N; i += 2)
structure.directPath.add (new IIR::Coefficients<FloatType> (static_cast<FloatType> (ai[i]),
0, 1, 1, 0, static_cast<FloatType> (ai[i])));
structure.delayedPath.add (new IIR::Coefficients<FloatType> (0, 1, 1, 0));
for (int i = 1; i < N; i += 2)
structure.delayedPath.add (new IIR::Coefficients<FloatType> (static_cast<FloatType> (ai[i]),
0, 1, 1, 0, static_cast<FloatType> (ai[i])));
structure.alpha.addArray (ai);
return structure;
}
template struct FilterDesign<float>;
template struct FilterDesign<double>;
} // namespace juce::dsp
juce::Array
Definition juce_Array.h:56
juce::Array::setUnchecked
void setUnchecked(int indexToChange, ParameterType newValue)
Definition juce_Array.h:568
juce::Array::addArray
void addArray(const Type *elementsToAdd, int numElementsToAdd)
Definition juce_Array.h:583
juce::Array::add
void add(const ElementType &newElement)
Definition juce_Array.h:418
juce::Array::resize
void resize(int targetNumItems)
Definition juce_Array.h:670
juce::Decibels::decibelsToGain
static Type decibelsToGain(Type decibels, Type minusInfinityDb=Type(defaultMinusInfinitydB))
Definition juce_Decibels.h:42
juce::ReferenceCountedArray
Definition juce_ReferenceCountedArray.h:51
juce::ReferenceCountedArray::add
ObjectClass * add(ObjectClass *newObject)
Definition juce_ReferenceCountedArray.h:354
juce::dsp::Matrix
Definition juce_Matrix.h:39
juce::dsp::Matrix::hankel
static Matrix hankel(const Matrix &vector, size_t size, size_t offset=0)
Definition juce_Matrix.cpp:59
juce::dsp::Matrix::identity
static Matrix identity(size_t size)
Definition juce_Matrix.cpp:30
juce::dsp::Matrix::toeplitz
static Matrix toeplitz(const Matrix &vector, size_t size)
Definition juce_Matrix.cpp:41
juce::dsp::WindowingFunction
Definition juce_Windowing.h:41
juce::dsp::WindowingFunction::multiplyWithWindowingTable
void multiplyWithWindowingTable(FloatType *samples, size_t size) const noexcept
Definition juce_Windowing.cpp:167
juce::MathConstants::halfPi
static constexpr FloatType halfPi
Definition juce_MathsFunctions.h:147
juce::MathConstants::twoPi
static constexpr FloatType twoPi
Definition juce_MathsFunctions.h:144
juce::MathConstants::pi
static constexpr FloatType pi
Definition juce_MathsFunctions.h:141
juce::dsp::FIR::Coefficients
Definition juce_FIRFilter.h:214
juce::dsp::FIR::Coefficients::Ptr
ReferenceCountedObjectPtr< Coefficients > Ptr
Definition juce_FIRFilter.h:233
juce::dsp::FilterDesign
Definition juce_FilterDesign.h:40
juce::dsp::FilterDesign::designIIRLowpassHighOrderButterworthMethod
static ReferenceCountedArray< IIRCoefficients > designIIRLowpassHighOrderButterworthMethod(FloatType frequency, double sampleRate, FloatType normalisedTransitionWidth, FloatType passbandAmplitudedB, FloatType stopbandAmplitudedB)
Definition juce_FilterDesign.cpp:337
juce::dsp::FilterDesign::designFIRLowpassLeastSquaresMethod
static FIRCoefficientsPtr designFIRLowpassLeastSquaresMethod(FloatType frequency, double sampleRate, size_t order, FloatType normalisedTransitionWidth, FloatType stopBandWeight)
Definition juce_FilterDesign.cpp:126
juce::dsp::FilterDesign::designFIRLowpassKaiserMethod
static FIRCoefficientsPtr designFIRLowpassKaiserMethod(FloatType frequency, double sampleRate, FloatType normalisedTransitionWidth, FloatType amplitudedB)
Definition juce_FilterDesign.cpp:65
juce::dsp::FilterDesign::designIIRLowpassHalfBandPolyphaseAllpassMethod
static IIRPolyphaseAllpassStructure designIIRLowpassHalfBandPolyphaseAllpassMethod(FloatType normalisedTransitionWidth, FloatType stopbandAmplitudedB)
Definition juce_FilterDesign.cpp:612
juce::dsp::FilterDesign::designIIRLowpassHighOrderChebyshev1Method
static ReferenceCountedArray< IIRCoefficients > designIIRLowpassHighOrderChebyshev1Method(FloatType frequency, double sampleRate, FloatType normalisedTransitionWidth, FloatType passbandAmplitudedB, FloatType stopbandAmplitudedB)
Definition juce_FilterDesign.cpp:348
juce::dsp::FilterDesign::designIIRHighpassHighOrderButterworthMethod
static ReferenceCountedArray< IIRCoefficients > designIIRHighpassHighOrderButterworthMethod(FloatType frequency, double sampleRate, int order)
Definition juce_FilterDesign.cpp:577
juce::dsp::FilterDesign::designFIRLowpassTransitionMethod
static FIRCoefficientsPtr designFIRLowpassTransitionMethod(FloatType frequency, double sampleRate, size_t order, FloatType normalisedTransitionWidth, FloatType spline)
Definition juce_FilterDesign.cpp:93
juce::dsp::FilterDesign::designIIRLowpassHighOrderChebyshev2Method
static ReferenceCountedArray< IIRCoefficients > designIIRLowpassHighOrderChebyshev2Method(FloatType frequency, double sampleRate, FloatType normalisedTransitionWidth, FloatType passbandAmplitudedB, FloatType stopbandAmplitudedB)
Definition juce_FilterDesign.cpp:359
juce::dsp::FilterDesign::designIIRLowpassHighOrderEllipticMethod
static ReferenceCountedArray< IIRCoefficients > designIIRLowpassHighOrderEllipticMethod(FloatType frequency, double sampleRate, FloatType normalisedTransitionWidth, FloatType passbandAmplitudedB, FloatType stopbandAmplitudedB)
Definition juce_FilterDesign.cpp:370
juce::dsp::FilterDesign::designFIRLowpassHalfBandEquirippleMethod
static FIRCoefficientsPtr designFIRLowpassHalfBandEquirippleMethod(FloatType normalisedTransitionWidth, FloatType amplitudedB)
Definition juce_FilterDesign.cpp:235
juce::dsp::FilterDesign::designFIRLowpassWindowMethod
static FIRCoefficientsPtr designFIRLowpassWindowMethod(FloatType frequency, double sampleRate, size_t order, WindowingMethod type, FloatType beta=static_cast< FloatType >(2))
Definition juce_FilterDesign.cpp:31
juce::dsp::FilterDesign::IIRPolyphaseAllpassStructure
Definition juce_FilterDesign.h:257
juce::dsp::IIR::Coefficients
Definition juce_dsp/processors/juce_IIRFilter.h:124
juce::dsp::IIR::Coefficients::makeFirstOrderHighPass
static Ptr makeFirstOrderHighPass(double sampleRate, NumericType frequency)
Definition juce_dsp/processors/juce_IIRFilter.cpp:314
juce::dsp::IIR::Coefficients::makeHighPass
static Ptr makeHighPass(double sampleRate, NumericType frequency)
Definition juce_dsp/processors/juce_IIRFilter.cpp:343
juce::dsp::IIR::Coefficients::makeLowPass
static Ptr makeLowPass(double sampleRate, NumericType frequency)
Definition juce_dsp/processors/juce_IIRFilter.cpp:328
juce::dsp::IIR::Coefficients::makeFirstOrderLowPass
static Ptr makeFirstOrderLowPass(double sampleRate, NumericType frequency)
Definition juce_dsp/processors/juce_IIRFilter.cpp:307
juce::dsp::SpecialFunctions::sne
static Complex< double > sne(Complex< double > u, double k) noexcept
Definition juce_SpecialFunctions.cpp:95
juce::dsp::SpecialFunctions::cde
static Complex< double > cde(Complex< double > u, double k) noexcept
Definition juce_SpecialFunctions.cpp:72
juce::dsp::SpecialFunctions::asne
static Complex< double > asne(Complex< double > w, double k) noexcept
Definition juce_SpecialFunctions.cpp:118
juce::dsp::SpecialFunctions::ellipticIntegralK
static void ellipticIntegralK(double k, double &K, double &Kp) noexcept
Definition juce_SpecialFunctions.cpp:49