#include <stdlib.h>
#include <stdio.h>
#include <math.h>

/*
 *  We want to distort an image to give it a logarithmic vertical axis instead
 * a linear one.
 *
 * Each column of the input if the result of an FFT whose elements represent
 * the energy in the sound at frequencies from 0 Hz to max_freq_in Hz
 *
 * We want to produce an output from 0 to max_y_out pixels height
 * where the values represent frequencies from min_freq_out to max_freq_out, 
 *
 * For each value y_in, where would that frequency end up on the y axis?
 * freq = y_in * max_freq_in / max_y_in	// the frequency it represents
 * log(1) = 0, log(base) = 1, so scale input to 1..base range
 * freq01 = (freq - min_freq_out)/(max_freq_out - min_freq_out)
 *    is 0..1 for the values on the output axis
 * base1 = (freq01 * (base - 1)) + 1
 *    is 1..base
 * log(base1) is 0..1 so
 * y_out = log(base1) * max_y_out
 *
 * To be able to generate the output image, we need to find
 * the reverse mepping function from output y coordinates to
 * input y coordinates. First, inline all the terms:
 *
 * y_out = max_y_out * log( 1 + (((y_in * max_freq_in / max_y_in) - min_freq_out)/(max_freq_out - min_freq_out)) * (base - 1) )
 *
 * and rearrange it until we have y_in = ...
 *
 * exp(y_out / max_y_out) - 1 = (((y_in * max_freq_in / max_y_in) - min_freq_out)/(max_freq_out - min_freq_out)) * (base - 1)
 * (exp(y_out / max_y_out) - 1) / (base - 1) = ((y_in * max_freq_in / max_y_in) - min_freq_out)/(max_freq_out - min_freq_out)
 * ((exp(y_out / max_y_out) - 1) / (base - 1)) * (max_freq_out - min_freq_out) = (y_in * max_freq_in / max_y_in) - min_freq_out
 * ((exp(y_out / max_y_out) - 1) / (base - 1)) * (max_freq_out - min_freq_out) + min_freq_out = (y_in * max_freq_in / max_y_in)
 * y_in = (((exp(y_out / max_y_out) - 1) / (base - 1)) * (max_freq_out - min_freq_out) + min_freq_out) * max_y_in / max_freq_in

 * Or... calculate the frequency value of the output pixel and
 * map that back to the coresponding input pixel
 * freq = 50 * (3150/50)^(y_out/max_y_out)
 * y_in = max_y_in * (freq / max_freq_in)
 */

main(int argc, char **argv)
{
  int max_freq_in = 22050;
  int max_y_in = 8191;
  double min_freq_out = 50;
  double max_freq_out = 3150;
  int max_y_out = 1023;
  double y_out;
  double y_in;
  double y_in2;
  double freq;
#define base M_E

  printf(" Out   In   In2  Freq\n");
  for (y_out = 0; y_out <= max_y_out; y_out += 93) {
    y_in = (((exp(y_out / max_y_out) - 1) / (base - 1)) * (max_freq_out - min_freq_out) + min_freq_out) * max_y_in / max_freq_in;
    freq = min_freq_out * pow((max_freq_out/min_freq_out), y_out / max_y_out);
    y_in2 = max_y_in * (freq / max_freq_in);

    freq = y_in * max_freq_in / max_y_in;
    printf("%4ld %4ld %4ld %5ld\n", lrint(y_out), lrint(y_in), lrint(y_in2), lrint(freq));
  }
  exit(0);
}