clear
hold off
format long e
N = 4096; %No. of FFT samples
f0=500;     %fundamental freq of input triangular wave
T0 = 1/f0;  %period 
T = T0/2; 50% duty cycle square wave
M=3; % In this example square wave has 2*M+1 cycles

fmax = N/(4*(M+0.5)*T0);
%fmax = 100;
sampling_rate = 2*fmax; %Nyquist rate 
tstep = 1/sampling_rate;
tmax = N*tstep/2;
tmin = -tmax;
tt = tmin:tstep:tmax-tstep;
%Square wave generation
gp = zeros(1,N);
for m= -M:M
    n=1;
    for t = tt
        gp(n) = gp(n) + rect((t-m*T0)/T);
        n=n+1;
    end
end


Hp1 = plot(tt,gp)
set(Hp1,'LineWidth',2)
Ha = gca;
set(Ha,'Fontsize',16)
title('input - time domain')
axis([tmin tmax 0 1.1])
pause
fmin = -fmax;
fstep = (fmax-fmin)/N;
freq = fmin:fstep:fmax-fstep;



Gpf1 = fft(gp);

Gpf = fftshift(Gpf1);
Hp1=plot(freq,abs(Gpf))
set(Hp1,'LineWidth',2)
Ha = gca;
set(Ha,'Fontsize',16)
Hx=xlabel('Frequency (Hz) ');
set(Hx,'FontWeight','bold','Fontsize',16)
Hx=ylabel('|Gp(f)|');
set(Hx,'FontWeight','bold','Fontsize',16)
title('magnitude spectrum Gp(f)');
axis([-10*f0 10*f0 0 2000])
pause
%hold on

MM=10;
for n = -MM:MM
    c(n+MM+1) = sinc(n/2)/2;
end
fvec = (-MM:MM)*f0;
plot(fvec,abs(c),'rx')
%plot(fvec,abs(c)*4000,'rx')