9 #define LANCZOS_RADIUS 10
10 #define RESOLUTION 256
11 #define WAVE_FREQ 44100
12 #define C64_FREQ 985248
14 // Cutoff frequency of low-pass filter (must be max. WAVE_FREQ / 2 or you
16 #define LPFILTER_FREQ 22050
18 #define NORMALIZED_LPFILTER_FREQ (float(LPFILTER_FREQ) / float(WAVE_FREQ))
19 #define LANCZOS_EFFECTIVE_RADIUS (LANCZOS_RADIUS / (NORMALIZED_LPFILTER_FREQ * 2.0))
21 #define SIZEOF_HALF_TABLE ((int)((LANCZOS_EFFECTIVE_RADIUS * RESOLUTION)))
22 static double *integrated_window;
27 return 1.0f - fabs(x);
33 double weight(double x)
35 if (fabs(x) > LANCZOS_EFFECTIVE_RADIUS) {
38 float t = 2.0 * M_PI * NORMALIZED_LPFILTER_FREQ * x;
39 return sinc(t) * sinc(t / LANCZOS_RADIUS);
44 integrated_window = new double[2 * SIZEOF_HALF_TABLE + 1];
47 for (int i = 0; i <= 2 * SIZEOF_HALF_TABLE; ++i) {
48 float t = (i - SIZEOF_HALF_TABLE) / (float)RESOLUTION;
49 integrated_window[i] = sum;
50 sum += weight(t) * NORMALIZED_LPFILTER_FREQ * 2.0 / (float)RESOLUTION;
51 //printf("%f %f %f\n", t, weight(t), sum);
56 // integral from -inf to window
57 double window_one_sided_integral(double to)
59 double array_pos = to * RESOLUTION + SIZEOF_HALF_TABLE;
61 int whole = int(floor(array_pos));
62 double frac = array_pos - whole;
66 if (whole >= 2 * SIZEOF_HALF_TABLE) {
69 return integrated_window[whole] + frac * (integrated_window[whole + 1] - integrated_window[whole]);
72 double window_integral(double from, double to)
74 return window_one_sided_integral(to) - window_one_sided_integral(from);
78 // all values in samples
83 int main(int argc, char **argv)
87 FILE *fp = fopen(argv[1], "rb");
88 fseek(fp, 14, SEEK_SET);
99 cycles = a | (b << 8) | (c << 16);
104 p.start = float(x) * WAVE_FREQ / C64_FREQ;
105 p.end = (float(x) + cycles * 0.5) * WAVE_FREQ / C64_FREQ;
111 int len_samples = int(ceil(float(len_cycles) * WAVE_FREQ / C64_FREQ));
113 fprintf(stderr, "%d pulses, total %.2f seconds (%d samples)\n", pulses.size(), len_cycles / float(C64_FREQ), len_samples);
117 for (int i = 0; i < len_samples; ++i) {
118 // for (int i = 50000000; i < 51000000; ++i) {
119 // double t = i * 0.01;
122 for (unsigned j = pulse_begin; j < pulses.size(); ++j) {
123 if (t - pulses[j].end > LANCZOS_EFFECTIVE_RADIUS) {
127 if (pulses[j].start - t > LANCZOS_EFFECTIVE_RADIUS) {
130 float contribution = window_integral(pulses[j].start - t, pulses[j].end - t);
131 sample += contribution;
133 printf("%f\n", sample);
134 // short s = lrintf((sample - 0.5f) * 16384.0f);
135 // fwrite(&s, 2, 1, stdout);