Convolution de 2 tableaux avec FFTW

Question

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

cs_louis14 · Answer

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

cs_louis14 · Answer

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

Convolution de 2 tableaux avec FFTW

2 réponses

Votre réponse

Discussions similaires