/*
     fhti.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)    \
 {/*long*/ double temp;             \
   temp = Pr;                       \
   Pr = Pr * Nr - Pi   * Ni + Pr;   \
   Pi = Pi * Nr + temp * Ni + Pi;   \
 }

#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;                  \
 }

#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 DiF FHT */
#define FHT_F1Butterfly(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[i2]);D1=D1-temp; \
  }                                          \
  Right[i2]=D1*(C)+D2*(S);                   \
 }

#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_T(double *Data, int Len)
/* Decimation in Time Hartley transform */
{
  int Step, Step2, Step4;
  double Sin0,Sin,Cos0,Cos;
  int a,b;
  int TrigIndex=1;
  double *Right,*Left;

  TrigIndex=1;
  Step = 1;
  while (Step < Len)
    {
      Step *= 2;
      Step2 = Step/2;
      Step4 = Step/4;
      Left  = &Data[0];
      Right = &Data[Step2];

      Sin0 = SineTable[TrigIndex++];
      Cos0 = SineTable[TrigIndex];Cos0 = -2.0 * Cos0 * Cos0;
      Sin=Sin0;Cos=1.0+Cos0;

      /* Do the special b=0 loop below */
      for (a = 0; a < Len; a+=Step)
         {
          FHT_T1Butterfly(a,a,1.0,0.0);
          if (Step4)
            FHT_T1Butterfly(a+Step4,a+Step4,1.0,0.0);
         }
      for (b = 1; b < Step4; b++)
        {int I1,I2;
         for (I1=b,I2=Step2-b; I1 < Len; I1+=Step,I2+=Step)
            FHT_T2Butterfly(I1,I2,Cos,Sin);
         NextTrigPow(Cos,Sin,Cos0,Sin0);
        }
    }
}

void
FHT_F(double *Data, int Len)
/* Decimation in frequency Hartley transform */
{
  int Step, Step2, Step4;
 /* long */ double Sin0,Sin,Cos0,Cos;
  int a,b;
  int TrigIndex;
  double *Right,*Left;

  TrigIndex=Log2(Len);
  Step = Len;
  while (Step > 1)
    {
      Step2 = Step/2;
      Step4 = Step/4;
      Left  = &Data[0];
      Right = &Data[Step2];

      Sin0 = SineTable[TrigIndex];
      Cos0 = SineTable[TrigIndex+1];Cos0 = -2.0 * Cos0 * Cos0;
      Sin=Sin0;Cos=1.0+Cos0;
      TrigIndex--;

      /* Do the special b=0 loop below */
      for (a = 0; a < Len; a+=Step)
        {
         FHT_F1Butterfly(a,a,1.0,0.0);
         if (Step4)
           FHT_F1Butterfly(a+Step4,a+Step4,1.0,0.0);
        }

      for (b = 1; b < Step4; b++)
        {int I1,I2;
         for (I1=b,I2=Step2-b; I1 < Len; I1+=Step,I2+=Step)
            FHT_F2Butterfly(I1,I2,Cos,Sin);
         NextTrigPow(Cos,Sin,Cos0,Sin0);
        }
      Step /= 2;
    }
}

#endif