FFT inverse, algorithme adapté ?

Version imprimable

06/04/2006, 21h13
progfou

FFT inverse, algorithme adapté ?

Bonsoir à tout le monde !
J'ai un algorithme de fft inverse, mais il ne me parait pas adapté à mon architecture (les multiplications mettent plus de temps que les additions, au moins 5 fois plus). Il semblerait que ce soit le radix (mais avec un programme pas commenté...). Je pex fournir un bout de code correspondant, si nécessaire, mais ce que je cherche c'est surtout tous les algos qui existent, qui font un calcul de fft inverse complexe, avec le moins de multiplications possibles...

Merci d'avance pour vos réponses !
06/04/2006, 21h46
j.p.mignot

j'uitilise la + part du temps l'algorithme de Cooley-Tukey qui est en n.log(n).
L'algorithme de Winograd [WFTA] permet de réduire sensiblement le nombre de multiplications par contre il necessite plus de mise en forme des data.

Voilà un bout de l'algo que j'ai :

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
 
        bg = fftn >> 3;                 /* 128 or 64 >> 3 */
        gp = 4;                         /* groups per pass */
 
        for (k = fftnlg2m3; k != 0; k--)
        {
            bfyrptr1 = fftrptr;
            bfyiptr1 = fftiptr;
            bfyrptr2 = fftrptr + bg;
            bfyiptr2 = fftiptr + bg;
            brxcosptr = brxcos;
            brxsinptr = brxsin;
 
            for (j=gp ; j!=0 ; j++) {
                cr = *brxcosptr;
                brxcosptr++;
                ci = *brxsinptr;
                brxsinptr++;
 
                for (i = bg; i !=0 ; i--)
                {
                    ar = *bfyrptr1;
                    ai = *bfyiptr1;
                    br = *bfyrptr2;
                    bi = *bfyiptr2;
 
                    rtemp = br*cr - bi*ci;
                    itemp = br*ci + bi*cr;
 
                    *bfyrptr1 = (ar - rtemp);
                    bfyrptr1++;
                    *bfyiptr1 = (ai - itemp);
                    bfyiptr1++;
                    *bfyrptr2 = (ar + rtemp);
                    bfyrptr2++;
                    *bfyiptr2 = (ai + itemp);
                    bfyiptr2++;
                }
                bfyrptr1 += bg;
                bfyiptr1 += bg;
                bfyrptr2 += bg;
                bfyiptr2 += bg;
            }
            bg >>= 1;
            gp <<= 1;
        }

En voyant bien que les brxcos/sin sont des pointeurs sur des tableaux de cos/sin précalculés...

Mon problème ne situe pas dans la compréhension du code (ce ne serait pas le bon forum :P) mais dans l'amélioration.
Je ne connais pas Winograd, mais je vais me renseigner. Les données ont vraiment besoin de beaucoup de mise en forme ?

Quelle architecture utilises-tu ? (Moi je connais que Intel)

Ton algorithme utilise de toute évidence un radix de 2.
A titre de comparaison, tu trouveras ci-après mon implémentation en C++ de la FFT à radix 2, facilement transposable ce C.

Un truc qui pourrait améliorer la vitesse de ton implémentation :
Faire des vecteurs de complexes (RIRIRI…) plutôt que 2 vecteurs de réels (RRR…, III…)
=> moins de pointeurs à gérer (2 au lieux de 4, ça peut faire une grande différence sur un Pentium), peut-être moins d’accès mémoire…

Pour des implémentations beaucoup plus performantes, il faut augmenter le radix, et utiliser les instructions parallèles (SSE, SSE2, SSE3 pour les Pentiums).
Mon implémentation perso:
=> www.ient.rwth-aachen.de/team/laurent/genial/genial.html
Mais FFTW est plus connue :
=> http://www.fftw.org
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 template<class G> Vector<G> &complex_radix2_fft_in_place(Vector<G> &X) { typedef Vector<G> array_type; typedef typename array_type::size_type size_type; typedef PROMOTE2(complex<float>,typename array_type::value_type) value_type; typedef PROMOTE2(complex<float>,typename value_type_traits<value_type>::value_type) complex_type; typedef typename complex_type::value_type real_type; size_type n = X.size(); assert(log2n(n)>=0); bitreverse_order(X); for (int dual=1; dual<n; dual*=2) { real_type theta = -PI/dual; real_type s = sin(theta); real_type s2 = real_type(2)*sqr(sin(real_type(0.5)*theta)); complex_type w(1,0); value_type wd; for (size_t a=0; a<dual; ++a, w+=s*complex_type(-w.imag(),w.real())-s2*w) for (size_t b=0; b<n; b+=2*dual) { value_type &Ti=X[b+a]; value_type &Tj=X[b+a+dual]; wd=w*Tj; Tj=Ti-wd; Ti+=wd; } } return X; }

07/04/2006, 22h42
progfou

A priori, recherches à l'appui, l'implémentation serait split-radix.
L'architecture est un Leon2 (SparcV8).
Je vois qu'il existe aussi une implémentation du split-radix à 3 mult 3 add au lieu de 4 mult 2 add...
07/04/2006, 23h29
Charlemagne

Ton "split-radix", c'est peut-être ce que j'appelle "avec radix mixte".
Ca consiste à décomposer la longueur N en un produit de facteurs.
Le radix-2 ne permet que des puissances de 2.
Alors qu'avec un radix mixte permet bien d'autres longueurs, du moment qu'elle se laisse décomposer.
FFTW et mon implémentation de la FFT utilisent des facteurs allant jusqu'à 64.

Deux avantages majeurs pour la vitesse:
1) moins d'opérations
2) meilleure utilisation des registres => moins d'acces mémoire (goulot d'étranglement, en tout cas pour les Pentiums)

Inconvénient: plus compliqué

J'ai pas réfléchi au problème, mais je doute qu'il soit possible de réduire le nombre de multiplications d'une autre façon (si vraiment les Sparcs sont lents par rapport à l'addition). Il existe une limite inférieure théorique dont je ne connais pas la formule.

Si tu recherches la vitesse, si tu n'est pas obligé d'implémenter, et si la licence te le permet alors FFTW se laisse très probablement compiler sur un Sparc.

Ca existe les instructions SIMD (calcul parallèle) sur les Sparcs? Très probablement.
09/04/2006, 11h21
progfou

Hors de question d'utiliser FFTW, question de licence comme tu l'as si justement soupçonné :wink: .
Je peux m'en inspirer, et c'est ce que je fais ^^.
Maintenant, pour la question de l'existence d'instructions SIMD sur le Sparc, je ne crois pas, mais je ne le maîtrise pas encore assez pour répondre...
Je suis en plein dans la doc ^^.
Pour l'algo qui est implémenté, il y a un radix 4 puis un radix 2.
L'idée de mettre les réels et imaginaires dans le même tableau est à creuser, merci pour cette piste ;).
09/04/2006, 17h08
progfou
J'ai trouvé dans le bouquin "The Art of Computer Programming" volume 2 une solution pour faire une multiplication complexe avec seulement 3 multiplications et 5 additions :
Code:

1 2 3 (a+bi)(c+di) =a(c+d)-(a+b)d+i[a(c+d)+(b-a)c]
Ca a le mérite d'être facile à montrer...