#*****************************************************
# programme d'analyse de Fourier
# Dominique Lefebvre  mars 2013
# www.tangenteX.com
#*****************************************************

# chargement des packages
with(plots):

# définition de la fonction f(t), ici un créneau
f := t -> (-1)^floor(t/Pi):

# définition de la période et de la pulsation de f
T := 2*Pi:
w := 2*Pi/T:
n := 20:

# définition des procédures de calcul des coefficients an et bn
an := proc(f::procedure,m)
 (2/T)*int(f(t)*cos(m*2*Pi*t/T), t= 0..T)
end:

bn := proc(f::procedure,m)
 (2/T)*int(f(t)*sin(m*2*Pi*t/T), t= 0..T)
end:
 
# procedure de calcul de la série de Fourier
SerieF := proc(f::procedure,w,m)
  local i:
  an(f,0)/2 + sum('an(f,i)*cos(i*w*t) + bn(f,i)*sin(i*w*t)', i=1..m):
end:

# calcul de la décomposition de f en série de Fourier
SerieF(f,w,n);

# calcul et affichage du spectre d'amplitude de f
points := [0,0], [0,an(f,0)/2], [0,0]:
for j from 1 to n do
  points := points, [j,0],[j,abs(an(f,j) - I*bn(f,j))],[j,0]:
od:
sp1 := plot([points],0..(n+1),color=red, numpoints=1000,title = "Spectre d'amplitude", 
           labels = ["harmonique","amplitude"]):
display([sp1]);

# calcul et affichage du spectre de puissance de f
points := [0,0], [0,an(f,0)/2], [0,0]:
for j from 1 to n do
  points := points, [j,0],[j,(abs(an(f,j) - I*bn(f,j)))^2],[j,0]:
od:
sp2 := plot([points],0..(n+1),color=red, numpoints=1000,title = "Spectre de puissance", 
           labels = ["harmonique","puissance"]):
display([sp2]);