Thanks for everyone's help. 

I've already noticed the possible normalization and FFT direction (+i,
-i) problem, and tried that before the 1st post.

Matteo's method works! I just switched SizeX and SizeZ for FFTW and it
is now fine.

In my understanding, FFTW will take the fastest changing dimension as
Z, and the slowest changing dimension as X. It took me some time to
realize this because in our area (medical image processing), many
people take it for granted that Z is the slowest changing dimension.

This is described in FFTW document for higher dimensions:
http://www.fftw.org/fftw3_doc/Complex-Multi-Dimensional-DFTs.html#Complex%20Multi-Dimensional%20DFTs.
If I read it careful enough, I should have assumed the same thing in
3D.

I would like to thank everybody for the patience to read such a long
post (and code).



Matteo Frigo <athena@fftw.org> wrote in message news:<81y8w5geaj.fsf@ab-initio.mit.edu>...
> hudu2k@yahoo.com (hudu2k) writes:
> 
> >     for( k = 0; k < szz; k++ ) {
> >         for( j = 0; j < szy; j++ ) {
> >             for( i = 0; i < szx; i++ ) {
> >
> >                 in[ k*szx*szy + i + j*szx ][0] = k*szx*szy + i*szy +
> > j;
> >                 in[ k*szx*szy + i + j*szx ][1] = 0;
> >             }
> >         }
> >     }
> 
> Your indexing is wrong.  A(i,j,k) maps to A[k + j * nk + i * nj * nk]
> in C.  (Your order would be correct in FORTRAN.)
> 
> You have three choices:
> 
> 1) fix the indexing to use C conventions
> 2) keep FORTRAN conventions, use the FFTW FORTRAN interface (probably
>    not a good idea if you are writing a C program)
> 3) keep FORTRAN conventions, and flip dimensions when calling
>    the FFTW planner.  E.g., use:
> 
>     p = fftw_plan_dft_3d(szz, szy, szx, in, out,
> 			 FFTW_FORWARD, FFTW_ESTIMATE);
> 
>    instead of
> 
>     p = fftw_plan_dft_3d(szx, szy, szz, in, out,
> 			 FFTW_FORWARD, FFTW_ESTIMATE);
> 
> Regards,
> Matteo Frigo

For what its worth I believe the respective formulae used to calculate
the FFT are:

MATLAB FFT:

                 N
    X(k) =      sum  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
                n=1

FFTW FFT:

                n-1
    Y(k) =      Sum  X(j)*exp(-2*pi*j*k*sqrt(-1)/n)      0<=k<=n-1
                j=0


Matlab arrays are 1-indexed, so we can express the FFTW formulae in
the same form and get:

                n-1
    Y(k) =      Sum  X(j)*exp(-2*pi*j*(k-1)*sqrt(-1)/n)      1<=k<=n
                j=0

and again:

                 N
    Y(k) =      Sum  X(n)*exp(-2*pi*(n-1)*(k-1)*sqrt(-1)/N)  1<=k<=N
                n=1


and rearranging, get back to the MATLAB form.

                 N
    Y(k) =      Sum  X(n)*exp(-j*2*pi*(k-1)*(n-1)/N)      1<=k<=N
                n=1


I.e. the forward transforms should yield the same results.
It is possible you have the same error in the both the FFTW FFT and
IFFT, and as such they cancel. i.e. you are performing an invertable
transform, it is just not an FFT.

Al

You will see Rune's post.
In addition, FFTW doesn't scale the results in *either* direction I
believe - leaving it up to you to scale the results using your favorite
convention.
I believe Matlab and many others scale in the inverse direction by 1/N.
So, as far as scaling goes, I believe that they should agree in the forward
transform results.

With FFTW, if you forward and inverse transform you will get a result that
is N times too big.  This suggests that your Matlab code has a problem
doesn't it if Forward followed by Inverse generates the same results for
either set of code.

