Trasformata di Fourier discreta: DFT

Scegliamo come numero di punti sull'intervallo [0, 2π) una potenza di 2:

nPti = 8;

Introduciamo la radice ottava primitiva dell'unità:

ω  =  Exp[2 π I /nPti]

Coefficienti approssimati (da -4 a 3):

c[n_] :=  
      Sum[f[2 k π/nPti] ω^(-k n), {k, 0, nPti - 1}]/
        nPti

Introduciamo una funzione di prova:

f[x_] :=  x

Ad esempio:

c[3] // ComplexExpand

Calcoliamo tutti i coefficienti:

coeff  =  Table[c[k], {k, -4, 3}] // ComplexExpand

Attenzione: le posizioni occupate dai termini nella lista appena calcolata vanno da 1 a 8, in corrispondenza degli indici che vanno da -4 a 3.

p[x_]  :=  
      Sum[ coeff[[ 
            k + 5]]  Exp[I k x], {k, -4, 3}]

p[x] // ComplexExpand

Poiché la funzione introdotta assume solo valori reali, la parte immaginaria del polinomio interpolante deve annullarsi nei punti d'interpolazione e la parte reale deve assumere gli stessi valori di f(x):

pReale  =  ((7 π)/8) - (1/4) π Cos[x] - 
          (1/4) π Cos[2 x] - (1/4) π Cos[3 x] - 
          (1/8) π Cos[4 x] - 
          (1/4) π Sin[x] - ((π Sin[x])/(2 (Sqrt[2]))) - 
          (1/4) π Sin[2 x] + 
          (1/4) π Sin[
              3 x] - ((π Sin[3 x])/(2 (Sqrt[2])));

gr1  =  
        Plot[{x, pReale}, {x, 0, 
            2  π},  
PlotStyle -> {{}, {RGBColor[
                  0, 0, 1],  
                Thickness[0.01]}}, 
AspectRatio -> 
            Automatic];

Aggiungiamo dei segmenti verticali in corrispondenza dei punti d'interpolazione:

gr2  =  
        Show[Graphics[{RGBColor[1, 0, 0], 
              Table[Line[
{{x, 0}, {x, f[x]}}],  {x, 0, 
                  2  π,  2  π/nPti}]}] ];

Show[gr1,  gr2];

pImmaginaria  =  (1/8) π Sin[4 x];

Plot[{x, pImmaginaria}, {x, 0, 2  π}, 
        AspectRatio -> Automatic,  PlotRange -> All];

Calcoliamo anche i coefficienti esatti:

coeffEsatti  =  
      Table[(1/(2  π)) Integrate[
              x Exp[-I k x], {x, 0, 
                2  π}], {k, -4, 3}] // ComplexExpand

Secondo esempio

Introduciamo un'onda quadra:

f[x_] := If[x ≤ π, 1, 0]

coeff  =  Table[c[k], {k, -4, 3}] // ComplexExpand

p[x_]  :=  
      Sum[ coeff[[ k + 5]]  Exp[I k x], {k, -4, 
            3}] // ComplexExpand

p[x]

pReale  =  
        (5/8) + (1/4) Cos[2 x] + (1/8) Cos[4 x] + Sin[x]/4 + 
          Sin[x]/(2 (Sqrt[2])) - (1/4) Sin[3 x] + 
          Sin[3 x]/(2 (Sqrt[2]));

gr1  =  
        Plot[{f[x], pReale}, {x, 0, 
            2  π},  
PlotStyle -> {{}, {RGBColor[
                  0, 0, 1],  Thickness[0.01]}}];

gr2  =  
        Show[Graphics[{RGBColor[1, 0, 0], 
              Table[Line[
{{x, 0}, {x, f[x]}}],  {x, 0, 
                  2  π,  2  π/nPti}]}] ];

Show[gr1, gr2];

Terzo esempio

Introduciamo un'onda triangolare (funzione pari):

f[x_] := If[x ≤ π, x, 2  π  -  x]

coeff  =  Table[c[k], {k, -4, 3}] // ComplexExpand

p[x_]  :=  
      Sum[ coeff[[ 
              k + 5]]  Exp[I k x], {k, -4, 3}] // ComplexExpand

p[x]

gr1  =  
        Plot[{f[x], Chop[p[x]]}, {x, 0, 
            2  π},  
PlotStyle -> {{}, {RGBColor[
                  0, 0, 1],  Thickness[0.01]}}];

gr2  =  
        Show[Graphics[{RGBColor[1, 0, 0], 
              Table[Line[
{{x, 0}, {x, f[x]}}],  {x, 0, 
                  2  π,  2  π/nPti}]}] ];

Show[gr1,  gr2];