13
pp. 495-507 495 Jean-Luc ZOLESIO ** Transformations de Fourier discr. tes i algorithmes rap des Analyse Cet article pr~sente une revue, non exhaustive, des algorithmes de transformations de Fourier discrbtes : les algorithmes d base premibre rdcursive, d base 2k, division non rdgulibre, les algorithmes courts, facteurs premiers et ceux de Winograd. Pour chacun d'eux, la technique de construction de l'algorithme est r~sumde en ddgageant les iddes essentielles ou dventuellement les structures algdbriques sous-jacentes. De plus les caractdristiques principales des algorithmes sont mentionndes en essayant de montrer que les critbres de choix d'un algorithme rapide d~pendent ~troitement de la far de le mettre en ceuvre. En particulier les charges de calcul respectives des algorithmes sont compardes. Mots cl~s : Transformation Fourier discr6te, Transfor- mation Fourier rapide, Algorithme, Etude comparative, Article synth6se, Complexit6 algorithme, Op6ration arithm&ique. DISCRET FOURIER TRANSFORMS ; FAST ALGORITHMS Abstract This paper reviews some discret Fourier transform algorithms : the algorithms with prime recursive radix, 2k radix, non regular division, short algo- rithms, prime factor and Winograd ones are consi- dered. For each of them, the technique to built the algorithm is resumed and algebraic properties on basic ideas used are pointed out. The most important proper- ties of the algorithms are shown underlying the links between criteria choice and implementation of algo- rithms. In particular, the computational load of the algorithms is analysed. Key words : Discrete Fourier transformation, Fast Fourier transformation, Algorithm, Comparative study, Review, Algo- rithm complexity, Arithmetic op6ration. Sommaire I. Introduction. II. La charge de calcul. III. Transformation de Fourier discrbte d base premibre rdcursive. IV. Les algorithmes de base 2~. V. Les algorithmes ~ division non rdgulibre. VI. Les algorithmes courts. VII. Les algorithmes ~ facteurs premiers. VIII. Conclusion. Bibliographic (17 rdf.). I. INTRODUCTION I1 y a une vingtaine d'ann6es, l'apparition des techniques num6riques en traitement du signal a suscit6 un int6r& croissant pour la transformation de Fourier discr6te (TFD). Suivant les applications, deux situations ont 6t6 rencontr6es : le calcul en temps r6el avec des cadences de renouvellement des donn6es tr~s 61ev6es, ou le calcul en temps diff6r6 et la simula- tion portant sur des taffies d'6chantillons tr6s grandes. * C.ette 6tude a pu 6tre men6e h bien en partie grace au soutien de la Direction des recherches et 6tudes techniques, I~16gation g6n6rale pour l'armement, clue nous remercions pour son autorisation de publication. ** Laboratoire central de t616communication, 18-20, rue Grange Dame-Rose, B.P. 40, 78141 Velizy-Villacoublay. 1/13 ANN. T~Li~COMMUN., 40, n ~ 9-10, 1985

Transformations de Fourier discrètes: algorithmes rapides

Embed Size (px)

Citation preview

Page 1: Transformations de Fourier discrètes: algorithmes rapides

pp. 495-507 495

Jean-Luc Z O L E S I O **

Transformations de Fourier discr. tes i algorithmes rap des

Analyse

Cet article pr~sente une revue, non exhaustive, des algorithmes de transformations de Fourier discrbtes : les algorithmes d base premibre rdcursive, d base 2 k,

division non rdgulibre, les algorithmes courts, facteurs premiers et ceux de Winograd. Pour chacun d'eux, la technique de construction de l'algorithme est r~sumde en ddgageant les iddes essentielles ou dventuellement les structures algdbriques sous-jacentes. De plus les caractdristiques principales des algorithmes sont mentionndes en essayant de montrer que les critbres de choix d'un algorithme rapide d~pendent ~troitement de la far de le mettre en ceuvre. En particulier les charges de calcul respectives des algorithmes sont compardes.

Mots cl~s : Transformation Fourier discr6te, Transfor- mation Fourier rapide, Algorithme, Etude comparative, Article synth6se, Complexit6 algorithme, Op6ration arithm&ique.

D I S C R E T F O U R I E R T R A N S F O R M S ; F A S T A L G O R I T H M S

Abstract

This paper reviews some discret Fourier transform algorithms : the algorithms with prime recursive radix, 2 k radix, non regular division, short algo- rithms, prime factor and Winograd ones are consi- dered. For each of them, the technique to built the algorithm is resumed and algebraic properties on basic ideas used are pointed out. The most important proper- ties of the algorithms are shown underlying the links

between criteria choice and implementation of algo- rithms. In particular, the computational load of the algorithms is analysed.

Key words : Discrete Fourier transformation, Fast Fourier transformation, Algorithm, Comparative study, Review, Algo- rithm complexity, Arithmetic op6ration.

Sommaire

I. Introduction.

II. La charge de calcul.

III . Transformation de Fourier discrbte d base premibre rdcursive.

IV. Les algorithmes de base 2 ~.

V. Les algorithmes ~ division non rdgulibre.

VI. Les algorithmes courts.

VII. Les algorithmes ~ facteurs premiers.

VIII . Conclusion.

Bibliographic (17 rdf.).

I. I N T R O D U C T I O N

I1 y a une vingtaine d'ann6es, l ' appar i t ion des techniques num6riques en traitement du signal a suscit6 un int6r& croissant pour la t ransformat ion de Fourier discr6te (TFD). Suivant les applications, deux situations ont 6t6 rencontr6es : le calcul en temps r6el avec des cadences de renouvellement des donn6es tr~s 61ev6es, ou le calcul en temps diff6r6 et la simula- t ion por tan t sur des taffies d '6chantil lons tr6s grandes.

* C.ette 6tude a pu 6tre men6e h bien en partie grace au soutien de la Direction des recherches et 6tudes techniques, I~16gation g6n6rale pour l'armement, clue nous remercions pour son autorisation de publication. ** Laboratoire central de t616communication, 18-20, rue Grange Dame-Rose, B.P. 40, 78141 Velizy-Villacoublay.

1/13 ANN. T~Li~COMMUN., 40, n ~ 9-10, 1985

Page 2: Transformations de Fourier discrètes: algorithmes rapides

496 J . -L . ZOLESIO. - TRANSFORMATIONS DE FOURIER DISCRETES

Dans les deux cas, on a cherch6 h r6duire le temps de calcul des TFD en construisant des algorithmes rapides dont la charge de caleul, c'est-/t-dire le volume d'op6rations arithm6tiques n6cessaires, est r6duite.

Le but de cet article est de pr6senter une synth6se de quelques algorithmes de TFD rapides en essayant d'en d6gager les id6es qui ont conduit/ t leur 61abora- tion et en 6valuant leur charge de calcul et leur effica- cit6 respective.

