Discussion :

[Heimdall]

Salut,

J'ai un petit probleme de compréhension de la doc de FFTW. j'aurai besoin de votre aide ou d'un petit exemple

Alors voila, j'ai un signal 2D S(x,y) dont j'aimerais bien faire la TF sur une et une seule dimension (x disons).

Mes données sont dans un tableau de nx+1 points suivant la direction x, par ny+1 points suivant la direction y. (nx = 128 et ny = 64 par exemple)

Le point S(x,y) est S_{i,j} sur la grille, ou l'indice du tableau (1D) correspondant est i + j*(nx+1).


Apparement FFTW permet de faire ce genre de chose, via l'interface avancée (advanced interface). Mais je ne comprends pas vraiment tous les parametres. En particulier, qu'elle est la différence entre n et inembed ?

je poste ici un petit code ou j'essaie de tester ça, mais apparement je me goure dans les tailles je pense :

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
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
#include <fftw3.h>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
 
 
int main(void)
{
fftw_plan pb, pf;
double *in, *back;
fftw_complex *out;
 
int nx, ny;
int i, j, ij;
double xm, ym;
 
 
int n[1];
int inembed[1];
int onembed[1];
int rank;
int howmany;
int istride, ostride;
int idist, odist;
double x,y;
 
 
FILE *fpin, *fpback;
 
nx = 128;
ny = 64;
xm = 20.0;
ym = 20.0;
 
in  = fftw_malloc((nx+1)*(ny+1)*sizeof(*in));
out = fftw_malloc((nx/2+1)*(ny+1)*sizeof(*out)); 
back  = fftw_malloc((nx+1)*(ny+1)*sizeof(*back));
 
rank = 1;
howmany = ny+1;
istride = 1;
ostride = 1;
idist = nx+1; 
odist = nx/2+1;
 
 
n[0] = nx;
 
inembed[0] = nx+1; 
onembed[0] = nx/2+1;
 
 
fpin = fopen("in.txt" , "w");
for (i = 0; i < nx+1; i++)
    {
    for (j = 0; j < ny+1;  j++)
        {
        ij = i + j*(nx+1);
        x = i*xm/(double)nx;
        y = j*ym/(double)ny;
        in[ij] = cos(2.0*M_PI*x/xm)*sin(2.0*M_PI*y/ym);
        fprintf(fpin, "%lg %lg %lg\n", x, y, in[ij]);
        }
    }
fclose(fpin);
 
 
 
pf =  fftw_plan_many_dft_r2c(rank, n, howmany, in, 
                             inembed, istride, idist, out,
                             onembed, ostride, odist, FFTW_MEASURE); 
 
fftw_execute(pf);
 
rank = 1;
howmany = ny+1;
istride = 1;
ostride = 1;
idist = nx/2+1; 
odist = nx+1;
 
n[0] = nx/2+1;
 
inembed[0] = nx/2+1; 
onembed[0] = nx+1;
 
 
 
pb =  fftw_plan_many_dft_c2r(rank, n, howmany, 
                             out, inembed, istride, 
                             idist, back, onembed,
                             ostride, odist, FFTW_MEASURE); 
 
fftw_execute(pb);
 
 
fpback = fopen("back.txt" , "w");
for (i = 0; i < nx+1; i++)
    {
    for (j = 0; j < ny+1;  j++)
        {
        ij = i + j*(nx+1);
        x = i*xm/(double)nx;
        y = j*ym/(double)ny;
        fprintf(fpback, "%lg %lg %lg\n", x, y, back[ij]);
        }
    }
fclose(fpback);
 
 
 
fftw_destroy_plan(pf);
fftw_destroy_plan(pb);
free(in);
free(out);
free(back);
 
return 0;
}

[mchk0123]

Si tu veux lancer la FFT uniquement sur la dimension x, je te conseil de faire plusieurs fftw_execute de taille nx, et ce, pour chaque ligne.

Donc déjà tu doit modifier :

Code :

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

1
2
3
4
5
in  = fftw_malloc((nx+1)*sizeof(*in));
out = fftw_malloc((nx/2+1)*sizeof(*out)); 
back  = fftw_malloc((nx+1)*sizeof(*back));
...

Le reste se déduit aisément.
Dans un premier temps, essaye de faire plusieurs fftw_execute plutôt que de te lancer dans l'utilisation de "_many" qui est plus compliquée. Si ça marche, tu pourras alors faire compliqué.

[Heimdall]

Si tu veux lancer la FFT uniquement sur la dimension x, je te conseil de faire plusieurs fftw_execute de taille nx, et ce, pour chaque ligne.

Donc déjà tu doit modifier :

Code :

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

1
2
3
4
5
in  = fftw_malloc((nx+1)*sizeof(*in));
out = fftw_malloc((nx/2+1)*sizeof(*out)); 
back  = fftw_malloc((nx+1)*sizeof(*back));
...

Le reste se déduit aisément.
Dans un premier temps, essaye de faire plusieurs fftw_execute plutôt que de te lancer dans l'utilisation de "_many" qui est plus compliquée. Si ça marche, tu pourras alors faire compliqué.

Ok merci de ta réponse. Dans la doc ils disent qu'utiliserl'interface avanceée avec les _many_dft est plus rapide que de fais les N tf "a la main".

J'ai "réussi" a m'en servir mais c'est un immense basard lorsuqe tu veux faire unr r2r dans l'autre direction.

[mchk0123]

Ok merci de ta réponse. Dans la doc ils disent qu'utiliserl'interface avanceée avec les _many_dft est plus rapide que de fais les N tf "a la main".

J'ai "réussi" a m'en servir mais c'est un immense basard lorsuqe tu veux faire unr r2r dans l'autre direction.

C'est vrai j'avait lu ça dans la documentation. Maintenant la FFTW est suffisament optimisée pour que l'utilisation des fftw_execute simples soient déjà acceptablement rapide.
Et à la reflexion, je me demande si l'utilisation des _many ne font pas uniquement l'économie des appels multiples de fonctions (les formules de calcul restent les mêmes) ; donc le gain en temps devrait être trés petit.