(I didn't bother to read your code).

Fred

hudu2k@yahoo.com (hudu2k) writes:

>     for( k = 0; k < szz; k++ ) {
>         for( j = 0; j < szy; j++ ) {
>             for( i = 0; i < szx; i++ ) {
>
>                 in[ k*szx*szy + i + j*szx ][0] = k*szx*szy + i*szy +
> j;
>                 in[ k*szx*szy + i + j*szx ][1] = 0;
>             }
>         }
>     }

Your indexing is wrong.  A(i,j,k) maps to A[k + j * nk + i * nj * nk]
in C.  (Your order would be correct in FORTRAN.)

You have three choices:

1) fix the indexing to use C conventions
2) keep FORTRAN conventions, use the FFTW FORTRAN interface (probably
   not a good idea if you are writing a C program)
3) keep FORTRAN conventions, and flip dimensions when calling
   the FFTW planner.  E.g., use:

    p = fftw_plan_dft_3d(szz, szy, szx, in, out,
			 FFTW_FORWARD, FFTW_ESTIMATE);

   instead of

    p = fftw_plan_dft_3d(szx, szy, szz, in, out,
			 FFTW_FORWARD, FFTW_ESTIMATE);

Regards,
Matteo Frigo

hudu2k@yahoo.com (hudu2k) wrote in message news:<7182cb2a.0309301855.48a58651@posting.google.com>...
> Hi, 
> 
> Just for testing purpose, I want to calculate the FFT for a 4x4x2 (3D)
> matrix. The matrix looks like this:
...
> They are aparently very different. The wierd thing is when I run
> inverse transform, both Matlab and FFTW returns the original result.
> 
> Can anyone tell me which result is correct and why the other is wrong?

Well, there are a couple of reasons why different FFT implementations 
give different results. None of the results need be wrong, though.

First, there is the matter of normalization. In each 1D FFT/IFFT
transform
pair you need to divide with N somewhere along the route, to reproduce
the
original sequence. Some implementations divide by N in the forward 
transform and leave the inverse transform, others divide by N in the 
inverse transform and leave the forward transform, and yet others
divide
by sqrt(N) in both.

