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