/*
     fht.c
*/
/*
** This file is placed into the public domain by its author,
** Carey Bloodworth (Carey@Bloodworth.org) on July 16, 2001
**
** This multiplication demo is not designed for high performance.
** It's a tutorial program designed to be used with the information
** on my web site at www.Bloodworth.org
*/
/*
** This file demonstrates several versions of the Fast Hartley
** Transform (FHT).  The FHT algorithm is no longer patented and
** is in the public domain, like the FFT algorithm.
**
** To compile this using GCC:
** gcc main.c fht.c -o fht.exe
**
** You can optionally include a few options
**
** gcc -DPURE_FHT -DRECURSIVE main.c fht.c -o fht.exe
** gcc -DPURE_FHT -DITERATIVE main.c fht.c -o fht.exe
** gcc -DREAL_FHT -DRECURSIVE main.c fht.c -o fht.exe
** gcc -DREAL_FHT -DITERATIVE main.c fht.c -o fht.exe
*/
#include 
#include 
#include 
#include 
#include 

/*
#define REAL_FHT
#define PURE_FHT
** Make FHT look like a real value FFT?  If not, just treat it as
** a regular scrambled FHT.
*/
#if !defined(PURE_FHT) && !defined(REAL_FHT)
#define PURE_FHT
#endif

/*
#define ITERATIVE
#define RECURSIVE
** Do you want to use iterative or recursive fht's?
*/
#if !defined(ITERATIVE) && !defined(RECURSIVE)
#define ITERATIVE
#endif

typedef short int Short;

extern double MaxFFTError;

static double *FHTNum1=NULL, *FHTNum2=NULL;
static double BASE;
static int BASE_DIG;
static int IBASE;
static double SineTable[80];

#define CalcFHTLen(_NumLen) ((((_NumLen)*1)*2)/1)


#define MY_PI      3.1415926535897932384626434
#define MY_SQRT_2  0.7071067811865475244008443621
#define MY_SQRT2   1.41421356237309504880

static int Log2(int Num)
{int x=-1;
if (Num==0) return 0;
while (Num) {x++;Num/=2;}
return x;
}

#ifdef ITERATIVE
#include "fhti.h"
#endif
#ifdef RECURSIVE
#include "fhtr.h"
#endif


#ifdef REAL_FHT

void
FHT_Scramble(double *Data,int Len)
{int k1,k2,k;
 for (k1=1,k2=0;k1>1; (!((k2^=k)&k)); k>>=1) ;
     if (k1>k2)
       {double Temp;
         Temp=Data[k1];
         Data[k1]=Data[k2];
         Data[k2]=Temp;
       }
    }
}

void
FwdRealFFT(double *Data,int Len)
{double a,b;
 int i,j,k;
 FHT_Scramble(Data,Len);FHT_T(Data,Len);
/* FHT_F(Data,Len);FHT_Scramble(Data,Len); */
 for (i=1,j=Len-1,k=Len/2;i 0; x--)
    {
      RawPyramid = FHTNum1[Len2 - x] * inv + Carry;
      Pyramid = floor(RawPyramid + 0.5);
      PyramidError = fabs(RawPyramid - Pyramid);
      if (PyramidError > MaxFFTError)
        {
          MaxFFTError = PyramidError;
/*        printf("New Max FFT Error=%f\n",MaxFFTError);*/
        }
      Carry = floor(Pyramid / BASE);
      prod[x - 1] = Pyramid - Carry * BASE;
    }
}


void
InitFHT(unsigned int Len)
/*
** If you need to use your own trig tables, you can calculate
** them here.
*/
{int x;unsigned int P;
SineTable[0]=-9.0;
x=1;P=1;
while (P<=Len*4)
  {
   SineTable[x]=sin(MY_PI/P);
   P*=2;
   x++;
  }
}

void
InitFFT(unsigned long int Len,int Base,int BaseDig)
{int Bytes;

BASE=Base;
IBASE=Base;
BASE_DIG=BaseDig;

Bytes=sizeof(double)*CalcFHTLen(Len);

if (BaseDig > 4)
  {
   printf("Error:  The fft is slightly hardwired for <= 4 digits in the base.\n");
   exit(0);
  }

FHTNum1=(double*)malloc(Bytes);
FHTNum2=(double*)malloc(Bytes);

if ((FHTNum1==NULL) || (FHTNum2==NULL))
   {
    printf("Unable to allocate memory for FHTNum.\n");
    printf("Len=%d Bytes=%d\n",(int)Len,(int)Bytes);
    exit(0);
   }
InitFHT(Len);
}

void DeInitFFT(unsigned long int Len)
{
free(FHTNum1);free(FHTNum2);
}