// system include files
//
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

// isip include files
//
#define __ISIP_INTEGRAL
#define __ISIP_INTEGRAL_OLD

#include "integral.h"

// define external functions
//
extern void find_fft(int_4 fft_len, float_4* xr, float_4* xi, 
		     float_4* fft_real, float_4* fft_imag);

// define constants
//
# define PI 3.1415927
# define SAMPLE_FREQ 8000
# define SCALING_FACTOR 3.051850e-05
# define FFT_ORDER 1024


/* main program
   ------------ */

main(int argc, char **argv)
{
  if (argc != 5)
    {
      printf ("Wrong number of arguments! Error!\n");
      exit (-1);
    }
  
  // get variables from stdin
  //
  float_4 preemphasis        = atof (argv[1]);
  int_4 window_length        = (int_4) (atof(argv[2]) * SAMPLE_FREQ);
  int_4 window_center_time   = (int_4) (atof(argv[3]) * SAMPLE_FREQ);
  int_4 lp_order             = (int_4) (atoi(argv[4]));
  
  // declare variables
  //
  int_2 i,j;
  int_4 offset =(int_4)((window_center_time - window_length * 0.5) * 
			sizeof(int_2));
  int_4 no_of_samples = window_length;

  // allocate memory
  //
  int_2 *signal_buf = new (int_2)[no_of_samples];
  float_4 *signal_r = new (float_4)[FFT_ORDER];
  float_4 *signal_i = new (float_4)[FFT_ORDER];

  //read in the required values from the input file
  //
  fseek(stdin, offset, SEEK_SET);
  fread(signal_buf, sizeof(int_2), (int_4)(no_of_samples), stdin);
  
  /* de-bias, pre-emphasize and (hamming) window the data
     ---------------------------------------------------- */
  
  float_4 average = 0;
  for(i = 0; i < no_of_samples; i++)
    {
      average += (float_4)(signal_buf[i]);
    }
  average /= no_of_samples;
  
  float_4 win_factor  = 2.0 * PI / (no_of_samples - 1);
  float_4 signal_temp_minus_1 = 0;
  float_4 signal_temp = 0;

  float_4 *signal = new (float_4)[FFT_ORDER];
  for(i = 0; i < no_of_samples; i++)
    {
      signal_temp = (float_4(signal_buf[i]) - average) * SCALING_FACTOR;
      signal[i]  = signal_temp - preemphasis * signal_temp_minus_1;
      signal_temp_minus_1 = signal_temp;
      signal[i] *= (0.54 - 0.46 * cos ((float_4)i * win_factor)); 
    }
  delete signal_buf;
  
  
  /* compute auto-correlation values
     ------------------------------- */
  
  float_4 *r = new (float_4)[lp_order + 1];
  for(i = 0; i < lp_order + 1; i++)
    {
      r[i] = 0;
      for(j = 0; j < no_of_samples - i; j++)
	r[i] += (signal[j] * signal[j+i]);
    }
  
    
  /* The Levinson-Durbin Recursion for computing the lp coefficients
     --------------------------------------------------------------- */
  
  // define memory space for LP variables
  //
  
  // predictor coefficients
  //
  float_4 *a = new float_4[lp_order + 1];
  
  // reflection coefficients
  //
  float_4 *k = new float_4[lp_order + 1];
  
  // error coefficients
  //
  float_4 *E = new float_4[lp_order + 1];
  
  // temporary coefficients
  //
  float_4 *temp = new float_4[lp_order + 1];
  
  E[0] = r[0];
  
  for(i = 1; i < lp_order + 1; i++)
    {
      // compute reflection coefficients
      //
      k[i] = r[i];
      for(j = 1; j < i; j++)
	k[i] -= temp[j] * r[i-j];
      k[i] /= E[i-1];
      
      // compute predictor coefficients
      //
      a[i] = k[i];      
      for(j = 1; j < i; j++)
	a[j] = temp[j] - k[i] * temp[i-j];
      
      // error 
      //
      E[i] = (1 - k[i] * k[i]) * E[i-1];
            
      for(j = 0; j <= i; j++)
	temp[j] = a[j];
    }
  
  delete temp;
  delete k;
  delete E;
  delete r;

  float_4 gain = 10 * log10(E[lp_order]);
  
  /* generate values for e^(2*PI/N) for computing the fft and dft
     ------------------------------------------------------------ */
  
  // allocate memory for the constant values
  //
  float_4 *fft_real_coeff = new float_4[FFT_ORDER + 1];
  float_4 *fft_imag_coeff = new float_4[FFT_ORDER + 1];
  
  //compute the constant values
  //  
  float pi_factor = 2 * PI / FFT_ORDER;
  for(int_2 z = 0; z < FFT_ORDER; z++)
    {
      fft_real_coeff[z] = cos (pi_factor * z);
      fft_imag_coeff[z] = sin (pi_factor * z);
    } 
  
  /* generate the LPC model's spectrum
     --------------------------------- */

  float_4 *lp_dft = new float_4[FFT_ORDER+1]; 
  for (i = 0; i < FFT_ORDER + 1; i++)
    {
      float_4 real_part = 1.0;
      float_4 imag_part = 0;
      
      for (j = 1; j < lp_order + 1; j++) 
        {  
          real_part -= a[j] * fft_real_coeff[(i*j) % FFT_ORDER];
          imag_part += a[j] * fft_imag_coeff[(i*j) % FFT_ORDER];
        }
      lp_dft[i] = 10 * log10(real_part * real_part + imag_part * imag_part);
    }
  
  float_4 const_factor = float_4((float_4)SAMPLE_FREQ/(float_4)FFT_ORDER); 
  for (int z=1; z < FFT_ORDER/2  ; z++)
    fprintf(stdout,"%f %f\n", const_factor * z, gain - lp_dft[z]);
  
  fprintf(stdout,"\n");
  
  delete lp_dft;  
  
  /* zero-pad the remaining values of signal[i] for computing the fft
     ---------------------------------------------------------------- */
  
  for(i = 0; i < no_of_samples; i++)
    {
      signal_r[i] = signal[i];
      signal_i[i] = 0;
    }
  for(i = no_of_samples; i < FFT_ORDER; i++)
    {
      signal_r[i] = 0;
      signal_i[i] = 0;
    }
  
  /* compute the fft values
     ---------------------- */

  float_4 scale_factor = (float_4)SAMPLE_FREQ / (float_4)FFT_ORDER ;
  find_fft(FFT_ORDER, signal_r, signal_i, fft_real_coeff, fft_imag_coeff);
  
  for( i = 0; i < FFT_ORDER / 2 + 1; i++)
    {
      float_4 magnitude = 10 * log10 (signal_r[i] * signal_r[i] + 
				      signal_i[i] * signal_i[i]);
      fprintf(stdout, "%f %f\n", scale_factor * i, magnitude);
    }
  
  // free memory space
  //
  delete a;
  delete fft_real_coeff;
  delete fft_imag_coeff;
  delete signal_r;
  delete signal_i;
  delete signal;
}

