/*****************************************************************
*  psd_fixpt_ccs.c - C program for computing FFT-based PSD		
******************************************************************
*  System configuration:
*
*          _____     ____     __________     __________
*         |     |   |    |   |          |   |          |
*  x(n)-->| Buf |-->| BR |-->| (1/N)*FFT|-->| |X(k)|^2 |--> P(k)
*         |_____|   |____|   |__________|   |__________|                         
*           
*  (a) Periodigram
******************************************************************
*  System simulation configuration:
*
*     x(n) is the input data from data file "in_int.dat"
*     out(n) is the power spectrum data to data file "out_int.dat"
*
*****************************************************************/

/* header files  */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>   
#include "def_complex_fixpt.h"	/* complex.h header file */

int in_buffer[512];
#define TRUE 1 

/*Function */
static void dataIO(void);

/* external functions  */
extern void ditr2fft_fixpt_ccs(complex *, unsigned int, complex *, unsigned int);
extern void ibit_reversal_ccs(complex *, unsigned int);


/*variables and constants */
#define N 512           /* Number of FFT points */
#define EXP 9           /* EXP=log2(N)          */
#define pi 3.1415926535897                         
complex X[N];           /* Declare input array  */
complex W[EXP];         /* Twiddle e^(-j2pi/N) table */    
lcomplex ltemp;
int spectrum[N];
int spectrum_db[N];


void main()
{
  unsigned int i,L,LE,LE1; 
  int *input = &in_buffer[0]; 
	
    /* Step 1: Create a twiddle factor table */

    for (L=1; L<=EXP; L++) /* Create twiddle factor table */
    {
        LE=1<<L;           /* LE=2^L=points of sub DFT */
        LE1=LE>>1;     	   /* Number of butterflies in sub-DFT */
	W[L-1].re = (int)((0x7fff*cos(pi/LE1))+0.5);
        W[L-1].im = -(int)((0x7fff*sin(pi/LE1))+0.5);
    }    
    
    /* Step 2: Enter input data to reference buffer */
    
    dataIO();
    for(i=0; i<N; i++)
    { 	 
	X[i].re = *input++;
	X[i].im = 0;	
    } 

	

    /* Step 3: Perform bit reversal, follow by FFT */

    /* Start FFT */
    ibit_reversal_ccs(X,EXP);        /* Arrange X[] in bit-reverse order */
    ditr2fft_fixpt_ccs(X,EXP,W,1);   /* Perform FFT with scaling of 0.5 in each stage */

    /* Step 4: Perform Magnitude-Square and display results in MATLAB */

    for (i=0; i<N; i++) /* Verify FFT result */
    {
 	ltemp.re = ((long)X[i].re*(long)X[i].re); 
  	ltemp.im = ((long)X[i].im*(long)X[i].im);   
 	spectrum[i] = (int)((ltemp.re+ltemp.im)>>9);/*Scaling is perform in fft*/       
        spectrum_db[i] = 10*log10(spectrum[i]+1);
    }
} 

/* Function for dataIO */

static void dataIO()
{
    return;
}