Move the tap_hdr definition into a separate header file.
[c64tapwav] / synth.cpp
1 #include <stdio.h>
2 #include <math.h>
3 #include <stdlib.h>
4 #include <vector>
5 #include <assert.h>
6
7 #include "common.h"
8
9 using namespace std;
10
11 #define LANCZOS_RADIUS 10
12 #define RESOLUTION 256
13
14 // Cutoff frequency of low-pass filter (must be max. WAVE_FREQ / 2 or you
15 // will get aliasing)
16 #define LPFILTER_FREQ 22050
17
18 // Cutoff frequency of high-pass filter
19 #define HPFILTER_FREQ 40
20
21 #define NORMALIZED_LPFILTER_FREQ (float(LPFILTER_FREQ) / float(WAVE_FREQ))
22 #define LANCZOS_EFFECTIVE_RADIUS (LANCZOS_RADIUS / (NORMALIZED_LPFILTER_FREQ * 2.0))
23
24 #define SIZEOF_HALF_TABLE ((int)((LANCZOS_EFFECTIVE_RADIUS * RESOLUTION)))
25 static double *integrated_window;
26
27 double sinc(double x)
28 {
29         if (fabs(x) < 1e-6) {
30                 return 1.0f - fabs(x);
31         } else {
32                 return sin(x) / x;
33         }
34 }
35
36 double weight(double x)
37 {
38         if (fabs(x) > LANCZOS_EFFECTIVE_RADIUS) {
39                 return 0.0f;
40         }
41         float t = 2.0 * M_PI * NORMALIZED_LPFILTER_FREQ * x;
42         return sinc(t) * sinc(t / LANCZOS_RADIUS);
43 }
44
45 void make_table()
46 {
47         integrated_window = new double[2 * SIZEOF_HALF_TABLE + 1];
48
49         double sum = 0.0;
50         for (int i = 0; i <= 2 * SIZEOF_HALF_TABLE; ++i) {
51                 float t = (i - SIZEOF_HALF_TABLE) / (float)RESOLUTION;
52                 integrated_window[i] = sum;
53                 sum += weight(t) * NORMALIZED_LPFILTER_FREQ * 2.0 / (float)RESOLUTION;
54                 //printf("%f %f %f\n", t, weight(t), sum);
55         }
56 //      exit(0);
57 }
58
59 // integral from -inf to window
60 double window_one_sided_integral(double to)
61 {
62         double array_pos = to * RESOLUTION + SIZEOF_HALF_TABLE;
63
64         int whole = int(floor(array_pos));
65         double frac = array_pos - whole;
66         if (whole < 0) {
67                 return 0.0;
68         }
69         if (whole >= 2 * SIZEOF_HALF_TABLE) {
70                 return 1.0;
71         }
72         return integrated_window[whole] + frac * (integrated_window[whole + 1] - integrated_window[whole]);
73 }
74
75 double window_integral(double from, double to)
76 {
77         return window_one_sided_integral(to) - window_one_sided_integral(from);
78 }
79
80 struct pulse {
81         // all values in samples
82         double start, end;
83 };
84 vector<pulse> pulses;
85
86 static float a1, a2, b0, b1, b2;
87 static float d0, d1;
88
89 static void filter_init(float cutoff_radians)
90 {
91         float resonance = 1.0f / sqrt(2.0f);
92         float sn, cs;
93         sincosf(cutoff_radians, &sn, &cs);
94
95         float alpha = float(sn / (2 * resonance));
96
97         // coefficients for highpass filter
98         float a0 = 1 + alpha;
99         b0 = (1 + cs) * 0.5f;
100         b1 = -(1 + cs);
101         b2 = b0;
102         a1 = -2 * cs;
103         a2 = 1 - alpha;
104
105         float invA0 = 1.0f / a0;
106         b0 *= invA0;
107         b1 *= invA0;
108         b2 *= invA0;
109         a1 *= invA0;
110         a2 *= invA0;
111
112         // reset filter delays
113         d0 = d1 = 0.0f;
114 }
115
116 static float filter_update(float in)
117 {
118         float out = b0*in + d0;
119         d0 = b1 * in - a1 * out + d1;
120         d1 = b2 * in - a2 * out;
121         return out;
122 }
123
124 vector<float> synth(const vector<pulse> &pulses)
125 {
126         make_table();
127
128         int len_samples = int(ceil(pulses.back().start + (pulses.back().end - pulses.back().start) * 2));
129         
130         fprintf(stderr, "%d pulses, total %.2f seconds (%d samples)\n", int(pulses.size()), len_samples / float(WAVE_FREQ), len_samples);
131
132         int pulse_begin = 0;
133
134         vector<float> samples;
135         samples.reserve(len_samples);
136
137         for (int i = 0; i < len_samples; ++i) {
138 //      for (int i = 50000000; i < 51000000; ++i) {
139 //              double t = i * 0.01;
140                 double t = i;
141                 double sample = -0.5;
142                 for (unsigned j = pulse_begin; j < pulses.size(); ++j) {
143                         if (t - pulses[j].end > LANCZOS_EFFECTIVE_RADIUS) {
144                                 ++pulse_begin;
145                                 continue;
146                         }
147                         if (pulses[j].start - t > LANCZOS_EFFECTIVE_RADIUS) {
148                                 break;
149                         }
150                         float contribution = window_integral(pulses[j].start - t, pulses[j].end - t);
151                         sample += contribution;
152                 }
153                 samples.push_back(sample);
154         }
155
156         // filter forwards, then backwards (perfect phase filtering)
157         vector<float> filtered_samples, refiltered_samples;
158         filtered_samples.resize(len_samples);
159         refiltered_samples.resize(len_samples);
160
161         filter_init(M_PI * HPFILTER_FREQ / WAVE_FREQ);
162         for (int i = 0; i < len_samples; ++i) {
163                 filtered_samples[i] = filter_update(samples[i]);
164         }
165         filter_init(M_PI * HPFILTER_FREQ / WAVE_FREQ);
166         for (int i = len_samples; i --> 0; ) {
167                 refiltered_samples[i] = filter_update(filtered_samples[i]);
168         }
169
170         return refiltered_samples;
171 }