//fastconvosim.c Overlap-add fast convolution demonstration program
//Uses FFT function from C31 DSK book. Not real-time, and therefore
//does not need McBSP I/O. Run with breakpoints inserted as indicated
//in source file and view graphs of

// start     acquisition   display    index      DSP
// address   buffer size   data size  increment  data type

// h         512           256        2          32-bit float
// samples   512           256        2          32-bit float
// iobuffer  128           128        1          32-bit float
// overlap   128           128        1          32-bit float

#include <c6x.h>
#include <math.h>
#include <stdio.h>
#include "coeffs.h"               //time domain FIR filter coefficients
#define PTS 256                   //number of points used in FFT
 
struct cmpx                       //complex data structure used by FFT
    {
    float real;
    float imag;
    };
typedef struct cmpx COMPLEX;
#define PI 3.14159265358979
short poweroftwo[] = {1,2,4,8,16,32,64,128,256,512,1024};  
short buffercount = 0;            //number of new input samples in iobuffer 
float iobuffer[PTS/2];            //primary input/output buffer
COMPLEX samples[PTS];             //primary working buffer used by FFT
float overlap[PTS/2];             //intermediate result buffer
COMPLEX w[PTS];                   //twiddle factors stored in w
COMPLEX h[PTS];                   //freq domain filter coeffs stored in h
short i;                          //general purpose index variable 
float a,b;                        //variables used in complex multiply
short NUMCOEFFS = sizeof(coeffs)/sizeof(float);

