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

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

/* external functions used by this program */
extern void ditr2fft_fixpt(complex *, unsigned int, complex *, unsigned int);
extern void ibit_reversal(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 xn;

void main()
{
  unsigned int i,L,LE,LE1; 
  /*****************************************************************
  *    Declare file pointers
  *****************************************************************/

  FILE *xn_in;                       /* file pointer of x(n)      */
  FILE *yn_out;                      /* file pointer of y(n)      */
  xn_in = fopen("in_int.dat","r");   /* open file for input x(n)  */
  yn_out = fopen("out_int.dat","w"); /* open file for output y(n) */

  /* Step 1: Create a twiddle factor table                        */

  for (L=1; L<=EXP; L++)
  {
    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                 */

  for(i=0; i<N; i++)
  {	         /* Input signal                                  */
	fscanf(xn_in,"%d",&xn) ;
	{	
	  X[i].re = xn;
	  X[i].im = 0;
	}
  }

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

  ibit_reversal(X,EXP);      /* arrange X[] in bit-reverse order  */
  ditr2fft_fixpt(X,EXP,W,1); /* FFT with scaling 0.5 in each stage*/


  /* Step 4: Perform magnitude-square                             */

  for (i=0; i<N; i++) /* verify FFT result                        */
  {
    ltemp.re = (long)X[i].re*X[i].re;
    ltemp.im = (long)X[i].im*X[i].im;        
    spectrum[i] = (int)((ltemp.re+ltemp.im)>>10); /* scaling      */
	fprintf(yn_out,"%d\n",spectrum[i]);
  }
  fcloseall();
}