close all hidden
clear
clc

set(groot, 'DefaultTextInterpreter', 'LaTeX');
set(groot, 'DefaultAxesTickLabelInterpreter', 'LaTeX');
set(groot, 'DefaultAxesFontName', 'LaTeX');
set(groot, 'DefaultAxesFontSize', 14);
set(groot, 'DefaultLegendInterpreter', 'LaTeX'); 
cmap = lines;

%**************************************************************************
% Basic paramters
%**************************************************************************
f_s_h = 100;    % Sampling freq for "continuous" signal in Hz
T_0 = 4;          % Fundamental period for normalization in s
f_1 = 1;           % Fundamental freq for sin signal in Hz
T_max = 8;      % Signal duration in seconds

%**************************************************************************
% Calculate support variables
%**************************************************************************
% Continuous signal --------------------------------------
N_h = T_max * f_s_h + 1;    % Duration in samples
t_h = linspace(0, T_max / T_0, T_max * f_s_h + 1); % Time scale in seconds

%**************************************************************************
% Signal selection (uncomment a single signal to use)
%**************************************************************************
% Single cos signal ---------------------------------------
v = cos(2*pi*f_1*t_h);

% Multi cos signal -----------------------------------------
% v = cos(2*pi*f_1*t_h) + cos(2*pi*2*f_1*t_h);

% Rect signal -----------------------------------------------
% v = ones(N_h, 1);
% v(ceil(N_h/2) : end) = v(ceil(N_h/2) : end) * -1;

%**************************************************************************
% Plots
%**************************************************************************
% Evaluate section by section to get plots after each other

figure
% set(gcf, 'Position', [2000 0 1400 1350]) % Set window size

subplot(2, 1, 1)
plot(t_h, v, 'LineWidth', 1)
hold on
grid on
xlabel('$\frac{t}{T_0}$', 'FontSize', 30);
ylabel('$v(t)$', 'FontSize', 20);
title('Continuous signal  $v(t)$', 'FontSize', 20)
%%
% Periodicity check ---------------------------------------
T_p = 1 * T_0;  

plot([0 T_p / T_0], [0.01 0.01], 'LineWidth', 3, 'Color', cmap(3, :))
legend('Cont. signal', 'Periodicity')
%%
% Highest frequency check -----------------------------
T_n = 1 * T_0;

plot([0 T_n / T_0], [-0.01 -0.01], 'LineWidth', 3, 'Color', cmap(4, :))
legend('Cont. signal', 'Periodicity', 'Highest freq. component')

%%
% Sampling -------------------------------------------------
T_s = 0.125 * T_0;
t_s = downsample(t_h, f_s_h * T_s);
u = downsample(v, f_s_h * T_s);
N_1 = length(t_s);
n_1 = 0 : N_1 - 1;

stem(t_s, u, 'LineWidth', 2, 'Color', cmap(2, :), 'LineWidth', 2)
legend('Cont. signal', 'Periodicity', 'Highest freq. component', 'Sampling points')

%%
% Reconstruction using SI interpolation ------------
y_1 = zeros(N_h, 1);

r = T_0 / T_s;
t_n = t_h * r;

for i = 0 : N_1 - 1
    y_1 = y_1 + sinc((t_n - i))' * u(i + 1);
end

subplot(2, 1, 2)
stem(n_1, u, 'LineWidth', 2, 'Color', cmap(2, :))
hold on
plot(t_n, y_1)
grid on
xlabel('$n$', 'FontSize', 20);
ylabel('$u(n)$', 'FontSize', 20);
title('Discrete signal  $u(n)$', 'FontSize', 20)
legend('Discrete signal', 'Reconstruction')