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