Utilisation de fftw

Bonjour tout le monde. Je poste ici car je suis vraiment au bout du rouleau et au bout de tous mes tests. Je suis en train de devenir fou. En effet, j'essaye d'utiliser FFTW pour calculer la transformer de Fourrier d'une image. A des fins de test j'utilise une image carrée de 4px par 4 px.
La transformation directe est un succès (même pour les plus grandes images)

Voici un bout de code issu de cette partie et me donnant par conséquent les valeurs de mon input:
Code:

1 2 3 4 for (i=0; i<nbdata; i++){ spatial_repr[i]=(double)gbuf[i]; printf("%lf\n", (double)spatial_repr[i]); }
valeurs affichées :
0,000000
0,000000
0,000000
0,000000
0,000000
255,000000
255,000000
0,000000
0,000000
255,000000
255,000000
0,000000
0,000000
0,000000
0,000000
0,000000

Cette autre partie de la transformée directe qui marche bien
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 ......... /* we run the Fourier Transform */ fftw_plan p = fftw_plan_dft_2d(width, height, spatial_repr, frequency_repr, FFTW_FORWARD, FFTW_ESTIMATE); fftw_execute(p); /* we compute magnitude and phase */ double max_m; double *magnitude = malloc(nbdata * sizeof(double)); double *phase = malloc(nbdata * sizeof(double)); fft_to_spectra(nbdata, frequency_repr, magnitude, phase, &max_m); printf("max=%lf\n", max_m); for(j=0;j<height;j++){ for(i=0;i<width;i++){ unsigned int x=i; unsigned int y=j; //~ change_symmetry(width, height, i, j, &x, &y); frequency_repr[j*width+i] = magnitude[y*width+x] + I*phase[y*width+x]; printf("%lf + I * %lf\n", creal(frequency_repr[j*width+i]), cimag(frequency_repr[j*width+i])); } } ......
Donne ces valeurs en sortie :

1020,000000 + I * 0,000000
721,248917 + I * -2,356194
0,000000 + I * 0,000000
721,248917 + I * 2,356194
721,248917 + I * -2,356194
510,000000 + I * 1,570796
0,000000 + I * 0,000000
510,000000 + I * 0,000000
0,000000 + I * 0,000000
0,000000 + I * 0,000000
0,000000 + I * 0,000000
0,000000 + I * 0,000000
721,248917 + I * 2,356194
510,000000 + I * 0,000000
0,000000 + I * 0,000000
510,000000 + I * -1,570796

Mon problème apparait quand j'essaye de faire la transformée inverse. Cela me donne de mauvais résultats, et ce, meme si je remplie frequency_repr des mêmes valeurs qu'il a en sorti de ma FFT Direct ..... Du coup, je ne comprend pas du tout :(.
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 ....... /* TEST */ frequency_repr[0] = 1020.000000 + I* (0.000000 ); frequency_repr[1] = 721.248917 + I* (-2.356194); frequency_repr[2] = 0.000000 + I* (0.000000 ); frequency_repr[3] = 721.248917 + I* (2.356194 ); frequency_repr[4] = 721.248917 + I* (-2.356194); frequency_repr[5] = 510.000000 + I* (1.570796 ); frequency_repr[6] = 0.000000 + I* (0.000000 ); frequency_repr[7] = 510.000000 + I* (0.000000 ); frequency_repr[8] = 0.000000 + I* (0.000000 ); frequency_repr[9] = 0.000000 + I* (0.000000 ); frequency_repr[10] = 0.000000 + I* (0.000000 ); frequency_repr[11] = 0.000000 + I* (0.000000 ); frequency_repr[12] = 721.248917 + I* (2.356194 ); frequency_repr[13] = 510.000000 + I* (0.000000 ); frequency_repr[14] = 0.000000 + I* (0.000000 ); frequency_repr[15] = 510.000000 + I* (-1.570796); /* */ fftw_plan p = fftw_plan_dft_2d(width, height, frequency_repr, spatial_repr, FFTW_BACKWARD, FFTW_ESTIMATE); fftw_execute(p); for (i=0; i<nbdata; i++){ printf("%lf\n", creal(spatial_repr[i]) / nbdata); double pxl = creal(spatial_repr[i]) / nbdata; gbuf[i] = (pxl < 0) ? 0 : ((pxl > USHRT_MAX) ? USHRT_MAX : (WORD)pxl); } ........
Si quelqu'un pouvait m'aider sur ce problème, je lui serais très reconnaissant. Il y'a quelque chose que je dois mal comprendre et mal faire dans l'utilisation de la librairie, c'est certain, mais cela m'échappe.

Je vous souhaite une bonne journée,
Bien cordialement

Bonjour,

Cette ligne est à reconsidérer ;)

Citation:

Envoyé par lock42

Code:

frequency_repr[j*width+i] = magnitude[y*width+x] + I*phase[y*width+x];

Bonjour,
Tout d'abord merci de t'être arrété sur mon problème.

Ensuite, je ne suis pas sur de comprendre pourquoi je dois reconsidérer cette ligne, désolé. J'ai de bon résultats pour ma FFT direct. Donc le code marche dans mon esprit.

La j'ai commenté la partie qui me change la symétrie de l'image, mais c'est assez utile quand je veux centrer le tout.

Voici ma fonction fft_to_spectra(int nbdata, fftw_complex *frequency_repr, double *as, double *ps, double *max_m).
Mais la encore, je le rappelle, ca semble bien marcher dans ce sens.
Code:

1 2 3 4 5 6 7 8 9 void fft_to_spectra(int nbdata, fftw_complex *frequency_repr, double *as, double *ps, double *max_m){ *max_m = 0.0; unsigned int i; for(i = 0; i < nbdata; i++){ as[i] = hypot(creal(frequency_repr[i]), cimag(frequency_repr[i])); ps[i] = atan2(cimag(frequency_repr[i]), creal(frequency_repr[i])); *max_m=max(*max_m, as[i]); } }

Citation:

Envoyé par Winjerome

Bonjour,

Cette ligne est à reconsidérer ;)

En effet je commence a comprendre. Je ne tenais pas en compte ce changement ds ma fft inverse.

Quand tu me demandes de la reconsidérer la ligne qu'entends tu par la ? As tu un conseil ?
En tout cas merci ca vient de m'aiguiller sur mon problème et ca c'est déjà énorme !!

EDIT : En fait je peux peut être virer cette ligne et intégrer ma normalisation dans ce bloc.
Bref la je suis sans ordi donc c'est dur a imaginer :)

Encore merci

Ton calcul est juste faux : tu n'obtiens pas la valeur complexe avec mag + i * phase, révise un peu tes cours de lycée :P

Erf.
Oui au temps pour moi. Si j'avais écrit le truc sur le papier avant ....

Je suis désolé, honteux et confus.
Merci.