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 // Nyquist frequency of low-pass filter (must be max. WAVE_FREQ / 2 or you
  15 // will get aliasing)
  16 #define LPFILTER_FREQ 8000
  17
  18 #define NORMALIZED_LPFILTER_FREQ (float(LPFILTER_FREQ) / float(WAVE_FREQ))
  19 #define LANCZOS_EFFECTIVE_RADIUS (LANCZOS_RADIUS / (NORMALIZED_LPFILTER_FREQ * 2.0))
  20
  21 #define SIZEOF_HALF_TABLE ((int)((LANCZOS_EFFECTIVE_RADIUS * RESOLUTION)))
  22 static double *integrated_window;
  23
  24 double sinc(double x)
  25 {
  26         if (fabs(x) < 1e-6) {
  27                 return 1.0f - fabs(x);
  28         } else {
  29                 return sin(x) / x;
  30         }
  31 }
  32
  33 double weight(double x)
  34 {
  35         if (fabs(x) > LANCZOS_EFFECTIVE_RADIUS) {
  36                 return 0.0f;
  37         }
  38         float t = 2.0 * M_PI * NORMALIZED_LPFILTER_FREQ * x;
  39         return sinc(t) * sinc(t / LANCZOS_RADIUS);
  40 }
  41
  42 void make_table()
  43 {
  44         integrated_window = new double[2 * SIZEOF_HALF_TABLE + 1];
  45
  46         double sum = 0.0;
  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);
  52         }
  53 //      exit(0);
  54 }
  55
  56 // integral from -inf to window
  57 double window_one_sided_integral(double to)
  58 {
  59         double array_pos = to * RESOLUTION + SIZEOF_HALF_TABLE;
  60
  61         int whole = int(floor(array_pos));
  62         double frac = array_pos - whole;
  63         if (whole < 0) {
  64                 return 0.0;
  65         }
  66         if (whole >= 2 * SIZEOF_HALF_TABLE) {
  67                 return 1.0;
  68         }
  69         return integrated_window[whole] + frac * (integrated_window[whole + 1] - integrated_window[whole]);
  70 }
  71
  72 double window_integral(double from, double to)
  73 {
  74         return window_one_sided_integral(to) - window_one_sided_integral(from);
  75 }
  76
  77 struct pulse {
  78         // all values in samples
  79         double start, end;
  80 };
  81 vector<pulse> pulses;
  82
  83 int main(int argc, char **argv)
  84 {
  85         make_table();
  86
  87         FILE *fp = fopen(argv[1], "rb");
  88         fseek(fp, 14, SEEK_SET);
  89
  90         int x = 0;
  91
  92         while (!feof(fp)) {
  93                 int len = getc(fp);
  94                 int cycles;
  95                 if (len == 0) {
  96                         int a = getc(fp);
  97                         int b = getc(fp);
  98                         int c = getc(fp);
  99                         cycles = a | (b << 8) | (c << 16);
 100                 } else {
 101                         cycles = len * 8;
 102                 }
 103                 pulse p;
 104                 p.start = float(x) * WAVE_FREQ / C64_FREQ;
 105                 p.end = (float(x) + cycles * 0.5) * WAVE_FREQ / C64_FREQ;
 106                 pulses.push_back(p);
 107                 x += cycles;
 108         }
 109
 110         int len_cycles = x;
 111         int len_samples = int(ceil(float(len_cycles) * WAVE_FREQ / C64_FREQ));
 112
 113         fprintf(stderr, "%d pulses, total %.2f seconds (%d samples)\n", pulses.size(), len_cycles / float(C64_FREQ), len_samples);
 114
 115         int pulse_begin = 0;
 116
 117         for (int i = 0; i < len_samples; ++i) {
 118 //      for (int i = 50000000; i < 51000000; ++i) {
 119 //              double t = i * 0.01;
 120                 double t = i;
 121                 double sample = 0.0;
 122                 for (unsigned j = pulse_begin; j < pulses.size(); ++j) {
 123                         if (t - pulses[j].end > LANCZOS_EFFECTIVE_RADIUS) {
 124                                 ++pulse_begin;
 125                                 continue;
 126                         }
 127                         if (pulses[j].start - t > LANCZOS_EFFECTIVE_RADIUS) {
 128                                 break;
 129                         }
 130                         float contribution = window_integral(pulses[j].start - t, pulses[j].end - t);
 131                         sample += contribution;
 132                 }
 133                 printf("%f\n", sample);
 134 //              short s = lrintf((sample - 0.5f) * 16384.0f);
 135 //              fwrite(&s, 2, 1, stdout);
 136         }
 137 }