core/bsinc_tables.cpp

   1
   2 #include "bsinc_tables.h"
   3
   4 #include <algorithm>
   5 #include <array>
   6 #include <cassert>
   7 #include <cmath>
   8 #include <cstddef>
   9 #include <limits>
  10 #include <stdexcept>
  11 #include <vector>
  12
  13 #include "alnumbers.h"
  14 #include "alnumeric.h"
  15 #include "alspan.h"
  16 #include "bsinc_defs.h"
  17 #include "resampler_limits.h"
  18
  19
  20 namespace {
  21
  22 using uint = unsigned int;
  23
  24 #if __cpp_lib_math_special_functions >= 201603L
  25 using std::cyl_bessel_i;
  26
  27 #else
  28
  29 /* The zero-order modified Bessel function of the first kind, used for the
  30  * Kaiser window.
  31  *
  32  *   I_0(x) = sum_{k=0}^inf (1 / k!)^2 (x / 2)^(2 k)
  33  *          = sum_{k=0}^inf ((x / 2)^k / k!)^2
  34  *
  35  * This implementation only handles nu = 0, and isn't the most precise (it
  36  * starts with the largest value and accumulates successively smaller values,
  37  * compounding the rounding and precision error), but it's good enough.
  38  */
  39 template<typename T, typename U>
  40 U cyl_bessel_i(T nu, U x)
  41 {
  42     if(nu != T{0})
  43         throw std::runtime_error{"cyl_bessel_i: nu != 0"};
  44
  45     /* Start at k=1 since k=0 is trivial. */
  46     const double x2{x/2.0};
  47     double term{1.0};
  48     double sum{1.0};
  49     int k{1};
  50
  51     /* Let the integration converge until the term of the sum is no longer
  52      * significant.
  53      */
  54     double last_sum{};
  55     do {
  56         const double y{x2 / k};
  57         ++k;
  58         last_sum = sum;
  59         term *= y * y;
  60         sum += term;
  61     } while(sum != last_sum);
  62     return static_cast<U>(sum);
  63 }
  64 #endif
  65
  66 /* This is the normalized cardinal sine (sinc) function.
  67  *
  68  *   sinc(x) = { 1,                   x = 0
  69  *             { sin(pi x) / (pi x),  otherwise.
  70  */
  71 constexpr double Sinc(const double x)
  72 {
  73     constexpr double epsilon{std::numeric_limits<double>::epsilon()};
  74     if(!(x > epsilon || x < -epsilon))
  75         return 1.0;
  76     return std::sin(al::numbers::pi*x) / (al::numbers::pi*x);
  77 }
  78
  79 /* Calculate a Kaiser window from the given beta value and a normalized k
  80  * [-1, 1].
  81  *
  82  *   w(k) = { I_0(B sqrt(1 - k^2)) / I_0(B),  -1 <= k <= 1
  83  *          { 0,                              elsewhere.
  84  *
  85  * Where k can be calculated as:
  86  *
  87  *   k = i / l,         where -l <= i <= l.
  88  *
  89  * or:
  90  *
  91  *   k = 2 i / M - 1,   where 0 <= i <= M.
  92  */
  93 constexpr double Kaiser(const double beta, const double k, const double besseli_0_beta)
  94 {
  95     if(!(k >= -1.0 && k <= 1.0))
  96         return 0.0;
  97     return cyl_bessel_i(0, beta * std::sqrt(1.0 - k*k)) / besseli_0_beta;
  98 }
  99
 100 /* Calculates the (normalized frequency) transition width of the Kaiser window.
 101  * Rejection is in dB.
 102  */
 103 constexpr double CalcKaiserWidth(const double rejection, const uint order) noexcept
 104 {
 105     if(rejection > 21.19)
 106         return (rejection - 7.95) / (2.285 * al::numbers::pi*2.0 * order);
 107     /* This enforces a minimum rejection of just above 21.18dB */
 108     return 5.79 / (al::numbers::pi*2.0 * order);
 109 }
 110
 111 /* Calculates the beta value of the Kaiser window. Rejection is in dB. */
 112 constexpr double CalcKaiserBeta(const double rejection)
 113 {
 114     if(rejection > 50.0)
 115         return 0.1102 * (rejection-8.7);
 116     if(rejection >= 21.0)
 117         return (0.5842 * std::pow(rejection-21.0, 0.4)) + (0.07886 * (rejection-21.0));
 118     return 0.0;
 119 }
 120
 121
 122 struct BSincHeader {
 123     double width{};
 124     double beta{};
 125     double scaleBase{};
 126
 127     std::array<uint,BSincScaleCount> a{};
 128     uint total_size{};
 129
 130     constexpr BSincHeader(uint Rejection, uint Order) noexcept
 131         : width{CalcKaiserWidth(Rejection, Order)}, beta{CalcKaiserBeta(Rejection)}
 132         , scaleBase{width / 2.0}
 133     {
 134         uint num_points{Order+1};
 135         for(uint si{0};si < BSincScaleCount;++si)
 136         {
 137             const double scale{lerpd(scaleBase, 1.0, (si+1) / double{BSincScaleCount})};
 138             const uint a_{std::min(static_cast<uint>(num_points / 2.0 / scale), num_points)};
 139             const uint m{2 * a_};
 140
 141             a[si] = a_;
 142             total_size += 4 * BSincPhaseCount * ((m+3) & ~3u);
 143         }
 144     }
 145 };
 146
 147 /* 11th and 23rd order filters (12 and 24-point respectively) with a 60dB drop
 148  * at nyquist. Each filter will scale up the order when downsampling, to 23rd
 149  * and 47th order respectively.
 150  */
 151 constexpr BSincHeader bsinc12_hdr{60, 11};
 152 constexpr BSincHeader bsinc24_hdr{60, 23};
 153
 154
 155 template<const BSincHeader &hdr>
 156 struct BSincFilterArray {
 157     alignas(16) std::array<float, hdr.total_size> mTable{};
 158
 159     BSincFilterArray()
 160     {
 161         static constexpr uint BSincPointsMax{(hdr.a[0]*2u + 3u) & ~3u};
 162         static_assert(BSincPointsMax <= MaxResamplerPadding, "MaxResamplerPadding is too small");
 163
 164         using filter_type = std::array<std::array<double,BSincPointsMax>,BSincPhaseCount>;
 165         auto filter = std::vector<filter_type>(BSincScaleCount);
 166
 167         const double besseli_0_beta{cyl_bessel_i(0, hdr.beta)};
 168
 169         /* Calculate the Kaiser-windowed Sinc filter coefficients for each
 170          * scale and phase index.
 171          */
 172         for(uint si{0};si < BSincScaleCount;++si)
 173         {
 174             const uint m{hdr.a[si] * 2};
 175             const size_t o{(BSincPointsMax-m) / 2};
 176             const double scale{lerpd(hdr.scaleBase, 1.0, (si+1) / double{BSincScaleCount})};
 177             const double cutoff{scale - (hdr.scaleBase * std::max(1.0, scale*2.0))};
 178             const auto a = static_cast<double>(hdr.a[si]);
 179             const double l{a - 1.0/BSincPhaseCount};
 180
 181             for(uint pi{0};pi < BSincPhaseCount;++pi)
 182             {
 183                 const double phase{std::floor(l) + (pi/double{BSincPhaseCount})};
 184
 185                 for(uint i{0};i < m;++i)
 186                 {
 187                     const double x{i - phase};
 188                     filter[si][pi][o+i] = Kaiser(hdr.beta, x/l, besseli_0_beta) * cutoff *
 189                         Sinc(cutoff*x);
 190                 }
 191             }
 192         }
 193
 194         size_t idx{0};
 195         for(size_t si{0};si < BSincScaleCount;++si)
 196         {
 197             const size_t m{((hdr.a[si]*2) + 3) & ~3u};
 198             const size_t o{(BSincPointsMax-m) / 2};
 199
 200             /* Write out each phase index's filter and phase delta for this
 201              * quality scale.
 202              */
 203             for(size_t pi{0};pi < BSincPhaseCount;++pi)
 204             {
 205                 for(size_t i{0};i < m;++i)
 206                     mTable[idx++] = static_cast<float>(filter[si][pi][o+i]);
 207
 208                 /* Linear interpolation between phases is simplified by pre-
 209                  * calculating the delta (b - a) in: x = a + f (b - a)
 210                  */
 211                 if(pi < BSincPhaseCount-1)
 212                 {
 213                     for(size_t i{0};i < m;++i)
 214                     {
 215                         const double phDelta{filter[si][pi+1][o+i] - filter[si][pi][o+i]};
 216                         mTable[idx++] = static_cast<float>(phDelta);
 217                     }
 218                 }
 219                 else
 220                 {
 221                     /* The delta target for the last phase index is the first
 222                      * phase index with the coefficients offset by one. The
 223                      * first delta targets 0, as it represents a coefficient
 224                      * for a sample that won't be part of the filter.
 225                      */
 226                     mTable[idx++] = static_cast<float>(0.0 - filter[si][pi][o]);
 227                     for(size_t i{1};i < m;++i)
 228                     {
 229                         const double phDelta{filter[si][0][o+i-1] - filter[si][pi][o+i]};
 230                         mTable[idx++] = static_cast<float>(phDelta);
 231                     }
 232                 }
 233             }
 234
 235             /* Now write out each phase index's scale and phase+scale deltas,
 236              * to complete the bilinear equation for the combination of phase
 237              * and scale.
 238              */
 239             if(si < BSincScaleCount-1)
 240             {
 241                 for(size_t pi{0};pi < BSincPhaseCount;++pi)
 242                 {
 243                     for(size_t i{0};i < m;++i)
 244                     {
 245                         const double scDelta{filter[si+1][pi][o+i] - filter[si][pi][o+i]};
 246                         mTable[idx++] = static_cast<float>(scDelta);
 247                     }
 248
 249                     if(pi < BSincPhaseCount-1)
 250                     {
 251                         for(size_t i{0};i < m;++i)
 252                         {
 253                             const double spDelta{(filter[si+1][pi+1][o+i]-filter[si+1][pi][o+i]) -
 254                                 (filter[si][pi+1][o+i]-filter[si][pi][o+i])};
 255                             mTable[idx++] = static_cast<float>(spDelta);
 256                         }
 257                     }
 258                     else
 259                     {
 260                         mTable[idx++] = static_cast<float>((0.0 - filter[si+1][pi][o]) -
 261                             (0.0 - filter[si][pi][o]));
 262                         for(size_t i{1};i < m;++i)
 263                         {
 264                             const double spDelta{(filter[si+1][0][o+i-1] - filter[si+1][pi][o+i]) -
 265                                 (filter[si][0][o+i-1] - filter[si][pi][o+i])};
 266                             mTable[idx++] = static_cast<float>(spDelta);
 267                         }
 268                     }
 269                 }
 270             }
 271             else
 272             {
 273                 /* The last scale index doesn't have scale-related deltas. */
 274                 for(size_t i{0};i < BSincPhaseCount*m*2;++i)
 275                     mTable[idx++] = 0.0f;
 276             }
 277         }
 278         assert(idx == hdr.total_size);
 279     }
 280
 281     [[nodiscard]] constexpr auto getHeader() const noexcept -> const BSincHeader& { return hdr; }
 282     [[nodiscard]] constexpr auto getTable() const noexcept { return al::span{mTable}; }
 283 };
 284
 285 const BSincFilterArray<bsinc12_hdr> bsinc12_filter{};
 286 const BSincFilterArray<bsinc24_hdr> bsinc24_filter{};
 287
 288 template<typename T>
 289 constexpr BSincTable GenerateBSincTable(const T &filter)
 290 {
 291     BSincTable ret{};
 292     const BSincHeader &hdr = filter.getHeader();
 293     ret.scaleBase = static_cast<float>(hdr.scaleBase);
 294     ret.scaleRange = static_cast<float>(1.0 / (1.0 - hdr.scaleBase));
 295     for(size_t i{0};i < BSincScaleCount;++i)
 296         ret.m[i] = ((hdr.a[i]*2) + 3) & ~3u;
 297     ret.filterOffset[0] = 0;
 298     for(size_t i{1};i < BSincScaleCount;++i)
 299         ret.filterOffset[i] = ret.filterOffset[i-1] + ret.m[i-1]*4*BSincPhaseCount;
 300     ret.Tab = filter.getTable();
 301     return ret;
 302 }
 303
 304 } // namespace
 305
 306 const BSincTable gBSinc12{GenerateBSincTable(bsinc12_filter)};
 307 const BSincTable gBSinc24{GenerateBSincTable(bsinc24_filter)};