#include <inttypes.h>
#include <avr/io.h>
#include <avr/interrupt.h>
#include <stdio.h>
#include "uart.h"
#include <math.h>

// UART file descriptor
// putchar and getchar are in uart.c
FILE uart_str = FDEV_SETUP_STREAM(uart_putchar, uart_getchar, _FDEV_SETUP_RW);
       

//I like these definitions
#define begin {
#define end   } 

#define int2fix(a)   (((int)(a))<<8) 
#define float2fix(a) ((int)((a)*256.0)) 
#define fix2float(a) ((float)(a)/256.0)  

#define N_WAVE          16    /* size of FFT */
#define LOG2_N_WAVE     4     /* log2(N_WAVE) */

int Sinewave[N_WAVE]; // a table of sines for the FFT

//======================================================== 

void FFTfix(int fr[], int fi[], int m)
//Adapted from code by:
//Tom Roberts 11/8/89 and Malcolm Slaney 12/15/94 malcolm@interval.com
//fr[n],fi[n] are real,imaginary arrays, INPUT AND RESULT.
//size of data = 2**m 
// This routine does foward transform only
begin 
    int mr,nn,i,j,L,k,istep, n;
    int qr,qi,tr,ti,wr,wi; 
    
    mr = 0;
    n = 1<<m;   //number of points
    nn = n - 1;

    /* decimation in time - re-order data */
    for(m=1; m<=nn; ++m) 
    begin
        L = n;
        do L >>= 1; while(mr+L > nn);
        mr = (mr & (L-1)) + L;
        if(mr <= m) continue;
        tr = fr[m];
        fr[m] = fr[mr];
        fr[mr] = tr;
        //ti = fi[m];   //for real inputs, don't need this
        //fi[m] = fi[mr];
        //fi[mr] = ti;
    end
    
    L = 1;
    k = LOG2_N_WAVE-1;
    while(L < n) 
    begin
        istep = L << 1;
        for(m=0; m<L; ++m) 
        begin
            j = m << k;
            wr =  Sinewave[j+N_WAVE/4];
            wi = -Sinewave[j];
            wr >>= 1;
            wi >>= 1;
            
            for(i=m; i<n; i+=istep) 
            begin
                j = i + L;
                tr = multfix(wr,fr[j]) - multfix(wi,fi[j]);
                ti = multfix(wr,fi[j]) + multfix(wi,fr[j]);
                qr = fr[i] >> 1;
                qi = fi[i] >> 1;     
                fr[j] = qr - tr;
                fi[j] = qi - ti;
                fr[i] = qr + tr;
                fi[i] = qi + ti;
            end
        end
        --k;
        L = istep;
    end
end

//========================================================         

unsigned char i; 

int fr[N_WAVE], fi[N_WAVE]; // sin[N_WAVE], cos[N_WAVE];

int main(void)
begin  
  //init the UART -- uart_init() is in uart.c
  uart_init();
  stdout = stdin = stderr = &uart_str;
  fprintf(stdout,"Starting...\n\r");  
  
  // one cycle sine table
  //  required for FFT
  for (i=0; i<N_WAVE; i++)
    Sinewave[i] = float2fix(sin(6.283*((float)i)/N_WAVE)); 
    
  while(1)
  begin
    
	// compute an input signal for FFT
    for (i=0; i<N_WAVE; i++)
    begin
        fr[i]=float2fix(sin(6.283*4*((float)i)/N_WAVE)); fi[i]=0;
    end
    
    FFTfix(fr, fi, LOG2_N_WAVE);
    
    for (i=0; i<N_WAVE; i++)
    begin
        fprintf(stdout,"%2.6f\t%2.6f\t%2.6f\r\n",
        fix2float(Sinewave[i]),  fix2float(fr[i]), fix2float(fi[i]));
    end
   
    
    //============================================
	// time the FFT function
    TCCR1B = 1 ;
    TCNT1 = 0;
    FFTfix(fr, fi, LOG2_N_WAVE);
    TCCR1B = 0;
    fprintf(stdout,"FFT cycles=%d\n\r",TCNT1)  ;
    
    //============================================
  while(1);   
  end
end