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 using uint = unsigned int;
  19
  20
  21 /* This is the normalized cardinal sine (sinc) function.
  22  *
  23  *   sinc(x) = { 1,                   x = 0
  24  *             { sin(pi x) / (pi x),  otherwise.
  25  */
  26 constexpr double Sinc(const double x)
  27 {
  28     if(!(x > 1e-15 || x < -1e-15))
  29         return 1.0;
  30     return std::sin(al::MathDefs<double>::Pi()*x) / (al::MathDefs<double>::Pi()*x);
  31 }
  32
  33 /* The zero-order modified Bessel function of the first kind, used for the
  34  * Kaiser window.
  35  *
  36  *   I_0(x) = sum_{k=0}^inf (1 / k!)^2 (x / 2)^(2 k)
  37  *          = sum_{k=0}^inf ((x / 2)^k / k!)^2
  38  */
  39 constexpr double BesselI_0(const double x)
  40 {
  41     /* Start at k=1 since k=0 is trivial. */
  42     const double x2{x / 2.0};
  43     double term{1.0};
  44     double sum{1.0};
  45     double last_sum{};
  46     int k{1};
  47
  48     /* Let the integration converge until the term of the sum is no longer
  49      * significant.
  50      */
  51     do {
  52         const double y{x2 / k};
  53         ++k;
  54         last_sum = sum;
  55         term *= y * y;
  56         sum += term;
  57     } while(sum != last_sum);
  58
  59     return sum;
  60 }
  61
  62 /* Calculate a Kaiser window from the given beta value and a normalized k
  63  * [-1, 1].
  64  *
  65  *   w(k) = { I_0(B sqrt(1 - k^2)) / I_0(B),  -1 <= k <= 1
  66  *          { 0,                              elsewhere.
  67  *
  68  * Where k can be calculated as:
  69  *
  70  *   k = i / l,         where -l <= i <= l.
  71  *
  72  * or:
  73  *
  74  *   k = 2 i / M - 1,   where 0 <= i <= M.
  75  */
  76 constexpr double Kaiser(const double beta, const double k, const double besseli_0_beta)
  77 {
  78     if(!(k >= -1.0 && k <= 1.0))
  79         return 0.0;
  80     return BesselI_0(beta * std::sqrt(1.0 - k*k)) / besseli_0_beta;
  81 }
  82
  83 /* Calculates the (normalized frequency) transition width of the Kaiser window.
  84  * Rejection is in dB.
  85  */
  86 constexpr double CalcKaiserWidth(const double rejection, const uint order)
  87 {
  88     if(rejection > 21.19)
  89         return (rejection - 7.95) / (order * 2.285 * al::MathDefs<double>::Tau());
  90     /* This enforces a minimum rejection of just above 21.18dB */
  91     return 5.79 / (order * al::MathDefs<double>::Tau());
  92 }
  93
  94 /* Calculates the beta value of the Kaiser window. Rejection is in dB. */
  95 constexpr double CalcKaiserBeta(const double rejection)
  96 {
  97     if(rejection > 50.0)
  98         return 0.1102 * (rejection-8.7);
  99     else if(rejection >= 21.0)
 100         return (0.5842 * std::pow(rejection-21.0, 0.4)) + (0.07886 * (rejection-21.0));
 101     return 0.0;
 102 }
 103
 104 /* NOTE: GCC 5 has an issue with BSincHeader objects being in an anonymous
 105  * namespace while also being used as non-type template parameters.
 106  */
 107 #if !defined(__clang__) && defined(__GNUC__) && __GNUC__ < 6
 108 } // namespace
 109 #endif
 110
 111 struct BSincHeader {
 112     double width{};
 113     double beta{};
 114     double scaleBase{};
 115     double scaleRange{};
 116     double besseli_0_beta{};
 117
 118     uint a[BSincScaleCount]{};
 119     uint total_size{};
 120
 121     constexpr BSincHeader(uint Rejection, uint Order) noexcept
 122     {
 123         width = CalcKaiserWidth(Rejection, Order);
 124         beta = CalcKaiserBeta(Rejection);
 125         scaleBase = width / 2.0;
 126         scaleRange = 1.0 - scaleBase;
 127         besseli_0_beta = BesselI_0(beta);
 128
 129         uint num_points{Order+1};
 130         for(uint si{0};si < BSincScaleCount;++si)
 131         {
 132             const double scale{scaleBase + (scaleRange * si / (BSincScaleCount-1))};
 133             const uint a_{std::min(static_cast<uint>(num_points / 2.0 / scale), num_points)};
 134             const uint m{2 * a_};
 135
 136             a[si] = a_;
 137             total_size += 4 * BSincPhaseCount * ((m+3) & ~3u);
 138         }
 139     }
 140 };
 141
 142 /* 11th and 23rd order filters (12 and 24-point respectively) with a 60dB drop
 143  * at nyquist. Each filter will scale up the order when downsampling, to 23rd
 144  * and 47th order respectively.
 145  */
 146 constexpr BSincHeader bsinc12_hdr{60, 11};
 147 constexpr BSincHeader bsinc24_hdr{60, 23};
 148
 149 #if !defined(__clang__) && defined(__GNUC__) && __GNUC__ < 6
 150 namespace {
 151 #endif
 152
 153 /* FIXME: This should be constexpr, but the temporary filter arrays are too
 154  * big. This requires using heap space, which is not allowed in a constexpr
 155  * function (maybe in C++20).
 156  */
 157 template<const BSincHeader &hdr>
 158 struct BSincFilterArray {
 159     alignas(16) std::array<float, hdr.total_size> mTable;
 160
 161     BSincFilterArray()
 162     {
 163         using filter_type = double[][BSincPhaseCount+1][BSincPointsMax];
 164         auto filter = std::make_unique<filter_type>(BSincScaleCount);
 165
 166         /* Calculate the Kaiser-windowed Sinc filter coefficients for each
 167          * scale and phase index.
 168          */
 169         for(uint si{0};si < BSincScaleCount;++si)
 170         {
 171             const uint m{hdr.a[si] * 2};
 172             const size_t o{(BSincPointsMax-m) / 2};
 173             const double scale{hdr.scaleBase + (hdr.scaleRange * si / (BSincScaleCount-1))};
 174             const double cutoff{scale - (hdr.scaleBase * std::max(0.5, scale) * 2.0)};
 175             const auto a = static_cast<double>(hdr.a[si]);
 176             const double l{a - 1.0};
 177
 178             /* Do one extra phase index so that the phase delta has a proper
 179              * target for its last index.
 180              */
 181             for(uint pi{0};pi <= BSincPhaseCount;++pi)
 182             {
 183                 const double phase{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/a, hdr.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-1;++si)
 196         {
 197             const size_t m{((hdr.a[si]*2) + 3) & ~3u};
 198             const size_t o{(BSincPointsMax-m) / 2};
 199
 200             for(size_t pi{0};pi < BSincPhaseCount;++pi)
 201             {
 202                 /* Write out the filter. Also calculate and write out the phase
 203                  * and scale deltas.
 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                 for(size_t i{0};i < m;++i)
 212                 {
 213                     const double phDelta{filter[si][pi+1][o+i] - filter[si][pi][o+i]};
 214                     mTable[idx++] = static_cast<float>(phDelta);
 215                 }
 216
 217                 /* Linear interpolation between scales is also simplified.
 218                  *
 219                  * Given a difference in points between scales, the destination
 220                  * points will be 0, thus: x = a + f (-a)
 221                  */
 222                 for(size_t i{0};i < m;++i)
 223                 {
 224                     const double scDelta{filter[si+1][pi][o+i] - filter[si][pi][o+i]};
 225                     mTable[idx++] = static_cast<float>(scDelta);
 226                 }
 227
 228                 /* This last simplification is done to complete the bilinear
 229                  * equation for the combination of phase and scale.
 230                  */
 231                 for(size_t i{0};i < m;++i)
 232                 {
 233                     const double spDelta{(filter[si+1][pi+1][o+i] - filter[si+1][pi][o+i]) -
 234                         (filter[si][pi+1][o+i] - filter[si][pi][o+i])};
 235                     mTable[idx++] = static_cast<float>(spDelta);
 236                 }
 237             }
 238         }
 239         {
 240             /* The last scale index doesn't have any scale or scale-phase
 241              * deltas.
 242              */
 243             constexpr size_t si{BSincScaleCount-1};
 244             const size_t m{((hdr.a[si]*2) + 3) & ~3u};
 245             const size_t o{(BSincPointsMax-m) / 2};
 246
 247             for(size_t pi{0};pi < BSincPhaseCount;++pi)
 248             {
 249                 for(size_t i{0};i < m;++i)
 250                     mTable[idx++] = static_cast<float>(filter[si][pi][o+i]);
 251                 for(size_t i{0};i < m;++i)
 252                 {
 253                     const double phDelta{filter[si][pi+1][o+i] - filter[si][pi][o+i]};
 254                     mTable[idx++] = static_cast<float>(phDelta);
 255                 }
 256                 for(size_t i{0};i < m;++i)
 257                     mTable[idx++] = 0.0f;
 258                 for(size_t i{0};i < m;++i)
 259                     mTable[idx++] = 0.0f;
 260             }
 261         }
 262         assert(idx == hdr.total_size);
 263     }
 264 };
 265
 266 /* FIXME: These can't be constexpr due to the calls reaching the compiler's
 267  * step limit.
 268  */
 269 const BSincFilterArray<bsinc12_hdr> bsinc12_filter{};
 270 const BSincFilterArray<bsinc24_hdr> bsinc24_filter{};
 271
 272
 273 constexpr BSincTable GenerateBSincTable(const BSincHeader &hdr, const float *tab)
 274 {
 275     BSincTable ret{};
 276     ret.scaleBase = static_cast<float>(hdr.scaleBase);
 277     ret.scaleRange = static_cast<float>(1.0 / hdr.scaleRange);
 278     for(size_t i{0};i < BSincScaleCount;++i)
 279         ret.m[i] = ((hdr.a[i]*2) + 3) & ~3u;
 280     ret.filterOffset[0] = 0;
 281     for(size_t i{1};i < BSincScaleCount;++i)
 282         ret.filterOffset[i] = ret.filterOffset[i-1] + ret.m[i-1]*4*BSincPhaseCount;
 283     ret.Tab = tab;
 284     return ret;
 285 }
 286
 287 } // namespace
 288
 289 const BSincTable bsinc12{GenerateBSincTable(bsinc12_hdr, &bsinc12_filter.mTable.front())};
 290 const BSincTable bsinc24{GenerateBSincTable(bsinc24_hdr, &bsinc24_filter.mTable.front())};