Well, according to the amount of E-mail I received, and as I cannot
answer to all of them, I am posting the 2D FFT program I have written in C.

	How to use it : first, it can be used only on matrix of the type
complexe, as defined at the begining of the program, with a size in
2^maxlevel (2 power maxlevel). One just have to call the program as following :

	InvFFT(matrix,maxlevel);

	As you have noticed, it is in fact an Inverse Fourier Transform
(sampling matrix in the frequency space -> matrix of heights)  which
is performed, but just an initialisation line at the begining of the
procedure  InvFFT has to be change to obtain the FFT :

	Omega = 2 * M_PI / dimension;      must be change into
	Omega = - 2 * M_PI / dimension;


	I am clearly aware that this program has not been optimised as it
could be, mainly in order to keep the the clarity of the algorithm, but I am
open to any suggestion that you could make in order to improve the process.

	I have tested it in a fractal landscape project and it seems to work
correctly (when there were bugs in it such as forgetting the multiplication
of omega by two at the recursive calls, I had funny landscapes looking like
a micro-chip seen through an electronic microscop).

	Well, I think that is all, for the moment. Here comes the program.

	Note : you will have probably to cut the signature at the end of the
file as well.

					Jean-Marc

----cut here -----cut here -----cut here -----cut here -----cut here -----

/* It is a pity I have to add this notice, but it must be done.           */

/*                         NOTICE                                         *
 * Copyright (c) 1990, Jean-Marc de Lauwereyns                            *
 *               for the procedures mirror, RecFFT, InvFFT.               *
 * These procedures can be used and given freely to anybody but can not   *
 * be sold or used in any commercial program or in any book without the   *
 * permission of the original author.                                     *
 * These procedures can be modified at the condition that the nature of   *
 * modification is mentioned in comments at the place of the modification *
 * with the original lines accompanied by the name of the original author.*
 * This notice itself can not be removed or modified except by the        *
 * original author himself.                                               */

/* you must have included the file /usr/include/math.h at the begining of *
 * your program								  */

typedef struct {
	double a,b;
	} complexe;

/****************************************************/
/* function mirror just transform the binary        */
/* representation of an integer into its mirror     */
/* binary representation on maxlevel bits           */
/* eg : with maxlevel = 6,  000101 becomes 101000   */
/****************************************************/
int mirror(index,maxlevel)
int index;
int maxlevel;
{
  register int i;
  int res;

  res = 0;
  for (i=0;i<maxlevel;i++)
    {
      res = (res << 1) | (index & (int)1);
      index = (index >> 1);
    }
  return(res);
}

/*****************************************************/
/* recursive procedure to calculate the FFT on 2^N   */
/*****************************************************/
void RecFFT(A,k,l,p,q,Omega)
complexe **A;
int k,l; /* indexes for the begin and the end on the lines taken */
int p,q; /* indexes for the begin and the end of the columns taken */
double Omega;
{
  int i,j,m;
  complexe a,b,c,d;
  double argument[3];

  if (k != l)
    {
      m = (l-k)/2;

      for (i = k;i < k + m;i++)
	{
	  argument[0] = Omega * (double)(i-k);
	  for (j = p;j < p + m;j++)
	    {
	      argument[1] = Omega * (double)(j-p);
	      argument[2] = Omega * (double)(i-k + j-p);

	      a.a = A[i][j].a;a.b = A[i][j].b;
	      b.a = A[i+m][j].a;b.b = A[i+m][j].b;
	      c.a = A[i][j+m].a;c.b = A[i][j+m].b;
	      d.a = A[i+m][j+m].a;d.b = A[i+m][j+m].b;

	      A[i][j].a = a.a + b.a + c.a + d.a;
	      A[i][j].b = a.b + b.b + c.b + d.b;

	      A[i+m][j].a = (a.a - b.a + c.a - d.a)*cos(argument[0]) -
		            (a.b - b.b + c.b - d.b)*sin(argument[0]);
	      A[i+m][j].b = (a.a - b.a + c.a - d.a)*sin(argument[0]) +
		            (a.b - b.b + c.b - d.b)*cos(argument[0]);

	      A[i][j+m].a = (a.a + b.a - c.a - d.a)*cos(argument[1]) -
		            (a.b + b.b - c.b - d.b)*sin(argument[1]);
	      A[i][j+m].b = (a.a + b.a - c.a - d.a)*sin(argument[1]) +
		            (a.b + b.b - c.b - d.b)*cos(argument[1]);

	      A[i+m][j+m].a = (a.a - b.a - c.a + d.a)*cos(argument[2]) -
	                      (a.b - b.b - c.b + d.b)*sin(argument[2]);
	      B[i+m][j+m].b = (a.a - b.a - c.a + d.a)*sin(argument[2]) +
		              (a.b - b.b - c.b + d.b)*cos(argument[2]);
	    }
	}
      if (m != 0)
	{
	  RecFFT(k,k+m,p,p+m,Omega*2);
	  RecFFT(k+m,l,p,p+m,Omega*2);
	  RecFFT(k,k+m,p+m,q,Omega*2);
	  RecFFT(k+m,l,p+m,q,Omega*2);
	}
    }
}

/********************************************************/
/* InvFFT is a procedure which produce the inverse of   */
/* Fourier Transformation.                              */
/********************************************************/
void InvFFT(A,interlevel)
complexe **A;
int interlevel;
{
  int dimension;
  int i,j,i1,j1;
  double coeff,Omega;

  dimension = (1 << interlevel);
  Omega = 2 * M_PI / dimension;
  RecFFT(A,0,dimension,0,dimension,Omega);

  for (i=0;i<=(dimension - 1);i++)
    /* swap every rows i with its mirror image */
    {
      i1 = mirror(i,interlevel);
      if (i1 > i)   /* if i1 <= i, the swaaping have been already made */
	for (j=0;j<dimension;j++)
	  {
	    coeff = A[i][j].a;
	    A[i][j].a = A[i1][j].a;
	    A[i1][j].a = coeff;

	    /* this may be remove if you are sure that the result is a  *
	     * matrix of real which means that the original matrix must *
	     * verify the following conditions :                        *
	     *               ______              ______                 *
	     *  A(N-i,N-j) = A(i,j) , A(0,N-j) = A(0,j)                 *
	     *             ______                                       *
	     *  A(N-i,0) = A(i,0) for i,j > 0    and A(0,0) is a pure   *
	     *  real                                                    */

	    coeff = A[i][j].b;
	    A[i][j].a = A[i1][j].a;
	    A[i1][j].a = coeff;
	  }
    }

  for (j=0;j<=(dimension - 1);j++) /* idem for the columns */
    {
      j1 = mirror(j,interlevel);
      if (j1 >j)
	for (i=0;i<dimension;i++)
	  {
	    coeff = A[i][j].a;
	    A[i][j].a = A[i][j1].a;
	    A[i][j1].a = coeff;

	    /* part which may be removed if (see above)   */
	    coeff = A[i][j].b;
	    A[i][j].b = A[i][j1].b;
	    A[i][j1].b = coeff;
	  }
    }
}

----cut here -----cut here -----cut here -----cut here -----cut here ----
Jean-Marc de Lauwereyns     |              ____  |     e-mail addresses :
Heriot-Watt University      |     |\  /|   |   \ | JANET: lauwer@uk.ac.hw.cs
Computer Science Department |     | \/ |   |___/ | ARPA.: lauwer@cs.hw.ac.uk
Edinburgh                   | \___/    |___|   \ |