This is C code for 2 dim FFTs for data length a power of two.
thanks to John W. Herman for the pointer to the code !
this is also avialable at
ftp://ftp.teleos.com/VISION-LIST-ARCHIVE/SHAREWARE/CODE/FFT/FFT.shar
from the file:
** forward_fft(array, rows, cols)
** COMPLEX *array;
** int rows, cols;
**
** inverse_fft(array, rows, cols)
** COMPLEX *array;
** int rows, cols;
**
** These entry points compute the forward and inverse DFT's, respectively,
** of a single-precision COMPLEX array.
**
** The result is a COMPLEX array of the same size, returned in
** the same space as the input array. That is, the original array is
** overwritten and destroyed.
**
** Rows and columns must each be an integral power of 2.
**
** These routines return integer value ERROR if an error was detected,
** NO_ERROR otherwise
**
** This implementation of the DFT uses the transform pair defined as follows.
**
** Let there be two COMPLEX arrays each with n rows and m columns
** Index them as
** f(x,y): 0 <= x <= m - 1, 0 <= y <= n - 1
** F(u,v): -m/2 <= u <= m/2 - 1, -n/2 <= v <= n/2 - 1
**
** Then the forward and inverse transforms are related as
**
** Forward:
**
** F(u,v) = \sum_{x=0}^{m-1} \sum_{y=0}^{n-1}
** f(x,y) \exp{-2\pi i (ux/m + vy/n)}
**
**
** Inverse:
**
** f(x,y) = 1/(mn) \sum_{u=-m/2}^{m/2-1} \sum_{v=-n/2}^{n/2-1}
** F(u,v) \exp{2\pi i (ux/m + vy/n)}
**
** Therefore, the transforms have these properties:
** 1. \sum_x \sum_y f(x,y) = F(0,0)
** 2. m n \sum_x \sum_y |f(x,y)|^2 = \sum_u \sum_v |F(u,v)|^2