Le point de vue adopt6 est relatif h la logique cabl6e, technique particuli6rement performante pour les calculs en temps r6el /t tr6s grandes cadences, ofa la charge de calcul, si elle n'est pas le seul crit6re de choix d 'un algorithme, d&ermine n6anmoins pour une part importante ce choix. (I1 est/L noter que lorsqu'il s'agit de sa programmation, la rapidit6 d 'un algorithme n'est pas une notion intrins~que et d6pend en fait plus du mat6riel utilis6 que de l'algorithme lui-mSme [1].)

Apr6s avoir d6fini pr6cis6ment ce que l 'on entend par charge de calcul et apr6s avoir pr6sent6 un algo- rithme rapide bien que pourtant de charge de calcul identique/l celle du calcul direct d'une TFD, on analyse les quatre families d'algorithmes suivantes :

- - l e s algorithmes de base 2 r,

- - les algorithmes /l division non r6guli~re,

les algorithmes courts (ou diagonaux),

les algorithmes /t facteurs premiers.

H. LA CHARGE DE CALCUL

L'eflicacit6 d 'un algorithme rapide d6pend de plu- sieurs crit6res dont l 'ordre d'importance r6sulte du contexte de sa mise en oeuvre.

Sans dresser une liste exhaustive, on peut 6num6rer les plus importants : la charge de calcul, l 'agencement temporel des donn6es tout au long du calcul, la nature des op6rations/t r6aliser, la masse de m6moire n6ees- saire aux stockages interm6diaires, la robustesse au bruit de calcul, la r6p6titivit6 et la modularit6, etc.

I1 serait done peu r6aliste d'6tablir une hi6rarchie entre les diff6rents algorithmes h partir seulement de leur charge de calcul respective bien que ce soit eependant un crit6re important de choix, particuli6re- ment en logique cabl6e.

Or, comme on le verra par la suite, la charge de calcul de certains algorithmes est essentiellement due /t des multiplications de variables complexes par des donn6es r6elies ou imaginaires pures, alors que darts d'autres, elle est due/ l des multiplications de nombres complexes entre eux.

Afin de pouvoir comparer objectivement ces diff6- rentes techniques sur ce point, on est done amen6 /t d6composer les multiplications entre nombres com- plexes en op6rations arithm6tiques entre hombres r6els,

op6rations qui seront alors appel6es 616mentaires.

Consid6rons le produit : P = (a + bj) (x + y j) oh a et b sont donn6s, d'apr~s une remarque due /t Winograd [2], on a :

P = [a(x + y ) - - ( a + b)y] + j[a(x + y ) - - ( a - - b ) x ] ,

trois multiplications r6elles et trois additions r6elles sont n6cessaires.

Dans le cas d 'un processeur en logique c~tbl6e, on peut retenir que la charge de calcul provoqu6e par une multiplication complexe est conforme ~t cette analyse (une multiplication par • j ne repr6sente qu'un efiblage diff6rent et pas une op6ration suppl6- mentaire).

HI. TRANSFORMATION DE FOURIER DISCRETE

A BASE PREMI]gRE RI~CURSIVE

Cet algorithme, mis en oeuvre par T. E. Curtis [3], est un exemple int6ressant mettant en 6vidence l 'impor- tance de la nature m~me des multiplications ~ r6aliser. En effet, bien que sa charge de calcul soit strictement identique/t celle r6sultant du calcul direct d'une TFD, l'algorithme /t base premi6re r6cursive, est cependant beaucoup plus rapide que la TFD.

En notant, x = (Xo, x l , ..., xN-1) r u n vecteur complexe de dimension N, l'indice T d6signant la transposition, le calcul direct de la TFO A du vecteur a, est le suivant :

N - 1

(1) AK = ~ a l W~ ~, K = 0 , . . . , N - - I / = 0

oh W~ = exp ( - -2~ j l /N ) , t ~ Z , j2 = __ 1.

(lorsque toute ambiguR6 sera lev6e, l'indice N de W~ sera omis).

L'expression (1) peut encore s'6crire sous forme matricielle comme suit �9

(2) A = wN a,

oh wN est la matrice de la TFD de dimension N. L'id6e originale de cet algorithme est d'6crire (1)

pour K non nul sous la forme �9

(3) A = or(a) w,

ofa e(a) est une matrice obtenue/ t partir de permuta- tions des donn6es et o/a w = (1, w, ..., wrY-x) r.

Cette 6criture n'est possible que lorsque N e s t premier. En effet, dans ce cas, dans l'expression (1), pour chaque valeur de K non nulle, toutes les compo- santes de W sont utilis6es pour le calcul de AK.

On peut done 6crire : N - 1

(4) A K = ~ alK* W ~; 1 ~ < K ~ < N - - 1 , l=O

avee K* d6fini par : KK* = 1 [mod N].

ANN. T~L~COMMUN., 40, n ~ 9-10, 1985 2/13

Page 3: Transformations de Fourier discrètes: algorithmes rapides

J.-L. ZOLESIO. - TRANSFORMATIONS DE FOURIER DISCRETES 497

En posant a~ = a~r* et en utilisant Ie sch6ma de Horner pour le calcul de (4), on a :

t (5) A k = ( . . . ( (a '~_lW+a'N_2) W + a ~ _ a ) . . . +ao),

pour K = 1 . . . . , N - - 1, et la m~me expression pour K nul dans laquelle Wes t remplac6 par 1.

On voit donc que la charge de calcul de cet algo- rithme est strictement la m~me que celle r6sultante du calcul de (1). Mais une mise en oeuvre c~bl6e astucieuse (voir [3]) de cet algorithme conduit h u n calcul de la TFD assez rapide pour certaines appli- cations. En effet, les multiplications intervenant dans (5) ont toutes le m~me multiplicateur W ; on peut donc calculer (5) de fa~on r6cursive en utilisant une m6moire contenant x W ~ l 'adresse x.

De cette fagon, le calcul de (5) en pratique ne fait pas intervenir d 'op6rateur de multiplication et se ram6ne ~t une succession de lectures m6moire et d 'additions ; ce qui conduit, de plus, ~t des mat6riels de bas prix.

IV. LES A L G O R I T H M E S D E B A S E 2 r

Ces algorithmes sont bas6s sur les sym6tries de rotations qui existent entre les bases de l'unit6 lorsque la dimension de la TFD est : N = 2 ~'.

Dans ce cas, on peut d6composer la TFD en 2 K TFD de dimension 2 v-r .

En effet, en posant M = 2 r, (1) peut aussi s'6crire :

M--1 N I M - I Ar = ~ ~ aM.+m W ~M~+m~r ; K = 0 . . . . . N - - l ,

m=O n = O

soit encore : M - - 1

(6) A K = ~ W m x T ~ ; K = 0 , . . . , N - - 1 , m = O

avec : N/M- 1

(7) r ~ = Z aMn+m WMnK, n=O

ce qui n 'est autre qu 'une TFD de dimension N [ M car :

W~I I M .

En r6it6rant l 'op6ration, on se ram6ne A ne calculer que des TFD de dimension 2 r.

I V . 1 . K = 1 : l ' a l g o r i t h m e de b a s e 2 .

Mis au point par J . W . Cooley et J . W . Tukey en 1965 [4], cet algorithme, connu sous le nom d'algo- rithme de transformation de Fourier rapide (FFT en anglais), a 6t6 le pr6curseur de l 'arsenal des m6thodes de calcul de TFD actuellement ~t notre disposition. I1 est fr6quemment utilis6 en raison de sa simplicit6 d'adressage.

En tenant compte des propri6t6s de sym6trie des

3/13

bases de l'unit6, (6) se d6compose dans ce cas comme suit :

l A t = T O + WKT1K,

(8) Autz+r T ~ - - W K T~,

oh les T~, d6finis par (7), sont des TFD de dimension NI2.

La charge de calcul de cet algorithme est la sui- vante :

- - nombre de multiplications 616mentaires :

- - nombre d'additions 616mentaires :

N ( ~ p - - 5 ) + 8 , of a N = 2".

La figure 1 pr6sente h titre d 'exemple le graphe de circulation d 'une transformation de Fourier rapide de dimension 8.

FIG. 1. - - Graphe de circulation des donn6es de l 'algorithme h base 2 pour N = 8.

Signal flow graph of the radix - 2 algorithm, N = 8 .

Un algorithme dquivalcnt h (8) peut ~tre obtenu en regroupant les termes comme suit :

NI2 -- 1 NI2 -- 1 AK"~'~ ~ an w n K 21 - ~ an+N~2 W (n+NI2)K,

(9) , = o , = o

K = 0, ..., N - - 1.

En tenant alors compte des 6galit6s :

W 'v = 1 et W m2 = - - 1,

on a encore :

i NI2 - 1 A2K = .=Z ~ (a. + a.+m2) W e"K,

(10) Nt2-~ A2K+I= ~ ( a . - - a . + n n ) W" W 2"r,

\ u = O

K = 0 . . . . , N I 2 - 1.

ANN. T[L~COMMON., 40, n ~ 9-10, 1985

Page 4: Transformations de Fourier discrètes: algorithmes rapides

498 J.-L. ZOLESIO. - TRANSFORMATIONS DE FOURIER DISCRETES

Comme dans l 'approche pr6c6dente, iI est possible de calculer une TFD par 6tapes successives, chacune d'elles permettant de diviser par deux la dimension des transformations.

IV.2. Algorithme de base 2 /t coefficients r6els.

Cet algorithme, du ~t C. M. Radar et N. M. Bren- ner [5] est d6riv6 de (10).

En consid&ant une suite auxiliaire d6finie par :

I x, = (a, - - a,+m2)/[2 cos(2rcn]N)], x , = 0 s i n t = 0 , N/4,

et en calculant la TFD de dimension N]2 de cette suite �9

NI2 - 1

XK = Z x. W 2"K, 11=0

on peut exprimer A2K+~ dans (10) en fonction de X~: ; on a :

A2K+I = XK + XK+t + Vo, pour K pair, (11)

A2K+I XK + XK+I + va, pour K impair,

avec vi = (ao - - a l v l 2 ) - - (-- 1) l (as14 - - aam,). Les N]2 multiplications par W" clans (10) ont 6t6

