From d44112514d84614af6c16ec48d9450706da3b335 Mon Sep 17 00:00:00 2001 From: Chris Robinson Date: Tue, 7 Nov 2023 15:04:51 -0800 Subject: Precompute a value used multiple times --- common/polyphase_resampler.cpp | 28 +++++++++------------------- 1 file changed, 9 insertions(+), 19 deletions(-) (limited to 'common/polyphase_resampler.cpp') diff --git a/common/polyphase_resampler.cpp b/common/polyphase_resampler.cpp index 14f7e40d..23102b9d 100644 --- a/common/polyphase_resampler.cpp +++ b/common/polyphase_resampler.cpp @@ -3,6 +3,7 @@ #include #include +#include #include "alnumbers.h" #include "opthelpers.h" @@ -67,23 +68,11 @@ constexpr double BesselI_0(const double x) * * k = 2 i / M - 1, where 0 <= i <= M. */ -double Kaiser(const double b, const double k) +double Kaiser(const double beta, const double k, const double besseli_0_beta) { if(!(k >= -1.0 && k <= 1.0)) return 0.0; - return BesselI_0(b * std::sqrt(1.0 - k*k)) / BesselI_0(b); -} - -// Calculates the greatest common divisor of a and b. -constexpr uint Gcd(uint x, uint y) -{ - while(y > 0) - { - const uint z{y}; - y = x % y; - x = z; - } - return x; + return BesselI_0(beta * std::sqrt(1.0 - k*k)) / besseli_0_beta; } /* Calculates the size (order) of the Kaiser window. Rejection is in dB and @@ -124,11 +113,11 @@ constexpr double CalcKaiserBeta(const double rejection) * p -- gain compensation factor when sampling * f_t -- normalized center frequency (or cutoff; 0.5 is nyquist) */ -double SincFilter(const uint l, const double b, const double gain, const double cutoff, - const uint i) +double SincFilter(const uint l, const double beta, const double besseli_0_beta, const double gain, + const double cutoff, const uint i) { const double x{static_cast(i) - l}; - return Kaiser(b, x / l) * 2.0 * gain * cutoff * Sinc(2.0 * cutoff * x); + return Kaiser(beta, x/l, besseli_0_beta) * 2.0 * gain * cutoff * Sinc(2.0 * cutoff * x); } } // namespace @@ -137,7 +126,7 @@ double SincFilter(const uint l, const double b, const double gain, const double // that's used to cut frequencies above the destination nyquist. void PPhaseResampler::init(const uint srcRate, const uint dstRate) { - const uint gcd{Gcd(srcRate, dstRate)}; + const uint gcd{std::gcd(srcRate, dstRate)}; mP = dstRate / gcd; mQ = srcRate / gcd; @@ -160,11 +149,12 @@ void PPhaseResampler::init(const uint srcRate, const uint dstRate) // calculating the left offset to avoid increasing the transition width. const uint l{(CalcKaiserOrder(180.0, width)+1) / 2}; const double beta{CalcKaiserBeta(180.0)}; + const double besseli_0_beta{BesselI_0(beta)}; mM = l*2 + 1; mL = l; mF.resize(mM); for(uint i{0};i < mM;i++) - mF[i] = SincFilter(l, beta, mP, cutoff, i); + mF[i] = SincFilter(l, beta, besseli_0_beta, mP, cutoff, i); } // Perform the upsample-filter-downsample resampling operation using a -- cgit v1.2.3 From 33cd77b81b1dce9a5b55ec7e215315c500f9d0bf Mon Sep 17 00:00:00 2001 From: Chris Robinson Date: Wed, 8 Nov 2023 13:01:14 -0800 Subject: Use the C++ standard's regular modified Bessel function --- common/polyphase_resampler.cpp | 31 ++----------------------------- core/bsinc_tables.cpp | 37 ++++--------------------------------- 2 files changed, 6 insertions(+), 62 deletions(-) (limited to 'common/polyphase_resampler.cpp') diff --git a/common/polyphase_resampler.cpp b/common/polyphase_resampler.cpp index 23102b9d..f1e9eaf3 100644 --- a/common/polyphase_resampler.cpp +++ b/common/polyphase_resampler.cpp @@ -27,33 +27,6 @@ double Sinc(const double x) return std::sin(al::numbers::pi*x) / (al::numbers::pi*x); } -/* The zero-order modified Bessel function of the first kind, used for the - * Kaiser window. - * - * I_0(x) = sum_{k=0}^inf (1 / k!)^2 (x / 2)^(2 k) - * = sum_{k=0}^inf ((x / 2)^k / k!)^2 - */ -constexpr double BesselI_0(const double x) -{ - // Start at k=1 since k=0 is trivial. - const double x2{x/2.0}; - double term{1.0}; - double sum{1.0}; - int k{1}; - - // Let the integration converge until the term of the sum is no longer - // significant. - double last_sum{}; - do { - const double y{x2 / k}; - ++k; - last_sum = sum; - term *= y * y; - sum += term; - } while(sum != last_sum); - return sum; -} - /* Calculate a Kaiser window from the given beta value and a normalized k * [-1, 1]. * @@ -72,7 +45,7 @@ double Kaiser(const double beta, const double k, const double besseli_0_beta) { if(!(k >= -1.0 && k <= 1.0)) return 0.0; - return BesselI_0(beta * std::sqrt(1.0 - k*k)) / besseli_0_beta; + return std::cyl_bessel_i(0, beta * std::sqrt(1.0 - k*k)) / besseli_0_beta; } /* Calculates the size (order) of the Kaiser window. Rejection is in dB and @@ -149,7 +122,7 @@ void PPhaseResampler::init(const uint srcRate, const uint dstRate) // calculating the left offset to avoid increasing the transition width. const uint l{(CalcKaiserOrder(180.0, width)+1) / 2}; const double beta{CalcKaiserBeta(180.0)}; - const double besseli_0_beta{BesselI_0(beta)}; + const double besseli_0_beta{std::cyl_bessel_i(0, beta)}; mM = l*2 + 1; mL = l; mF.resize(mM); diff --git a/core/bsinc_tables.cpp b/core/bsinc_tables.cpp index c734b0e9..74684fa5 100644 --- a/core/bsinc_tables.cpp +++ b/core/bsinc_tables.cpp @@ -33,35 +33,6 @@ constexpr double Sinc(const double x) return std::sin(al::numbers::pi*x) / (al::numbers::pi*x); } -/* The zero-order modified Bessel function of the first kind, used for the - * Kaiser window. - * - * I_0(x) = sum_{k=0}^inf (1 / k!)^2 (x / 2)^(2 k) - * = sum_{k=0}^inf ((x / 2)^k / k!)^2 - */ -constexpr double BesselI_0(const double x) noexcept -{ - /* Start at k=1 since k=0 is trivial. */ - const double x2{x / 2.0}; - double term{1.0}; - double sum{1.0}; - double last_sum{}; - int k{1}; - - /* Let the integration converge until the term of the sum is no longer - * significant. - */ - do { - const double y{x2 / k}; - ++k; - last_sum = sum; - term *= y * y; - sum += term; - } while(sum != last_sum); - - return sum; -} - /* Calculate a Kaiser window from the given beta value and a normalized k * [-1, 1]. * @@ -80,7 +51,7 @@ constexpr double Kaiser(const double beta, const double k, const double besseli_ { if(!(k >= -1.0 && k <= 1.0)) return 0.0; - return BesselI_0(beta * std::sqrt(1.0 - k*k)) / besseli_0_beta; + return std::cyl_bessel_i(0, beta * std::sqrt(1.0 - k*k)) / besseli_0_beta; } /* Calculates the (normalized frequency) transition width of the Kaiser window. @@ -109,7 +80,6 @@ struct BSincHeader { double width{}; double beta{}; double scaleBase{}; - double besseli_0_beta{}; uint a[BSincScaleCount]{}; uint total_size{}; @@ -119,7 +89,6 @@ struct BSincHeader { width = CalcKaiserWidth(Rejection, Order); beta = CalcKaiserBeta(Rejection); scaleBase = width / 2.0; - besseli_0_beta = BesselI_0(beta); uint num_points{Order+1}; for(uint si{0};si < BSincScaleCount;++si) @@ -154,6 +123,8 @@ struct BSincFilterArray { using filter_type = double[BSincPhaseCount+1][BSincPointsMax]; auto filter = std::make_unique(BSincScaleCount); + const double besseli_0_beta{std::cyl_bessel_i(0, hdr.beta)}; + /* Calculate the Kaiser-windowed Sinc filter coefficients for each * scale and phase index. */ @@ -176,7 +147,7 @@ struct BSincFilterArray { for(uint i{0};i < m;++i) { const double x{i - phase}; - filter[si][pi][o+i] = Kaiser(hdr.beta, x/l, hdr.besseli_0_beta) * cutoff * + filter[si][pi][o+i] = Kaiser(hdr.beta, x/l, besseli_0_beta) * cutoff * Sinc(cutoff*x); } } -- cgit v1.2.3 From d6d572df66eb2b5beeb3da7f3b11a328500c33c3 Mon Sep 17 00:00:00 2001 From: Chris Robinson Date: Sat, 18 Nov 2023 22:11:58 -0800 Subject: Handle systems that don't support std::cyl_bessel_i --- common/polyphase_resampler.cpp | 47 ++++++++++++++++++++++++++++++++++++++++-- core/bsinc_tables.cpp | 46 +++++++++++++++++++++++++++++++++++++++-- 2 files changed, 89 insertions(+), 4 deletions(-) (limited to 'common/polyphase_resampler.cpp') diff --git a/common/polyphase_resampler.cpp b/common/polyphase_resampler.cpp index f1e9eaf3..9c569b3a 100644 --- a/common/polyphase_resampler.cpp +++ b/common/polyphase_resampler.cpp @@ -4,6 +4,7 @@ #include #include #include +#include #include "alnumbers.h" #include "opthelpers.h" @@ -15,6 +16,48 @@ constexpr double Epsilon{1e-9}; using uint = unsigned int; +#if __cpp_lib_math_special_functions >= 201603L +using std::cyl_bessel_i; + +#else + +/* The zero-order modified Bessel function of the first kind, used for the + * Kaiser window. + * + * I_0(x) = sum_{k=0}^inf (1 / k!)^2 (x / 2)^(2 k) + * = sum_{k=0}^inf ((x / 2)^k / k!)^2 + * + * This implementation only handles nu = 0, and isn't the most precise (it + * starts with the largest value and accumulates successively smaller values, + * compounding the rounding and precision error), but it's good enough. + */ +template +U cyl_bessel_i(T nu, U x) +{ + if(nu != T{0}) + throw std::runtime_error{"cyl_bessel_i: nu != 0"}; + + /* Start at k=1 since k=0 is trivial. */ + const double x2{x/2.0}; + double term{1.0}; + double sum{1.0}; + int k{1}; + + /* Let the integration converge until the term of the sum is no longer + * significant. + */ + double last_sum{}; + do { + const double y{x2 / k}; + ++k; + last_sum = sum; + term *= y * y; + sum += term; + } while(sum != last_sum); + return static_cast(sum); +} +#endif + /* This is the normalized cardinal sine (sinc) function. * * sinc(x) = { 1, x = 0 @@ -45,7 +88,7 @@ double Kaiser(const double beta, const double k, const double besseli_0_beta) { if(!(k >= -1.0 && k <= 1.0)) return 0.0; - return std::cyl_bessel_i(0, beta * std::sqrt(1.0 - k*k)) / besseli_0_beta; + return cyl_bessel_i(0, beta * std::sqrt(1.0 - k*k)) / besseli_0_beta; } /* Calculates the size (order) of the Kaiser window. Rejection is in dB and @@ -122,7 +165,7 @@ void PPhaseResampler::init(const uint srcRate, const uint dstRate) // calculating the left offset to avoid increasing the transition width. const uint l{(CalcKaiserOrder(180.0, width)+1) / 2}; const double beta{CalcKaiserBeta(180.0)}; - const double besseli_0_beta{std::cyl_bessel_i(0, beta)}; + const double besseli_0_beta{cyl_bessel_i(0, beta)}; mM = l*2 + 1; mL = l; mF.resize(mM); diff --git a/core/bsinc_tables.cpp b/core/bsinc_tables.cpp index 74684fa5..41102e9a 100644 --- a/core/bsinc_tables.cpp +++ b/core/bsinc_tables.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include "alnumbers.h" #include "alnumeric.h" @@ -19,6 +20,47 @@ namespace { using uint = unsigned int; +#if __cpp_lib_math_special_functions >= 201603L +using std::cyl_bessel_i; + +#else + +/* The zero-order modified Bessel function of the first kind, used for the + * Kaiser window. + * + * I_0(x) = sum_{k=0}^inf (1 / k!)^2 (x / 2)^(2 k) + * = sum_{k=0}^inf ((x / 2)^k / k!)^2 + * + * This implementation only handles nu = 0, and isn't the most precise (it + * starts with the largest value and accumulates successively smaller values, + * compounding the rounding and precision error), but it's good enough. + */ +template +U cyl_bessel_i(T nu, U x) +{ + if(nu != T{0}) + throw std::runtime_error{"cyl_bessel_i: nu != 0"}; + + /* Start at k=1 since k=0 is trivial. */ + const double x2{x/2.0}; + double term{1.0}; + double sum{1.0}; + int k{1}; + + /* Let the integration converge until the term of the sum is no longer + * significant. + */ + double last_sum{}; + do { + const double y{x2 / k}; + ++k; + last_sum = sum; + term *= y * y; + sum += term; + } while(sum != last_sum); + return static_cast(sum); +} +#endif /* This is the normalized cardinal sine (sinc) function. * @@ -51,7 +93,7 @@ constexpr double Kaiser(const double beta, const double k, const double besseli_ { if(!(k >= -1.0 && k <= 1.0)) return 0.0; - return std::cyl_bessel_i(0, beta * std::sqrt(1.0 - k*k)) / besseli_0_beta; + return cyl_bessel_i(0, beta * std::sqrt(1.0 - k*k)) / besseli_0_beta; } /* Calculates the (normalized frequency) transition width of the Kaiser window. @@ -123,7 +165,7 @@ struct BSincFilterArray { using filter_type = double[BSincPhaseCount+1][BSincPointsMax]; auto filter = std::make_unique(BSincScaleCount); - const double besseli_0_beta{std::cyl_bessel_i(0, hdr.beta)}; + const double besseli_0_beta{cyl_bessel_i(0, hdr.beta)}; /* Calculate the Kaiser-windowed Sinc filter coefficients for each * scale and phase index. -- cgit v1.2.3