/* fft128.c  complex data stored  re, im, re, im, re, ... im */

void fft128(float Z[256]) /* input data points and output  [0] to [255] */
{
  static float W[256];    /* scratch vector, used many times, */
  static float E[130] =   /* constants for FFT algorithm */
  {1.0, 0.0, 0.998796, 0.0490677,
   0.995185, 0.0980171, 0.989177, 0.14673,
   0.980785, 0.19509, 0.970031, 0.24298,
   0.95694, 0.290285, 0.941544, 0.33689,
   0.92388, 0.382683, 0.903989, 0.427555,
   0.881921, 0.471397, 0.857729, 0.514103,
   0.83147, 0.55557, 0.803208, 0.595699,
   0.77301, 0.634393, 0.740951, 0.671559,
   0.707107, 0.707107, 0.671559, 0.740951,
   0.634393, 0.77301, 0.595699, 0.803208,
   0.55557, 0.83147, 0.514103, 0.857729,
   0.471397, 0.881921, 0.427555, 0.903989,
   0.382683, 0.92388, 0.33689, 0.941544,
   0.290285, 0.95694, 0.24298, 0.970031,
   0.19509, 0.980785, 0.14673, 0.989177,
   0.0980171, 0.995185, 0.0490677, 0.998796,
   0.0, 1.0, -0.0490677, 0.998796,
   -0.0980171, 0.995185, -0.14673, 0.989177,
   -0.19509, 0.980785, -0.24298, 0.970031,
   -0.290285, 0.95694, -0.33689, 0.941544,
   -0.382683, 0.92388, -0.427555, 0.903989,
   -0.471397, 0.881921, -0.514103, 0.857729,
   -0.55557, 0.83147, -0.595699, 0.803208,
   -0.634393, 0.77301, -0.671559, 0.740951,
   -0.707107, 0.707107, -0.740951, 0.671559,
   -0.77301, 0.634393, -0.803208, 0.595699,
   -0.83147, 0.55557, -0.857729, 0.514103,
   -0.881921, 0.471397, -0.903989, 0.427555,
   -0.92388, 0.382683, -0.941544, 0.33689,
   -0.95694, 0.290285, -0.970031, 0.24298,
   -0.980785, 0.19509, -0.989177, 0.14673,
   -0.995185, 0.0980171, -0.998796, 0.0490677,
   -1.0, 0.0};
  float Tre, Tim;
  int i, j, k, l, m;

   m = 64;
   l = 1;
   while(1)
   { 
      k = 0;
      j = l;
      i = 0;
      while(1)
      { 
         while(1)
	 {
	   /* W[i+k] = Z[i] + Z[m+i]; complex */
           W[2*(i+k)]   = Z[2*i]   + Z[2*(m+i)];
           W[2*(i+k)+1] = Z[2*i+1] + Z[2*(m+i)+1];

           /* W[i+j] = E[k] * (Z[i] - Z[m+i]); complex */
           Tre = Z[2*i]   - Z[2*(m+i)];
           Tim = Z[2*i+1] - Z[2*(m+i)+1];
           W[2*(i+j)]   = E[2*k] * Tre - E[2*k+1] * Tim;
           W[2*(i+j)+1] = E[2*k] * Tim + E[2*k+1] * Tre; 
           i++;
           if(i >= j) break;
	 }
         k = j;
         j = k+l;
         if(j > m) break;
      }
      l = l+l;
      if(l > m)
      {
        for(i=0; i<4*m; i++) Z[i] = W[i]; // odd power finish
        return;
      }
                  /* work back other way without copying */
      k = 0;
      j = l;
      i = 0;
      while(1)
      {
        while(1)
	{
          /* Z[i+k] = W[i] + W[m+i]; complex */
          Z[2*(i+k)]   = W[2*i]   + W[2*(m+i)];
          Z[2*(i+k)+1] = W[2*i+1] + W[2*(m+i)+1];

          /* Z[i+j] = E[k] * (W[i] - W[m+i]); complex */
          Tre = W[2*i]   - W[2*(m+i)];
          Tim = W[2*i+1] - W[2*(m+i)+1];
          Z[2*(i+j)]   = E[2*k] * Tre - E[2*k+1] * Tim;
          Z[2*(i+j)+1] = E[2*k] * Tim + E[2*k+1] * Tre;
          i++;
          if(i >= j) break;
	}
        k = j;
        j = k+l;
        if(j > m) break;
      }
      l = l+l;
      if(l > m) break; // result is in Z
   }
} /* end fft128 */