% Dr. Durant, EE3221, W1D3, 12/4/2020, revised 12/2/2021
% Demonstrate aliased sampling of sine.

% aliased_sampling

clearvars
close all

fs = 40; % sampling frequency, Hz
f = 30; % frequency of input sine (no antialias filter; AA filter would wipe out any component above fs/2)
N = 10;  % number of samples to illustrate
n = 0:N-1; % sample number
Ts = 1/fs; % sampling interval, s

assert(f <= 3*fs/2, 'This code only considers the n=±1 images, so it requires f <= 3*(fs/2)')

t_samp = n*Ts; % sampling time
y_dig = sin(2*pi*f*t_samp); % samples of sine

t = 0 : Ts/100 : t_samp(end); % approximate continuous time; 100 is arbitrary but should be large
y_ana = sin(2*pi*f*t); % values at highly sampled, continuous approximation

f_al = mod(f + fs/2, fs) - fs/2; % in Hz, find the image in the Nyquist range
% (which will be the original frequency if |f| < fs/2)

y_al_ana = sin(2*pi*f_al*t); % reconstructed (reconstruction filter limits to 0-fs/2 Hz, or 0-pi in digital frequency)

figure

% Plot all analog ("continuous" time) and digital (sampled) signals vs. time
subplot(211)
h = plot(t     ,y_ana    ,'b-',...
         t_samp,y_dig    ,'bo',...
         t     ,y_al_ana,'r--');
% h is array of handles (like pointers) to the 3 traces (series) in the plot
set(h, 'LineWidth', 4) % Make all lines thicker; can also modify select lines, e.g., h([1 3])
set(h(2), 'MarkerSize', 15) % Make the circles on the 2nd series bigger
legend('Analog source', 'Sampled', 'Analog reconstructed (may be aliased)')
xlabel('Time (s)'), ylabel('Signal level')
title('Aliasing Example')
set(gca, 'FontSize', 14)

% Plot frequency domain view
subplot(212)
h = plot(    [-f f], [1 1], 'ro', ...
          fs+[-f f], [1 1], 'g^', ...
         -fs+[-f f], [1 1], 'bv');
set(h,'LineWidth',4,'MarkerSize',15)
hold on
h = plot([-fs/2 -fs/2], [0 1], 'k-', ...
         [ fs/2  fs/2], [0 1], 'k-');
set(h,'LineWidth',10)
legend('Original', 'Periodic copy at fs (2\pi rad/sample)', 'Periodic copy at -fs (-2\pi rad/sample)', '-fs/2', 'fs/2')
xlabel('Frequency (Hz)'), ylabel('Signal level')
title('Frequency domain view - vertical bars indicate |f| < fs/2')
set(gca, 'FontSize', 14)