git.sesse.net Git - c64tapwav/blob - interpolate.h

   1 #ifndef _INTERPOLATE_H
   2 #define _INTERPOLATE_H 1
   3
   4 #include <math.h>
   5
   6 #include <algorithm>
   7 #include <vector>
   8
   9 #define LANCZOS_RADIUS 30
  10 #define LANCZOS_RESOLUTION 256
  11
  12 inline double sinc(double x)
  13 {
  14         if (fabs(x) < 1e-6) {
  15                 return 1.0f - fabs(x);
  16         } else {
  17                 return sin(x) / x;
  18         }
  19 }
  20
  21 inline double lanczos_weight(double x)
  22 {
  23         if (fabs(x) > LANCZOS_RADIUS) {
  24                 return 0.0f;
  25         }
  26         return sinc(M_PI * x) * sinc(M_PI * x / LANCZOS_RADIUS);
  27 }
  28
  29 extern double lanczos_table[LANCZOS_RADIUS * LANCZOS_RESOLUTION];
  30 void make_lanczos_weight_table();
  31
  32 inline double lanczos_weight_table(float x)
  33 {
  34         int table_id = lrintf(fabsf(x) * LANCZOS_RESOLUTION);
  35         if (table_id >= LANCZOS_RADIUS * LANCZOS_RESOLUTION) {
  36                 return 0.0;
  37         }
  38         return lanczos_table[table_id];
  39 }
  40
  41 template<class T>
  42 inline double lanczos_interpolate(const std::vector<T> &pcm, double i)
  43 {
  44         int lower = std::max<int>(ceil(i - LANCZOS_RADIUS), 0);
  45         int upper = std::min<int>(floor(i + LANCZOS_RADIUS), pcm.size() - 1);
  46         double sum = 0.0f;
  47
  48         for (int x = lower; x <= upper; ++x) {
  49                 sum += pcm[x] * lanczos_weight_table(i - x);
  50         }
  51         return sum;
  52 }
  53
  54 template<class T>
  55 inline double lanczos_interpolate_right(const std::vector<T> &pcm, double i)
  56 {
  57         int lower = std::max<int>(ceil(i - LANCZOS_RADIUS), 0);
  58         int upper = std::min<int>(floor(i + LANCZOS_RADIUS), pcm.size() - 1);
  59         double sum = 0.0f;
  60
  61         for (int x = lower; x <= upper; ++x) {
  62                 sum += pcm[x].right * lanczos_weight_table(i - x);
  63         }
  64         return sum;
  65 }
  66
  67 template<class T>
  68 inline double linear_interpolate(const std::vector<T> &pcm, double i)
  69 {
  70         int ii = int(i);
  71         if (ii < 0 || ii >= int(pcm.size() - 1)) {
  72                 return 0.0;
  73         }
  74         double frac = i - ii;
  75
  76         return pcm[ii] + frac * (pcm[ii + 1] - pcm[ii]);
  77 }
  78
  79 template<class T>
  80 inline double linear_interpolate_right(const std::vector<T> &pcm, double i)
  81 {
  82         int ii = int(i);
  83         if (ii < 0 || ii >= int(pcm.size() - 1)) {
  84                 return 0.0;
  85         }
  86         double frac = i - ii;
  87
  88         return pcm[ii].right + frac * (pcm[ii + 1].right - pcm[ii].right);
  89 }
  90
  91 #endif  // !defined(_INTERPOLATE_H)