1 files changed, 41 insertions, 35 deletions

diff --git a/common/polyphase_resampler.cpp b/common/polyphase_resampler.cpp
index 14f7e40d..9c569b3a 100644
--- a/common/polyphase_resampler.cpp
+++ b/common/polyphase_resampler.cpp
@@ -3,6 +3,8 @@
 
 #include <algorithm>
 #include <cmath>
+#include <numeric>
+#include <stdexcept>
 
 #include "alnumbers.h"
 #include "opthelpers.h"
@@ -14,34 +16,36 @@ constexpr double Epsilon{1e-9};
 
 using uint = unsigned int;
 
-/* This is the normalized cardinal sine (sinc) function.
- *
- *   sinc(x) = { 1,                   x = 0
- *             { sin(pi x) / (pi x),  otherwise.
- */
-double Sinc(const double x)
-{
-    if(std::abs(x) < Epsilon) UNLIKELY
-        return 1.0;
-    return std::sin(al::numbers::pi*x) / (al::numbers::pi*x);
-}
+#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.
  */
-constexpr double BesselI_0(const double x)
+template<typename T, typename U>
+U cyl_bessel_i(T nu, U x)
 {
-    // Start at k=1 since k=0 is trivial.
+    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.
+    /* Let the integration converge until the term of the sum is no longer
+     * significant.
+     */
     double last_sum{};
     do {
         const double y{x2 / k};
@@ -50,7 +54,20 @@ constexpr double BesselI_0(const double x)
         term *= y * y;
         sum += term;
     } while(sum != last_sum);
-    return sum;
+    return static_cast<U>(sum);
+}
+#endif
+
+/* This is the normalized cardinal sine (sinc) function.
+ *
+ *   sinc(x) = { 1,                   x = 0
+ *             { sin(pi x) / (pi x),  otherwise.
+ */
+double Sinc(const double x)
+{
+    if(std::abs(x) < Epsilon) UNLIKELY
+        return 1.0;
+    return std::sin(al::numbers::pi*x) / (al::numbers::pi*x);
 }
 
 /* Calculate a Kaiser window from the given beta value and a normalized k
@@ -67,23 +84,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 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
@@ -124,11 +129,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<double>(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 +142,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 +165,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{cyl_bessel_i(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