Obtenez des résultats égaux avec fft2 (matlab) et fftw (C)

Question

J'essaie d'implémenter la fonction Matlab fft2() en C en utilisant la bibliothèque FFTW3.

Cependant, j'ai des résultats différents. Considérant une matrice [1 2; 3 4], le résultat avec Matlab fft2 est [10 -2; -4 0] et avec FFTW3 [10 -2; (3.05455e-314 + 2.122e-314i) (9.31763e-315 + 9.32558e-315i)]

J'ai utilisé le code suivant:

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); 

Qu'est-ce que je fais mal?

Merci d'avance.

Answer 1

Il semble fonctionner OK pour moi:

// 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; 
} 

Compilez et exécutez:

$ 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 } 
$ 

Je me demande si c'est juste que le code que vous utilisez pour afficher les résultats est incorrect ?

---

Pour le cas 3x3 Je suis aussi obtenir correspondant:

// 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; 
} 

Compilez et exécutez:

$ 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 } 
$ 

Comparer avec Octave (clone Matlab):

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

Votre code a l'air correct. Je l'ai essayé avec GCC gcc-branche révision 4_7 189773 sur une machine 64 bits, et il revient

10 -2 
-4 0 

comme il se doit.

Également essayé différentes techniques de sortie avec printf formats (% f,% g,% lg), et l'accès aux données avec creal pour faire bonne mesure - il renvoie toujours les résultats attendus. Utilisez-vous GCC avec quad_precision activé? C'est l'une des choses que je ne peux pas tester maintenant.