remplac6es par (N[2 - - 2) multiplications par les r6els :

[2 cos (2 ~ I N ) ] - '.

Ceci permet de ramener le nombre de multiplications 616mentaires h : N (p - - 3) + 4, ce qui repr6sente un gain asymptotique de 33 % par rapport b, l'algo- rithme de base 2 pour un nombre d'additions supdrieur d'environ 10 ~o.

Le tableau I permet d'appr6cier le gain en charge de calcul pour des petites dimensions par rapport h l'algorithme de base 2.

TABL. I. - - Gain en charge de calcul de l'algorithme de base 2 ~t coefficients r~els.

Dimension de la TFD

Gain en multipli- cations (%)

Gain en additions (%)

[ 8 t 16 32 64

i ]

0__0_[ 16 22 25

0 [ - - 2 - - 4 - - 7 - - 9 - - 1 0

Cependant, si cette r6duction de la charge de calcul est int6ressante, il faut rioter que cet algorithme g6n6re un bruit de calcul important lorsque N e s t grand. En effet, dans ce cas cos(2= n[N) prend des valeurs tr6s petites pour certaines vateurs de n e t les multiplications par [2 cos(2 rc n/N)]- i introduisent alors des erreurs importantes.

IV.3. K = 2 : l 'algorithme de base 4.

En reprenant (6) dans ce cas particulier et en tenant compte des 6galit6s suivantes :

Win(K+ INI4) = [ (__ j)l]m w m K ,

m = 0 , . . . , 3 et l = 0 , . . . , 3 ,

on obtient :

t Ax+lNt, = (T ~ + (-- 1)' W:KT~ + (12) ( - - j ) ' (WKTI + (-- 1) ~ W ar T~),

K = 0 . . . . , N [ 4 - - 1 et l = 0 , . . . , 3 .

Comme pr6c6demment, on peut r6it6rer l'op6ration jusqu'~t obtenir une d6composition en TFD de dimen- sion 4 si N &ait une puissance de 4. Lorsque ce n'est pas le cas, h l'avant-derni6re 6tape, on se trouve en pr6sence de TFD de dimension 8 que l 'on d6compose eft deux TFD de dimension 4 suivant l'algorithme de racine 2.

Cet algorithme est trSs utilis6, car il est d'adressage relativement simple, tr6s r6p&itif et poss6de une charge de calcul sensiblement du m~me ordre que celui de base 2 ~t coefficients r6els, sans en avoir les inconv6nients.

La charge de calcul est la suivante :

- - nombre de multiplications 616mentaires :

N(918 p - b) + 16/3,

- - n o m b r e d'additions 616mentaires :

N(25/8 p - - b) + 16/3,

off b e s t 6gal respectivement ~t 43/12 ou fi 85/24 suivant que p e s t pair ou non.

Le gain asymptotique par rapport fi l'algorithme de base 2 du paragraphe IV.1. est le suivant :

�9 25 % de multiplications 616mentaires,

�9 10% d'additions 616mentaires.

(On pourra remarquer que l'algorithme de base 2 ~t coefficients r6els repr6sentait un << gain >> relatif en additions de - - 10% contre + 10% pour cet algorithme.)

Le tableau II permet d'appr6cier le gain en charge de calcul pour des petites dimensions par rapport ~t l'algorithme de base 2.

TABL. II. -- Gain en charge de calcul de l'algorithme de base 4.

Dimension de la TFD [ 81 161 321 64 11281 256 I

Gain en multiplications (%)

Gain en additions (%)

Un exemple de graphe de circulation des donn6es de cet algorithme est donn6 par la figure 2 alors que

ANN. T~LigCOMMtJN., 40, n ~ 9-10, 1985 4/13

Page 5: Transformations de Fourier discrètes: algorithmes rapides

J.-L. ZOLESIO. - TRANSFORMATIONS DE FOURIER DISCRETES 499

M o d u l e a

13~_ .~, \ / \ _ ~ - ~ ,*.

,, X X \ / / / . \ ; ~ , " [ J " ~ , - . - j ~,~ ,, \ \ \ V / / / ~ - - A~

- ~ \ \ X X / / / / X X X _ w~ V ~ " ,~ \ X X X X O ( / X / / ~ " ~ =~ X X - -A,0

_- _- - A14

~ / -- \ / -A, ~, XOO<XXXX- \ / w' ~<- =A, , , o / / " X ) ( ~ \ " < \ \ / / w~ 50( -

, A~a

, , , / / / / ~ \ \ " < -J X X X ~ . A. , . 1 / / \ \ " ~ - J / X X ~ ' ~ w" [ . ~ / X . A,,

, , , / , % - ~ / % v~ I /%_-iX_. ' " : A1s

Moduh~ b

FIG. 2. - - Graphe de circulation des donn6es de l'algorithme ~t base 4 pour N = 16.

Le module a est h double base et le module best h base 4.

Signal flow graph o f the radix - 4 algorithm, N = 16.

aa a12

a o ~ Aa --All A7 A~ 5 (11

ao I . . . . . . . I a~ I : .~,/" a~ I ~ • I Aa iX--J X I A~ a~ I . . . . [ Aa

Fie. 3. - - Structure modulaire de l'algorithme h base 4 pour N = 16.

Modularity o f the radix - 4 algorithm, N ~ 16.

la figure 3 mont re la s tructure tr6s r6p6titive et modu- laire, caract6ristiques des algori thmes de base 2 k.

Ces propri6t6s permet ten t des mises en oeuvres en cascade (<< pipe-line >> en anglais) tr6s ais&s con- duisant ~t des mat6ricls per formants .

A titre d 'exemple , Ia figure 4 mont re un analyseur de spectre effcctuant la pondf ra t ion , la TFD calcul6e par l ' a lgor i thme de base 4 et le calcul du module . Ce processeur cgtbl6, d6vclopp6 cn technologie stan- dard pour un 6quipemcnt radar par le Labora to i re central de t616communication en 1980, pe rmet d 'effec- tuer une TFD c o m p k x e de dimension 16 en 4,5 ~s, les temps d 'entr6e et de sortie des donn6es 6tant compris . Sa dynamique maximale est sup6rieure ~t 80 dB et sa capacit6 de t ra i tement cst de 5 • 10 5 spcctres pa r seconde.

Ces algori thmes, don t deux ~xcmplcs K = 1 el K = 2 ont 6t6 d6taill6s (pour K = 3 voir [6]), ont l ' avan tage d 'a l l ier h une charge de calcul tr6s int6- ressante, une structure modulai re ct tr6s r@6titive. Les cas K = 1 et K = 2 sont les plus utilis6s car ils utilisent d . s m o d u k s de base (cncore app~16s des << papil lons >>) simples ne faisant in t : rv .n i r aucunc mult ipl icat ion (voir Fig. 3).

V. L E S A L G O R I T H M E S A D I V I S I O N N O N R I ~ G U L I E R E

Les algori thmes consid6r6s au paragraphe IV con- duisent ~t chaque 6tape h diviser une TFD , n TFO de m~me dim'_nsion 2 k. On peut c :pend:mt ~nvisagcr

FIG. 4. - - Analyseur de spectre du Laboratoire central de t616communications.

A << Laboratoire central de t616communication >> spectrum analyzer.

5/13 ANN. TI~LI~COMMUN., 40, n ~ 9-10, 1985

Page 6: Transformations de Fourier discrètes: algorithmes rapides

5 0 0 J.-L. ZOLESIO. - TRANSFORMATIONS DE FOURIER DISCRETES

de diviser une TFD en plusieurs parties, chacune d'elles &ant d 'un calcul ais6 et ayant sa propre dimension. C'est le cas des algorithmes pr6sent& ci-apr6s qui conduisent tousles deux A la m~me charge de calcul.

V.1. TFD /t double base.

Cet algorithme a 6t6 d6velopp6 en 1984 par Duhamel et Hollmann [7], la remarque de base &ant qu'une base 4 est plus int6ressante pour les termes impairs de (10).

En effet, en reprenant l'expression de A2K+I de (10) et en posant K = 2 p, on a :

N I 2 - 1 np A4o+l = 2~ ( a n - a.+Nn) W~ Wfv/4,

n = 0

o r s i n = l + N ] 4 a v e c l = 0 . . . . , N ] 4 - - 1, o n a : nO w; w;., j w~ ' "

d'ofi finalement :

N I 4 - 1

A,,~+I = Z [ ( a l - - al+N/2) + 1 = O

j (a ,+m, - - a,+am,)] W~ Wk/,'P,

p = O . . . . , N ] 4 - - 1 ,

de m~me en prenant les termes tels que K = 2 p + 1, on obtient pour p = 0, ..., NI4 - - 1 :

N I 4 - 1

A,,+3 = 5] [ (a l - - al+Nn) - - j(al+NI,--a~+~tq,) • /=0

W3NI lp W]~14 �9

On peut donc ramener le calcul de (I0) au calcul de :

�9 A2K ; K = O, ..., N / 2 - 1,

�9 A4v+a et A,v+a ; p = 0 . . . . . N / 4 - - 1,

soit le calcul d 'une TFD de dimension N]2 et celui de deux TFD de dimension N]4. La charge de calcul en fonction de N = 2 ves t alors la suivante :

- - n o m b r e de multiplications 616mentaires :

N(p - - 3) + 4,

- - n o m b r e d'additions 616mentaires :

3 N ( p - - 1) + 4 .

cette charge de calcul est 16g~rement inf6rieure celle relative h l'algorithme de base 4.

Les gains asymptotiques par rapport respectivement aux algorithmes de bases 2 et 4 sont les suivants :

�9 33 % et 11 ~/o de multiplications 616mentaires,

�9 14~o et 40/o d'additions 616mentaires.

Cet algorithme est de charge de calcul inf6rieure h celle de tous les autres algorithmes connus h ce jour, pour des dimensions 2 ~, ex aequo avec l'algo-

rithme d6crit au paragraphe V.2. et l'algorithme de Martens [8] qui est de m~me structure mais de cons- truction moins 6vidente.

Pour la dimension N = 16, le graphe de circulation des donn6es est identique ~t celui de l'algorithme base 4 de la figure 2 ~ la position pr6s des multipli- cations. Le module (ou << papillon >>) de base est enca- dr6 en haut de la figure 2. Sa dissym&rie implique d'utiliser des modules ~. deux points ~ la derni6re 6tape ce qui restreint 16g6rement le caract6re r6p6titif de cet algorithme et peut &re g~nant dans des r6alisa- tions en logique c~bl6e ou microprogramm6e.

V.2. TFD par transformations trigonom6triques dis- cr6tes.

Mis au point par Vetterli et Nussbaumer en 1984 [9], cet algorithme utilise essentiellement la d6composition naturelle en parties r6elles et imaginaires de la TFD ainsi que la formule trigonom&rique du cosinus de la somme de deux angles.

En effet, en consid6rant la partie r&lle de (I) et en s6parant les termes d'indices respectivement pairs et impairs, on a :

X a . cos 2zc = Z a2vcos + ,=o v=o \ N/2 /

m*-I (2 n ( 2 p + l ) K ) X (a~p+~ + a~_(~p+.) cos

v= o 4 N[4 '

le calcul de la partie r6elle d 'une TFD de dimension N peut donc se ramener h celui de la partie r6elle d 'une TFD de dimension N]2 et t~ celui d'une transformation appel6e transformation en cosinus discrete (TCD) de dimension N[4.

De m~me pour la partie imaginaire, on a :

Y~ a. sin 2 n = ~ a2v sin + . = o o=o \ N / 2 /

Y. (a2v+l- -an_(2 ,§ sin ~_(2p + 1)K v=o 4 N ] 4

en tenant compte de l'6galit6 suivante :

s i n ( 2 ~ ( 2 p + l ) K )

= ( - - 1) ~ cos - .

On se ram6ne encore au calcul de la partie imaginaire d'une TFD de dimension N[2 et d'une TCD de dimension N]4 d'indice de sortie N]4 - - K, ce qui repr&ente une permutation des r6sultats. Or le calcul d'une TCD de dimension N]4 se ram6ne ~. celui des parties r6elle et imaginaire d'une TFD de dimension N]4, en effet, en consid6rant le changement de variable suivant :

t t xl = x2z et xu-z-1 = x21+l , l = 0, ..., N]8 - - 1,

o n a :

ANN. T~L~CO~IUN., 40, n ~ 9-10, 1985 6/13

Page 7: Transformations de Fourier discrètes: algorithmes rapides

J.-L. ZOLESIO. - TRANSFORMATIONS DE FOURIER DISCRETES 501

u/,~-i / 2 ~ ( 2 l + 1 ) K \ N ,~-~ , / 2 ~ ( 4 l + 1 ) K \ _ - ,cos\ ) ~g x, cos~ 4N]4 ) z=o

/ = 0

= c o s Z x', c o s - - 1=o \ N/4 /

sin ~] x~ sin . ~=o \ N[4 /

A chaque 6tape de l'algorithme le calcul des parties respectivement r6elle et imaginaire d'une TFO de dimension N e s t remplac6 par le calcul des parties r6elle et imaginaire d'une TFO de dimension N]2 et de la partie r6elle et imaginaire d'une TFO de dimen- sion N]4. Cet algorithme a la m~me charge de calcul que l'algorithme h double base du paragraphe pr6c6dent, mais poss~de une structure beaucoup moins r6guli6re ~ cause des permutations n6cessaires d chaque ~tape de calcul.

