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

#define BLOCK_SIZE 44100
#define abs(n) sqrt((n)[0]*(n)[0]+(n)[1]*(n)[1])

typedef signed short int16;

const char *temp_dat = "temp.dat";
const char *title0 = "E[|X(f)|]";
const char *title1 = "E[|X(f)|]dB";
const char *gnuplot =
   "set terminal png notransparent enhanced size 640,480\n"
   "set output \"%s\"\n"
   "%s set logscale x\n"
   "plot \"%s\" ti \"%s\" with lines\n";
const char *help =
   "%s input output [x_axis y_axis]\n"
   "\tinput\t\twav input\n"
   "\toutput\t\tpng output\n"
   "\tx_axis\t\t0 -> f; 1 -> log(f)\n"
   "\ty_axis\t\t0 -> E[|X(f))|]; 1 -> E[|X(f)|]dB\n";

int main(int argc, char **argv) {
   FILE *fp;
   int read, i, n = 0, x_axis = 0, y_axis = 0;
   fftw_complex *input, *transform;
   fftw_plan p;
   int16 audio_data[BLOCK_SIZE];
   float abs_sum[BLOCK_SIZE];
   
   if(argc < 3) {
      printf(help, argv[0]);
      return -1;
   }
   if(argc >= 4)
      x_axis = atoi(argv[3]);
   if(argc >= 5)
      y_axis = atoi(argv[4]);
   
   fp = fopen(argv[1], "rb");
   fseek(fp, 44, SEEK_SET);
   
   input = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * BLOCK_SIZE);
   transform = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * BLOCK_SIZE);
   
   for(i = 0; i < BLOCK_SIZE; i++) {
      abs_sum[i] = 0;
   }
   
   read = fread(audio_data, sizeof(int16), BLOCK_SIZE, fp);
   while(read == BLOCK_SIZE) {
      for(i = 0; i < BLOCK_SIZE; i++) {
         input[i][0] = audio_data[i];
         input[i][1] = 0;
      }
      p = fftw_plan_dft_1d(BLOCK_SIZE, input, transform, FFTW_FORWARD, FFTW_ESTIMATE);
      fftw_execute(p);
      for(i = 20; i < 20000 /*BLOCK_SIZE*/; i++) {
         abs_sum[i] += abs(transform[i]);
      }
      fftw_destroy_plan(p);
      read = fread(audio_data, sizeof(int16), BLOCK_SIZE, fp);
      n++;
   }
   fftw_free(input);
   fftw_free(transform);
   fclose(fp);
   
   fp = fopen(temp_dat, "w");
   for(i = 20; i < 20000/*BLOCK_SIZE*/; i++) {
      abs_sum[i] /= n;
      fprintf(fp, "%d %.5f\n", i, y_axis ? 20*log10(abs_sum[i]) : abs_sum[i]);
   }
   fclose(fp);
   
   fp = popen("gnuplot", "w");
   fprintf(fp, gnuplot, argv[2], x_axis ? "" : "#", temp_dat, y_axis ? title1 : title0);
   pclose(fp);
   
   remove(temp_dat);
   
   return 0;
}