Then there is the matter of element ordering. Some implementations
provide
the DC component in the first entry in the spectrum (assuming the
spectrum
to be valid in the frequency range [0,Fs>, Fs being sampling
frequency),
while other implementations provide the DC component as the middle
entry
in the spectrum, assuming a spectrum in the range [-Fs/2, Fs/2>. 

Last, there is the matter of transform conventions. Some people
(electrics
engineers) assume the forward transform to be

   X(w) = integral x(t) exp(-jwt) dt

while others (physicists, mathematicians) assume the forward transform
to be

   X(w) = integral x(t) exp(iwt) dt.

Note the different "-j" and "i". While some people claim that the two 
transforms are equal, since

   i =  sqrt(-1) 
   j = -sqrt(-1)

the important point is that different FFT implementations produce
spectra
that are complex conjugated from each other. Analyzing the data under
one transform convention while unknowingly using a FFT implementation
that
assumes the other,  is not fun. Believe me, it literally takes forever
to
find out what's going on.

Rune

Hi, 

Just for testing purpose, I want to calculate the FFT for a 4x4x2 (3D)
matrix. The matrix looks like this:

in(:,:,1) =
     0     4     8    12
     1     5     9    13
     2     6    10    14
     3     7    11    15
in(:,:,2) =
    16    20    24    28
    17    21    25    29
    18    22    26    30
    19    23    27    31

Matlab Result:

out(:,:,1) =
   1.0e+02 *
   4.9600            -0.6400 + 0.6400i  -0.6400            -0.6400 -
0.6400i
  -0.1600 + 0.1600i        0                  0                  0
  -0.1600                  0                  0                  0
  -0.1600 - 0.1600i        0                  0                  0
out(:,:,2) =
  -256     0     0     0
     0     0     0     0
     0     0     0     0
     0     0     0     0

FFTW Result: (x,y) means x+j*y

FFT
k=0
(496.00, 0.00)  (-64.00, 0.00)  (-8.00, 8.00)   (0.00, 0.00)
(-128.00, 0.00) (0.00, 0.00)    (-8.00, -8.00)  (0.00, 0.00)
(-128.00, 128.00)       (0.00, 0.00)    (0.00, 0.00)    (0.00, 0.00)
(0.00, 0.00)    (0.00, 0.00)    (0.00, 0.00)    (0.00, 0.00)
k=1
(-32.00, 0.00)  (0.00, 0.00)    (0.00, 0.00)    (0.00, 0.00)
(0.00, 0.00)    (0.00, 0.00)    (0.00, 0.00)    (0.00, 0.00)
(-128.00, -128.00)      (0.00, 0.00)    (0.00, 0.00)    (0.00, 0.00)
(0.00, 0.00)    (0.00, 0.00)    (0.00, 0.00)    (0.00, 0.00)

They are aparently very different. The wierd thing is when I run
inverse transform, both Matlab and FFTW returns the original result.

Can anyone tell me which result is correct and why the other is wrong?

Following is my source code.

Matlab Source:

in = zeros(4,4,2); in(:) = 0:31
out = fftn(in)

C Source: (require FFTW)

int test3d_r2c() {
    
    fftw_complex *in, *out;
    int szx, szy, szz;
    fftw_plan p,p2;
    int i, j, k;
    
    szx = 4; szy = 4; szz = 2; 

    in = fftw_malloc(sizeof(fftw_complex) * szx * szy * szz );
    out = fftw_malloc(sizeof(fftw_complex) * szx * szy * szz );

    p = fftw_plan_dft_3d( szx, szy, szz,in, out, FFTW_FORWARD,
FFTW_ESTIMATE );
    p2 = fftw_plan_dft_3d( szx, szy, szz,out, in, FFTW_BACKWARD ,
FFTW_ESTIMATE );

    for( k = 0; k < szz; k++ ) {
        for( j = 0; j < szy; j++ ) {
            for( i = 0; i < szx; i++ ) {

                in[ k*szx*szy + i + j*szx ][0] = k*szx*szy + i*szy +
j;
                in[ k*szx*szy + i + j*szx ][1] = 0;
            }
        }
    }

    fftw_execute(p);
    
    
    printf("Original\n");
    for( k = 0; k < szz; k++ ) {
    printf("k=%d\n",k); 
        for( j = 0; j < szy; j++ ) {
            for( i = 0; i < szx; i++ ) {

                printf("(%.2f, %.2f)\t", 
                    in[k*szx*szy + i + j * szx][0], 
                    in[k*szx*szy + i + j * szx][1]);

            }
            printf("\n");
        }
    
    }

    printf("\nFFT\n");
    for( k = 0; k < szz; k++ ) {
    printf("k=%d\n",k);  
        for( j = 0; j < szy; j++ ) {
            for( i = 0; i < szx; i++ ) {

                printf("(%.2f, %.2f)\t", (double)out[k*szx*szy + i + j
* szx][0], (double)out[k*szx*szy + i + j * szx][1]);

            }
            printf("\n");
        }
    
    }

    fftw_execute(p2);
    
    printf("\nInverse\n");
    for( k = 0; k < szz; k++ ) {
    printf("k=%d\n",k); 
        for( j = 0; j < szy; j++ ) {
            for( i = 0; i < szx; i++ ) {

                printf("(%.2f, %.2f)\t", 
                    in[k*szx*szy + i + j * szx][0]/szx/szy/szz, 
                    in[k*szx*szy + i + j * szx][1]/szx/szy/szz);

            }
            printf("\n");
        }
    
    }

    fftw_destroy_plan(p);
    fftw_destroy_plan(p2);
    fftw_free(in); fftw_free(out);
    return 0;
}