Enfin, pour des donn6es r6elles, la charge de calcul de ces deux algorithmes est exactement r6duite de moiti6.

VI. LES ALGORITHMES COURTS

L'id6e, g6n6ralement attribu6e h Rader [10], d'6va- luer une TFO ~t partir d'une convolution circulaire conjugude avec de nombreux travaux notamment de Winograd [2] concernant le calcul rapide de telles convolutions, a permis de construire des algorithmes tr6s performants pour des TFD de petites dimensions : les algorithmes ~ courts ~).

Afin d'en comprendre Ia construction, iI est n6ces- saire d'6tudier les diff6rentes 6tapes n6cessaires h leur 61aboration ; c'est l'objet des paragraphes V.1. et V.2. (on pourra voir [11] pour un expos6 plus d6taill6).

VI.1. TFD et convolution.

Dans certains cas, une TFD peut s'exprimer ~t partir d'une convolution circulaire gr~tce h des permutations effectu6es sur une partie de la matrice de la TFD. En effet, en permutant des lignes et des colonnes dans une sous-matrice de celle de la TFO, le calcul de celle-ci revient essentiellement au calcul d'une convolution circulaire ; h une renum6rotation pr6s des 6chantillons de sortie.

Pour Npremier, consid6rons la matrice du syst6me :

