/*
 *     Program: REALFFT.C
 *      Author: Philip VanBaren
 *        Date: 2 September 1993
 *
 * Description: These routines perform an FFT on real data.
 *              On a 486/33 compiled using Borland C++ 3.1 with full
 *              speed optimization and a small memory model, a 1024 point 
 *              FFT takes about 16ms.
 *              This code is for integer data, but could be converted
 *              to float or double simply by changing the data types
 *              and getting rid of the bit-shifting necessary to prevent
 *              overflow/underflow in fixed-point calculations.
 *
 *  Note: Output is BIT-REVERSED! so you must use the BitReversed to
 *        get legible output, (i.e. Real_i = buffer[ BitReversed[i] ]
 *                                  Imag_i = buffer[ BitReversed[i]+1 ] )
 *        Input is in normal order.
 */

#include <math.h>
#include <stdlib.h>
#include <stdio.h>

int *BitReversed;
int *SinTable;
int Points = 0;

/*
 *  Initialize the Sine table and Twiddle pointers (bit-reversed pointers)
 *  for the FFT routine.
 */
void InitializeFFT(int fftlen)
{
   int i;
   int temp;
   int mask;

   /*
    *  FFT size is only half the number of data points
    *  The full FFT output can be reconstructed from this FFT's output.
    *  (This optimization can be made since the data is real.)
    */
   Points = fftlen/2;

   if((SinTable=(int *)malloc(2*Points*sizeof(int)))==NULL)
   {
      puts("Error allocating memory for Sine table.");
      exit(1);
   }
   if((BitReversed=(int *)malloc(Points*sizeof(int)))==NULL)
   {
      puts("Error allocating memory for BitReversed.");
      exit(1);
   }

   for(i=0;i<Points;i++)
   {
      temp=0;
      for(mask=Points/2;mask>0;mask >>= 1)
         temp=(temp >> 1) + (i&mask ? Points : 0);

      BitReversed[i]=temp;
   }

   for(i=0;i<Points;i++)
   {
      SinTable[BitReversed[i]  ]=floor(-32768*sin(2*M_PI*i/(2*Points))+0.5);
      SinTable[BitReversed[i]+1]=floor(-32768*cos(2*M_PI*i/(2*Points))+0.5);
   }
}

/*
 *  Free up the memory allotted for Sin table and Twiddle Pointers
 */
void EndFFT()
{
   free(BitReversed);
   free(SinTable);
   Points=0;
}

int *A,*B;
int *sptr;
int *endptr1,*endptr2;
int *br1,*br2;
long HRplus,HRminus,HIplus,HIminus;

/*
 *  Actual FFT routine.  Must call InitializeFFT(fftlen) first!
 */
void realfft(int *buffer)
{
   int ButterfliesPerGroup=Points/2;

   endptr1=buffer+Points*2;

   /*
    *  Butterfly:
    *     Ain-----Aout
    *         \ /
    *         / \
    *     Bin-----Bout
    */

   while(ButterfliesPerGroup>0)
   {
      A=buffer;
      B=buffer+ButterfliesPerGroup*2;
      sptr=SinTable;

      while(A<endptr1)
      {
         register int sin=*sptr;
         register int cos=*(sptr+1);
         endptr2=B;
         while(A<endptr2)
         {
            long v1=((long)*B*cos + (long)*(B+1)*sin) >> 15;
            long v2=((long)*B*sin - (long)*(B+1)*cos) >> 15;
            *(A++)=(*(B++)=(*A+v1)>>1)-v1;
            *(A++)=(*(B++)=(*A-v2)>>1)+v2;
         }
         A=B;
         B+=ButterfliesPerGroup*2;
         sptr+=2;
      }
      ButterfliesPerGroup >>= 1;
   }
   /*
    *  Massage output to get the output for a real input sequence.
    */
   br1=BitReversed+1;
   br2=BitReversed+Points-1;

   while(br1<br2)
   {
      register long temp1;
      register long temp2;
      int sin=SinTable[*br1];
      int cos=SinTable[*br1+1];
      A=buffer+*br1;
      B=buffer+*br2;
      HRplus = (HRminus = *A     - *B    ) + (*B << 1);
      HIplus = (HIminus = *(A+1) - *(B+1)) + (*(B+1) << 1);
      temp1  = ((long)sin*HRminus - (long)cos*HIplus) >> 15;
      temp2  = ((long)cos*HRminus + (long)sin*HIplus) >> 15;
      *B     = (*A     = (HRplus  + temp1) >> 1) - temp1;
      *(B+1) = (*(A+1) = (HIminus + temp2) >> 1) - HIminus;

      br1++;
      br2--;
   }
   /*
    *  Handle DC bin separately
    */
   buffer[0]+=buffer[1];
   buffer[1]=0;
}