Soit N un entier.
On considère une suite x périodique de période N : elle est
entièrement déterminée par la liste x=[x0,x1,...xN−1].
La transformée de Fourier discréte est une transformation FN qui a une
suite x périodique de période N fait correspondre une suite y,
périodique de période N, définie pour k=0..N−1 par :
(FN(x))k=yk=∑j=0N−1xjωN−k· j
avec ωN racine N-ième de l’unité.
Le but est de calculer pour k=0..N−1:
| xjωN−k· j |
avec ωN=exp(2iπ/N)
La méthode de la transformée de Fourier rapide permet de calculer
rapidement ces sommes si N est une puissance de 2.
La transformée de Fourier discréte FN est une transformation bijective.
On a :
FN−1=1/N FN
c’est à dire :
(FN−1(x))k=1/N∑j=0N−1xjωNk· j
Notation
Avec Xcas on a les fonctions fft et ifft qui sont :
fft(x)=FN(x) et
ifft(x)=FN−1(x)
Définiton
Soient deux suites x et y périodiques de période N.
On définit :
- le produit de Hadamard, noté ·, par :
(x · y)k=xkyk
- le produit de convolution, noté *, par :
(x * y)k=∑j=0N−1xjyk−j
On a alors :
N*FN(x · y)=FN(x) * FN(y)
FN(x * y)=FN(x) · FN(y)
^
2 et
w:=i^
-1)=-1-i,^
-2)=0,^
-3)=-1+i.
^
2 et
w:=i^
1)=-1+i,^
2)=0,^
3)=-1-i.
On calcule l’intégrale donnant cn de manière approchée par la
méthode des trapèzes. Ici Romberg
ne sert a rien, car le développement d’Euler Mac Laurin a ses coefficients
déjà nuls puisque la fonction que l’on intègre est périodique et
donc toutes ses dérivées sont égales en 0 et en 2π.
Si cn≃ est la valeur approchée de cn obtenue par la
méthode des trapèzes, on a pour
−N/2 ≤ n<N/2:
cn≃=1/2π2π/N∑k=0N−1ykexp(−2inkπ/N)
puisque xk=2kπ/N et f(xk)=yk on a
f(xk)exp(−inxk)=ykexp(−2inkπ/N) et
f(0)exp(0)=f(2π)exp(−2inNπ/N)=y0=yN
On a donc :
[c0≃,..cN/2−1≃,c−N/2≃,..c−1≃]=1/NFN([y0,y1...y(N−1)])
car
si n≥0, on a cn≃=yn et
si n<0, on a cn≃=yn+N
si ωN=exp(2iπ/N), ωNn=ωNn+N puisque
ωNN=1
Propriétés
On retrouve les coefficients du polynôme interpolateur de f :
pn=cn≃ pour −N/2 ≤ n<N/2.
Ce qui veut dire que le polynôme trigonométrique
∑n=−N/2N/2−1cn≃exp(inx)
interpole f(x) aux points x=2kπ/N.
Donc si f est un polynôme trigonométrique P de degré m≤ N/2 :
f(t)=P(t)=∑k=−mm−1ckexp(2ikπ t),
le polynôme trigonométrique qui interpole f=P est P lui-même, les
coefficients approchés sont exacts (cn≃=cn).
On peut, plus généralement calculer l’erreur cn≃−cn.
On suppose que f est égale à sa série de Fourier, c’est à dire que
l’on a :
f(t)=∫m=−∞+∞cmexp(2iπ mt) avec
∑m=−∞+∞|cm|<∞
On a donc :
f(xk)=f(2kπ/N)=yk=∑m=−∞+∞cmωNkm et
cn≃=1/N∑k=0N−1ykωN−kn
En remplaçant yk on obtient :
cn≃=1/N∑k=0N−1∑m=−∞+∞cmωNkmωN−kn
Or :
ωNkmωN−kn=ωNk(m−n) et
puisque ωNm−n est une racine N-ième de l’unité, on a :
ωN(m−n)N=1 et ∑k=0N−1ωN(m−n)k=0.
Donc :
si m−n est un multiple de N (m=n+l· N) on a ∑k=0N−1ωNk(m−n)=N et
sinon ∑k=0N−1ωNk(m−n)=0
En intervertisant les deux sommes on a :
cn≃=1/N∑m=−∞+∞cm∑k=0N−1ωNk(m−n)
cn≃=∑l=−∞+∞c(n+l· N)
c’est à dire :
cn≃=...cn−2· N+cn−N+cncn+Ncn+2· N+.....
Exemple :
f(t):=cos(t)+cos(2*t)
x:=f(2*k*pi/8)$(k=0..7)
On obtient :
x={2,(√2)/2,-1,-((√2)/2),0,-((√2)/2),-1,(√2)/2}
fft(x)=[0.0,4.0,4.0,0.0,0.0,0.0,4.0,4.0]
donc en divisant par N=8 :
c0=0,c1=4.0/8,c2=4.0/2,c3=0.0,
c−4=0,c−3=0,c−2=4.0/8,=c−1=4.0/8
donc retrouve bien :
bk=0 et ak=c−k+ck vaut 1 si k=1,2 et 0 sinon.