%==========================================
% Gibbs phenomenon
%==========================================
clear,close all
t1=10;dt=0.001;
t=-t1:dt:t1; 
T=2;T1=T/4;
x1=u(t+T1)-u(t-T1);x=0;
for i=-t1:t1;
    x=x+u(t-i*T+T1)-u(t-i*T-T1);
end
N=input('Type in the number of harmonic components N=');
%=========================================================
% Computation of the Fourier series coefficients
%=========================================================
L=2*N+1;
k=-N:N;
an=(2*T1/T)*sinc(2*k*T1/T);
%=====================
% Synthesis of signal
%=====================
y=0;
for q=1:L;
      y=y+an(q)*exp(j*(-N+q-1)*2*pi*t/T);
end;
%=================================
subplot(221)
plot(t,x)
axis([-T,T,-0.2,1.2])
title('The original signal x(t)')
subplot(212)
plot(t,y)
axis([-T,T,-0.2,1.2])
title('The synthesis signal y(t)')
xlabel('Time t')
subplot(222)
k=-N:N;
stem(k,an,'r.')
title('The amplitude ak of x(t)')
axis([-N,N,-0.5*max(abs(an)),1.2*max(abs(an))])