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 // Cutoff frequency of high-pass filter
19 #define HPFILTER_FREQ 40
21 #define NORMALIZED_LPFILTER_FREQ (float(LPFILTER_FREQ) / float(WAVE_FREQ))
22 #define LANCZOS_EFFECTIVE_RADIUS (LANCZOS_RADIUS / (NORMALIZED_LPFILTER_FREQ * 2.0))
24 #define SIZEOF_HALF_TABLE ((int)((LANCZOS_EFFECTIVE_RADIUS * RESOLUTION)))
25 static double *integrated_window;
30 return 1.0f - fabs(x);
36 double weight(double x)
38 if (fabs(x) > LANCZOS_EFFECTIVE_RADIUS) {
41 float t = 2.0 * M_PI * NORMALIZED_LPFILTER_FREQ * x;
42 return sinc(t) * sinc(t / LANCZOS_RADIUS);
47 integrated_window = new double[2 * SIZEOF_HALF_TABLE + 1];
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);
59 // integral from -inf to window
60 double window_one_sided_integral(double to)
62 double array_pos = to * RESOLUTION + SIZEOF_HALF_TABLE;
64 int whole = int(floor(array_pos));
65 double frac = array_pos - whole;
69 if (whole >= 2 * SIZEOF_HALF_TABLE) {
72 return integrated_window[whole] + frac * (integrated_window[whole + 1] - integrated_window[whole]);
75 double window_integral(double from, double to)
77 return window_one_sided_integral(to) - window_one_sided_integral(from);
81 // all values in samples
86 static float a1, a2, b0, b1, b2;
89 static void filter_init(float cutoff_radians)
91 float resonance = 1.0f / sqrt(2.0f);
93 sincosf(cutoff_radians, &sn, &cs);
95 float alpha = float(sn / (2 * resonance));
97 // coefficients for highpass filter
105 float invA0 = 1.0f / a0;
112 // reset filter delays
116 static float filter_update(float in)
118 float out = b0*in + d0;
119 d0 = b1 * in - a1 * out + d1;
120 d1 = b2 * in - a2 * out;
124 int main(int argc, char **argv)
128 FILE *fp = fopen(argv[1], "rb");
129 fseek(fp, 14, SEEK_SET);
140 cycles = a | (b << 8) | (c << 16);
145 p.start = float(x) * WAVE_FREQ / C64_FREQ;
146 p.end = (float(x) + cycles * 0.5) * WAVE_FREQ / C64_FREQ;
152 int len_samples = int(ceil(float(len_cycles) * WAVE_FREQ / C64_FREQ));
154 fprintf(stderr, "%d pulses, total %.2f seconds (%d samples)\n", pulses.size(), len_cycles / float(C64_FREQ), len_samples);
158 vector<float> samples;
159 samples.reserve(len_samples);
161 for (int i = 0; i < len_samples; ++i) {
162 // for (int i = 50000000; i < 51000000; ++i) {
163 // double t = i * 0.01;
165 double sample = -0.5;
166 for (unsigned j = pulse_begin; j < pulses.size(); ++j) {
167 if (t - pulses[j].end > LANCZOS_EFFECTIVE_RADIUS) {
171 if (pulses[j].start - t > LANCZOS_EFFECTIVE_RADIUS) {
174 float contribution = window_integral(pulses[j].start - t, pulses[j].end - t);
175 sample += contribution;
177 samples.push_back(sample);
180 // filter forwards, then backwards (perfect phase filtering)
181 vector<float> filtered_samples, refiltered_samples;
182 filtered_samples.resize(len_samples);
183 refiltered_samples.resize(len_samples);
185 filter_init(M_PI * HPFILTER_FREQ / WAVE_FREQ);
186 for (int i = 0; i < len_samples; ++i) {
187 filtered_samples[i] = filter_update(samples[i]);
189 filter_init(M_PI * HPFILTER_FREQ / WAVE_FREQ);
190 for (int i = len_samples; i --> 0; ) {
191 refiltered_samples[i] = filter_update(filtered_samples[i]);
193 for (int i = 0; i < len_samples; ++i) {
194 //printf("%f %f\n", samples[i], refiltered_samples[i]);
195 short s = lrintf(refiltered_samples[i] * 16384.0f);
196 fwrite(&s, 2, 1, stdout);