Get equal results with fft2 (matlab) and fftw (C)

Question

I'm trying to implement the Matlab fft2() function in C using the FFTW3 library.

However, I've got different results.

Considering a matrix [1 2; 3 4], the result with Matlab fft2 is [10 -2; -4 0] and with FFTW3 [10 -2; (3.05455e-314 + 2.122e-314i) (9.31763e-315 + 9.32558e-315i)]

I used the following code:

fftw_plan planG;
fftw_complex *inG, *outG;

inG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * 2 * 2);
outG = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * 2 * 2);

inG[0][0]=1.0;
inG[1][0]=2.0;
inG[2][0]=3.0;
inG[3][0]=4.0;

inG[0][1]=0.0;
inG[1][1]=0.0;
inG[2][1]=0.0;
inG[3][1]=0.0;

planG = fftw_plan_dft_2d(2, 2, inG, outG, FFTW_FORWARD, FFTW_ESTIMATE);

fftw_execute(planG);

What am I doing wrong?

Thanks in advance.

Answer 1

It seems to work OK for me:

// fftw_test.c

#include <stdio.h>
#include <fftw3.h>

int main(int argc, char * argv[])
{
    fftw_plan planG;
    fftw_complex *inG, *outG;

    inG = fftw_malloc(sizeof(fftw_complex) * 2 * 2);
    outG = fftw_malloc(sizeof(fftw_complex) * 2 * 2);

    inG[0][0] = 1.0; // real input
    inG[1][0] = 2.0;
    inG[2][0] = 3.0;
    inG[3][0] = 4.0;

    inG[0][1] = 0.0; // imag input
    inG[1][1] = 0.0;
    inG[2][1] = 0.0;
    inG[3][1] = 0.0;

    planG = fftw_plan_dft_2d(2, 2, inG, outG, FFTW_FORWARD, FFTW_ESTIMATE);

    fftw_execute(planG);

    printf("outg[0] = { %g, %g }\n", outG[0][0], outG[0][1]);
    printf("outg[1] = { %g, %g }\n", outG[1][0], outG[1][1]);
    printf("outg[2] = { %g, %g }\n", outG[2][0], outG[2][1]);
    printf("outg[3] = { %g, %g }\n", outG[3][0], outG[3][1]);

    return 0;
}

Compile and run:

$ gcc -Wall -lfftw3 fftw_test.c -o fftw_test
$ ./fftw_test 
outg[0] = { 10, 0 }
outg[1] = { -2, 0 }
outg[2] = { -4, 0 }
outg[3] = { 0, 0 }
$ 

I'm wondering if it's just that the code you're using to display the results is incorrect ?

---

For the 3x3 case I'm also getting matching:

// fftw_test_3_x_3.c

#include <stdio.h>
#include <fftw3.h>

#define M 3
#define N 3

int main(int argc, char * argv[])
{
    fftw_plan planG;
    fftw_complex *inG, *outG;
    int i;

    inG = fftw_malloc(sizeof(fftw_complex) * M * N);
    outG = fftw_malloc(sizeof(fftw_complex) * M * N);

    for (i = 0; i < M * N; ++i)
    {
        inG[i][0] = (double)(i + 1);
        inG[i][1] = 0.0; // imag input
    }

    planG = fftw_plan_dft_2d(M, N, inG, outG, FFTW_FORWARD, FFTW_ESTIMATE);

    fftw_execute(planG);

    for (i = 0; i < M * N; ++i)
    {
        printf("outg[%d] = { %g, %g }\n", i, outG[i][0], outG[i][1]);
    }

    return 0;
}

Compile and run:

$ gcc -Wall -lfftw3 fftw_test_3_x_3.c -o fftw_test_3_x_3
$ ./fftw_test_3_x_3 
outg[0] = { 45, 0 }
outg[1] = { -4.5, 2.59808 }
outg[2] = { -4.5, -2.59808 }
outg[3] = { -13.5, 7.79423 }
outg[4] = { 0, 0 }
outg[5] = { 0, 0 }
outg[6] = { -13.5, -7.79423 }
outg[7] = { 0, 0 }
outg[8] = { 0, 0 }
$

Compare with Octave (MATLAB clone):

octave:1> a = [1 2 3 ; 4 5 6 ; 7 8 9]
a =

   1   2   3
   4   5   6
   7   8   9

octave:2> b = fft2(a, 3, 3)
b =

   45.00000 +  0.00000i   -4.50000 +  2.59808i   -4.50000 -  2.59808i
  -13.50000 +  7.79423i    0.00000 +  0.00000i    0.00000 +  0.00000i
  -13.50000 -  7.79423i    0.00000 -  0.00000i    0.00000 -  0.00000i

octave:2> 

Answer 2

Your code looks correct. I tried it with GCC gcc-4_7-branch revision 189773 on a 64bit machine, and it is returning

10 -2
-4 0

as it should.

Also tried different output techniques with printf formats (%f, %g, %lg), and accessing data with creal for good measure -- it is always returning the expected results.

Are you using GCC with quad_precision enabled? That's one of the things I can't test now.