Convolution de 2 tableaux unidimensionnels avec FFTW C++

Convolution de 2 tableaux unidimensionnels avec FFTW C++ - C++ - Programmation

Marsh Posté le 05-10-2011 à 16:49:41    

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 51 se trouve la fonction qui permet de faire la multiplication des transformées de fourrier
 

Code :
  1. void confftW (fftw_complex* A, fftw_complex* B, int M, int N) {
  3.     fftw_complex* Apadding, * Bpadding, ApaddingFFT,BpaddingFFT*R, *IR;
  4.     fftw_plan     plan_a, plan_b, plan_R;
  5.     int           i;
  7.     /* Allocate input & output array */
  8.     Apadding = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * M+N-1);
  9.     Bpadding = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * M+N-1);
  12.     ApaddingFFT = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * M+N-1);
  13.     BpaddingFFT = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * M+N-1);
  14.     IR = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * M+N-1);
  17. //0-padding de A et B
  18. for (i = 0; i < M+N-1; i++)
  19.     {
  20.        //Pour Simplifier on considère que A et B sont de même taille
  21.        //M=N
  22.        if (i<M)
  23.        {
  24.        Apadding [i][0] = A[i][0];
  25.        Bpadding [i][0] = B[i][0];
  26.        }
  27.        else
  28.        {
  29.        Apadding [i][0] = 0.0;
  30.        Bpadding [i][0] = 0.0; 
  31.        }
  33.        Apadding [i][1] = 0.0;
  34.        Bpadding [i][1] = 0.0; 
  35.     }
  37.     /* Create plans */
  38.     plan_a = fftw_plan_dft_1d(M+N-1, Apadding , ApaddingFFT, FFTW_FORWARD,  FFTW_ESTIMATE);
  40.     plan_b = fftw_plan_dft_1d(M+N-1, Bpadding , BpaddingFFT, FFTW_FORWARD,  FFTW_ESTIMATE);
  44.     plan_R = fftw_plan_dft_1d(M+N-1, R, IR, FFTW_BACKWARD, FFTW_ESTIMATE);
  47.     /*Exécution des TF de A et B*/
  48.     fftw_execute(plan_a);
  49.     fftw_execute(plan_b);
  51.     R = Multiply(BpaddingFFT,ApaddingFFT); // FONCTION QUE JE VEUX IMPLEMENTER
  53.     //Exécution de la transformée de fourier inverse de R qui va stocker
  54.     // le resultat dans IR
  55.     fftw_execute(plan_R);
  58.     /* Free memory */
  59.     fftw_destroy_plan(plan_a);
  60.     fftw_destroy_plan(plan_b);
  61.     fftw_destroy_plan(plan_R);
  62.     fftw_free(Apadding );
  63.     fftw_free(Bpadding );
  64.     fftw_free(ApaddingFFT );
  65.     fftw_free(BpaddingFFT );
  66.     fftw_free(IR);
  68. }


 
Je vous remercie d'avance

Reply

Marsh Posté le 05-10-2011 à 16:49:41   

Reply

Marsh Posté le 06-10-2011 à 17:30:09    

ya rien à faire, c'est une multiplication élément par élément

Reply

Sujets relatifs:

Leave a Replay

Make sure you enter the(*)required information where indicate.HTML code is not allowed