N - 1

(13) Ar = ~ a, W ~K, K = 1, ..., N - - 1 . t 1 = 1

Soit G le groupe multiplicatif des entiers (1, 2 . . . . . N - 1) muni de la multiplication modulo N ; G est un groupe cyclique. Soit alors une racine pri- mitive de l'unit6 d'ordre N - 1 de G, c'est-~t-dire

~ G , ~K ~ 1 p o u r 0 < K < N - - 1 et ~N-1 = 1. I1 existe un isomorphisme de G sur le groupe additif

ZN-1 d6fini comme suit :

(14) G -+ Zs_ 1 ; n -+ m tel que ~m = n.

En effectuant le changement de variables (14) dans (13), on a :

N - - 2

A,s = Y, am W *t~+m~, s = 0, . . . , N - - 2. m = O

En posant alors :

a" = a~_,. et b~_,. = W ~*-"~,

(ce changement de variable revient ~t inverser l'ordre des 6chantillons d'entr6e hormis ao, le vecteur d'entr6e sera donc: a~o, a t r N - - 2 , acrN-3 , ...,a~2 , aal), on a :

N - - 2 I A ~ = Y, am b,-m , S = O, ..., N - - 2.

r a = O

Done ~t une permutation pr6s, le calcul de -4K est celui d'une convolution circulaire.

Or, comme on a :

Ao = (ao + a~ + ... + aN-0,

A K = a o + A K , K = 1,..., N - - 1 .

Le calcul de la TFD peut s'effectuer par le calcul d'une convolution circulaire de dimension N J 1, d'oia l'int6r~t de savoir 6valuer rapidement de telles convo- lutions.

VI.2. Convolution circulaire rapide.

Le calcul d'une TFD pouvant s'effectuer ~ partir de celui d 'une convolution circulaire, on est amen6 ~t construire des algorithmes rapides pour le calcul de telles convolutions. Pour cela, on exprime les convolutions circulaires comme produit de poly- n6mes que l 'on calcule ~t partir de plusieurs sous- probl6mes plus simples.

En effet, consid6rons la convolution circulaire suivante :

M - 1

(15) Ys = ~ arab . . . . s = 0 . . . . . M - - 1 . r n = O

Un calcul simple prouve que les coefficients Ys sont les coefficients du polyn6me Y(Z) d6fini par :

Y(Z) = P(Z) Q(Z) [mod (Z M - - 1)], off :

M - 1 M - - 1

P(Z) = ~ amZ m et Q(Z) = ~ b~.Z m. n t = O m = O

En fait dans (15) on consid6re la multiplication darts : C(Z)[(Z M J 1). La m&hode de convolution rapide est bas6e sur l'id6e suivante : trouver K anneaux, A I , A2, ..., AK et T u n isomorphisme tel que :

C(Z)](Z ~ - 1 ) r> A1 • A2 • ... • AK.

Dans ce cas, le calcul de (15) se ram6nerait h calculer : T (P(Z)), T (Q(Z)), effectuer des multiplications

7/13 ANN. TI~L~COMMUN., 40, n ~ 9-10, 1985

Page 8: Transformations de Fourier discrètes: algorithmes rapides

502 J . -L . ZOLESIO. - TRANSFORMATIONS DE FOURIER DISCRf~TES

plus simples dans chaque anneau A~ puis calculer : T -1 (T(P(Z)) , T(Q(Z))) .

Cette m6thode 6tant int6r.ssante dans la mzsure off l'image et l'image inwrse d 'un polyn6me sont rapides /t calculer. Pour r6aliser une telle situation, on utilise le r6sultat d'alg6bre suivant [12].

Thdor~me 1.

Soit A un anneau principal et ] r J deux id6aux premiers de A, alors il existe un isomorphisme :

Al l J ~ AI i x Af t .

Le probl6me 6tant alors de d6composer Z M - - 1 en K facteurs premiers, Winograd [13] a d6montr6 que le nombre minimal de multiplications n6cessaires b. l'6valuation de (15) est : 2 M - - K.

Ce r6sultat ne prend pas en eompte les multiplica- tions off le multiplicateur est un scalaire du corps sur lequel on a d6compos6 Z M - - 1 en facteurs irr6ductibles.

Pour diminuer le nombre 2 M - - K on est done tent6 de consid6rer la d6composition de Z u - - 1 sur C en utilisant une racine primitive de l'unit6, K est alors maximal et on a : 2 M - - K = 2 M - - M = M (dans ce cas, T e s t la transform6e de Fourier). Mais un tel choix ne serait pas objectif puisqu'il conduirait ~t ne pas tenir compte des multiplications par des scalaires complexes.

Pour pallier ce probl6me, on choisit le corps Q pour d6composer Z u - - 1, les 6ventuelles multipli- cations par j n'6tant alors pas compt6es.

Pour d6composer Z rt - - 1 sur Q, on utilise les polyn6mes cyclotomiques dont on ne rappelle ici que la d6finition.

D dfinition.

On appclle polyn6me cyclotomique d'indice n, n e N, le polyn6me suivant :

F , (Z) = 1--[ ( Z - co), r

off P, d6signe l'ensemble des racines primitives d 'ordre n de l'unit6 sur C, c'est-h-dire :

P, = {exp ( - - 2 r~ j K/n), K ~ N, l ~< K ~< n - - 1, K et n premiers}.

Ces polyn6mes permettent en effct de d6composer Z ~ - - 1 en facteurs irr6ductibles sur Q, la d6composi- tion 6tant :

K

Z ~ - - 1 = I-I Fd,(Z),

off Fd(Z) est le polyn6me cyclotomique d'indice d et off d~ divise M.

Le probl6me (15) est alors scind6 en K sous- probl6mes :

Y,(Z) = P(Z) Q(Z) [mod Fd,(Z)], i = 1 . . . . , K.

I1 est ~t noter que ces r6ductions modulo un polynSme ne n6cessitent aucune multiplication darts les cas

consid6r6s au paragraphe VI.3. De plus, chaque sous-probl6me pouvant se simplifier,

sans multiplications suppl6mentaires, comme suit �9

Y,(Z) = P(Z) QfZ) [mod (Fd,(Z))]

= (P(Z) [mod (Fd,(Z))]) (Q(Z) [modFd,(Z)]),

[modFd,(Z)],

on est amen6 pour 6valuer Y~(Z) h n'effectuer des multi- plications que sur des polynSmes de degr6 beaueoup plus b a s q u e ceux de P(Z) et Q(Z).

Y(Z) 6tant ensuite reeonstruit de la fa~on suivante :

K

(16) Y(Z) = Y~ S~(Z) Y,(Z) [mod (Z M - - 1)], i = l

off S~(Z) est tel que : Sz(Z) = 8~j [mod F~j(Z)]. Cette m6thode r6pond bien au probl6me puisqu'elle

ne n6cessite que des multiplications entre des poly- nSmes simples, puisque r6duits modulo Fdj(Z), ee qui correspond bien h effectuer les multiplications dans des anneaux plus simples ~t savoir :

C(Z)/Fa~(Z).

D'autre part, T ( P ( Z ) ) s e caleule en 6valuant des r6ductions modulo Fdj(Z) ne n6cessitant, comme on l 'a vu, aucune multiplication.

Les seules multiplications n6cessaires sont done celles utiles pour exprimer les diff6rents coefficients de Y~(Z) pour i = 1 . . . . , K (aux multiplications par les coefficients de S~(Z) pr6s, coefficients rationnels rappelons-le).

VI.3. Aigorithmes courts ou ~ diagonaux )).

