alc/bsinc_tables.cpp

   1
   2 #include "bsinc_tables.h"
   3
   4 #include <algorithm>
   5 #include <array>
   6 #include <cassert>
   7 #include <cmath>
   8 #include <limits>
   9 #include <memory>
  10 #include <stdexcept>
  11
  12 #include "math_defs.h"
  13 #include "vector.h"
  14
  15
  16 namespace {
  17
  18 /* The max points includes the doubling for downsampling, so the maximum number
  19  * of base sample points is 24, which is 23rd order.
  20  */
  21 constexpr int BSincPointsMax{BSINC_POINTS_MAX};
  22 constexpr int BSincPointsHalf{BSincPointsMax / 2};
  23
  24 constexpr int BSincPhaseCount{BSINC_PHASE_COUNT};
  25 constexpr int BSincScaleCount{BSINC_SCALE_COUNT};
  26
  27
  28 template<typename T>
  29 constexpr std::enable_if_t<std::is_floating_point<T>::value,T> sqrt(T x)
  30 {
  31     if(!(x >= 0 && x < std::numeric_limits<double>::infinity()))
  32         throw std::domain_error{"Invalid sqrt value"};
  33
  34     T cur{x}, prev{0};
  35     while(cur != prev)
  36     {
  37         prev = cur;
  38         cur = 0.5f*(cur + x/cur);
  39     }
  40     return cur;
  41 }
  42
  43 template<typename T>
  44 constexpr std::enable_if_t<std::is_floating_point<T>::value,T> sin(T x)
  45 {
  46     if(x >= al::MathDefs<T>::Tau())
  47     {
  48         if(!(x < 65536))
  49             throw std::domain_error{"Invalid sin value"};
  50         do {
  51             x -= al::MathDefs<T>::Tau();
  52         } while(x >= al::MathDefs<T>::Tau());
  53     }
  54     else if(x < 0)
  55     {
  56         if(!(x > -65536))
  57             throw std::domain_error{"Invalid sin value"};
  58         do {
  59             x += al::MathDefs<T>::Tau();
  60         } while(x < 0);
  61     }
  62
  63     T prev{x}, n{6};
  64     int i{4}, s{-1};
  65     const T xx{x*x};
  66     T t{xx*x};
  67
  68     T cur{prev + t*s/n};
  69     while(prev != cur)
  70     {
  71         prev = cur;
  72         n *= i*(i+1);
  73         i += 2;
  74         s = -s;
  75         t *= xx;
  76
  77         cur += t*s/n;
  78     }
  79     return cur;
  80 }
  81
  82
  83 /* This is the normalized cardinal sine (sinc) function.
  84  *
  85  *   sinc(x) = { 1,                   x = 0
  86  *             { sin(pi x) / (pi x),  otherwise.
  87  */
  88 constexpr double Sinc(const double x)
  89 {
  90     if(!(x > 1e-15 || x < -1e-15))
  91         return 1.0;
  92     return sin(al::MathDefs<double>::Pi()*x) / (al::MathDefs<double>::Pi()*x);
  93 }
  94
  95 /* The zero-order modified Bessel function of the first kind, used for the
  96  * Kaiser window.
  97  *
  98  *   I_0(x) = sum_{k=0}^inf (1 / k!)^2 (x / 2)^(2 k)
  99  *          = sum_{k=0}^inf ((x / 2)^k / k!)^2
 100  */
 101 constexpr double BesselI_0(const double x)
 102 {
 103     /* Start at k=1 since k=0 is trivial. */
 104     const double x2{x / 2.0};
 105     double term{1.0};
 106     double sum{1.0};
 107     double last_sum{};
 108     int k{1};
 109
 110     /* Let the integration converge until the term of the sum is no longer
 111      * significant.
 112      */
 113     do {
 114         const double y{x2 / k};
 115         ++k;
 116         last_sum = sum;
 117         term *= y * y;
 118         sum += term;
 119     } while(sum != last_sum);
 120
 121     return sum;
 122 }
 123
 124 /* Calculate a Kaiser window from the given beta value and a normalized k
 125  * [-1, 1].
 126  *
 127  *   w(k) = { I_0(B sqrt(1 - k^2)) / I_0(B),  -1 <= k <= 1
 128  *          { 0,                              elsewhere.
 129  *
 130  * Where k can be calculated as:
 131  *
 132  *   k = i / l,         where -l <= i <= l.
 133  *
 134  * or:
 135  *
 136  *   k = 2 i / M - 1,   where 0 <= i <= M.
 137  */
 138 constexpr double Kaiser(const double beta, const double k, const double besseli_0_beta)
 139 {
 140     if(!(k >= -1.0 && k <= 1.0))
 141         return 0.0;
 142     return BesselI_0(beta * sqrt(1.0 - k*k)) / besseli_0_beta;
 143 }
 144
 145 /* Calculates the (normalized frequency) transition width of the Kaiser window.
 146  * Rejection is in dB.
 147  */
 148 constexpr double CalcKaiserWidth(const double rejection, const int order)
 149 {
 150     if(rejection > 21.19)
 151        return (rejection - 7.95) / (order * 2.285 * al::MathDefs<double>::Tau());
 152     /* This enforces a minimum rejection of just above 21.18dB */
 153     return 5.79 / (order * al::MathDefs<double>::Tau());
 154 }
 155
 156 /* Calculates the beta value of the Kaiser window. Rejection is in dB. */
 157 constexpr double CalcKaiserBeta(const double rejection)
 158 {
 159     if(rejection > 50.0)
 160        return 0.1102 * (rejection-8.7);
 161     else if(rejection >= 21.0)
 162        return (0.5842 * std::pow(rejection-21.0, 0.4)) + (0.07886 * (rejection-21.0));
 163     return 0.0;
 164 }
 165
 166
 167 struct BSincHeader {
 168     double width;
 169     double beta;
 170     double scaleBase;
 171     double scaleRange;
 172     double besseli_0_beta;
 173
 174     int a[BSINC_SCALE_COUNT];
 175     int total_size;
 176 };
 177
 178 constexpr BSincHeader GenerateBSincHeader(int Rejection, int Order)
 179 {
 180     BSincHeader ret{};
 181     ret.width = CalcKaiserWidth(Rejection, Order);
 182     ret.beta = CalcKaiserBeta(Rejection);
 183     ret.scaleBase = ret.width / 2.0;
 184     ret.scaleRange = 1.0 - ret.scaleBase;
 185     ret.besseli_0_beta = BesselI_0(ret.beta);
 186
 187     int num_points{Order+1};
 188     for(int si{0};si < BSincScaleCount;++si)
 189     {
 190         const double scale{ret.scaleBase + (ret.scaleRange * si / (BSincScaleCount - 1))};
 191         const int a{std::min(static_cast<int>(num_points / 2.0 / scale), num_points)};
 192         const int m{2 * a};
 193
 194         ret.a[si] = a;
 195         ret.total_size += 4 * BSincPhaseCount * ((m+3) & ~3);
 196     }
 197
 198     return ret;
 199 }
 200
 201 /* 11th and 23rd order filters (12 and 24-point respectively) with a 60dB drop
 202  * at nyquist. Each filter will scale up the order when downsampling, to 23rd
 203  * and 47th order respectively.
 204  */
 205 constexpr BSincHeader bsinc12_hdr{GenerateBSincHeader(60, 11)};
 206 constexpr BSincHeader bsinc24_hdr{GenerateBSincHeader(60, 23)};
 207
 208
 209 /* FIXME: This should be constexpr, but the temporary filter arrays are too
 210  * big. This requires using heap space, which is not allowed in a constexpr
 211  * function (maybe in C++20).
 212  */
 213 template<size_t total_size>
 214 std::array<float,total_size> GenerateBSincCoeffs(const BSincHeader &hdr)
 215 {
 216     auto filter = std::make_unique<double[][BSincPhaseCount+1][BSincPointsMax]>(BSincScaleCount);
 217
 218     /* Calculate the Kaiser-windowed Sinc filter coefficients for each scale
 219      * and phase index.
 220      */
 221     for(unsigned int si{0};si < BSincScaleCount;++si)
 222     {
 223         const int m{hdr.a[si] * 2};
 224         const int o{BSincPointsHalf - (m/2)};
 225         const int l{hdr.a[si] - 1};
 226         const int a{hdr.a[si]};
 227         const double scale{hdr.scaleBase + (hdr.scaleRange * si / (BSincScaleCount - 1))};
 228         const double cutoff{scale - (hdr.scaleBase * std::max(0.5, scale) * 2.0)};
 229
 230         /* Do one extra phase index so that the phase delta has a proper target
 231          * for its last index.
 232          */
 233         for(int pi{0};pi <= BSincPhaseCount;++pi)
 234         {
 235             const double phase{l + (pi/double{BSincPhaseCount})};
 236
 237             for(int i{0};i < m;++i)
 238             {
 239                 const double x{i - phase};
 240                 filter[si][pi][o+i] = Kaiser(hdr.beta, x/a, hdr.besseli_0_beta) * cutoff *
 241                     Sinc(cutoff*x);
 242             }
 243         }
 244     }
 245
 246     auto ret = std::make_unique<std::array<float,total_size>>();
 247     size_t idx{0};
 248
 249     for(unsigned int si{0};si < BSincScaleCount-1;++si)
 250     {
 251         const int m{((hdr.a[si]*2) + 3) & ~3};
 252         const int o{BSincPointsHalf - (m/2)};
 253
 254         for(int pi{0};pi < BSincPhaseCount;++pi)
 255         {
 256             /* Write out the filter. Also calculate and write out the phase and
 257              * scale deltas.
 258              */
 259             for(int i{0};i < m;++i)
 260                 (*ret)[idx++] = static_cast<float>(filter[si][pi][o+i]);
 261
 262             /* Linear interpolation between phases is simplified by pre-
 263              * calculating the delta (b - a) in: x = a + f (b - a)
 264              */
 265             for(int i{0};i < m;++i)
 266             {
 267                 const double phDelta{filter[si][pi+1][o+i] - filter[si][pi][o+i]};
 268                 (*ret)[idx++] = static_cast<float>(phDelta);
 269             }
 270
 271             /* Linear interpolation between scales is also simplified.
 272              *
 273              * Given a difference in points between scales, the destination
 274              * points will be 0, thus: x = a + f (-a)
 275              */
 276             for(int i{0};i < m;++i)
 277             {
 278                 const double scDelta{filter[si+1][pi][o+i] - filter[si][pi][o+i]};
 279                 (*ret)[idx++] = static_cast<float>(scDelta);
 280             }
 281
 282             /* This last simplification is done to complete the bilinear
 283              * equation for the combination of phase and scale.
 284              */
 285             for(int i{0};i < m;++i)
 286             {
 287                 const double spDelta{(filter[si+1][pi+1][o+i] - filter[si+1][pi][o+i]) -
 288                     (filter[si][pi+1][o+i] - filter[si][pi][o+i])};
 289                 (*ret)[idx++] = static_cast<float>(spDelta);
 290             }
 291         }
 292     }
 293     {
 294         /* The last scale index doesn't have any scale or scale-phase deltas. */
 295         const unsigned int si{BSincScaleCount - 1};
 296         const int m{((hdr.a[si]*2) + 3) & ~3};
 297         const int o{BSincPointsHalf - (m/2)};
 298
 299         for(int pi{0};pi < BSincPhaseCount;++pi)
 300         {
 301             for(int i{0};i < m;++i)
 302                 (*ret)[idx++] = static_cast<float>(filter[si][pi][o+i]);
 303             for(int i{0};i < m;++i)
 304             {
 305                 const double phDelta{filter[si][pi+1][o+i] - filter[si][pi][o+i]};
 306                 (*ret)[idx++] = static_cast<float>(phDelta);
 307             }
 308             for(int i{0};i < m;++i)
 309                 (*ret)[idx++] = 0.0f;
 310             for(int i{0};i < m;++i)
 311                 (*ret)[idx++] = 0.0f;
 312         }
 313     }
 314     assert(idx == total_size);
 315
 316     return *ret;
 317 }
 318
 319 /* FIXME: These can't be constexpr due to the calls reaching the compiler's
 320  * step limit.
 321  */
 322 alignas(16) const auto bsinc12_table = GenerateBSincCoeffs<bsinc12_hdr.total_size>(bsinc12_hdr);
 323 alignas(16) const auto bsinc24_table = GenerateBSincCoeffs<bsinc24_hdr.total_size>(bsinc24_hdr);
 324
 325
 326 constexpr BSincTable GenerateBSincTable(const BSincHeader &hdr, const float *tab)
 327 {
 328     BSincTable ret{};
 329     ret.scaleBase = static_cast<float>(hdr.scaleBase);
 330     ret.scaleRange = static_cast<float>(1.0 / hdr.scaleRange);
 331     for(int i{0};i < BSincScaleCount;++i)
 332         ret.m[i] = static_cast<unsigned int>(((hdr.a[i]*2) + 3) & ~3);
 333     ret.filterOffset[0] = 0;
 334     for(int i{1};i < BSincScaleCount;++i)
 335         ret.filterOffset[i] = ret.filterOffset[i-1] + ret.m[i-1]*4*BSincPhaseCount;
 336     ret.Tab = tab;
 337     return ret;
 338 }
 339
 340 } // namespace
 341
 342 const BSincTable bsinc12{GenerateBSincTable(bsinc12_hdr, &bsinc12_table.front())};
 343 const BSincTable bsinc24{GenerateBSincTable(bsinc24_hdr, &bsinc24_table.front())};