type = input('segnale: 1-impulso sinusoidale; 2-sin^2; 3-onda quadra; 4-triangolo di durata tau: ');
t0=4;
tstep=0.01;
t=-t0:tstep:t0; T=1;
nmax=input('termini sviluppo in serie: ');
n=1:nmax;
x=zeros(length(n),length(t));
a=zeros(1,length(n)); b=zeros(1,length(n));
switch type
    case 1
        %sin(2*pi*t/T)
        a0=1/pi;
        b(1)=-1/4;
        a(rem(n,2)==0)=-1./(pi*((n(rem(n,2)==0)).^2-1));
    case 2
        %(sin(2*pi*t/T))^2
        a0=1/4;
        a(2)=-1/8;
        b(rem(n,2)==1)=2./(pi*n(rem(n,2)==1).*((n(rem(n,2)==1)).^2-4));
    case 3
        %onda quadra
        a0=1/2;
        b(rem(n,2)==1)=-1./(pi*n(rem(n,2)==1));
    case 4
        %triangolo di durata tau
        tau=input('durata triangolo (<=1): ')*T;
        a0=tau/(2*T);
        a=sin(2*pi*n*tau/T)./(n*2*pi)-T/tau*(1-cos(2*pi*n*tau/T))./((n*2*pi).^2);
        b=cos(2*pi*n*tau/T)./(n*2*pi)-T/tau*sin(2*pi*n*tau/T)./((n*2*pi).^2);
end
if(1)
    for i=1:length(n), x(i,:)=2*a(i)*cos(2*pi*n(i)*t/T)-2*b(i)*sin(2*pi*n(i)*t/T); end
    figure(1)
    plot(t,a0+x(1,:),'k');
    hold on
    plot(t,a0+sum(x(1:2,:)),'r');
    plot(t,a0+sum(x(1:5,:)),'b');
    keyboard
    hold off
    figure(2)
    plot(t,a0+sum(x(1:length(n),:)),'k');
    keyboard
end
figure(3)
n1=[0 n];
b1=[0 b];
a1=[a0 a];
stem(n1,(b1.^2+a1.^2).^.5,'k')
hold on
set(gca,'YScale','log');