DSP BLIT synthesis in C++

dsp_blit.cpp

// This code generates a band limited impulse train, then uses the table
// to generate an antialiased square wave at two different frequencies.
// The computational cost of using BLIT is far less than a sinc resampler
// because points are only evaluated on sample changes instead of every
// sampled point.

const static lc_u32 AUDIO_RATE = 44100;

void agb_test_audio_blit() {
  const int FRAME_AUDIO_SIZE = (AUDIO_RATE + 59) / 60;
  lc_i16 wav[2 * FRAME_AUDIO_SIZE];

  CoInitializeEx(nullptr, COINIT_MULTITHREADED);
  AudioDev aud;
  aud.Init(5, AUDIO_RATE);

  const int blit_points = 16;
  const int blit_points_mask = blit_points - 1;
  const int blit_resolution = 32;
  double blit_lut[blit_resolution][blit_points];

  const double cutoff = 0.8;
  for (int n = 0; n < blit_points; n++) {
    for (int m = 0; m < blit_resolution; m++) {
      double t = (m + 0.5) / double(blit_resolution);

      double x0 = M_PI * (t - n + blit_points / 2);
      double sinc = sin(cutoff * x0) / x0;

      double x1 = 2 * M_PI * (n + t) / blit_points;
      double blackman = 0.42659 - 0.49656 * cos(x1) + 0.076849 * cos(2 * x1);

      blit_lut[m][n] = sinc * blackman;
    }
  }

  for (int m = 0; m < blit_resolution; m++) {
    double kernel_sum = 0.0;
    for (int n = 0; n < blit_points; n++) {
      kernel_sum += blit_lut[m][n];
    }

    for (int n = 0; n < blit_points; n++) {
      blit_lut[m][n] /= kernel_sum;
    }

    for (int n = 1; n < blit_points; n++) {
      blit_lut[m][n] += blit_lut[m][n - 1];
    }
  }

  for (int m = 0; m < blit_resolution; m++) {
    double error = 1.0;
    double prev = 0.0;

    for (int n = 0; n < blit_points; n++) {
      double cur = blit_lut[m][n];

      double delta = cur - prev;
      error = error - delta;
      prev = cur;
      blit_lut[m][n] = delta;
    }

    blit_lut[m][blit_points / 2 - 1] += error * 0.5;
    blit_lut[m][blit_points / 2] += error * 0.5;
  }

  const double sub_dt = 1.0 / (AUDIO_RATE * blit_resolution);
  const lc_i16 vol = 1 << 13;

  int tap_p = 0, wav_p = 0;
  double taps[blit_points] = { 0 };
  double running_sum = 0, clock = 0, prev_input = 0;
  for (int second = 0; second < 2; ++second) {
    for (int s = 0; s < AUDIO_RATE; s++) {
      const static double duty_cycles[4] = { 0.125, 0.25, 0.5, 0.75 };
      const double duty = duty_cycles[2];
      const double freq = 131072. / (2048 - (second < 1 ? 0x7A0 : 0x3A0));
      const double delay = 1.0 / freq;

      for (int x = 0; x < blit_resolution; x++) {
        double phase = clock / delay;
        phase -= int(phase);
        clock = phase * delay + sub_dt;
        double input = phase < duty ? -1 : 1;

        if (input != prev_input) {
          double add = input - prev_input;
          for (int n = 0; n < blit_points; n++) {
            taps[(tap_p + n) & blit_points_mask] += add * blit_lut[x][n];
          }

          prev_input = input;
        }
      }

      running_sum = 0.997 * running_sum + taps[tap_p];
      taps[tap_p] = 0.0;
      tap_p = (1 + tap_p) & blit_points_mask;

      wav[2 * wav_p + 0] = wav[2 * wav_p + 1] = running_sum * vol;
      ++wav_p;
      if (wav_p == FRAME_AUDIO_SIZE) {
        aud.Write(wav);
        wav_p = 0;
      }
    }
  }

  aud.Release();
}