void FFT(COMPLEX *Y, int N)       //input sample array, # of points
{
  COMPLEX temp1,temp2;            //temporary storage variables
  int i,j,k;                      //loop counter variables
  int upper_leg, lower_leg;       //index of upper/lower butterfly leg
  int leg_diff;                   //difference between upper/lower leg
  int num_stages=0;               //number of FFT stages, or iterations
  int index, step;                //index and step between twiddle factor
  i=1;                            //log(base 2) of # of points = # of stages
  do
  {
    num_stages+=1;
    i=i*2;
  } while (i!=N);

  leg_diff=N/2;                   //starting difference between upper & lower legs
  step=512/N;                     //step between values in twiddle.h              
  for (i=0;i<num_stages;i++)      //for N-point FFT                 
  {
    index=0;
    for (j=0;j<leg_diff;j++)
    {
      for (upper_leg=j;upper_leg<N;upper_leg+=(2*leg_diff))
      {
        lower_leg=upper_leg+leg_diff;
        temp1.real=(Y[upper_leg]).real + (Y[lower_leg]).real;
        temp1.imag=(Y[upper_leg]).imag + (Y[lower_leg]).imag;
        temp2.real=(Y[upper_leg]).real - (Y[lower_leg]).real;
        temp2.imag=(Y[upper_leg]).imag - (Y[lower_leg]).imag;
        (Y[lower_leg]).real=temp2.real*(w[index]).real-temp2.imag*(w[index]).imag;
        (Y[lower_leg]).imag=temp2.real*(w[index]).imag+temp2.imag*(w[index]).real;
        (Y[upper_leg]).real=temp1.real;
        (Y[upper_leg]).imag=temp1.imag;
      }
      index+=step;
    }
    leg_diff=leg_diff/2;
    step*=2;
  }
  j=0;
  for (i=1;i<(N-1);i++)           //bit reversal for resequencing data*/
  {
    k=N/2;
    while (k<=j)
    {
      j=j-k;
      k=k/2;
    }
    j=j+k;
    if (i<j)
    {
      temp1.real=(Y[j]).real;
      temp1.imag=(Y[j]).imag;
      (Y[j]).real=(Y[i]).real;
      (Y[j]).imag=(Y[i]).imag;
      (Y[i]).real=temp1.real;
      (Y[i]).imag=temp1.imag;
    }
  }
  return;
}                                 //end of FFT()

     
main()
{

  for (i=0 ; i<PTS ; i++)         //set up twiddle factors in array w
  {
    w[i].real = cos(PI*(i)/PTS);
    w[i].imag = -sin(PI*(i)/PTS);
  }
printf("twiddle factors set up ...\n");
printf("\n");                     //insert breakpoint
  for (i=0 ; i<PTS ; i++)         //set up complex time domain filter coeffs
  {
    h[i].real = 0.0;
    h[i].imag = 0.0;
  }
  for (i=0 ; i<NUMCOEFFS ; i++)
  {
    h[i].real = coeffs[i];
  }
printf("time-domain filter coefficients loaded into h ...\n");
printf("\n");                     //insert breakpoint 
  FFT(h,PTS);                     //transform filter coeffs to freq domain 
printf("and transformed to frequency domain\n");
printf("\n");                     //insert breakpoint
  for (i=0 ; i<PTS/2 ; i++)
    {
    samples[i+PTS/2].real = 0.0;
    overlap[i] = 0.0; 
    }
 for (i=0 ; i<PTS/2 ; i++)        //initialise iobuffer contents
    {
    iobuffer[i] = (float)(cos(2*PI*i/110)) + 0.25*cos(2*PI*i/3); 
    }
printf("iobuffer filled\n");
printf("\n");                     //insert breakpoint
  for (i=0 ; i<PTS/2 ; i++)
    {
    samples[i].real = iobuffer[i];
    }
  for (i=0 ; i<PTS/2 ; i++)
    {
   overlap[i] = samples[i+PTS/2].real;
    samples[i+PTS/2].real = 0.0;
    }
  for (i=0 ; i<PTS ; i++)
    {
    samples[i].imag = 0.0;
    }
printf("iobuffer contents copied to samples ...\n");
printf("\n");                     //insert breakpoint
  FFT(samples,PTS);
printf("and converted to frequency domain\n");
printf("\n");                     //insert breakpoint
  for (i=0 ; i<PTS ; i++)
    {
    a = samples[i].real;
    b = samples[i].imag;
    samples[i].real = h[i].real*a - h[i].imag*b;
    samples[i].imag = h[i].real*b + h[i].imag*a;
    }
printf("samples multiplied by filter coefficients\n");
printf("\n");                     //insert breakpoint
    for (i=0 ; i<PTS ; i++)       //inverse transform samples back to time domain
    {
      samples[i].imag = -samples[i].imag;
    }
    FFT(samples,PTS);
    for (i=0 ; i<PTS ; i++)
    {
      samples[i].real = samples[i].real/PTS;
      samples[i].imag = -samples[i].imag/PTS;
    }
printf("samples transformed back to time domain ...\n");
printf("\n");                     //insert breakpoint
  for (i=0 ; i<PTS/2 ; i++)
    {
    overlap[i] += samples[i].real;
    }
printf("and summed with (overlapping) previous result into overlap\n");
printf("\n");                     //insert breakpoint
  for (i=0 ; i<PTS/2 ; i++)
    {
    iobuffer[i] = (float)(cos(2*PI*(i+(PTS/2))/110)) + 0.25*cos(2*PI*(i+(PTS/2))/3); 
    }
printf("iobuffer filled\n");
printf("\n");                     //insert breakpoint
  for (i=0 ; i<PTS/2 ; i++)
    {
    samples[i].real = iobuffer[i];
    }
  for (i=0 ; i<PTS/2 ; i++)
    {
    overlap[i] = samples[i+PTS/2].real;
    samples[i+PTS/2].real = 0.0;
    }
  for (i=0 ; i<PTS ; i++)
    {
    samples[i].imag = 0.0;
    }
printf("iobuffer contents copied to samples ...\n");
printf("\n");                     //insert breakpoint
  FFT(samples,PTS);
printf("and converted to frequency domain\n");
printf("\n");                     //insert breakpoint
  for (i=0 ; i<PTS ; i++)
    {
    a = samples[i].real;
    b = samples[i].imag;
    samples[i].real = h[i].real*a - h[i].imag*b;
    samples[i].imag = h[i].real*b + h[i].imag*a;
    }
printf("samples multiplied by filter coefficients\n");
printf("\n");                     //insert breakpoint
    for (i=0 ; i<PTS ; i++)       //inverse transform samples back to time domain
    {
      samples[i].imag = -samples[i].imag;
    }
    FFT(samples,PTS);
    for (i=0 ; i<PTS ; i++)
    {
      samples[i].real = samples[i].real/PTS;
      samples[i].imag = -samples[i].imag/PTS;
    }
printf("samples transformed back to time domain ...\n");
printf("\n");                     //insert breakpoint
  for (i=0 ; i<PTS/2 ; i++)
    {
    overlap[i] += samples[i].real;
    }
printf("and summed with (overlapping) previous result\n");
printf("\n end of program\n");    //insert breakpoint
}                                 //end of main()