close all
clear all 
clc

Tp = 5; % period
d = .5; % duty cycle
om0 = 2*pi/Tp; % omega0
T = 0.001; % sampling spacing
N = 100; % number of Fourier coefficients

t = 0:T:Tp/2;
s = square_wave(t,Tp,d);
figure
plot(t,s)
grid
axis([1 1.35 .85 1.15])
xlabel('t')
ylabel('s(t)')
hold on

% true coefficients
tfs = zeros(size(t)); % truncated Fourier series
for k = -N:N % cicle on coefficients
    ak = d*sinc(k*d); % Fourier coefficient        
    tfs = tfs +  ak*exp(1i*k*om0*t);
end
tfs = real(tfs); % prevent numerical errors
plot(t,tfs)

% numerical coefficients
for M = [200,500,1000]
    tfs = zeros(size(t)); % truncated Fourier series
    nT = (0:M-1)*Tp/M;
    snT = square_wave(nT,Tp,d);
    for k = -N:N % cicle on coefficients
        bk = sum(snT.*exp(-1i*k*om0*nT))/M;
        tfs = tfs +  bk*exp(1i*k*om0*t);
    end
    tfs = real(tfs); % prevent numerical errors
    plot(t,tfs)
end
legend('signal','true coefs','M=200','M=500','M=1000')
title('approx truncated Fourier series')

% printing figure in png format
set(gcf,'PaperUnits','inches','PaperPosition',[0 0 6 4.5])
print -dpng exercise2.png -r100

function s = square_wave(t,Tp,d)
t1 = mod(t/Tp,1); 
s = rect(t1/d) + rect((t1-1)/d);
end

function s = rect(t)
s = (abs(t)<.5)+.5*(abs(t)==.5);
end