BULLSHIT

/* This code is licensed under The Anti-Viral License! */
// License: http://www.ecstaticlyrics.com/secret/The%20Anti-Viral%20License
// Information:  http://www.ecstaticlyrics.com/pinnwand/messages/the_anti_viral_license.html

// This function reads 128 ints from the pointer, un-filters the signal,
// then returns the unfiltered signal by writing it to the same pointer.
// Note that there is one second of lag, thus the data returned by the
// first call to this function is garbage.

void unfilter_magic(int *signal) {
  static double input[256] = {};
  static double output[256] = {};

  // First, copy signal to a buffer so we can work in floating point.

  for (int i = 0; i < 128; i++) {
    input[i + 128] = signal[i];
  };

  // Use a sliding window and average together the results.

  #define DIVISIONS 4
  #define STEP_SIZE (128 / DIVISIONS)
  //#define WINDOW 1.0
  #define WINDOW (2.0 * (1.0 - fabs(i - 64) / 64))
  //#define WINDOW (1.0 - cos(M_PI * (i - 64) / 64.0))
  //#define WINDOW (M_PI * sin(M_PI * i / 128.0) / 2.0)

  for (int d = 0; d < DIVISIONS; d++) {
    for (int f = 1; f <= 63; f++) {

      // For frequencies from 1 to 63 Hz, find their amplitude and phase.

      double x = 0, y = 0;
      for (int i = 0; i < 128; i++) {
        x += WINDOW * input[i] * sin(2 * M_PI * f * i / 128.0);
        y += WINDOW * input[i] * cos(2 * M_PI * f * i / 128.0);
      };
      x = 2 * x / 128.0;
      y = 2 * y / 128.0;
      double amplitude = sqrt(x * x + y * y);
      double phase = atan2(y, x);

      // Scale amplitude to remove measured attenuation from filters.

      amplitude /= measurements[f-1][0];

      // Adjust phase to reverse measured phase shift of filters.

      phase -= measurements[f-1][1];

      // Also, filter out shit over 20 Hz.  No one cares about that shit.

      if (f >= 40) {
        amplitude = 0;
      } else if (f >= 20) {
        amplitude *= (f - 40.0) / -20.0;
      };

      // Recreate waveform from the adjusted frequencies and phases.

      for (int i = 0; i < 128; i++) {
        output[i + 128] += WINDOW * amplitude * sin(2 * M_PI * f * i / 128.0 + phase);
      };

    };

    // Slide input and output windows...

    for (int i = 0; i < 256 - STEP_SIZE; i++) {
      input[i] = input[i + STEP_SIZE];
      output[i] = output[i + STEP_SIZE];
    };

    // Zero the addition to the end of the output window.

    for (int i = 256 - STEP_SIZE; i < 256; i++) {
      output[i] = 0;
    };

    // Return data...

    for (int i = 0; i < 128; i++) {
      signal[i] = output[i] / DIVISIONS;
    };

  };

};

static double measurements[63][2] = {
  {0.359, 0.000},
  {0.683, 0.657},
  {0.818, 1.600},
  {0.879, 2.734},
  {0.911, -2.379},
  {0.928, -1.133},
  {0.938, 0.098},
  {0.944, 1.396},
  {0.948, 2.633},
  {0.949, -2.340},
  {0.950, -1.055},
  {0.949, 0.188},
  {0.947, 1.506},
  {0.945, 2.806},
  {0.943, -2.204},
  {0.940, -0.899},
  {0.937, 0.361},
  {0.933, 1.717},
  {0.929, 3.064},
  {0.926, -2.017},
  {0.921, -0.713},
  {0.917, 0.646},
  {0.913, 2.013},
  {0.908, -3.054},
  {0.903, -1.755},
  {0.897, -0.440},
  {0.890, 0.902},
  {0.884, 2.186},
  {0.876, -2.834},
  {0.868, -1.529},
  {0.859, -0.178},
  {0.849, 1.067},
  {0.839, 2.359},
  {0.828, -2.617},
  {0.817, -1.210},
  {0.806, 0.023},
  {0.794, 1.363},
  {0.783, 2.626},
  {0.771, -2.343},
  {0.760, -0.924},
  {0.748, 0.366},
  {0.736, 1.585},
  {0.725, 2.903},
  {0.713, -2.048},
  {0.701, -0.796},
  {0.687, 0.480},
  {0.673, 1.671},
  {0.659, -3.094},
  {0.643, -1.860},
  {0.625, -0.571},
  {0.604, 0.597},
  {0.586, 2.153},
  {0.562, -2.970},
  {0.534, -1.682},
  {0.506, -0.382},
  {0.476, 1.066},
  {0.453, 2.456},
  {0.415, -2.546},
  {0.380, -1.320},
  {0.386, -0.032},
  {0.302, 1.195},
  {0.282, 2.749},
  {0.247, -2.418},
};