Previous Up Next

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=[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ωNk· j avec ωN racine N-ième de l’unité.
Le but est de calculer pour k=0..N−1:

N−1
j=0
xjωNk· 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.

Les propriétes

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/Nj=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−1xjykj
On a alors :
N*FN(x · y)=FN(x) * FN(y)
FN(x * y)=FN(x) · FN(y)

Où interviennent ces sommes ?

  1. Valeur d’un polynôme
    Soit un polynôme P(x)=∑j=0N−1cjxj donné par la liste de ces coefficients [c0,c1,..cN−1] complété par des zéros pour que N soit une puissance de 2.
  2. Interpolation trigonométrique
    Soit f une fonction 2π-périodique dont on connait les valeurs en xk=2kπ/N, f(xk)=f(2kπ/N)=fk 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)=a0/2+∑n=0N/2−1ancos(nx)+bnsin(nx)+1/2 aN/2cos(Nx/2)
    on suppose donc b0=bN/2=0 et on choisit cette notation pour simplifier les calculs ultérieurs.
    On a aussi :
    p(x)=1/2 pN/2exp(−iNx/2)+∑n=−N/2+1N/2−1pnexp(inx)+1/2 pN/2exp(iNx/2).
    avec pour k=0..N/2
    pk=1/2(akibk)
    pk=1/2(ak+ibk)
    comme bN/2=0 on a pN/2=pN/2
    On veut déterminer les pk pour avoir :
    p(xk)=fk=
    1/2 pN/2exp(−ikπ)+∑n=−N/2+1N/2−1pnexp(in2kπ/N)+1/2 pN/2exp(ikπ)
    puisque exp(−ikπ)=exp(ikπ)=−1 et pN/2=pN/2 on a
    p(xk)=fk=−pN/2+∑n=−N/2+1N/2−1pnexp(in2kπ/N)
    On transforme ∑n=−N/2+1−1 en posant j=n+N on obtient :
    n=−N/2−1 pn exp(inx)=∑j=N/2N−1pjNexp(i(jN)x)
    puisque pour k=0..(N−1), exp(i(jN)xk)=exp(ijxk),
    on a :
    p(xk)=∑j=N/2N−1pjNexp(ijxk)+∑n=0N/2−1pnexp(inxk).
    ou encore :
    p(xk)=∑n=0N/2−1pnexp(inxk)+∑n=N/2N−1pnNexp(inxk).
    On pose :
    pour n=0..(N/2−1), qn=pn et
    pour n=N/2..(N−1), qn=pnN.
    Donc pour k=0..(N−1), p(xk)=∑n=0N−1qnexp(inxk)=fk.
    Il suffit donc de résoudre ce système d’inconnues qn.
    On a :
    fk=f(2kπ/N)=p(xk)=N*FN−1(q)=N*ifft(q)
    soit :
    q=[p0,..pN/2−1,pN/2,..,p−1]=1/NFN([f0,..f(N−1)])=
    1/Nfft([f0,...f(N-1)])
    Remarque
    Si la fonction f est réelle on a :
    pk=pk pour k=1..(N/2−1) et p0 et pN/2 sont réels.
    Donc si la fonction f est réelle on a :
    p(x)=p0+2*re(∑k=0N/2−1pkexp(ikx)+1/2pN/2*exp(−iNx/2))
    En résumé :
    pour 0 ≤ n<N/2, pn est le n-ième élément de 1/NFN([f0,..f(N−1)]) et,
    pour −N/2 ≤ n<0, pn est le (n+N)-ième élément de 1/NFN([f0,..f(N−1)]).
  3. Série de Fourier
    Soit f une fonction 2π-périodique dont on connait les valeurs en xk=2kπ/N, f(xk)=f(2kπ/N)=yk pour k=0..(N−1).
    On suppose que sa série de Fourier converge simplement vers f, et on veut connaître les coefficients cn pour −N/2 ≤ n<N/2.
    Il fait donc calculer de façon approchée :
    cn=1/2π ∫0f(t)exp(−int)dt
    Remarque
    Lorsque la fonction périodique f est égale à sa série de Fourier, on a :
    f(x)=a0/2+∑k=1akcos(kx)+bksin(kx) avec ck=akibk/2.

    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π/Nk=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,cN/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), ωNnNn+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−1cnexp(inx) interpole f(x) aux points x=2kπ/N.
    Donc si f est un polynôme trigonométrique P de degré mN/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 cncn.
    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/Nk=0N−1ykωNkn
    En remplaçant yk on obtient :
    cn=1/Nk=0N−1m=−∞+∞cmωNkmωNkn
    Or :
    ωNkmωNknNk(mn) et
    puisque ωNmn est une racine N-ième de l’unité, on a :
    ωN(mn)N=1 et ∑k=0N−1ωN(mn)k=0.
    Donc :
    si mn est un multiple de N (m=n+l· N) on a ∑k=0N−1ωNk(mn)=N et
    sinon ∑k=0N−1ωNk(mn)=0
    En intervertisant les deux sommes on a :
    cn=1/Nm=−∞+∞cmk=0N−1ωNk(mn)
    cn=∑l=−∞+∞c(n+l· N)
    c’est à dire :
    cn=...cn−2· N+cnN+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=ck+ck vaut 1 si k=1,2 et 0 sinon.

  4. Produit de convolution
    Soient deux polynômes P(x)=∑j=0n−1ajxj et Q(x)=∑j=0m−1bjxj donné par la liste de leurs coefficients a=[a0,a1,..an−1] et b=[b0,b1,..bm−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=[a0,a1,..an−1,0..0] et b=[b0,b1,..bm−1,0..0]:
    P(x)Q(x)=∑j=0n+m−1(a*b)jxj
    On calcule donc :
    Fn(a), Fn(b), Fn−1(Fn(aFn(b)) et puisque
    nFn(x · y)=Fn(x) * Fn(y)
    Fn(x * y)=Fn(x) · Fn(y)
    on a :
    a*b=Fn−1(Fn(aFn(b))

Previous Up Next