clc
clear all
close all

% Number of coefficients to plot
N = 5;
% Number of coefficients to use for `fs`
M = 500;
% Period length
T = 2;
% Fundamental frreqency
w0 = 2*pi/T;
% Number of sample points for plotting
num_p = 1000*T;
% Number of periods to display
P = 5;

% Time sampling points for one period
t = linspace(-T/2, T/2, num_p);
% Sampling points for P periods
tp = linspace(-P*T/2, P*T/2, P*num_p);
% Original sawtooth function that is odd
% Created by replicating the `t` vector P times into the single row
f = repmat(t, 1, P);

% Storage of coefficients
bk = zeros(M,1);
% Coefficient `b_0` is zero, hence the start of the loop at index 2.
% Remember that Matlab starts indexing at 1 and not at 0 ...
for i = 2:M
    % Index `k` goes from 0 to M-1
    k = i-1;
    bk(i) = (-1)^(k+1)*T/(k*pi);
end

% Plot first N components of the sin/cos sum.
% As our sawtooth is an odd function, we have only sin() components
figure(1);
for i = 1:N
    % Index `k` goes from 0 to M-1
    k = i-1;
    bf = bk(i)*sin(k*w0*t);
    subplot(N,1,i);
    plot(t, bf);
    legend(sprintf('b_%d', k));
    ylim([-0.7, 0.7]);
end

% Sum all M components together
fs = zeros(1, P*num_p);
% Store intermediate results for 5,15,25, and 55 components
f05 = [];
f15 = [];
f25 = [];
f55 = [];
for i = 1:M
    % Index `k` goes from 0 to M-1
    k = i-1;
    fs = fs + bk(i)*sin(k*w0*tp);
    % Store the intermediate results
    if i == 5
        f05 = fs;
    elseif i == 15
        f15 = fs;
    elseif i == 25
        f25 = fs;
    elseif i == 55
        f55 = fs;
    end
end

% Plot the partial sum of M components over P periods
figure(2);
plot(tp, f, tp, fs);
title(sprintf('Fourier series representation with %d components', M));
legend('original', 'Fourier series');

% Plot the intermediate results over P periods
figure(3);
plot(tp, f, tp, f05, tp, f55);
title(sprintf('Fourier series representation with %d components', M));
legend('original', '5 element partial sum', '55 element partial sum');
