#include <stdio.h>

 #include <stdlib.h>

 #include <math.h>

 #include <assert.h>

 #include <string.h>

 #include <math.h>

 #include "SndfileWavUtil.h"

 

 const  int Taps = 9;

 const  int Points = 8;

 const  int log_Points = 3;

 static int newest;

 static short int x[Taps] = {0,0,0,0,0,0,0,0,0};

 static short int fir(short int);

 

 // foreward declaration. This is defined in DFT.cpp

 int DFT(int dir,int m,double *x1,double *y1);

 

 

 static int x_to_y(int x, int y)

 {

   assert(y>0);

   int rv = x;

   while(y-->1)

     x *= x;

   return rv;

 }

 

 // function implemented in LinearWindowingFunction.bsv

 static double windowing_fn(int index)

 {

   double pi = 3.14159265;

   double i = (double)index;

   return sin(0.5*pi* (sin((i + 0.5)/Points * pi) * sin((i + 0.5)/Points * pi)));

 }

 

 // lifted from the bluespec code

 static double h[Taps] =  {-0.0124,

 			  0.0,

 			  -0.0133,

 			  0.0,

 			  0.8181,

 			  0.0,

 			  -0.0133,

 			  0.0,

 			  -0.0124};

 

 

 short int

 fir(short int sample)

 {

   x[newest] = sample; 

   double y = 0;

   for (int k = 0; k < Taps; k++) { 

     short int bits = x[(newest+k)%Taps];

     double x = (double)bits;

     y += h[k] * x; 

   }

   newest = (newest == 0) ? (Taps-1) : (newest-1);

   short int rv = ((short int)y)&0x0000FFFF;

   return rv;

 }

 

 

 int

 main (int argc, char * argv [])

 {

   const char* inputWavFileName;

   const char* outputWavFileName;

 

   FILE *inputPcmFile;

   FILE *firPcmFile;

   FILE *ifftPcmFile;

   FILE *outputPcmFile;

 

   short int sample;

   short int samples[Points];

   const unsigned int short_sz = sizeof(short int);

 

   inputWavFileName = argv[1];

   outputWavFileName = argv[2];

 

   // Convert input wav to pcm

   generate_pcm(inputWavFileName, "input.pcm");

 

   inputPcmFile  = fopen("input.pcm", "r");

   firPcmFile    = fopen("fir.pcm", "w");

 

   assert(inputPcmFile);

   assert(firPcmFile);

 

   newest = 0;

   memset(x, 0, sizeof(x));

 

   while(fread(&sample, short_sz, 1, inputPcmFile)) {

     sample = fir(sample);

     assert(fwrite(&sample,short_sz, 1, firPcmFile));

   }

 

   fclose(inputPcmFile);

   fclose(firPcmFile);

 

   firPcmFile  = fopen("fir.pcm", "r");

   ifftPcmFile = fopen("ifft.pcm", "w");

 

   assert(firPcmFile);

   assert(ifftPcmFile);

 

   int read = 0;

 

   // read in half a frame

   for(int i = 0; i < Points/2; i++)

     read += fread(&samples[(Points/2)+i], short_sz, 1, firPcmFile);

 

   if(read > 0)

     assert(read==Points/2);

 

   // we will do an 'Points' point fft, and then undo it

   while(true) {

 

     read = 0;

 

     // shift last half frame to head

     for(int i = 0; i < Points/2; i++)

       samples[i] = samples[(Points/2)+i];

 

     // read in another half frame

     for(int i = 0; i < Points/2; i++)

       read += fread(&samples[(Points/2)+i], short_sz, 1, firPcmFile);

 

     if(read == 0)

       break;

 

     // pad out the 

     if(read < Points/2){

       for(int i = read; i < Points/2; i++)

 	samples[Points/2+i] = 0;

     }

 

     double dsamples[Points];

     double dimag[Points];

 

     // this shift is performed in the Bluespec

     for(int i = 0; i < Points; i++){

       dsamples[i] = (double)(samples[i]>>log_Points);

       dimag[i] = 0.0;

     }

 

 

     //for(int i = 0; i < Points; i++)

     //  printf("%d ", (int)dsamples[i]);

     //printf("\n");

 

     DFT(1,Points,dsamples,dimag);

 

     // shift up by 'shift' (in)

     int shift = 2;

 

     for(int i = 0; i < Points/2; i++){

 

       dsamples[Points/2-1-i] = dsamples[Points/2-1-i-shift];

       dimag[Points/2-1-i] = dimag[Points/2-1-i-shift];

 

 

       dsamples[Points/2+i] = dsamples[Points/2+i+shift];

       dimag[Points/2+i] = dimag[Points/2+i+shift];

 

     }

 

     // fill in the lower points

     for(int i = 0; i < shift; i++){

 

       dsamples[i] = 0.0;

       dimag[i] = 0.0;

 

       dsamples[Points-1-i] = 0.0;

       dimag[Points-1-i] = 0.0;

     }

 

     DFT(-1,Points,dsamples,dimag);

 

     short int ifftResult[Points];

     for(int i = 0; i < Points; i++)

       ifftResult[i] = (int)dsamples[i];

 

     // strip off the padding

     int write = (read < Points/2) ? read+(Points/2) : Points;

     for(int i = 0; i < write; i++)

       assert(fwrite(&ifftResult[i],short_sz,1,ifftPcmFile));

 

     // break out if we're at the end

     if(read < Points/2){

        break;

     }

   }

 

   fclose(firPcmFile);

   fclose(ifftPcmFile);

 

   // do the windowing

 

   ifftPcmFile   = fopen("ifft.pcm", "r");

   outputPcmFile = fopen("output.pcm", "w");

 

   assert(ifftPcmFile);

   assert(outputPcmFile);

 

   short int samplesA[Points];

   short int samplesB[Points];

 

   bool valid_stream = true;

   // write out the first half frame

   for(int i = 0; i < Points/2; i++){

     valid_stream &= fread(samplesA, short_sz, 1,ifftPcmFile);

     valid_stream &= fwrite(samplesA, short_sz,1,outputPcmFile);

   }

 

   while(valid_stream){

 

     for(int i = 0; i < Points/2; i++)

       valid_stream &= fread(&samplesA[i], short_sz, 1,ifftPcmFile);

 

     for(int i = 0; i < Points; i++)

       valid_stream &= fread(&samplesB[i], short_sz, 1,ifftPcmFile);

 

     for(int i = 0; i < Points/2; i++)

       valid_stream &= fread(&samplesA[(Points/2)+i], short_sz, 1,ifftPcmFile); 

 

 

     // this isn't quite right

     if(!valid_stream)

       break;

 

     for(int i = 0; i < Points; i++){

       double window = windowing_fn(i);

       short int rv = (short int)(((double)samplesA[i])*(1.0-window) + window*((double)samplesB[i]));

       assert(fwrite(&rv,short_sz,1,outputPcmFile));    

     }

 

   }

 

   fclose(ifftPcmFile);

   fclose(outputPcmFile);

 

   generate_wav("output.pcm", inputWavFileName, outputWavFileName);

 

 }