from http://astronomy.swin.edu.au/pbourke/analysis/fft2d/
/*-------------------------------------------------------------------------
Perform a 2D FFT inplace given a complex 2D array
The direction dir, 1 for forward, -1 for reverse
The size of the array (nx,ny)
Return false if there are memory problems or
the dimensions are not powers of 2
*/
int FFT2D(COMPLEX **c,int nx,int ny,int dir)
{
int i,j;
int m,twopm;
double *real,*imag;
/* Transform the rows */
real = (double *)malloc(nx * sizeof(double));
imag = (double *)malloc(nx * sizeof(double));
if (real == NULL || imag == NULL)
return(FALSE);
if (!Powerof2(nx,&m,&twopm) || twopm != nx)
return(FALSE);
for (j=0;j<ny;j++) {
for (i=0;i<nx;i++) {
real[i] = c[i][j].real;
imag[i] = c[i][j].imag;
}
FFT(dir,m,real,imag);
for (i=0;i<nx;i++) {
c[i][j].real = real[i];
c[i][j].imag = imag[i];
}
}
free(real);
free(imag);
/* Transform the columns */
real = (double *)malloc(ny * sizeof(double));
imag = (double *)malloc(ny * sizeof(double));
if (real == NULL || imag == NULL)
return(FALSE);
if (!Powerof2(ny,&m,&twopm) || twopm != ny)
return(FALSE);
for (i=0;i<nx;i++) {
for (j=0;j<ny;j++) {
real[j] = c[i][j].real;
imag[j] = c[i][j].imag;
}
FFT(dir,m,real,imag);
for (j=0;j<ny;j++) {
c[i][j].real = real[j];
c[i][j].imag = imag[j];
}
}
free(real);
free(imag);
return(TRUE);
}
/*-------------------------------------------------------------------------
This computes an in-place complex-to-complex FFT
x and y are the real and imaginary arrays of 2^m points.
dir = 1 gives forward transform
dir = -1 gives reverse transform
Formula: forward
N-1
---
1 \ - j k 2 pi n / N
X(n) = --- > x(k) e = forward transform
N / n=0..N-1
---
k=0
Formula: reverse
N-1
---
\ j k 2 pi n / N
X(n) = > x(k) e = forward transform
/ n=0..N-1
---
k=0
*/
int FFT(int dir,int m,double *x,double *y)
{
long nn,i,i1,j,k,i2,l,l1,l2;
double c1,c2,tx,ty,t1,t2,u1,u2,z;
/* Calculate the number of points */
nn = 1;
for (i=0;i<m;i++)
nn *= 2;
/* Do the bit reversal */
i2 = nn >> 1;
j = 0;
for (i=0;i<nn-1;i++) {
if (i < j) {
tx = x[i];
ty = y[i];
x[i] = x[j];
y[i] = y[j];
x[j] = tx;
y[j] = ty;
}
k = i2;
while (k <= j) {
j -= k;
k >>= 1;
}
j += k;
}
/* Compute the FFT */
c1 = -1.0;
c2 = 0.0;
l2 = 1;
for (l=0;l<m;l++) {
l1 = l2;
l2 <<= 1;
u1 = 1.0;
u2 = 0.0;
for (j=0;j<l1;j++) {
for (i=j;i<nn;i+=l2) {
i1 = i + l1;
t1 = u1 * x[i1] - u2 * y[i1];
t2 = u1 * y[i1] + u2 * x[i1];
x[i1] = x[i] - t1;
y[i1] = y[i] - t2;
x[i] += t1;
y[i] += t2;
}
z = u1 * c1 - u2 * c2;
u2 = u1 * c2 + u2 * c1;
u1 = z;
}
c2 = sqrt((1.0 - c1) / 2.0);
if (dir == 1)
c2 = -c2;
c1 = sqrt((1.0 + c1) / 2.0);
}
/* Scaling for forward transform */
if (dir == 1) {
for (i=0;i<nn;i++) {
x[i] /= (double)nn;
y[i] /= (double)nn;
}
}
return(TRUE);
}
>I just tried to use your FFT algorithm. It compiles and run fine, >but the result is not what expected. What should be in the center >of the image is split between the four corners of the image. That is correct, if you want DC to be in the center you need to rearrange the grid. Have a look at how my 1D transform works, that will describe how the positive and negative frequencies are arranged. While having DC in the center is convenient for dispay purposes, it isn't the most convenient for typical subsequent processing.
> Could send me a sample program to show me how to use the FFT2D > function...... Try this http://astronomy.swin.edu.au/pbourke/analysis/fft2d/2dfft.c