Convolution de 2 tableaux unidimensionnels FFTW
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:
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
