Discussion :

[pachalcs]

Bonjour,

Je viens tout juste de télécharger la bibliotheque fftw3 et j'ai lu le tutorial.
En fait, mon but est de faire la convolution de deux énormes tableaux de double, et je veux utiliser la transformée de fourrier pour des critères de rapidité.

Je sais que pour faire la convolution, il faut suivre cet algorithme ci dessous:

Etape 1
TF(A)=FFT(A) avec A un de nos tableaux de départ de taille M
TF(B)=FFT(B) avec B un de nos tableaux de départ de taille N

Etape 2
Puis faire TF(A)*TF(B)

Etape 3
Et finir par faire l'inversion en obtenant Conv(A,B) = IFFT( TF(A)*TF(B) ) qui aura une taille égale à M+N-1.


Mon problème se situe à l'étape 2, je ne sais vraiment pas comment implémenter cette étape.
Je sais qu'il faut faire un zéro padding sur les vecteurs A ET B pour les avoir avec une taille égale à M+N-1 avant de faire leur transformée de fourrier

Aidez moi, c'est la seule étape qui reste pour que je finisse un projet.
Je vous mets mon code à la suite et à la ligne 53 se trouve la fonction qui permet de faire la multiplication des transformées de fourrier

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
 
 
void confftW (fftw_complex* A, fftw_complex* B, int M, int N) {
 
    fftw_complex* Apadding, * Bpadding, ApaddingFFT,BpaddingFFT*R, *IR;
    fftw_plan     plan_a, plan_b, plan_R;
    int           i;
 
    /* Allocate input & output array */
    Apadding = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * M+N-1);
    Bpadding = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * M+N-1);
 
 
    ApaddingFFT = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * M+N-1);
    BpaddingFFT = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * M+N-1);
    IR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * M+N-1);
 
 
//0-padding de A et B
 for (i = 0; i < M+N-1; i++)
    {
       //Pour Simplifier on considère que A et B sont de même taille
       //M=N
       if (i<M)
       {
       Apadding [i][0] = A[i][0];
       Bpadding [i][0] = B[i][0];
       }
       else
       {
       Apadding [i][0] = 0.0;
       Bpadding [i][0] = 0.0;  
       }
 
       Apadding [i][1] = 0.0;
       Bpadding [i][1] = 0.0;  
    }
 
    /* Create plans */
    plan_a = fftw_plan_dft_1d(M+N-1, Apadding , ApaddingFFT, FFTW_FORWARD,  FFTW_ESTIMATE);
 
    plan_b = fftw_plan_dft_1d(M+N-1, Bpadding , BpaddingFFT, FFTW_FORWARD,  FFTW_ESTIMATE);
 
 
 
    plan_R = fftw_plan_dft_1d(M+N-1, R, IR, FFTW_BACKWARD, FFTW_ESTIMATE);
 
 
    /*Exécution des TF de A et B*/
    fftw_execute(plan_a);
    fftw_execute(plan_b);
 
    R = Multiply(BpaddingFFT,ApaddingFFT); // FONCTION QUE JE VEUX IMPLEMENTER
 
    //Exécution de la transformée de fourier inverse de R qui va stocker
    // le resultat dans IR
    fftw_execute(plan_R);
 
 
    /* Free memory */
    fftw_destroy_plan(plan_a);
    fftw_destroy_plan(plan_b);
    fftw_destroy_plan(plan_R);
    fftw_free(Apadding );
    fftw_free(Bpadding );
    fftw_free(ApaddingFFT );
    fftw_free(BpaddingFFT );
    fftw_free(IR);
 
}

Je vous remercie d'avance

[Wolfen]

Un truc comme ça ?

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
// Multiplication complexe
for (i = 0; i < (M+N-1)/2; i++)
{   R[i][0] = ApaddingFFT[i][0]*BpaddingFFT[i][0] - ApaddingFFT[i][1]*BpaddingFFT[i][1];
    R[i][1] = ApaddingFFT[i][0]*BpaddingFFT[i][1] + ApaddingFFT[i][1]*BpaddingFFT[i][0];
}
 
fftw_execute(plan_R);
 
// Ne pas oublier de normaliser la IFFT
for (i = 0; i < (M+N-1); i++)
{   IR[i] = IR[i] / (float) (M+N-1);
}