2.5.1

// Campiono la funzione fx nell'intervallo da Tmin a Tmax in N campioni

// ed effettuo la trasformata di Fourier discreta (DFT)

//

// 09/10/2019 saurischio@gmail.com

//---------------------------------------------------------------------

clf;clc;clear;

// INPUT

function f=fx(t) // funzione fx

    f=20*sin(50*t)+10*cos(60*t)

endfunction

Tmin=-1; // tempo iniziale

Tmax=1; // tempo finale

N=100; // numero di campioni

maxm=40 // le prime maxm pulsazioni che ho deciso di graficare

// ESECUZIONE

T=Tmax-Tmin; // intervallo di tempi

Dt=T/N; // passo temporale

Dw=2*%pi/T; // passo delle pulsazioni

maxw=2*%pi/Dt; // pulsazione massima

// funzione fx campionata, che chiamo x

x=zeros(1,N);

t=zeros(1,N);

for n=1:N

    x(n)=fx(Tmin+(n-1)*Dt)

    t(n)=Tmin+(n-1)*Dt

end

// DFT della funzione fx, che chiamo X

X=zeros(1,N);

w=zeros(1,N);

for m=1:N

    for n=1:N

        X(m)=X(m)+x(n)*exp(-2*%pi*%i*(m-1)*(n-1)/N)

    end

    w(m)=(m-1)*Dw

end

// grafico in alto

subplot(2,1,1)

ts=linspace(Tmin,Tmax,5*N)

plot(t,x,"om") // campioni x (pallini magenta)

plot(ts,fx(ts),"-k") // funzione fx (linea nera)

xlabel("time [s]")

legend(["function fx","sampled x"])

// grafico in basso

subplot(2,1,2)

wshort=w(1:maxm); // mi limito alle prime maxm pulsazioni

Xshort=2*X(1:maxm)/N; // moltiplico la trasformata per 2/N per avere i coefficienti delle armoniche

plot2d3(wshort,abs(Xshort)) // barre nere verticali per i moduli della trasformata

xlabel("angular velocity [1/s]")

title("2*abs(X) / N     with     X=DFT(fx)")

plot(wshort,abs(real(Xshort)),"-+r") // rosso per il modulo della parte reale della trasformata

plot(wshort,abs(imag(Xshort)),"-+b") // blu per il modulo della parte immaginaria della trasformata

legend(["2*abs(X) / N","Re(2X/N)","Im(2X/N)"])

// tabella su finestra di console

absXshort=abs(Xshort);

ReXshort=real(Xshort);

ImXshort=imag(Xshort);

ordered=gsort(absXshort);

for i=1:maxm // ordino la tabella mettendo in alto le pulsazioni che hanno più peso

    nn(i)=find(absXshort==ordered(i))(1);

    wwshort(i)=wshort(nn(i));

    ReXshort(i)=ReXshort(nn(i));

    ImXshort(i)=ImXshort(nn(i));

end

titolo="Coefficienti della DFT"

nomi=["  pulsazione","2abs(X)/N","2Re(X)/N","  2Im(X)/N"]

dati=cat(2,wwshort,ordered',ReXshort',ImXshort')

disp(nomi)

disp(dati)