git.sesse.net Git - c64tapwav/blob - synth.cpp

   1 #include <stdio.h>
   2 #include <math.h>
   3 #include <stdlib.h>
   4 #include <vector>
   5 #include <assert.h>
   6
   7 using namespace std;
   8
   9 #define LANCZOS_RADIUS 10
  10 #define RESOLUTION 256
  11 #define WAVE_FREQ 44100
  12 #define C64_FREQ 985248
  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 int main(int argc, char **argv)
 125 {
 126         make_table();
 127
 128         FILE *fp = fopen(argv[1], "rb");
 129         fseek(fp, 14, SEEK_SET);
 130
 131         int x = 0;
 132
 133         while (!feof(fp)) {
 134                 int len = getc(fp);
 135                 int cycles;
 136                 if (len == 0) {
 137                         int a = getc(fp);
 138                         int b = getc(fp);
 139                         int c = getc(fp);
 140                         cycles = a | (b << 8) | (c << 16);
 141                 } else {
 142                         cycles = len * 8;
 143                 }
 144                 pulse p;
 145                 p.start = float(x) * WAVE_FREQ / C64_FREQ;
 146                 p.end = (float(x) + cycles * 0.5) * WAVE_FREQ / C64_FREQ;
 147                 pulses.push_back(p);
 148                 x += cycles;
 149         }
 150
 151         int len_cycles = x;
 152         int len_samples = int(ceil(float(len_cycles) * WAVE_FREQ / C64_FREQ));
 153
 154         fprintf(stderr, "%d pulses, total %.2f seconds (%d samples)\n", pulses.size(), len_cycles / float(C64_FREQ), len_samples);
 155
 156         int pulse_begin = 0;
 157
 158         vector<float> samples;
 159         samples.reserve(len_samples);
 160
 161         for (int i = 0; i < len_samples; ++i) {
 162 //      for (int i = 50000000; i < 51000000; ++i) {
 163 //              double t = i * 0.01;
 164                 double t = i;
 165                 double sample = -0.5;
 166                 for (unsigned j = pulse_begin; j < pulses.size(); ++j) {
 167                         if (t - pulses[j].end > LANCZOS_EFFECTIVE_RADIUS) {
 168                                 ++pulse_begin;
 169                                 continue;
 170                         }
 171                         if (pulses[j].start - t > LANCZOS_EFFECTIVE_RADIUS) {
 172                                 break;
 173                         }
 174                         float contribution = window_integral(pulses[j].start - t, pulses[j].end - t);
 175                         sample += contribution;
 176                 }
 177                 samples.push_back(sample);
 178         }
 179
 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);
 184
 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]);
 188         }
 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]);
 192         }
 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);
 197         }
 198 }