Convolution de 2 tableaux avec FFTW

cs_pachalcs Messages postés 5 Date d'inscription samedi 9 mai 2009 Statut Membre Dernière intervention 5 octobre 2011 - 5 oct. 2011 à 16:59
cs_louis14 Messages postés 793 Date d'inscription mardi 8 juillet 2003 Statut Membre Dernière intervention 10 février 2021 - 6 oct. 2011 à 09:23
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 (R = Multiply(BpaddingFFT,ApaddingFFT);) se trouve la fonction qui permet de faire la multiplication des transformées de fourrier.

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

2 réponses

cs_louis14 Messages postés 793 Date d'inscription mardi 8 juillet 2003 Statut Membre Dernière intervention 10 février 2021 8
6 oct. 2011 à 09:18
Bonjour,
Après m'avoir fait confirmer mes souvenirs qui sont un peu vagues, il semble bon de faire du padding et en ce qui concerne la multiplication il faut la faire élément par élément.
Maintenant pour la convolution as-tu regardé du côté de CLAPACK, il y a une fonction pour la faire. De plus tu peux trouver des dll optimisées pour le processeur de ton PC (MKL pour Intel qui contient aussi une foncton fft et ACML pour AMD).
Tu trouveras une discussion sur ton sujet à cette adresse ave un morcau de code qui te servira :http://cboard.cprogramming.com/general-discussions/94954-fft-convolution.html


louis
0
cs_louis14 Messages postés 793 Date d'inscription mardi 8 juillet 2003 Statut Membre Dernière intervention 10 février 2021 8
6 oct. 2011 à 09:23
Bonjour,
Après m'être fait confirmer mes souvenirs , je te conseille de faire le padding et pour la multiplication il faut le faire élément par élément. Tu trouveras une discussion sur ce sujet à cette adresse avec un morceau de code qui peut te servir :http://cboard.cprogramming.com/general-discussions/94954-fft-convolution.html

Mantenat pour la convolution, as-tu regardé du côté de Lapack ( CLapack pour le c++) qui possède une fonction de convolution. De plus il y a des bibliothèque optimisées en foncton du proceseur du PC ( MKL pour intel qui contient auss une fft et ACML pour AMD)

Bon codage.

louis
0