358
9 Polynˆ omes orthogonaux en th´ eorie de l’approximation
9.9.1 La transformation de Fourier rapide
Comme on l’a signal´ e ` a la section pr´ ec´ edente, le calcul de la transformation
de Fourier discr` ete (DFT) ou de son inverse (IDFT) par un produit matricevecteur, n´ ecessiterait N
2 op´ erations. Dans cette section nous illustrons les
´ etapes de base de l’algorithme de Cooley-Tukey [CT65], commun´ ement appel´ e
transformation de Fourier rapide (ou FFT pour Fast Fourier Transform). Le
calcul d’une DFT d’ordre N est d´ ecompos´ e en DFT d’ordre p 0 , . . ., p m , o` u les
p i sont les facteurs premiers de N . Si N est une puissance de 2, le coˆ ut du
calcul est de l’ordre de N log 2 N flops.
Voici un algorithme r´ ecursif pour calculer la DFT quand N est une puissance de 2. Soit f = (f i )
T , i = 0, . . . , N − 1 et soit p(x) =
1
N
N−1
j=0 f j x
j .
Alors, le calcul de la DFT du vecteur f revient ` a ´ evaluer p(W
k−
N
2
N
) pour
k = 0, . . . , N − 1. Introduisons les polynˆ omes
p e (x) =
1
N
f 0 + f 2 x + . . . + f N−2 x
N
2 −1
,
p o (x) =
1
N
f 1 + f 3 x + . . . + f N−1 x
N
2 −1
.
Remarquer que
p(x) = p e (x
2 ) + xp o (x
2 ),
d’o` u on d´ eduit que le calcul de la DFT de f peut ˆ etre effectu´ e en ´ evaluant les
polynˆ omes p e et p o aux points W
2(k−
N
2 )
N
, k = 0, . . ., N − 1. Puisque
W
2(k−
N
2 )
N
= W
2k−N
N
= exp
−i
2πk
N/2
exp(i2π) = W
k
N/2 ,
on doit ´ evaluer p e et p o aux racines principales de l’unit´ e d’ordre N/2. De
cette mani` ere, la DFT d’ordre N est r´ ecrite en termes de deux DFT d’ordre
N/2 ; naturellement, on peut appliquer r´ ecursivement ce proc´ ed´ e pour p o et
p e . Le processus s’ach` eve quand le degr´ e des derniers polynˆ omes construits
est ´ egal ` a un.
Dans le Programme 78, nous proposons une impl´ ementation r´ ecursive simple
de la FFT. Les param` etres d’entr´ ee sont f et NN, o` u f est un vecteur contenant
NN valeurs f k , et o` u NN est une puissance de 2.
Programme 78 - fftrec : Algorithme de FFT r´ ecursif
function [fftv]=fftrec(f,NN)
%FFTREC Algorithme de FFT r´ ecursif.
N = length(f); w = exp(-2*pi*sqrt(-1)/N);
if N == 2
fftv = f(1)+w.ˆ[-NN/2:NN-1-NN/2]*f(2);
9 Polynˆ omes orthogonaux en th´ eorie de l’approximation
9.9.1 La transformation de Fourier rapide
Comme on l’a signal´ e ` a la section pr´ ec´ edente, le calcul de la transformation
de Fourier discr` ete (DFT) ou de son inverse (IDFT) par un produit matricevecteur, n´ ecessiterait N
2 op´ erations. Dans cette section nous illustrons les
´ etapes de base de l’algorithme de Cooley-Tukey [CT65], commun´ ement appel´ e
transformation de Fourier rapide (ou FFT pour Fast Fourier Transform). Le
calcul d’une DFT d’ordre N est d´ ecompos´ e en DFT d’ordre p 0 , . . ., p m , o` u les
p i sont les facteurs premiers de N . Si N est une puissance de 2, le coˆ ut du
calcul est de l’ordre de N log 2 N flops.
Voici un algorithme r´ ecursif pour calculer la DFT quand N est une puissance de 2. Soit f = (f i )
T , i = 0, . . . , N − 1 et soit p(x) =
1
N
N−1
j=0 f j x
j .
Alors, le calcul de la DFT du vecteur f revient ` a ´ evaluer p(W
k−
N
2
N
) pour
k = 0, . . . , N − 1. Introduisons les polynˆ omes
p e (x) =
1
N
f 0 + f 2 x + . . . + f N−2 x
N
2 −1
,
p o (x) =
1
N
f 1 + f 3 x + . . . + f N−1 x
N
2 −1
.
Remarquer que
p(x) = p e (x
2 ) + xp o (x
2 ),
d’o` u on d´ eduit que le calcul de la DFT de f peut ˆ etre effectu´ e en ´ evaluant les
polynˆ omes p e et p o aux points W
2(k−
N
2 )
N
, k = 0, . . ., N − 1. Puisque
W
2(k−
N
2 )
N
= W
2k−N
N
= exp
−i
2πk
N/2
exp(i2π) = W
k
N/2 ,
on doit ´ evaluer p e et p o aux racines principales de l’unit´ e d’ordre N/2. De
cette mani` ere, la DFT d’ordre N est r´ ecrite en termes de deux DFT d’ordre
N/2 ; naturellement, on peut appliquer r´ ecursivement ce proc´ ed´ e pour p o et
p e . Le processus s’ach` eve quand le degr´ e des derniers polynˆ omes construits
est ´ egal ` a un.
Dans le Programme 78, nous proposons une impl´ ementation r´ ecursive simple
de la FFT. Les param` etres d’entr´ ee sont f et NN, o` u f est un vecteur contenant
NN valeurs f k , et o` u NN est une puissance de 2.
Programme 78 - fftrec : Algorithme de FFT r´ ecursif
function [fftv]=fftrec(f,NN)
%FFTREC Algorithme de FFT r´ ecursif.
N = length(f); w = exp(-2*pi*sqrt(-1)/N);
if N == 2
fftv = f(1)+w.ˆ[-NN/2:NN-1-NN/2]*f(2);