Ces algorithmes sont obtenus h partir des m&hodes expos6es en VI.1. et VI.2.

G6n6ralement attribu6s /t Winograd [2], ces algo- rithmes permettent le calcul de la TFD de petites dimensions 3, 5, 7, 8, 9, 16 d'une mani6re tr6s per- formante. Leurs caract6ristiques sont essentiellement de deux ordres, hormis la charge de calcul tr6s faible qu'ils pr6sentent. D'une part, ceux-ci ne comportent que des multiplications par des r6els ou par des imaginaires puts, cette caract6ristique 6tant obtenue grace aux propri6t6s de sym6trie de la suite W ~ du paragraphe VI. 1. et de la d6composition de (Z u - - 1) en (Z u/2 - - 1) (Z u/2 + 1).

D'autre part, ces algorithmes peuvent 8tre syntM- tis6s comrne 6tant le produit de trois matrices (voir [141, [15]) :

A = (I D J),

off I e t J sont des matrices d'incidences (i.e. dont les coefficients sont - - l, 0, q- l) et oh D est une matrice diagonale, I e t J 6tant en g6n6ral rectangu- laires.

Les coefficients de D correspondent aux multiplica- tions, J h la formation des diffSrents coefficients de Y~(Z), i = 1, ..., K et I / t la reconstitution de A grace

(16).

ANN. TI~LI~COMMUN., 40, n ~ 9-10, 1985 8/13

Page 9: Transformations de Fourier discrètes: algorithmes rapides

J.-L. ZOLESIO. - TRANSFORMATIONS DE FOURIER DISCRETES 503

(L'existence d'une d6composition de la forme I D J ne pose aucun probl~me, le probl~me 6tant de r6duire au maximum la taille de la matrice D.)

Cet aspect << diagonal>> de ces algorithmes est comme on le verra fondamental pour les algorithmes de Winograd et rend leur utilisation tr~s performante.

Les figures 5 ~t 10 pr6sentent les graphes de circula- tions des donn6es pour des TFD de dimension 3, 4, 5, 7, 8, 9 et 16 points.

ao_- 1,o ~ _- Ao

2 t FIG. 5. -- Graphe de circulation des donn6es d'un algorithme

court h 3 points.

Signal flow graph of the 3 points short algorithm.

no. _ / 1 . 0 _ ~ a o

al _ ~ _ 1 , 5 3 - 1 ' 2 5 - - ~ 7 " ~ j ~ A4 a4 8841769-: ~ A1

aa _ ~ 1 5 - ' 9 0 1 6 9 9 4 4 - - J / - ~ - ~ ~ ~ , Aa

a2 0,363271264 -" _ =i A3

3,5877852523 -" , J FiG. 6. - - Graphe de circulat ion des donn6es d 'un algorithme

court h 5 points.

Signal flow graph of the 5 points short algorithm.

On notera que les multiplications par des r~els ou des imaginaires purs n'engendrent aucune addition, leur charge de calcul respective est r6sum6e par le tableau III.

TABL. IIL -- Charges de calcul des algorithmes courts exprim~es en op6rations 616mentaires.

Dimension dela TrD 3 4 5 7 8 9

