Clean up interpolation a bit.
[c64tapwav] / 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)