Version imprimable
La première étape calcule r[2i] et r[2i+1] en fonction des anciens r[2i] et r[2i+1].
La deuxième r[i] et r[i+2] en fonction de r[i] et r[i+2].
La troisième r[i] et r[i+4].
Bon bah en fait tu peux coder tout le parcours de l'arbre avec 2 boucles imbriquées, une globale j variant de 0 à 2 et tu calcules pour chaque j r[i] et r[i+2^j]. Non en fait, c'est un peu plus compliqué car i ne doit pas toujours varier de 0 à 4.
Bon je commence à comprendre que tu fais.
Bon bah vu l'heure bye bye.
Otlore,
Je m'en sors avec deux boucles imbriquées qui varie de 0 à 2 puis de 0 à 3. A la réflexion, c'est comme tes trois boucles.
Bon, maintenant le plus curieux, c'est que mes valeurs paires sont juste mais pas les impaires.
Ca paraît assez étonnant car le calcul des valeurs paires s'appuie sur des valeurs impaires intermédiaires. Seules les valeurs impaires finales seraient alors fausses.
Pourrais-tu vérifier avec autre chose que Mathematica ? Il faut paramétrer la fonction fft et on se trompe peut-être dans ces paramètres. J'utilise trop peu souvent Mathematica pour te dire. Il existe différentes implémentations de FFT sur le Net. Il est toujours possible de vérifier les résultats avec une autre.
www.fftw.org
Est-ce qu'un modérateur pourrait déplacer cette discussion sur le forum Algorithmes parce que ça n'a plus grand rapport avec le C ?
Pour Mathematica il donne :"In Mathematica, the discrete Fourier transform v[s] of a list u[ r] of length n is by default defined to be 1/Sqrt[n]*(Somme(u[r]e^(2*I*Pi*(r-1)(s-1)/n),r=1..n-1) "
Avec les bons paramètres {-1,-1} j'ai la bonne définition...
Et on trouve encore que les termes paires sont justes...
ça y est j'ai vu où était l'erreur!!!!!!
le wN dont on prend les puissances successives n'est pas le même à chaque étape!
dans les formules de passage de l'article (là où il y à l'arbre qui est écrit explicitement)
pour l'étape avec 4 calculs indépendants le wN=exp(2*I*Pi/N) est en fait wN^(2^(p-1))
pour l'étape avec 2 calculs indépendants le wN=exp(2*I*Pi/N) est en fait wN^(2^(p-2))
pour l'étape avec 1 calculs indépendants le wN=exp(2*I*Pi/N) est en fait wN^(2^(p-3)), pour cette étape, en remplaçant N par 8 et p par 3
on réobtient wN...
en effet les formules de passage sont du au fait que l'on trouve qu'a chaque étape les P(k) ou I(k) (qui sont différent à chaque étape et qui ne correspondent pas à une même lonqueur de liste, divisé par 2 à chaque fois) sont eux même des transformés de Fourier d'une autre liste mais d'ordre inferieur (c'est la transformé de fourier qui est d'ordre inferieur) donc les formules de passages change pour passé de liste de longueur m à 2*m
on a Y(k)=1/2*(P(k)+w2m^(-k)*I(k))
Y(k+m)=1/2*(P(k)-w2m^(-k)*I(k))
suis-je clair?
Voilà pour le code final!
Code:
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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119#include <stdio.h> #include <stdlib.h> #include <math.h> #include <complex.h> void *FFT(double complex *l,int p) { int N,i1; double complex *w,r1inter,r2inter,*r; N=2 << (p-1); r=(double complex *)malloc(N*sizeof(double complex)); w=(double complex *)malloc(N*sizeof(double complex)); w[0]=1+0*I; /*mise en place de la liste des puissances de w*/ double complex c; c=cos(2.*M_PI/N)-sin(2.*M_PI/N)*I; for(i1=0;i1<=(N-2);i1++) {w[i1+1]=w[i1]*c; } int k,w1,d,k1,z,s; /*inversions des bits*/ for(k=0;k<=(N-1) ;k++) {w1=0; z=k; s=0; d=N/2; for(k1=0; k1<=(p-1);k1++) { if (z==0){ break; } else if(z!=0) {w1=z%2; s=s+w1*d; d=d/2; } z=(z-w1)/2; } r[k]=l[s]; } /*la partie calcul proprement dite*/ { int d1,d2,i,j,b,f; d1=N/2; d2=1; for (j=1; j<=p; j++) { for(i=0; i<=(d1-1); i++) { for(b=0; b<=(d2-1); b++) { r1inter=r[i*2*d2+b]; r2inter=r[i*2*d2+b+d2]; r[i*2*d2+b]=(r1inter+r2inter*w[d1*b])*0.5; r[i*2*d2+b+d2]=(r1inter-r2inter*w[d1*b])*0.5; } } d1=d1/2; d2=d2*2; } } free(w); return(r); } int main(void) { int lambda; double complex *test; test=(double complex *)malloc(16*sizeof(double complex)); test[0]=40.; test[1]=23.; test[2]=1.; test[3]=54.; test[4]=1.; test[5]=8.; test[6]=1.; test[7]=14.; test[8]=34.; test[9]=89.; test[10]=3.; test[11]=4.; test[12]=5.; test[13]=27.; test[14]=13.; test[15]=5.; double complex *result; result=(double complex *)malloc(16*sizeof(double complex)); result=FFT(test,4); free(FFT(test,4)); for(lambda=0; lambda<=15; lambda++) {printf("%lg + %lg*i \n",creal(result[lambda]),cimag(result[lambda]));} system("pause"); free(test); free(result); return(0); }
en rouge ce qu'il fallait "tout simplement" rajouter! :bug: :lol:
J'ai regarde cet algorithme qui ameliore le calcul de la transformee de fourrier
Je vous paraphrase pardonnez moi
Les calculs au niveau i sont des sommes de calcul sur 2 sous ensemble de niveau i-1.
par exemple :
au niveau 3
(0 1 2 3 4 5 6 7)
sont calcules depuis le niveau 2 de
(0 2 4 6) et (1 3 5 7)
calcules depuis 1 de
(0 4), (2 6) et (1 5), (3 7)
A) Il faut d'un cote ... partir du niveau le plus haut ..Pour connaitre l'odre des paires
Partie permuttation
On scinde notre ensemble en 2 en gardant paires et impaires dans 2 sous ensembles disctincts
On reitere sur les ensembles obtenus
En descendant ainsi ... On obtient les permuttations de niveau 1 avec N/2 paires
Notez qu'une permutation et decoupage en 2 sous ensembles
admettent une fonction inverse celle du melange de 2 ensembles successif
Par exemple (0 4) et (2 6) niveau 1 donnent (0 2 4 6) au niveau 2
Connaissant les valeurs de (0 4), (2 6) au niveau 1 ... on calcule (0 2 4 6) au niveau 2
B) calculer depuis le niveau i pour connaitre le niveau i+1
Debut du calcul
... on calcule alors le niveau 1 sur les N/2 paires.. j'ai pas verifié... mais il doit etre tres simple.
et on remonte ..
On reassemble nos ensembles de niveau i et on calcule le resultat pour obtenir alors
le resultat sur le niveau i+1
C'est clair ... c'est pas le plus simple a programmer
Mais pour verifier qu'on est en phase .. Je ferais aussi le calcul normal en parallèle
Il est plutot simple. Cela validerait pas mal de choses.
Mais je pense que vous avez trouvé
Bravo
Felicitations... Implementer la FFT pour la premiere fois n'est jamais facile. De plus, partir directement sur la version iterative ca ne simplifie pas les choses.Citation:
Envoyé par otlore
Pour info, un petit papier explicatif avec le code C qui va bien:
http://cermics.enpc.fr/polys/info1/main/node97.html
Bon bah merci à tous pour votre aide...