#include <stdio.h>
#include <string.h>
#include <math.h>
#include <errno.h>
#include <stdlib.h>
#include <time.h>

#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp
#define LIW 15
#define MYDEBUG 1
#define CSWAP(a,b) ctemp=(a);(a)=(b);(b)=ctemp
#define Switches (volatile char *) 0x01003020 
#define LEDs (char *) 0x01003030

#pragma off
/*-------------------------------------------------------------------------*/
/*  C obfuncation test V1.0 using FFT_NIOS.c from the book DSP with FPGAs  */
/*  CUSTOM INSTRUCTION in software for NIOS system without BFPROC          */
/*  This file is copyrighted 2011 by Dr. Uwe Meyer-Baese                   */
/*-------------------------------------------------------------------------*/
#pragma on

typedef struct { long int r, i; } icomplex;
char            str[50];
int             LOOP, S,N,T=0;
int           Cy[15], LMAX, Ly[15];
long int xrIN1, xrIN2, xiIN1, xiIN2, c, s;
time_t          t1, t2;

typedef unsigned long DWORD;
long nm_bfp3_pfx(pfx, a, b)
{
  long tr, ti;
  long xrOUT1, xiOUT1, xrOUT2, xiOUT2;

  if (pfx==1) {xrIN1=a; xiIN1=b;};
  if (pfx==2) {xrIN2=a; xiIN2=b;};
  if (pfx==3) {c=a; s=b;};


        tr = xrIN1 - xrIN2;
        xrOUT1 = xrIN2 + xrIN1;
        ti = xiIN1 - xiIN2;
        xiOUT1 = xiIN2 + xiIN1;
        xrOUT2 = (c * tr - s * ti) >> LIW;
        xiOUT2 = (c * ti + s * tr) >> LIW;

        if (pfx==4) return xrOUT1; else
        if (pfx==5) return xiOUT1; else
        if (pfx==6) return xrOUT2; else
        if (pfx==7) return xiOUT2; else
                return 0;

}


