core/cubic_tables.cpp

   1
   2 #include "cubic_tables.h"
   3
   4 #include <array>
   5 #include <cmath>
   6 #include <cstddef>
   7
   8 #include "alnumbers.h"
   9 #include "alnumeric.h"
  10 #include "cubic_defs.h"
  11
  12 /* These gaussian filter tables are inspired by the gaussian-like filter found
  13  * in the SNES. This is based on the public domain code developed by Near, with
  14  * the help of Ryphecha and nocash, from the nesdev.org forums.
  15  *
  16  * <https://forums.nesdev.org/viewtopic.php?p=251534#p251534>
  17  *
  18  * Additional changes were made here, the most obvious being that is has full
  19  * floating-point precision instead of 11-bit fixed-point, but also an offset
  20  * adjustment for the coefficients to better preserve phase.
  21  */
  22 namespace {
  23
  24 [[nodiscard]]
  25 auto GetCoeff(double idx) noexcept -> double
  26 {
  27     const double k{0.5 + idx};
  28     if(k > 512.0) return 0.0;
  29     const double s{ std::sin(al::numbers::pi*1.280/1024 * k)};
  30     const double t{(std::cos(al::numbers::pi*2.000/1023 * k) - 1.0) * 0.50};
  31     const double u{(std::cos(al::numbers::pi*4.000/1023 * k) - 1.0) * 0.08};
  32     return s * (t + u + 1.0) / k;
  33 }
  34
  35 } // namespace
  36
  37 GaussianTable::GaussianTable()
  38 {
  39     static constexpr double IndexScale{512.0 / double{CubicPhaseCount*2}};
  40     /* Fill in the main coefficients. */
  41     for(std::size_t pi{0};pi < CubicPhaseCount;++pi)
  42     {
  43         const double coeff0{GetCoeff(static_cast<double>(CubicPhaseCount + pi)*IndexScale)};
  44         const double coeff1{GetCoeff(static_cast<double>(pi)*IndexScale)};
  45         const double coeff2{GetCoeff(static_cast<double>(CubicPhaseCount - pi)*IndexScale)};
  46         const double coeff3{GetCoeff(static_cast<double>(CubicPhaseCount*2_uz-pi)*IndexScale)};
  47
  48         const double scale{1.0 / (coeff0 + coeff1 + coeff2 + coeff3)};
  49         mTable[pi].mCoeffs[0] = static_cast<float>(coeff0 * scale);
  50         mTable[pi].mCoeffs[1] = static_cast<float>(coeff1 * scale);
  51         mTable[pi].mCoeffs[2] = static_cast<float>(coeff2 * scale);
  52         mTable[pi].mCoeffs[3] = static_cast<float>(coeff3 * scale);
  53     }
  54
  55     /* Fill in the coefficient deltas. */
  56     for(std::size_t pi{0};pi < CubicPhaseCount-1;++pi)
  57     {
  58         mTable[pi].mDeltas[0] = mTable[pi+1].mCoeffs[0] - mTable[pi].mCoeffs[0];
  59         mTable[pi].mDeltas[1] = mTable[pi+1].mCoeffs[1] - mTable[pi].mCoeffs[1];
  60         mTable[pi].mDeltas[2] = mTable[pi+1].mCoeffs[2] - mTable[pi].mCoeffs[2];
  61         mTable[pi].mDeltas[3] = mTable[pi+1].mCoeffs[3] - mTable[pi].mCoeffs[3];
  62     }
  63
  64     const std::size_t pi{CubicPhaseCount - 1};
  65     mTable[pi].mDeltas[0] =                 0.0f - mTable[pi].mCoeffs[0];
  66     mTable[pi].mDeltas[1] = mTable[0].mCoeffs[0] - mTable[pi].mCoeffs[1];
  67     mTable[pi].mDeltas[2] = mTable[0].mCoeffs[1] - mTable[pi].mCoeffs[2];
  68     mTable[pi].mDeltas[3] = mTable[0].mCoeffs[2] - mTable[pi].mCoeffs[3];
  69 }
  70
  71 SplineTable::SplineTable()
  72 {
  73     /* This filter table is based on a Catmull-Rom spline. It retains more of
  74      * the original high-frequency content, at the cost of increased harmonics.
  75      */
  76     for(std::size_t pi{0};pi < CubicPhaseCount;++pi)
  77     {
  78         const double mu{static_cast<double>(pi) / double{CubicPhaseCount}};
  79         const double mu2{mu*mu}, mu3{mu2*mu};
  80         mTable[pi].mCoeffs[0] = static_cast<float>(-0.5*mu3 +      mu2 + -0.5*mu);
  81         mTable[pi].mCoeffs[1] = static_cast<float>( 1.5*mu3 + -2.5*mu2           + 1.0);
  82         mTable[pi].mCoeffs[2] = static_cast<float>(-1.5*mu3 +  2.0*mu2 +  0.5*mu);
  83         mTable[pi].mCoeffs[3] = static_cast<float>( 0.5*mu3 + -0.5*mu2);
  84     }
  85
  86     for(std::size_t pi{0};pi < CubicPhaseCount-1;++pi)
  87     {
  88         mTable[pi].mDeltas[0] = mTable[pi+1].mCoeffs[0] - mTable[pi].mCoeffs[0];
  89         mTable[pi].mDeltas[1] = mTable[pi+1].mCoeffs[1] - mTable[pi].mCoeffs[1];
  90         mTable[pi].mDeltas[2] = mTable[pi+1].mCoeffs[2] - mTable[pi].mCoeffs[2];
  91         mTable[pi].mDeltas[3] = mTable[pi+1].mCoeffs[3] - mTable[pi].mCoeffs[3];
  92     }
  93
  94     const std::size_t pi{CubicPhaseCount - 1};
  95     mTable[pi].mDeltas[0] =                 0.0f - mTable[pi].mCoeffs[0];
  96     mTable[pi].mDeltas[1] = mTable[0].mCoeffs[0] - mTable[pi].mCoeffs[1];
  97     mTable[pi].mDeltas[2] = mTable[0].mCoeffs[1] - mTable[pi].mCoeffs[2];
  98     mTable[pi].mDeltas[3] = mTable[0].mCoeffs[2] - mTable[pi].mCoeffs[3];
  99 }
 100
 101
 102 CubicFilter::CubicFilter()
 103 {
 104     static constexpr double IndexScale{512.0 / double{sTableSteps*2}};
 105     /* Only half the coefficients need to be iterated here, since Coeff2 and
 106      * Coeff3 are just Coeff1 and Coeff0 in reverse respectively.
 107      */
 108     for(size_t i{0};i < sTableSteps/2 + 1;++i)
 109     {
 110         const double coeff0{GetCoeff(static_cast<double>(sTableSteps + i)*IndexScale)};
 111         const double coeff1{GetCoeff(static_cast<double>(i)*IndexScale)};
 112         const double coeff2{GetCoeff(static_cast<double>(sTableSteps - i)*IndexScale)};
 113         const double coeff3{GetCoeff(static_cast<double>(sTableSteps*2_uz - i)*IndexScale)};
 114
 115         const double scale{1.0 / (coeff0 + coeff1 + coeff2 + coeff3)};
 116         mFilter[sTableSteps + i] = static_cast<float>(coeff0 * scale);
 117         mFilter[i] = static_cast<float>(coeff1 * scale);
 118         mFilter[sTableSteps - i] = static_cast<float>(coeff2 * scale);
 119         mFilter[sTableSteps*2 - i] = static_cast<float>(coeff3 * scale);
 120     }
 121 }