clear all
close all

%definizione funzione periodica. es: onda quadra

T=xx;                                                                %periodo                  
Fmax=xx;                                                           %modulo della forza
N=1000;                                                             %nr. di punti per periodo (PARI)
dt=T/N;                
vett_t=dt:dt:T;                                                     %vettore dei tempi (1 periodo)
% vett_F=[Fmax*ones(1,N/4) -Fmax*ones(1,N/2) Fmax*ones(1,N/4)];       %vettore della forza
vett_F=[Fmax*ones(1,N/2) -Fmax*ones(1,N/2) ];       %vettore della forza
plot(vett_t,vett_F,'LineWidth',2)
grid
axis([0 T -Fmax*1.05 Fmax*1.05])

%analisi armonica con la funzione fft di matlab

grf=fft(vett_F,N);
nmax=N/2;
delta_freq=1/T;                             %passo in frequenza 
freqF=delta_freq*(0:(nmax-1));              %vettore delle frequenze 
modF(1)=grf(1)/N;                           %modulo armonica a 0 Hz 
modF(2:nmax)=2/N*abs(grf(2:nmax));          %modulo delle armoniche con f>0
immF=imag(grf(2:nmax));                     %parte immaginaria
reaF=real(grf(2:nmax));                     %parte immaginaria
faseF(2:nmax)=atan2(immF,reaF);             %fase in radianti

figure
subplot(211)
bar(freqF,modF,0.1)
subplot(212)
bar(freqF,faseF,0.1)

%modulo/fase prime tre armoniche diverse da zero

mod1=modF(1)
mod2=modF(2)
mod3=modF(3)

fas1=faseF(1)
fas2=faseF(2)
fas3=faseF(3)