Быстрое преобразование Хартли

//=====================================
/*Program
FHT.C*********************************************************************
/*
/* This program assumes a maximum array length of 2^M = N where */
/* M=9 and N=512. */
/* See Line 52 if the array length is increased. */
# include <stdio.h>
# include <math.h>
#include <conio.h>
# define M 3
# define N 8
float* myFht ( );
float* myFht()
    float X;
    int n, m;
    {
      int i, j, k, kk, _1, _10, _11, _12, _13, _14, _15, m1, n1, n2, NN, s;
      int diff=0, diff2, gamma, gamma2=2, n2_2, n2_4, n_2, n_4, n_8, n_16;
      int itemp, ntemp, phi, theta_by_2;
      float ee, temp1, temp2, xtemp1, xtemp2;
      float h_sec_b, x0, x1, x2, x3, x4, x5, xtemp;
      double cc1, cc2, ss1, ss2;
      double sine[256];
/***********************************************************************/
/* Digit reverse counter. */
/***********************************************************************/
      int powers_of_2[16], seed[256];
      int firstj, log2_n, log2_seed_size;
      int group_no, nn, offset;
      log2_n=m >> 1;
      nn=2<<(log2_n - 1);
      if ( (m % 2) ==1 )
      log2_n = log2_n + 1;
      seed[0] = 0; seed[1] = 1;
      for (log2_seed_size = 2; log2_seed_size <= log2_n; log2_seed_size++)
        {
          for ( i = 0; i <2 << (log2_seed_size - 2); i++)
            {
              seed[i] = 2*seed[i];
                for (k = 1; k < 2; k++)
              seed[ i + k * (2 << (log2_seed_size - 1)>>1) ] = seed[i];
            }
        }
      for (offset = 1; offset < nn; offset++)
        {
           {firstj = nn*seed[offset];
           i = offset; j=firstj;
           xtemp = x[i];
           x[i] = x[j];
           x[j] = xtemp;
              for ( group_no = 1; group_no < seed[offset]; group_no++)
                {
                  i = i + nn; j = firstj + seed[group_no];
                  xtemp = x[i];
                  x[i] = x[j];
                  x[j] = xtemp;
                }
        }
      j = 0;
      n1 = n - 1;
      n_16 = n >> 4;
      n_8 = n >> 3;
      n_4 = n >> 2;
      n_2 = n >> 1;
/************************************************************************/
/* Start the transform computation with 2-point butterflies. */
/************************************************************************/
      for (i = 0; i < n; i += 2)
        {
          s = i+1;
          xtemp = x[i];
          x[i] += x[s];
          x[s] = xtemp - x[s];
        }
/************************************************************************/
/* Now, the 4-point butterflies. */
/************************************************************************/
      for ( i = 0; i < n; i += 4)
        {
          xtemp = x[i];
          x[i] += x[i+2];
          x[i+2] = xtemp - x[i+2];
          xtemp = x[i+1];
          x[i+1] += x[i+3];
          x[i+3] = xtemp - x[i+3];
        }
/*******************************************************************/
/* Sine table initialization. */
/*******************************************************************/
      NN = n_4;
      sine[0] = 0;
      sine[n_16] = 0.382683432;
      sine[n_8] = 0.707106781;
      sine[3*n_16] = 0.923879533;
      sine[n_4] = 1.000000000;
      h_sec_b = 0.509795579;
      diff = n_16;
      theta_by_2 = n_4 >> 3;
      j = 0;
      while (theta_by_2 >= 1)
        {
          for ( i = 0; i <= n_4; i += diff)
            {
              sine[j + theta_by_2] = h_sec_b * (sine[j] + sine[j + diff] );
              j = j + diff;
            }
          j = 0;
          diff = diff >> 1;
          theta_by_2 = theta_by_2 >> 1;
          h_sec_b = 1 / sqrt(2 + 1/h_sec_b);
        }
/***********************************************************************/
/* Other butterflies. */
/***********************************************************************/
      for ( i = 3; i <= m; i++ )
        {
          diff = 1; gamma = 0;
          ntemp = 0; phi = 2 << (m-i) >> 1;
          ss1 = sine[phi];
          cc1 = sine[n_4 - phi];
          n2 = 2 << (i-1);
          n2_2 = n2 >>1;
          n2_4 = n2 >> 2;
          gamma2 = n2_4;
          diff2 = gamma2 + gamma2 - 1;
          itemp = n2_4;
          k=0;
/************************************************************************/
/* Initial section of stages 3, 4, …for which sines & cosines are */
/* not required. */
/************************************************************************/
             for (k = 0; k < (2 << (m-i)>>1); k++)
               {
                 _10 = gamma;
                 _11 = _10 + n2_2;
                 _13 = gamma2;
                 _14 = gamma2 + n2_2;
                 _15 = _11 + itemp;
                 x0 = x[_10];
                 x1 = x[_11];
                 x3 = x[_13];
                 x5 = x[_15];
                 x1 = x0 + x1;
                 x[_11] = x0 - x1;
                 x[_13] = x3 + x5;
                 x[_14] = x3 - x5;
                 gamma = gamma + n2;
                 gamma2 = gamma2 + n2;
               }
          gamma = diff;
          gamma2 = diff2;
/************************************************************************/
/* Next sections of stages 3, 4, … */
/************************************************************************/
               for ( j = 1; j < 2 << (i-3); j++ )
                 {
                   for ( k = 0; k < (2 << (m-i) >> 1); k++)
                     {
                       _10 = gamma;
                       _11 = _10 + n2_2;
                       _13 = gamma2;
                       _14 = _13 + n2_2;
                       x0 = x[_10];
                       x1 = x[_11];
                       x3 = x[_13];
                       x4 = x[_14];
                       x[_10] = x0 + x1 * cc1 + x4 * ss1;
                       x[_11] = x0 - x1 * cc1 - x4 * ss1;
                       x[_13] = x3 - x4 * cc1 + x1 * ss1;
                       x[_14] = x3 + x4 * cc1 - x1 * ss1;
                       gamma = gamma + n2;
                       gamma2 = gamma2 + n2;
                      }
                   itemp = 0;
                   phi = phi + ( 2 << (m-i) >> 1 );
                   ntemp=(phi < n_4) ? 0 : n_4;
                   ss1=sine[phi - ntemp];
                   cc1=sine[n_4 - phi + ntemp];
                   diff++; diff2--;
                   gamma=diff;
                   gamma2=diff2;
                 }
        }
     }
//
//--------- Main -------------
main(X,N,M)
{
/* Read the integer values 1, …, N into the vector X[N]. */
int i;
float X[N];
for (i = 0; i < N; i++)
  {
    X[i] = i+1;
      for (i = 0; i < N; i++)
      {
        printf ("%f\n", X[i]);
        myFht(X, N, M);
        printf ("\n");
          for (i = 0; i < N; i++)
            {
               printf ("%d: %f\n", i, X[i]/N);
            }
        }
    }
              /* It is assumed that the user divides by the integer N. */
//-------- Stop Program !!!! -------------
    getch();
    printf("For Stop Program press qey [q]: \n")
    return 0;
}
//==========================================