/* plot_fft.c  part of series, also .m  .py  */
/*             uses  execl  to run gnuplot   */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <unistd.h>
#include "fft16.h"
#include "ifft16.h"

int main(int argc, char * argv[])
{
  int nfft = 16; /* number of real values, number of complex values */
  int i, j, stat;
  float b16[16], ub16[16];
  float cxa16[32], cxb16[32], cxc16[32]; /* 16 complex values */
  double a, b;
  float a16[] = {0.0, 0.3826, 0.7071, 0.9238, 1.0, 0.9238, 0.7071, 0.3826,
                 0.0,-0.3826,-0.7071,-0.9238,-1.0,-0.9238,-0.7071,-0.3826};
  float aa[128], bb[128], ubb[128];
  float cxbb[256]; /* 128 complex values */


  FILE * outp, * inp;
  /* OK, this could have been just a file, generated here. */
  char line1[] = {"set title \"16 sample wave amplitude\" "};
  char line2[] = {"set xrange[0:15]"};
  char line3[] = {"set xlabel \"sample number, time\" "};
  char line4[] = {"set ylabel \"wave amplitude\" "};
  char line5[] = {"set term jpeg "};
  char line6[] = {"set output \"plot_fft_c0.jpg\" "};
  char line7[] = {"plot \"temp.dat\" with lines"};
  char line8[] = {"clear"};
  char line9[]  = {"set title \"FFT abs magnitude \" "};
  char line10[] = {"set xrange[0:15]"};
  char line11[] = {"set xlabel \"frequency bin\" "};
  char line12[] = {"set ylabel \"abs spectrum\" "};
  char line13[] = {"set term jpeg "};
  char line14[] = {"set output \"plot_fft_c1.jpg\" "};
  char line15[] = {"plot \"temp1.dat\" with lines"};
  char line16[] = {"replot"};
  char line17[] = {"clear"};
  char line18[] = {"set title \"FFT magnitude, unambiguous \" "};
  char line19[] = {"set xrange[0:15]"};
  char line20[] = {"set xlabel \"frequency bin\" "};
  char line21[] = {"set ylabel \"abs folded spectrum\" "};
  char line22[] = {"set term jpeg "};
  char line23[] = {"set output \"plot_fft_c2.jpg\" "};
  char line24[] = {"plot \"temp2.dat\" with lines"};
  char line25[] = {"replot"};
  char line26[] = {"clear"};
  char line27[] = {"set title \"tdoa_cwfm1.dat\" "};
  char line28[] = {"set xrange[0:127]"};
  char line29[] = {"set xlabel \"sample number, time\" "};
  char line30[] = {"set ylabel \"wave amplitude\" "};
  char line31[] = {"set term jpeg "};
  char line32[] = {"set output \"plot_fft_c3.jpg\" "};
  char line33[] = {"plot \"temp3.dat\" with lines"};
  char line34[] = {"replot"};
  char line35[] = {"clear"};
  char line36[] = {"set title \"FFT tdoa abs magnitude \" "};
  char line37[] = {"set xrange[0:127]"};
  char line38[] = {"set xlabel \"frequency bin\" "};
  char line39[] = {"set ylabel \"abs spectrum\" "};
  char line40[] = {"set term jpeg "};
  char line41[] = {"set output \"plot_fft_c4.jpg\" "};
  char line42[] = {"plot \"temp4.dat\" with lines"};
  char line43[] = {"replot"};
  char line44[] = {"clear"};
  char line45[] = {"set title \"FFT tdoa magnitude, unambiguous \" "};
  char line46[] = {"set xrange[0:127]"};
  char line47[] = {"set xlabel \"frequency bin\" "};
  char line48[] = {"set ylabel \"abs folded spectrum\" "};
  char line49[] = {"set term jpeg "};
  char line50[] = {"set output \"plot_fft_c5.jpg\" "};
  char line51[] = {"plot \"temp5.dat\" with lines"};
  char line52[] = {"replot"};

  printf("plot_fft.c running \n");

  printf("a16=\n");
  for(i=0; i<4; i++) printf("%f, %f, %f, %f \n",
                     a16[4*i], a16[4*i+1], a16[4*i+2], a16[4*i+3]);
  outp = fopen("temp.dat", "w");
  for(i=0; i<16; i++) fprintf(outp, "%f \n", a16[i]);
  fclose(outp);

  for(i=0; i<nfft; i++)
  {
    cxb16[2*i] = a16[i];
    cxb16[2*i+1] = 0.0;
  }
  fft16(cxb16); /* complex FFT in place */
  printf("FFT of a16, complex cxb16=\n");
  for(i=0; i<4; i++) printf("(%f,%f), (%f,%f), (%f,%f), (%f,%f) \n", 
                            cxb16[8*i], cxb16[8*i+1], cxb16[8*i+2], cxb16[8*i+3],
                            cxb16[8*i+4], cxb16[8*i+5], cxb16[8*i+6], cxb16[8*i+7]);

  for(i=0; i<nfft; i++)
  {
    cxc16[2*i] = cxb16[2*i];
    cxc16[2*i+1] = cxb16[2*i+1];
    b16[i] = sqrt(cxb16[2*i]*cxb16[2*i]+cxb16[2*i+1]*cxb16[2*i+1])/
      (double)nfft;
  }
  printf("absolute value of FFT, b16= \n");
  for(i=0; i<4; i++) printf("%f, %f, %f, %f \n",
                     b16[4*i], b16[4*i+1], b16[4*i+2], b16[4*i+3]);
  outp = fopen("temp1.dat", "w");
  for(i=0; i<16; i++) fprintf(outp, "%f \n", b16[i]);
  fclose(outp);


  ifft16(cxc16);
  printf("inverse FFT of cxb16, complex cxc16=\n");
  for(i=0; i<4; i++) printf("(%f,%f), (%f,%f), (%f,%f), (%f,%f) \n", 
                            cxc16[8*i], cxc16[8*i+1], cxc16[8*i+2], cxc16[8*i+3],
                            cxc16[8*i+4], cxc16[8*i+5], cxc16[8*i+6], cxc16[8*i+7]);

  printf("just real part of cxc16, should be a16=\n");
  for(i=0; i<4; i++) printf("%f, %f, %f, %f \n",
                     cxc16[8*i], cxc16[8*i+2], cxc16[8*i+4], cxc16[8*i+6]);

  ub16[0] = 0.0;
  ub16[nfft/2] = 0.0;
  for(i=1; i<nfft/2; i++)
  {
    a = cxb16[2*i]+cxb16[2*(nfft-i)];
    b = cxb16[2*i+1]-cxb16[2*(nfft-i)+1];
    ub16[i] = sqrt(a*a+b*b)/(double)nfft;
    ub16[nfft/2+i] = 0.0;
  }
  outp = fopen("temp2.dat", "w");
  for(i=0; i<nfft; i++) fprintf(outp, "%f \n", ub16[i]);
  fclose(outp);

  nfft = 128;
  inp = fopen("tdoa_cwfm1.dat","r");
  for(i=0; i<nfft; i++)
  {
    stat = fscanf(inp, "%f", &aa[i]);
  }
  fclose(inp);
  outp = fopen("temp3.dat", "w");
  for(i=0; i<nfft; i++) fprintf(outp, "%f \n", aa[i]);
  fclose(outp);

  for(i=0; i<nfft; i++)
  {
    cxbb[2*i] = aa[i];
    cxbb[2*i+1] = 0.0;
  }
  fft128(cxbb); /* FFT in place */

  for(i=0; i<nfft; i++)
  {
    bb[i] = sqrt(cxbb[2*i]*cxbb[2*i]+cxbb[2*i+1]*cxbb[2*i+1])/
      (double)nfft;
  }
  outp = fopen("temp4.dat", "w");
  for(i=0; i<nfft; i++) fprintf(outp, "%f \n", bb[i]);
  fclose(outp);

  ubb[0] = 0.0;
  ubb[nfft/2] = 0.0;
  for(i=1; i<nfft/2; i++)
  {
    a = cxbb[2*i]+cxbb[2*(nfft-i)];
    b = cxbb[2*i+1]-cxbb[2*(nfft-i)+1];
    ubb[i] = sqrt(a*a+b*b)/(double)nfft;
    ubb[nfft/2+i] = 0.0;
  }
  outp = fopen("temp5.dat", "w");
  for(i=0; i<nfft; i++) fprintf(outp, "%f \n", ubb[i]);
  fclose(outp);

  outp = fopen("temp.plot","w");
  fprintf(outp, "%s \n", line1);
  fprintf(outp, "%s \n", line2);
  fprintf(outp, "%s \n", line3);
  fprintf(outp, "%s \n", line4);
  fprintf(outp, "%s \n", line5);
  fprintf(outp, "%s \n", line6);
  fprintf(outp, "%s \n", line7);
  fprintf(outp, "%s \n", line8);
  fprintf(outp, "%s \n", line9);
  fprintf(outp, "%s \n", line10);
  fprintf(outp, "%s \n", line11);
  fprintf(outp, "%s \n", line12);
  fprintf(outp, "%s \n", line13);
  fprintf(outp, "%s \n", line14);
  fprintf(outp, "%s \n", line15);
  fprintf(outp, "%s \n", line16);
  fprintf(outp, "%s \n", line17);
  fprintf(outp, "%s \n", line18);
  fprintf(outp, "%s \n", line19);
  fprintf(outp, "%s \n", line20);
  fprintf(outp, "%s \n", line21);
  fprintf(outp, "%s \n", line22);
  fprintf(outp, "%s \n", line23);
  fprintf(outp, "%s \n", line24);
  fprintf(outp, "%s \n", line25);
  fprintf(outp, "%s \n", line26);
  fprintf(outp, "%s \n", line27);
  fprintf(outp, "%s \n", line28);
  fprintf(outp, "%s \n", line29);
  fprintf(outp, "%s \n", line30);
  fprintf(outp, "%s \n", line31);
  fprintf(outp, "%s \n", line32);
  fprintf(outp, "%s \n", line33);
  fprintf(outp, "%s \n", line34);
  fprintf(outp, "%s \n", line35);
  fprintf(outp, "%s \n", line36);
  fprintf(outp, "%s \n", line37);
  fprintf(outp, "%s \n", line38);
  fprintf(outp, "%s \n", line39);
  fprintf(outp, "%s \n", line40);
  fprintf(outp, "%s \n", line41);
  fprintf(outp, "%s \n", line42);
  fprintf(outp, "%s \n", line43);
  fprintf(outp, "%s \n", line44);
  fprintf(outp, "%s \n", line45);
  fprintf(outp, "%s \n", line46);
  fprintf(outp, "%s \n", line47);
  fprintf(outp, "%s \n", line48);
  fprintf(outp, "%s \n", line49);
  fprintf(outp, "%s \n", line50);
  fprintf(outp, "%s \n", line51);
  fprintf(outp, "%s \n", line51);
  fclose(outp);
  printf("plot_fft.c finished, should have .jpg plots \n");
  fflush(stdout);
  fclose(stdout);

  stat = execl("/usr/bin/gnuplot","gnuplot","-persist","temp.plot", NULL);
  /* only returns if error */
  printf("plot_fft.c finished, should have plots \n");
} /* end plot_fft.c */