/*****************************************************************
*  psd_bar_floatpt.c - C program for FFT-based PSD using
*                      Bartlett periodogram
******************************************************************
*  System configuration:
*
*          _____     ____     __________     __________     ____    
*         |     |   |    |   |          |   |          |   |    |
*  x(n)-->| Buf |-->| BR |-->| (1/N)*FFT|-->| |X(k)|^2 |-->| Av |-> Pi(k)
*         |_____|   |____|   |__________|   |__________|   |____|  
*                   
*           
*  (b) Bartlett Periodigram
******************************************************************
*  System simulation configuration:
*
*     x(n) is the input data from data file "in.dat"
*     (assuming that the length of x(n) has been zero padded to 
*      the right length for averaged periodogram)
*     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.h"    /* complex.h header file */

/* External functions                                */

extern void ditr2fft(complex *, unsigned int, complex *, unsigned int);
extern void bit_reversal(complex *, unsigned int);

/* variables and constants                           */

#define N 256           /* number of FFT points      */
#define EXP 8           /* EXP=log2(N)               */
#define pi 3.1415926535897  
#define K 2             /* K = int(length(x)/N)      */ 
                       
complex X[N];           /* declare input array       */
complex W[EXP];         /* twiddle e^(-j2pi/N) table */    
complex temp;
float spectrum[K][N];
float spectrum_av[N];
float xn;

void main()
{
  unsigned int i,j,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.dat","r");    /* open file for input x(n)   */
  yn_out = fopen("out.dat","w");  /* open file for output y(n)  */


  /* 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 = cos(pi/LE1);
    W[L-1].im = -sin(pi/LE1);
  }    
  for(i=0; i<N; i++)
  {
    spectrum_av[i] = 0.0;
  }

  for (j=0; j<K; j++)
  {
    /* Step 2: Enter input data to reference buffer */

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

    /* Step 3: Perform bit reversal, follow by FFT */
    /* Start FFT */
    bit_reversal(X,EXP);        /* arrange X[] in bit-reverse order              */
    ditr2fft(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                             */
    {
      temp.re = X[i].re*X[i].re;
      temp.im = X[i].im*X[i].im;        
      spectrum[j][i] = (temp.re+temp.im);  /* scaling is perform in fft          */	
    }
  }

  /* Step 5: Perform averaging... */
  for (j=0; j<K; j++)
  {
    for (i=0; i<N; i++)
    {
      spectrum_av[i] = spectrum[j][i]+spectrum_av[i];
    }
  }

  for (i=0;i<N; i++)
  {
    spectrum_av[i]=spectrum_av[i]/K;
    fprintf(yn_out,"%f\n",spectrum_av[i]);
  }
  fcloseall();
}