Discussion :

[lock42]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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

[Invité]

Bonjour,

Cette ligne est à reconsidérer

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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

[lock42]

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.

[lock42]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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

[lock42]

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

[Invité]

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

[lock42]

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

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