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
  11 inline double sinc(double x)
  12 {
  13         if (fabs(x) < 1e-6) {
  14                 return 1.0f - fabs(x);
  15         } else {
  16                 return sin(x) / x;
  17         }
  18 }
  19
  20 inline double lanczos_weight(double x)
  21 {
  22         if (fabs(x) > LANCZOS_RADIUS) {
  23                 return 0.0f;
  24         }
  25         return sinc(M_PI * x) * sinc(M_PI * x / LANCZOS_RADIUS);
  26 }
  27
  28 template<class T>
  29 inline double lanczos_interpolate(const std::vector<T> &pcm, double i)
  30 {
  31         int lower = std::max<int>(ceil(i - LANCZOS_RADIUS), 0);
  32         int upper = std::min<int>(floor(i + LANCZOS_RADIUS), pcm.size() - 1);
  33         double sum = 0.0f;
  34
  35         for (int x = lower; x <= upper; ++x) {
  36                 sum += pcm[x] * lanczos_weight(i - x);
  37         }
  38         return sum;
  39 }
  40
  41 template<class T>
  42 inline double lanczos_interpolate_right(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].right * lanczos_weight(i - x);
  50         }
  51         return sum;
  52 }
  53
  54 template<class T>
  55 inline double linear_interpolate(const std::vector<T> &pcm, double i)
  56 {
  57         int ii = int(i);
  58         if (ii < 0 || ii >= int(pcm.size() - 1)) {
  59                 return 0.0;
  60         }
  61         double frac = i - ii;
  62
  63         return pcm[ii] + frac * (pcm[ii + 1] - pcm[ii]);
  64 }
  65
  66 template<class T>
  67 inline double linear_interpolate_right(const std::vector<T> &pcm, double i)
  68 {
  69         int ii = int(i);
  70         if (ii < 0 || ii >= int(pcm.size() - 1)) {
  71                 return 0.0;
  72         }
  73         double frac = i - ii;
  74
  75         return pcm[ii].right + frac * (pcm[ii + 1].right - pcm[ii].right);
  76 }
  77
  78 #endif  // !defined(_INTERPOLATE_H)