6.26.2  Transformée de Fourier discrète

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=[x₀,x₁,...x_N−1].
La transformée de Fourier discréte est une transformation F_N 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 :
(F_N(x))_k=y_k=∑_j=0^N−1x_jω_N^{−k· j} avec ω_N racine N-ième de l’unité.
Le but est de calculer pour k=0..N−1:

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.

Les propriétes

La transformée de Fourier discréte F_N est une transformation bijective.
On a :
F_N⁻¹=1/N F_N
c’est à dire :
(F_N⁻¹(x))_k=1/N∑_j=0^N−1x_jω_N^{k· j}
Notation
Avec Xcas on a les fonctions fft et ifft qui sont :
fft(x)=F_N(x) et
ifft(x)=F_N⁻¹(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=x_ky_k
- le produit de convolution, noté *, par :
(x * y)_k=∑_j=0^N−1x_jy_k−j
On a alors :
N*F_N(x · y)=F_N(x) * F_N(y)
F_N(x * y)=F_N(x) · F_N(y)

Où interviennent ces sommes ?

1. Valeur d’un polynôme
   Soit un polynôme P(x)=∑_j=0^N−1c_jx^j donné par la liste de ces coefficients [c₀,c₁,..c_N−1] complété par des zéros pour que N soit une puissance de 2.
   • On veut calculer les valeurs prises par P(x) aux points a_k=ω_N^−k=exp(−2ikπ/N) pour k=0..N−1, c’est à dire les valeurs ∑_j=0^N−1c_j(ω_N^−k)^j pour k=0..N−1.
     On a :
     [P(a₀), P(a₁),..P(a_N−1)]=F_N([c₀,c₁,..c_N−1])
     Par exemple on tape :
     P(x):=x+x^2 et w:=i
     Les coefficients de P sont [0,1,1,0], ω=w=i, N=4 et exp(2iπ/4)=i.
     On a :
     fft([0,1,1,0])=[2,-1-i,0,-1+i] et donc
     P(0)=2,
     P(-i)=P(w^-1)=-1-i,
     P(-1)=P(w^-2)=0,
     P(i)=P(w^-3)=-1+i.
   • On veut calculer les valeurs prises par P(x) aux points :
     b_k=ω_N^k=exp(2ikπ/N) pour k=0..(N−1), c’est à dire calculer;les valeurs :
     ∑_j=0^N−1c_j(ω_N^k)^j pour k=0..(N−1).
     On a :
     [P(a₀), P(a₁),..P(a_(N−1))]=N*F_N⁻¹([c₀,c₁,..c_(N−1)])
     Par exemple on tape :
     P(x):=x+x^2 et w:=i
     Les coefficients de P sont [0,1,1,0], ω=w=i, N=4 et exp(2iπ/4)=i.
     On a :
     4*ifft([0,1,1,0])=[2,-1+i,0,-1-i] et donc
     P(0)=2,
     P(i)=P(w^1)=-1+i,
     P(-1)=P(w^2)=0,
     P(-i)=P(w^3)=-1-i.
2. Interpolation trigonométrique
   Soit f une fonction 2π-périodique dont on connait les valeurs en x_k=2kπ/N, f(x_k)=f(2kπ/N)=f_k pour k=0..(N−1) avec N=2*m.
   On veut déterminer le polynôme p trigonométrique interpolateur de f aux points 2kπ/N pour k=0..(N−1), sous la forme :
   p(x)=a₀/2+∑_n=0^N/2−1a_ncos(nx)+b_nsin(nx)+1/2 a_N/2cos(Nx/2)
   on suppose donc b₀=b_N/2=0 et on choisit cette notation pour simplifier les calculs ultérieurs.
   On a aussi :
   p(x)=1/2 p_−N/2exp(−iNx/2)+∑_n=−N/2+1^N/2−1p_nexp(inx)+1/2 p_N/2exp(iNx/2).
   avec pour k=0..N/2
   p_k=1/2(a_k−ib_k)
   p_−k=1/2(a_k+ib_k)
   comme b_N/2=0 on a p_−N/2=p_N/2
   On veut déterminer les p_k pour avoir :
   p(x_k)=f_k=
   1/2 p_−N/2exp(−ikπ)+∑_n=−N/2+1^N/2−1p_nexp(in2kπ/N)+1/2 p_N/2exp(ikπ)
   puisque exp(−ikπ)=exp(ikπ)=−1 et p_−N/2=p_N/2 on a
   p(x_k)=f_k=−p_−N/2+∑_n=−N/2+1^N/2−1p_nexp(in2kπ/N)
   On transforme ∑_n=−N/2+1⁻¹ en posant j=n+N on obtient :
   ∑_n=−N/2⁻¹ p_n exp(inx)=∑_j=N/2^N−1p_j−Nexp(i(j−N)x)
   puisque pour k=0..(N−1), exp(i(j−N)x_k)=exp(ijx_k),
   on a :
   p(x_k)=∑_j=N/2^N−1p_j−Nexp(ijx_k)+∑_n=0^N/2−1p_nexp(inx_k).
   ou encore :
   p(x_k)=∑_n=0^N/2−1p_nexp(inx_k)+∑_n=N/2^N−1p_n−Nexp(inx_k).
   On pose :
   pour n=0..(N/2−1), q_n=p_n et
   pour n=N/2..(N−1), q_n=p_n−N.
   Donc pour k=0..(N−1), p(x_k)=∑_n=0^N−1q_nexp(inx_k)=f_k.
   Il suffit donc de résoudre ce système d’inconnues q_n.
   On a :
   f_k=f(2kπ/N)=p(x_k)=N*F_N⁻¹(q)=N*ifft(q)
   soit :
   q=[p₀,..p_N/2−1,p_−N/2,..,p₋₁]=1/NF_N([f₀,..f_(N−1)])=
   1/Nfft([f₀,...f_(N-1)])
   Remarque
   Si la fonction f est réelle on a :
   p_−k=p_k pour k=1..(N/2−1) et p₀ et p_−N/2 sont réels.
   Donc si la fonction f est réelle on a :
   p(x)=p₀+2*re(∑_k=0^N/2−1p_kexp(ikx)+1/2p_−N/2*exp(−iNx/2))
   En résumé :
   pour 0 ≤ n<N/2, p_n est le n-ième élément de 1/NF_N([f₀,..f_(N−1)]) et,
   pour −N/2 ≤ n<0, p_n est le (n+N)-ième élément de 1/NF_N([f₀,..f_(N−1)]).
3. Série de Fourier
   Soit f une fonction 2π-périodique dont on connait les valeurs en x_k=2kπ/N, f(x_k)=f(2kπ/N)=y_k pour k=0..(N−1).
   On suppose que sa série de Fourier converge simplement vers f, et on veut connaître les coefficients c_n pour −N/2 ≤ n<N/2.
   Il fait donc calculer de façon approchée :
   c_n=1/2π ∫₀^2πf(t)exp(−int)dt
   Remarque
   Lorsque la fonction périodique f est égale à sa série de Fourier, on a :
   f(x)=a₀/2+∑_k=1^∞a_kcos(kx)+b_ksin(kx) avec c_k=a_k−ib_k/2.

   On calcule l’intégrale donnant c_n 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 c_n^≃ est la valeur approchée de c_n obtenue par la méthode des trapèzes, on a pour −N/2 ≤ n<N/2:
   c_n^≃=1/2π2π/N∑_k=0^N−1y_kexp(−2inkπ/N)
   puisque x_k=2kπ/N et f(x_k)=y_k on a
   f(x_k)exp(−inx_k)=y_kexp(−2inkπ/N) et
   f(0)exp(0)=f(2π)exp(−2inNπ/N)=y₀=y_N
   On a donc :
   [c₀^≃,..c_N/2−1^≃,c_−N/2^≃,..c₋₁^≃]=1/NF_N([y₀,y₁...y_(N−1)]) car
   si n≥0, on a c_n^≃=y_n et
   si n<0, on a c_n^≃=y_n+N si ω_N=exp(2iπ/N), ω_Nⁿ=ω_N^n+N puisque ω_N^N=1

   Propriétés
   On retrouve les coefficients du polynôme interpolateur de f :
   p_n=c_n^≃ pour −N/2 ≤ n<N/2.
   Ce qui veut dire que le polynôme trigonométrique ∑_n=−N/2^N/2−1c_n^≃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=−m^m−1c_kexp(2ikπ t), le polynôme trigonométrique qui interpole f=P est P lui-même, les coefficients approchés sont exacts (c_n^≃=c_n).
   On peut, plus généralement calculer l’erreur c_n^≃−c_n.
   On suppose que f est égale à sa série de Fourier, c’est à dire que l’on a :
   f(t)=∫_m=−∞^+∞c_mexp(2iπ mt) avec ∑_m=−∞^+∞|c_m|<∞
   On a donc :
   f(x_k)=f(2kπ/N)=y_k=∑_m=−∞^+∞c_mω_N^km et
   c_n^≃=1/N∑_k=0^N−1y_kω_N^−kn
   En remplaçant y_k on obtient :
   c_n^≃=1/N∑_k=0^N−1∑_m=−∞^+∞c_mω_N^kmω_N^−kn
   Or :
   ω_N^kmω_N^−kn=ω_N^k(m−n) et
   puisque ω_N^m−n est une racine N-ième de l’unité, on a :
   ω_N^(m−n)N=1 et ∑_k=0^N−1ω_N^(m−n)k=0.
   Donc :
   si m−n est un multiple de N (m=n+l· N) on a ∑_k=0^N−1ω_N^k(m−n)=N et
   sinon ∑_k=0^N−1ω_N^k(m−n)=0
   En intervertisant les deux sommes on a :
   c_n^≃=1/N∑_m=−∞^+∞c_m∑_k=0^N−1ω_N^k(m−n)
   c_n^≃=∑_l=−∞^+∞c_{(n+l· N)}
   c’est à dire :
   c_n^≃=...c_{n−2· N}+c_n−N+c_nc_n+Nc_{n+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 :
   c₀=0,c₁=4.0/8,c₂=4.0/2,c₃=0.0,
   c₋₄=0,c₋₃=0,c₋₂=4.0/8,=c₋₁=4.0/8
   donc retrouve bien :
   b_k=0 et a_k=c_−k+c_k vaut 1 si k=1,2 et 0 sinon.

4. Produit de convolution
   Soient deux polynômes P(x)=∑_j=0ⁿ⁻¹a_jx^j et Q(x)=∑_j=0^m−1b_jx^j donné par la liste de leurs coefficients a=[a₀,a₁,..a_n−1] et b=[b₀,b₁,..b_m−1].
   On veut calculer le produit de ces deux polynômes c’est à dire le produit de convolution de leurs coefficients si on se place dans l’ensemble des suites périodiques de période supérieure ou égale à (n+m) en complétant a (resp b) par m+? (resp n+?) zéros (en pratique on choisit ? pour que n+m+? soit une puissance de 2).
   On alors pour a=[a₀,a₁,..a_n−1,0..0] et b=[b₀,b₁,..b_m−1,0..0]:
   P(x)Q(x)=∑_j=0^n+m−1(a*b)_jx^j
   On calcule donc :
   F_n(a), F_n(b), F_n⁻¹(F_n(a)· F_n(b)) et puisque
   nF_n(x · y)=F_n(x) * F_n(y)
   F_n(x * y)=F_n(x) · F_n(y)
   on a :
   a*b=F_n⁻¹(F_n(a)· F_n(b))