Nombre de multiplications 4 - [ ~ 0 ~6 ~___ 2 _ ~ ~

Nombre d'additions L12[16 32178~52 841148 [

Enfin, on peut remarquer que ces algorithmes pr6sentent le handicap de comporter des permutations.

VII. LES ALGORITHMES A FACTEURS PREMIERS

1,0

- ~ 0,7343022012~

a 4 1 1 ~ _ _ ~ i ~0790156485 = Y ~

I

Ao

FI~. 7. -- Graphe de circulation des donn6es d'un algorithme court de avo ~ 7 points.

Signal flow graph of the 7 points short algorithm.

ao ~,o ~ ~ Ao

al 1 ,o ~ = - A4

a2 1,o - - = , ~ , ~ Ae

a3 1 .o J _ Aa

a4 1,o - - ~ : ~ / A 7

as v'2/2 ~ A3 J a6 1 ,o As

ar _ _ v~?2 AI

FIG. 8. -- Graphe de circulation des donn~es d'un algorithme court de xro ~t 8 points.

Signal flow graph of the 8 points short algorithm.

80 o ~ 1,0 /

a3 X -1,5 a6 _ - v"~-

2 0,5

a7 o,766o444431

al -0.1736481777

a4 0,9396926208 a2 ~ -0,6427876097 a8 / ~ -0,984807753 as _ 0,3420201433

-v3" 2

A. 2

A3 1 ----o Ae

---~ A5

- - - -o~ Aa

A1 A4

Ar

FIG. 9. -- Graphe de circulation des donn6es d'un algorithme court de XFD h 9 points.

Signal flow graph of the 9 points short algorithm.

L'id6e de base de ces algorithmes, introduite par Good [16] puis am61ior6e par Winograd, consiste ~t remarquer que lorsque la dimension N de la TrD

est le produit de deux nombres premiers N1 et N2, on peut se ramener au calcul d'une TrD bidimension- nelle.

9/13 ANN. Tt~LI~COMMUN., 40, n ~ 9-10, 1985

Page 10: Transformations de Fourier discrètes: algorithmes rapides

504 J .-L. ZOLESIO. - TRANSFORMATIONS DE FOURIER DISCRETES

1 1 -1 i :

-I

A

-I -8

V'~-/2 ~ : 0 -

C = - O0

D B=cos 0+cos 3 0 D=sin 0

FIG. 10. ~ Graphe de circulation des donn6es d'un algorithme court de rFD A 16 points.

Signal .[tow ,graph of the 16 points short algorithm.

`% ,% A~2 A~ A~4 ,% A~o A2 A~s

\ / A~

`% , ~ . ~ A~

`%

0=2=/16

En effet, d'apr~s le thdor~me 1 vu prdcddemment, on a un isomorphisme :

ZN --> ZN~ • Z~ 2 ,

n ---> (r/l,/'/2).

Avec n = n~ N2 + n2 N~ [mod. N], on a alors :

NI--1 N2-1 (17) AK,N,+K~N2 = Z Z a.~N2+.2u~Wg~ "'r' W~ "r2,

~1=0 n2=O

(I7) est une transform6e bidimcnsionnelle mais qui a subi une permutation par rapport ~ l 'ordre lexico- graphique habituel, ~ savoir :

a(nl,n2) W~v~ gi W~v~ K2,

comme NI et N2 sont premiers, les indices {N2nlK1} sont les m~mes que {n~K~) ~ une permutation pr6s due au modulo N.

Pour les remettre dans l'ordre habituel, il faut permuter les A (K1, K2) en changeant de variables :

t K~ -+ KTK~, i ~ j, avec N~K'~ = 1 [mod N~], i, j = t, 2,

on a alors :

NI--1 N2--1

A~,,~) = Z Z a~.,.~) W~ ~' W~ ~, nl=O n2=O

(18)

avec :

et A(KxK2) = AKIK*N2+K2K*NI,

a(nl,n2 ) = anlN2+n2N 1 �9

Bien entendu, cette technique peut se g~n~raliser aux cas multi-dimensionnels, ce qui permet le calcul de TFD de grande dimension ~ partir de TFD de petites dimensions.

En effet (18) peut aussi s'6crire :

N I - 1 N2-1 (19) A(Ka,K2) = Z [ Z a(.,,. 2) W~ K'] W~. 2K2,

rtl=O It2=O

ce qui n'est autre qu'une TFD de dimension N1 pour chaque valeur de n2, suivie d'une TFD de dimension N2 pour chaque valeur de K1.

La figure 11 illustre l'exemple off N~ = 3 et N2 = 5.

Cependant, cette technique peut ~tre am~lior6e en mettant A profit la structure diagonale des algo- rithmes courts.

VII.1. Les algorithmes de Winograd.

La construction de ces algorithmes, due A Wino- grad [2], repose sur la stabilit6 de la structure dia- gonale des algorithmes par produit tensoriel.

En effet, (18) peut s'exprimer sous la forme d 'un produit tensoriel :

(20) [,,1] = (w,,, | wN2) [a],

off [x] d~igne le tenseur d6duit de x.

Or, comme la matrice d 'un algorithme court peut se d6composer sous la forme :

wN, = IlDidi,

on a dans (20) en tenant compte des propri&~s du produit tensoriel :

(21) wu t @ wN, = (I, | (DI Q D2)(J~ @d2),

et ~ nouveau, il s'agit d 'un produit d 'une matrice diagonale par deux matrices d'incidences.

La cons6quence de (21) est de concentrer toutes

ANN. T~LI~COMMUN., 40, n ~ 9-10, 1985 10113

Page 11: Transformations de Fourier discrètes: algorithmes rapides

J.-L. ZOLESIO. - TRANSFORMATIONS DE FOURIER DISCRETES 505

an a(n~, ha)

ao a (O,O)j a~o - - at1,0) ~ i ,

a~ a(2 ,0 ) ~ % - - J ' ~ / ' ' '" - ' . a,o:, / \ a , a ( 1 , 1 ) ~ ~ .

a . a ( 2 , 1 ) " z ~ . .

a, a(1,4) - ~

a+ - - = : : 7 - - " 7 - - - -

a;3 - - a ( 1 , 3 )

a a a ( 2 , 3 ) L f 1 r : ~ t ]

a r - - a ( 1 , 2 )

a~ . ~. a ( 2 , 2 ) �9 "

A(K;, K~) A.

I I A ( 2 , 0 1 - - A~o

A(" As A(0,4) - - A~

A ( 2 , 4 ) - Ar A ( , 4 ) A~

A(o:i ,% AI !,1) A~a

A ( 1 , 1 ) - ,% - - A(o,2) . . . . . . ,%

A(2+2) - - A+ �9 A(1,2) A.

/ ~ ( 0 , 3 ) - - A9

A ( 2 , 3 ) ,%

A ( 1 , 3 ) - A~+

l t 2 5 mu l t i p l i ca teu r s

So i t 50 mu l t i p l i ca t i ons O l e m e n t a i r e s

FIG. 11. - - Graphe de circulation des donn~es de l'algorithme ~ facteurs premiers 3 et 5.

Signal flow graph of the prime factor algorithm with'N = 3 • 5.

an a(nl, n=)

a, - - a<o.o}, Z / i / - ~ a,o - - a(1,0) , t

el - - a(1,1) a1~ - - a(2,1) , - - ae a(0,4) , / ~ ' - , ,-IVX, I I - - r . a~ a(1,4) ' / / - - - I A V I ~-: a,+ - - a(2,4) as ~ a(O~3) ~ ~ y~ [ '

aa a(2,31 ,l==-- �9

. . . . . . ~(0~) ~ - /'I V \ I - -

. , - a(~,m . - - - - - - - 7 " / ~ 1 A \ I - L-.-- ._. a= a (2 2) - ' ~.r /.;

"--__.J ,,

VXI

A(KL K9 A.

. ~ 7 ,: A(o,o) - - Ao : / ~ 7 . / A(2,O) . ~ A,o / - ~ f x / 1 i," - . , - A ( . , o ) - - ~ As

J. X(o,4~ I / / A ( 2 , I - A,A" I . " ~ AO,4) - - A=

" - - 7 ~ - A m : ) - I / , / " Aiz:) - - A,, ; / AO,.I) - - .%

- " 7 - - - - - - - - ~(o,2) _ _ ~, I / / ~,(2,')j_ A, ' " ~ A ( 1 , 2 ) - - A. 7 '(0.3, , . jt,(2,3) . _ ~,

A(1,3) - - A~,=

17 multip~icateurs /

Soit 34 multiplications 616mentaires

FIG. 12. - - Graphe de circulation des donn~es de l'algorithme de Winograd de dimension 15.

Signal flow graph o f the Winograd Fourier transform algorithm with N = 15.

les multiplications au mt~me endroit dans le graphe de circulation des donn6es comme le montre l 'exemple de la figure 12.

Le nombre de multiplications 616mentaires est consid6rablement r6duit puisqu'il est dans ce cas de :

2(d~d2 - - (dl - - M~) (d2 - - M2)),

au lieu de : 2(M1N2 + M2N1), dans le cas pr6c6dent, oO d+ d6sigue la dimension de la matrice Dt et o~ Mt d6signe le nombre de coefficients complexes non triviaux de D~ pour i = 1, 2.

La figure 12 montre le graphe de circulations des

donn6es de l 'algorithme de Winograd dans le cas NI = 3 e t N2 = 5.

On pourra remarquer que la charge de calcul, en ce qui concerne les additions, n 'est pas ind~pendante de l 'ordre dans lequel on proc6de quant ~ N~ et N2 �9 Pour avoir la charge de calcul minimale, on a toujours int6r& ~t regrouper le plus au centre du graphe de circulation des donn6es la TFD dont la matrice dia- gonale a la dimension la plus grande.

En tenant compte de cette remarque, la charge de calcul dans le cas de TFD de dimension de la forme N~N2 est donn6e par le tableau IV.

11/13 ANN. T~L~CO~tUN., 40, n ~ 9-I0, 1985

Page 12: Transformations de Fourier discrètes: algorithmes rapides

506

TABL. I V . - - Charge de calcul des algorithmes de Winograd pour de petites dimensions exprim6es en op6rations 616mentaires.

Dimension de la TFD 6 10 1 2 14 15 18 20 21 24

Nombre de multiplications

Nombre d'additions

8 20 16 32 ' 34 40 40 48

36 84 96 '184 172 204 208 318

36

252

281 3oi 35 36 40 45 48

---~! 68 106 80 84 130 9 2 - - - - - -

424 372 614 480 516 708 636

56 63 72 80

132 196 164 200

988 1 290 1 140 1 252

Comme les probl6mes/t r6soudre imposent g6n6rale- merit plus un ordre de dimension des TFD plut6t qu'une valeur pr6cise, et les algorithmes de Winograd ne pouvant pas s'appliquer aux dimensions de la forme 2 ~, on a compar6 la charge de calcul avee l'algorithme de base 4 en encadrant chaque dimension 2 ~ par deux dimensions produit de nombres premiers.

Le gain en charge de calcul est donn6 par le tableau V.

TABL. V. - - Gain en charge de calcul des algorithmes de Winograd.

Dimension TFD base 4

Dimension TFD Winograd

Multiplications gain en 70

Additions gain en %

16

15 18

70 - - 100

1 6 - 38

32

30 36

5,5 ~ II

5 22

60

34,6

. . . . . 11,5

64

72

21

- - 1 6 , 8

Dimension TFD base 4

Dimension rFD Winograd

Multiplications gain en %

Additions gain en %

120

50,7

140

24,3

256

504

53

13,8

128

28,2

252 280

43,7 38,8

--12,4 --22,3 - - 8

512

560

42,9

19,2

J.-L. ZOLES10. - TRANSFORMATIONS DE FOURIER DISCR~TES

I1 ressort de cette analyse que les algorithmes de Winograd s'ils ont d6j/L l'int6r& de permettre le choix d'une dimension autre que 2 ~, ont une charge de calcul tr6s int6ressante quand la dimension croit.

De plus, les multiplications qui interviennent dans le calcul &ant r6elles ou imaginaires pures, le nombre de multiplications n6cessaires est approximativement r6duit de moiti6 pour des donn6es r6elles.

Notons encore que leur macro-modularit6 peut offrir l'int6r&, ~ partir de quelques modules de base, de construire des algorithmes de TFD de dimensions tr6s diff6rentes.

Par contre, ces algorithmes poss6dent les deux inconv6nients d'engendrer un bruit de calcul plus important que celui des algorithmes de base 2 K et de n6cessiter des permutations des entr6es et des sorties, ce qui complique le diagramme des temps des processeurs.

VIH. CONCLUSION

Le choix d'un algorithme de TFD r6sulte, en g6n6ral, d'un compromis entre plusieurs crit6res dent l'ordre d'importance d6pend de l'application particuli6re, de la dimension de la TFD et de la technologie utilis6e. Aussi, n'est-il pas possible d'6tablir une hi6rarchie absolue entre les diff6rents algorithmes, les crit6res de choix 6voluant avec le temps, m~me pour une application donn6e.

En effet, pendant longtemps, la charge de calcul a 6t6 sinon le plus important crit6re, tout au moins un des crit6res privil6gi6s de ce choix. Avec les progr6s permanents de la technologie, la charge de calcul ne rev& pas la meme importance que par le pass6 alors qu'un int6r& croissant appara/t quant/L la modu- larit6 et la r6p6titivit6 des algorithmes (voir [17]).

Toutefois, on peut penser aujourd'hui, que primer uniquement des crit6res relatifs ~t la facilit6 de mise en ceuvre des algorithmes au d&riment total de route tentative pour r6duire la charge de calcul ne serait pas plus astucieux que l'attitude inverse, la meilleure solution &ant bien entendu d'essayer de concitier ces deux points de rue.

Manuscrit reeu le 25 avril 1985,

acceptd le 8 juillet 1985.

BIBLIOGRAPHIE

[1] VAN DER STEEN (g. J.). Timing of some fast Dr'r-algorithms. Signal Processing, NL (nov. 1983), 5, n ~ 6.

[2] WINOGRAD (S.). On computing the Discrete Fourier Transform. Prec. Nat. Acad. Sc., USA (avr. 1976), 73, n ~ 4, pp. 1005-1006.

[3l CURTIS (T. E.). A hardware-based Fourier Transform

Algorithm. Int. conf. on Digital Signal Processing, Florence (2-5 sep. 1981), pp. 613-626.

[4] COOLEY (J. W.), TUKEY (J. W.). An algorithm for the machine calculation of complex Fourier series. Math. Comput., USA (avr. 1965), 19, pp. 297-301.

[5] I~ADER (C. M.), BRENNER (N. M.). A new principle for

ANN. T~L~COMMU~., 40, n ~ 9-10, 1985 12/13

Page 13: Transformations de Fourier discrètes: algorithmes rapides

J.-L. ZOLESIO. - TRANSFORMATIONS DE FOURIER DISCRETES 507

fast Fourier Transform computation. IEEE Trans. ASSP, USA (1976), 24, pp. 264-265.

[6] BERGLAND (G. D.). A Fast Transform Algorithm using base 8 iterations. IEEE Trans. ASSP, USA (juin 1968), 17, n ~ 2, pp. 4-108.

[7] DUHAMEL (P.), HOLLMANN (H.). Split-radix FFT algorithm. Electronics Lett., GB (jan. 1984), 20, n ~ 1, pp. 14-16.

[8] MARTENS O-B.). Recursive cyclotomic factorization. A new algorithm for calculating the Discrete Fourier Trans- form. IEEE Trans. ASSP, USA (avr. 1984), 32, n ~ 2, pp. 750-761.

[9] VETTERLI (M.), NUSSBAUMER (H. J.). Simple FFT and DCT algorithms with reduced number of operations. Signal processing, NL (juil. 1984), 6, n ~ 4, pp. 267-278.

[10] RADEg (C. M.). Discrete Fourier Transform when the number of data samples is prime. Proc. IEEE, USA (juin 1968), 56, pp. 1107-1108.

[11] ZOLESlO (J. L.). Analyse de quelques m6thodes de calcul

de TFD. Revue du CETHEDEC, Fr (1 er trim. 1983), n ~ 74. [12] MCLANE (S.), BmKHOFF (G.). Alg6bre. Structure Fonda-

mentale, t. I. Gauthier-Villars, Fr (1971). [13] WlNOGRAD (S.). Some bilinear forms whose multiplicative

complexity depends on the field of constants. IBM Res. Rep, USA (oct. 1975), R e 5669.

[14] SILVERMAN (H. F.). An introduction to programming the Winograd Fourier Transform Algorithm. IEEE Trans. ASSP, USA (avr. 1977), 25, pp. 152-165.

[15] KOLBA (n.), PARKS (T.). A prime factor FFT algorithm using high-speed convolution. IEEE Trans. ASSP, USA (1977), 25, pp. 90-103.

[16] GOOD ([. J.). The interaction algorithm and practical Fourier analysis. J. Royal Statist. Soe., GB (1958), 20, s6rie B, pp. 361-372.

[17] Lt (G. J.), WAH (B. W.). The design of optimal systolic algorithms. IEEE. Computer software and applications conf., Chicago (7-11 nov. 1983), pp. 310-319.

13/13 ANN. T ~ c o ~ . , 40, n ~ 9-10, 1985