74 lines
2.9 KiB
C++
74 lines
2.9 KiB
C++
#include "resampler.h"
|
|
using namespace Xybrid::NodeLib;
|
|
|
|
#include <QDebug>
|
|
|
|
#include <iostream>
|
|
#include <array>
|
|
#include <cmath>
|
|
|
|
#ifdef WITH_BOOST
|
|
#include <boost/math/special_functions/bessel.hpp>
|
|
using boost::math::cyl_bessel_i;
|
|
#else
|
|
using std::cyl_bessel_i;
|
|
#endif
|
|
|
|
namespace {
|
|
const constexpr double PI = 3.141592653589793238462643383279502884197169399375105820974;
|
|
|
|
// 5.5 for downrate; example gave 7.5
|
|
const constexpr double KAISER_ALPHA = 5.5;
|
|
const constexpr double KAISER_BETA = PI * KAISER_ALPHA;
|
|
|
|
inline constexpr double sinc(double x) {
|
|
if (x == 0) return 1;
|
|
double px = x * PI/1;
|
|
return std::sin(px) / px;
|
|
}
|
|
|
|
#if __cplusplus >= 202002L
|
|
using std::lerp;
|
|
#else
|
|
inline constexpr double lerp(double a, double b, double t) { return (1.0 - t) * a + t * b; }
|
|
#endif
|
|
}
|
|
|
|
|
|
|
|
// generate
|
|
const std::array<std::array<std::array<double, LUT_TAPS>, LUT_LEVELS>, LUT_STEPS> Xybrid::NodeLib::resamplerLUT = [] {
|
|
|
|
double denom = cyl_bessel_i(0, KAISER_BETA);
|
|
|
|
std::array<std::array<std::array<double, LUT_TAPS>, LUT_LEVELS>, LUT_STEPS> t;
|
|
//t[0] = {0, 0, 0, 1, 0, 0, 0, 0}; // we already know the ideal integer step
|
|
for (size_t step = 0; step < LUT_STEPS; step++) {
|
|
double sv = static_cast<double>(step) / LUT_STEPS; // subvalue (offset of tap position)
|
|
for (size_t tap = 0; tap < LUT_TAPS; tap++) {
|
|
double x = static_cast<double>(tap) - sv; // x position of tap;
|
|
double sx = x-LUT_HTAPS;
|
|
double kaiser = cyl_bessel_i(0, KAISER_BETA * std::sqrt(1.0 - std::pow( ( (2.0*(x+1))/(LUT_TAPS) ) - 1.0, 2 ) ) ) / denom; // original kaiser window generation
|
|
//double kaiser = cyl_bessel_i(0, KAISER_BETA * std::sqrt(1.0 - std::pow( (2.0*x)/(LUT_TAPS-2) - 1.0, 2 ) ) ) / denom; // by-the-book kaiser window of length LUT_TAPS-1
|
|
//double idl = (2.0*PI)/(LUT_TAPS-1);
|
|
//double kaiser = 0.40243 - 0.49804 * std::cos(idl * x) + 0.09831 * std::cos(2.0 * idl * x) - 0.00122 * std::cos(3.0 * idl * x); // approximate
|
|
//kaiser = std::max(kaiser, 0.0);
|
|
|
|
for (size_t i = 0; i < LUT_LEVELS; i++) {
|
|
double m = 1.0/std::max(static_cast<double>(i), 1.0); // sinc function "expands" the higher we pitch things
|
|
double om = lerp(0.2, 1.0, m) // we need to compensate slightly for amplitude loss at higher levels
|
|
/* */ * lerp(1.0, kaiser, std::pow(m, 0.333)); // apply kaiser proportionally; we want it less the higher we go
|
|
t[step][i][tap] = sinc(sx*m) * om;
|
|
if (t[step][i][tap] != t[step][i][tap]) t[step][i][tap] = 0; // NaN guard
|
|
}
|
|
|
|
if (t[step][0][tap] != t[step][0][tap]) t[step][0][tap] = 0; // NaN guard
|
|
}
|
|
}
|
|
/*t[0] = {};
|
|
t[0][LUT_HTAPS] = 1;*/
|
|
/*for (auto v : t[0]) std::cout << v << ", ";
|
|
std::cout << std::endl;*/
|
|
return t;
|
|
}();
|