/*
     fhtr.h
*/
/*
** This file is placed into the public domain by its author,
** Carey Bloodworth (Carey@Bloodworth.org) on July 16, 2001
*/
#ifndef FHT_H_
#define FHT_H_ 1

#define NextTrigPow(Pr,Pi,Nr,Ni)    \
 {double temp;                      \
   temp = Pr;                       \
   Pr = Pr * Nr - Pi   * Ni + Pr;   \
   Pi = Pi * Nr + temp * Ni + Pi;   \
 }

/* Macro for the DiF FHT */
#define FHT_F2Butterfly(N1,N2,C,S)           \
 {double D1,D2;                              \
  int i1=N1, i2=N2;                          \
  D1=Left[i1];D2=Left[i2];                   \
  {double temp;                              \
   Left[i1] =D1+(temp=Right[i1]);D1=D1-temp; \
   Left[i2] =D2+(temp=Right[i2]);D2=D2-temp; \
  }                                          \
  Right[i1]=D1*(C)+D2*(S);                   \
  Right[i2]=D1*(S)-D2*(C);                   \
 }

void FHT_F(double *Data,int Len)
/*
** Recursive Decimation in Frequency style Fast Hartley Transform
*/
{int x,Len2,Len4;
 double Sin0,Cos0;
 double Sin,Cos;
 double *Left,*Right;

Len/=2;Left=&Data[0];Right=&Data[Len];
if (Len==2)
  {double d0=Data[0]; double d1=Data[1];
   double d2=Data[2]; double d3=Data[3];
   {double d02=d0+d2; double d13=d1+d3;
    Data[0]=d02+d13; Data[1]=d02-d13;
   }
   {double d02=d0-d2; double d13=d1-d3;
    Data[2]=d02+d13; Data[3]=d02-d13;
   }
   return;
  }

{double t1,t2;
 t1=Left[0];t2=Right[0];
 Left[0]=t1+t2;Right[0]=t1-t2;
 t1=Left[Len/2];t2=Right[Len/2];
 Left[Len/2]=t1+t2;Right[Len/2]=t1-t2;
}

x=Log2(Len)+1;
Sin0=SineTable[x];
Cos0=SineTable[x+1];Cos0=-2.0*Cos0*Cos0;
Sin=Sin0;Cos=1.0+Cos0;

Len2=Len/2;
Len4=Len/4;
for (x=1;x=2) FHT_F(Left, Len);
if (Len>=2) FHT_F(Right,Len);
}


/* Macro for the DiT FHT */
#define FHT_T2Butterfly(N1,N2,C,S) \
 {double Rx,Ri;                  \
  int i1=N1,i2=N2;               \
  Rx=Right[i1];Ri=Right[i2];     \
  {double cas1,Lx;               \
   cas1=Rx*(C)+Ri*(S);           \
   Lx=Left[i1];                  \
   Left[i1]  = Lx+cas1;          \
   Right[i1] = Lx-cas1;          \
  }                              \
  {double cas2,Li;               \
   cas2=Rx*(S)-Ri*(C);           \
   Li=Left[i2];                  \
   Left[i2]  = Li+cas2;          \
   Right[i2] = Li-cas2;          \
  }                              \
 }

/* Macro for the DiT FHT */
#define FHT_T1Butterfly(N1,N2,C,S)         \
 {int i1=N1,i2=N2;                         \
  double cas1=Right[i1]*(C)+Right[i2]*(S); \
  double temp=Left[i1];                    \
  Left[i1] = temp + cas1;                  \
  Right[i2]= temp - cas1;                  \
 }

void FHT_T(double *Data,int Len)
/*
** recursive Decimation in Time style Fast Hartley Transform
*/
{int x,Len2,Len4;
 double Sin0,Cos0;
 double Sin,Cos;
 double *Left,*Right;

Len/=2;Right=&Data[Len];Left=&Data[0];
if (Len==4)
  {double d45,d67,sd0123,dd0123;
   {double ss0123,ds0123,ss4567,ds4567;
    {double s01,s23,d01,d23;
     d01 = Data[0] - Data[1];
     s01 = Data[0] + Data[1];
     d23 = Data[2] - Data[3];
     s23 = Data[2] + Data[3];
     ds0123 = (s01 - s23);
     ss0123 = (s01 + s23);
     dd0123 = (d01 - d23);
     sd0123 = (d01 + d23);
    }
    {double s45,s67;
     s45 = Data[4] + Data[5];
     s67 = Data[6] + Data[7];
     d45 = Data[4] - Data[5];
     d67 = Data[6] - Data[7];
     ds4567 = (s45 - s67);
     ss4567 = (s45 + s67);
    }
    Data[4] = ss0123 - ss4567;
    Data[0] = ss0123 + ss4567;
    Data[6] = ds0123 - ds4567;
    Data[2] = ds0123 + ds4567;
   }
   d45 *= MY_SQRT2;
   d67 *= MY_SQRT2;
   Data[5] = sd0123 - d45;
   Data[1] = sd0123 + d45;
   Data[7] = dd0123 - d67;
   Data[3] = dd0123 + d67;
   return;
  }

FHT_T(&Left[0], Len);
FHT_T(&Right[0],Len);

/* Do the special x=0 loop below. */
FHT_T1Butterfly(0,0,1.0,0.0);

x=Log2(Len)+1;
Sin0=SineTable[x];
Cos0=SineTable[x+1];Cos0=-2.0*Cos0*Cos0;
Sin=Sin0;Cos=1.0+Cos0;

Len2=Len/2;
Len4=Len/4;
for (x=1;x