int main( void )
{
  int j, n, t;

  float f, PI = 3.14159265;

  int i1, i2, w, dw, k1, k2, l, k, d, e;
  long int  Wr, Wi, tr, ti;
  icomplex *coef, *x, *y, **z;
  icomplex ctemp;
  int start_time, finish_time, total_time;

  *LEDs= 255;
   printf("FFT beginning\n");
 
  for (S=3; S<=12; S++)  /* Stage Loop */
  {
  *LEDs= S;
    N = 1<<S;
    /* Allocate Menory  Note: DFT need N all other N/2 coeffs */
    coef = (icomplex *) malloc(sizeof(icomplex) * N/2);
        if (coef == NULL) {
      printf(" coef - Memory Allocation Failure\n");
      exit(1)     ; }
        x = (icomplex *) malloc(sizeof(icomplex) * N);
    if (x == NULL) {
      printf(" x - Memory Allocation Failure\n");
      exit(1)     ; }
        y = (icomplex *) malloc(sizeof(icomplex) * N);
    if (y == NULL) {
      printf(" y - Memory Allocation Failure\n");
      exit(1)     ; }


      /* Calculate Twiddle Factor Length */

    f = (2 * PI)/N;  /* Calculation of wnr and wni values */

    for (n=0; n<(N/2); n++)
    {
            coef[n].r = ( long int)  ((1<<LIW) * cos(f * n));
            coef[n].i = (long int )  ((1<<LIW) * sin(f * n));

    }

    if (MYDEBUG == 2) {
      printf("**** UMB Fast SIN/COS coeffs for length %d\n",N);
      for (n=0; n<(N/2); n++)
        printf("(%ld)+i(%ld)\n", coef[n].r,coef[n].i);
    }

    /* Beginning 3. FFT Algorithm using Icomplex (compact) data type    */
    for (n=0; n<8; n++) x[n].r=(n+1)*20;    /* Specify Real Input Values  */
 for (n=8; n<N; n++)  x[n].r = 0;
    for (n=0; n<N; n++)  x[n].i = 0;

    /** Display Radix-2 FFT Output **/
    if (MYDEBUG == 2) {
      printf("**** First 4 RADIX-2 FFT INPUT for length %d\n",N);
      for (n=0; n<4; n++)
        printf("(%ld)+i(%ld)\n", x[n].r,x[n].i);
    }

    LMAX=(int) 1.0e6 / ( N*log(N) );
    printf("**** ICOMPLEX FFT Start: LENGTH = %d LOOPS= %d\n",N,(int) LMAX);
    /* Record Start Time */
    start_time = alt_nticks(); 
    for (LOOP=1;LOOP<LMAX;LOOP++){
    k2 = N; dw = 1;
    for (l = 1; l <= S; l++) {  /* Loop Stages */
      k1 = k2; k2 >>= 1;
      w = 0;                    /* Angle count */
      for (k = 0; k < k2; k++) { /* Group loop */
        Wr = coef[w].r; Wi = coef[w].i;         /* Load precomputed twiddle factors */
        nm_bfp3_pfx(3,Wr,Wi);            /* Read COS+SIN      */
        w += dw;
        for (i1 = k; i1 < N; i1 += k1) { /* Butterfly loop */
          i2 = i1 + k2;       /* Butterfly calculation */
          /* Custom Instructions */
          tr = x[i1].r; ti = x[i1].i;
          nm_bfp3_pfx(1,tr,ti);            /* Read Are and Aim          */
          tr = x[i2].r; ti = x[i2].i;
          nm_bfp3_pfx(2,tr,ti);            /* Read Bre and Bim                    */
          x[i1].r = nm_bfp3_pfx(4,0,0);   /* Write the real value */
          x[i1].i = nm_bfp3_pfx(5,0,0);   /* Write the imaginary value */
          x[i2].r = nm_bfp3_pfx(6,0,0);   /* Write the real value */
          x[i2].i = nm_bfp3_pfx(7,0,0);   /* Write the imaginary value */
        }
      }
      dw <<= 1;
    }
    if (LOOP==1) for (j=0;j<N;j++) { y[j].r=x[j].r; y[j].i=x[j].i;}
    }

 /* Record End Time */
  finish_time = alt_nticks();
  total_time = ((finish_time - start_time) * 1000) / alt_ticks_per_second();
  printf("Total time: %d ms\n", total_time);
  printf("Start  ticks: %d \n", start_time);
  printf("Finish ticks: %d \n", finish_time);
    
  Cy[S]= total_time;
  Ly[S]= LMAX;
  
    T=4;
    printf("ICOMPLEX N=%d",N);

    /* Final Bit Reversal Process */
    d = N / 2;
    for (n=1; n < (N-1); n++)
    {
      if (n < d)
      {
        CSWAP(y[d], y[n]);      }

      e = N / 2;

      while (e <= d)
      {
        d = d - e;
        e = e / 2;
      }

      d = d + e;

    }



    /* Display Radix-2 FFT Output */
    if (MYDEBUG >= 1) {
      printf("**** UMB Fast ICOMPLEX first 4 RADIX-2 FFT OUTPUT for length %d\n",N);
      for (n=0; n<4; n++)
        printf("(%ld)+i(%ld)\n", y[n].r,y[n].i);
    }

    /*Free Up the Dynamic Allocation Space     */
    free (coef);
        free (x);
        free (y);

  }
  printf("\nFFT cycles using ICOMPLEX:\n");
    for (S=3;S<=12;S++) {
      N=1<<S;t=Cy[S];f=5.0 * N*log(N)/log(2.0) * Ly[S] /t;
      printf("%4d) T=%d ms L=%d 5N*log2(N)/T (us)=%d\n",N,t,Ly[S],(int) f );
  }
  printf("Ticks pers sec: %u\n", alt_ticks_per_second());

#pragma off
/* No Warranty! What so ever ! copyright 2011 Uwe Meyer-Baese */
#pragma on 

  return 0;
}