% This script file similates the application of the Nyquist Sampling
% Theorem to the following analog signal:
% x(t) = cos(100*pi*t) + cos(200*pi*t) + sin(500*pi*t) [1]
% for -t_max <=t <= t_max [2] with dt = 0.0001 [3].
% We assume that dt = 0.0001 approximates the continuoys time interval
% define by Equation 2.

% Mind that in the most abstract level a signal may be represented by the
% following mathematical form:
%                              S(t) = A * sin(2*pi*f*t+θ) [4]
%
% A: denotes the amplitude of the oscillation.
% F: denotes the frequency of the oscillation.
% t: denotes the independent time variable.
% Θ: denotes the phase of the oscillation.

% Initialy we need to determine the maximum frequency of oscillation
% present within the given analog signal. To this end, we need to express
% the given signal as a superposition of fundamnental sinusoidal signals
% given by Equation 4. Thus, we may write that:
%        ----
%        \
% x(t) = /       Sj(t) [5] such that Sj(t) = Aj * sin(2*pi*Fj*t+Θj) [6]                  
%        ----
%    j in {1,2,3}    

% According to the previous definitions we have for the various phases that:
%            pi
% Θ1 = Θ2 = ---- [7] and Θ3 = 0 [8]
%            2
% Moreover, for the various frenquencies we have that:
% F1 = 50Hz F2 = 100Hz and F3 = 250Hz [9]
% Therefore, we may write that Fmax = max{F1,F2,F3} = 250Hz [10]

% Thus, the minimum sampling frequency may be determined by:
%                                      1  
% Fs >= 2 * Fmax ==> Fs >= 500 Hz ==> ----- >= 500 ==> Ts <= 0.002 sec
%                                     Ts

% Original Signal Representation.
dt = 0.0001;
t_max = 0.1;
t = -t_max:dt:t_max;
x = cos(100*pi*t) + cos(200*pi*t) + sin(500*pi*t);
figure('Name','Original Signal Representation');
plot(t,x,'-r','LineWidth',1.3);
xlabel('-0.1 \leq t \leq +0.1');
ylabel('x(t)');
grid on

% Sampled Signal Representation with sampling period Ts = 0.002 witch is
% equal to the sampling period according to the Nyquist theorem.
% Therefore, signal can be effeciently reconstructed.
% x(t_max) = x(Nmax*Ts) ==> t_max = Nmax * Ts ==> Nmax = t_max / Ts
Ts = 0.002;
Nmax = t_max / Ts;
n1 = -Nmax:1:Nmax;
% Generate vector storing the samples of the original signal.
xs = cos(100*pi*n1*Ts) + cos(200*pi*n1*Ts) + sin(500*pi*n1*Ts);
% Initialize vector for storing the values of the reconstructed signal.
x1 = zeros(1,length(t));
% Generate the values of the reconstructed signal for each time instance
% within the initial time interval.
for k = 1:1:length(t)
    x1(k) = xs * sinc((t(k)-n1*Ts)/Ts)';
end
% Plot the reconstructed signal along with the original signal.
figure('Name','Reconstructed Signal for Ts = 0.002');
hold on
plot(t,x,'-r','LineWidth',1.5);
plot(t,x1,'*b','LineWidth',1.2);
xlabel('-0.1 \leq t \leq -0.1');
ylabel('x(t) ans x1(t)');
grid on


% Sampled Signal Representation with sampling period Ts = 0.001 witch is
% less than the maximum sampling period according to the Nyquist theorem.
% Therefore, signal can be effeciently reconstructed.
Ts = 0.001;
Nmax = t_max / Ts;
n2 = -Nmax:1:Nmax;
% Generate vector storing the samples of the original signal.
xs = cos(100*pi*n2*Ts) + cos(200*pi*n2*Ts) + sin(500*pi*n2*Ts);
% Initialize vector for storing the values of the reconstructed signal.
x2 = zeros(1,length(t));
% Generate the values of the reconstructed signal for each time instance
% within the initial time interval.
for k = 1:1:length(t)
    x2(k) = xs * sinc((t(k)-n2*Ts)/Ts)';
end
% Plot the reconstructed signal along with the original signal.
figure('Name','Reconstructed Signal for Ts = 0.001');
hold on
plot(t,x,'-r','LineWidth',1.5);
plot(t,x2,'*b','LineWidth',1.2);
xlabel('-0.1 \leq t \leq +0.1');
ylabel('x(t) and x2(t)');
grid on

% Sampled Signal Representation with sampling period Ts = 0.001 witch is
% greater than the maximum sampling period according to the Nyquist theorem.
% Therefore, signal cannot be effeciently reconstructed.
Ts = 0.01;
Nmax = t_max / Ts;
n3 = -Nmax:1:Nmax;
% Generate vector storing the samples of the original signal.
xs = cos(100*pi*n3*Ts) + cos(200*pi*n3*Ts) + sin(500*pi*n3*Ts);
% Initialize vector for storing the values of the reconstructed signal.
x3 = zeros(1,length(t));
% Generate the values of the reconstructed signal for each time instance
% within the initial time interval.
for k = 1:1:length(t)
    x3(k) = xs * sinc((t(k)-n3*Ts)/Ts)';
end
% Plot the reconstructed signal along with the original signal.
figure('Name','Reconstructed Signal for Ts = 0.01');
hold on
plot(t,x,'-r','LineWidth',1.5);
plot(t,x3,'*b','LineWidth',1.2);
xlabel('-0.1 \leq t \leq -0.1');
ylabel('x(t) ans x3(t)');
grid on