/*****************************************************************
*  psd_wel_floatpt.c - C program computes FFT-based PSD using
*                      Welch periodogram
******************************************************************
*  System configuration:
*
*          _____     ____     __________     __________     ____    
*         |     |   |    |   |          |   |          |   |    |
*  x(n)-->| Buf |-->| BR |-->| (1/N)*FFT|-->| |X(k)|^2 |-->| Av |-> Pi(k)
*         |_____|   |____|   |__________|   |__________|   |____|  
*                   
******************************************************************
*  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"
*
*     Number of segments have been precomputed
*****************************************************************/

#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 overlap 128     /* overlap = 50%                        */
#define K 2             /* K = int[(length(x)-N)/(N-overlap)+1] */ 

complex X[N];           /* declare input array                  */
float X_old[overlap];
float X_new[N-overlap];
complex W[EXP];         /* twiddle e^(-j2pi/N) table            */    
complex temp;
float spectrum[K][N];
float spectrum_av[N];
float win[N];
float xn;
float U=0.0;

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);
  }    

  /* Create Hamming window and init spectrum_av                   */
  for (i=0;i<N; i++)
  {
    win[i] = 0.54f-0.46f*cos(2*pi*i/N);
    spectrum_av[i] = 0.0;
  }

  /* Save old overlapped samples                                  */
  for (i=0; i<overlap; i++)
  {
    fscanf(xn_in,"%f",&xn);
    X_old[i] = xn;
  }
  for (i=0; i<(N-overlap); i++)
  {
    fscanf(xn_in,"%f",&xn);
    X_new[i] = xn;
  }

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

    if (j > 0)
    {
      for(i=0; i<(N-overlap); i++)
      {	/* new input signal */
	fscanf(xn_in,"%f",&xn) ;
	X_new[i] = xn;
      }
    }

    for(i=0;i<overlap; i++)
    {
      X[i].re = X_old[i]*win[i];
      X[i].im = 0.0;
    }
    for(i=overlap; i<N; i++)
    {
      X[i].re = X_new[i-overlap]*win[i];
      X[i].im = 0.0;
    }
	
    if (N-(2*overlap) >= 0) 
    {				/* for overlap =< 50% */
      for (i=0; i<overlap; i++)
      {
	X_old[i] = X_new[N-(2*overlap)+i];
      }
    }
    else 
    {				/* for overlap > 50% */
      for (i=0; i<((2*overlap)-N);i++)
      {
	X_old[i] = X_old[i+(N-overlap)];
      }
      for (i=0; i<(N-overlap); i++)
      {
	X_old[i+(N-overlap)] = X_new[i];
      }
    }

    /* 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++)		/* compute window normalizing factor         */
  { 
	U = win[i]*win[i] + U;
  }
  U = U/N;

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