Merge branch 'master' of /srv/git.sesse.net/www/c64tapwav
[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 #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)