/* Decimation-in-frequency FFT algorithm
   ------------------------------------- */

void find_fft(int_4 fft_len, float_4* xr, float_4* xi, float_4* fft_real,
	      float_4* fft_imag)
{
  int_4 m,l,i1,i2;
  int_4 n1,n2;
  float_4 real_twid,imag_twid;
  float_4 tempr, tempi; 
  
  m = (int)(log10(fft_len)/log10(2));
  n2 = fft_len;
  for(int i=1; i<m+1; i++)
    {
      n1 = n2;
      n2 = (int_4)(n2/2);
      i1 = 1;
      i2 = (int_4)((fft_len)/n1);
      for(int j=1; j<n2+1; j++)
	{
	  real_twid = fft_real[i1];
	  imag_twid = fft_imag[i1];
	  
	  for(int k = j; k <= fft_len; k = k + n1)
	    {
	      l = k + n2;
	      tempr = xr[k] - xr[l];
	      xr[k] = xr[k] + xr[l];
	      tempi = xi[k] - xi[l];
	      xi[k] = xi[k] + xi[l];
	      xr[l] = real_twid*tempr + imag_twid*tempi;
	      xi[l] = real_twid*tempi - imag_twid*tempr;
	    }
	  i1= i1 + i2;
	}
    }
  
  //procedure for bit reversal
  //
  float_4 temp = 0;
  int j = 1;
  for (int i = 1; i < fft_len; i++)
    {
      if (i < j)
	{
	  temp = xr[j];
	  xr[j] = xr[i];
	  xr[i] = temp;
	  temp = xi[j];
	  xi[j] = xi[i];
	  xi[i] = temp;
	}
      int_4 k = fft_len/2;
      while (k < j)
	{
	  j = j - k;
	  k = k/2;
	}
      j = j + k;
    }
}