%% JVK 301 - Fall 2020 - Week 7 - Fourier Series of Pulse Trains

% See Oppenheim, Willsky, Young Chapter 3

close all 
clear
clc

%% Sinc Envelope
% We derived in class the generic FS coefficients and found them to be 
% samples of the sinc envelope. Here we plot the sinc envelope and explore
% its behavior.

% Below omega represents the "continuous variable", of course it must be
% sampled to plot. Therefore, we make its resolution higher than any other
% set of samples that might interest us.

epsilon = 0.001;
N = 5;
omega = [-N*pi:epsilon:N*pi];          %Note: value zero will
                                       %never be found in omega.
                                       
sinc = sin(omega) ./ omega;

figure(1)
plot(omega,sinc)
title('The Sinc Envelope')
xlabel('\omega (radians)')
grid on

saveas(gcf,'SincPlotV1','jpg')

% Since pi is an important value for the sinc function, it is often better
% to label the independent axis in units of pi. We can do this in the plot
% by a simple scaling.

figure(2)
plot(omega/pi,sinc)
title('The Sinc Envelope')
xlabel('\omega (in multiples of \pi radians)')
grid on

saveas(gcf,'SincPlotV2','jpg')

%% Plot the FS Coefficients and Experiment with Duty Cycle 

% The pulse train duty cycle delta = Delta/T is the ratio of the pulse
% width to the period of the pulse train. The duty cycle is a number
% between 0 and 1.
delta = 0.25;

khigh = floor(5/delta);
ksamples = [-khigh:1:khigh];
omegasamples = ksamples*(pi*delta);

% Compute the sinc function. Note that there will be a divide by zero,
% resulting in an NaN value at the point where the frequency variable
% equals zero. The limiting value of the sinc at that spot equals 1.0.
X = sin(omegasamples) ./ omegasamples;
X(khigh+1) = 1.0;                      %replace NaN by 1.0
X = delta*X;                           %FS coefficients are proportional to
                                       %duty cycle.

figure(3)
stem(ksamples,X)
title('FS Coeffs of Square Pulse Train for Duty Cycle = 0.25')
xlabel('k')
grid on

saveas(gcf,'FSCoeffs_delta=0.25_V1','jpg')

figure(4)
plot(omega,delta*sinc)
hold on
stem(omegasamples,X)
grid on
hold off
xlabel('\omega')
title('FS Coeffs of Square Pulse Train for Duty Cycle = 0.25, with Envelope')
saveas(gcf,'FSCoeffs_delta=0.25_V2','jpg')


%% Reproduce OW Figure 3.7 on page 195

% Hold the envelope constant and sample finer.
m = 1;

for d = [0.5, 0.25, 0.125]
    
    delta = d;
    
    khigh = floor(5/delta);
    ksamples = [-khigh:1:khigh];
    omegasamples = ksamples*(pi*delta);

    X = sin(omegasamples) ./ omegasamples;
    X(khigh+1) = 1.0;                      %replace NaN by 1.0
    X = delta*X;
    
    figure(5)
    subplot(3,1,m)
    stem(ksamples,X)
    xlabel('k')
    title(strcat('Fix Envelope - Sample Denser - FS Coeffs for \delta = ', num2str(delta)))

    grid on
    
    m = m+1;

end

saveas(gcf,'FixEnvelopeSampleDenser','jpg')

% Hold the number of FS coeffs to compute constant
m = 1;

khigh = floor(10/0.5);
ksamples = [-khigh:1:khigh];

for d = [0.5, 0.25, 0.125]
    
    delta = d;
    

    omegasamples = ksamples*(pi*delta);

    X = sin(omegasamples) ./ omegasamples;
    X(khigh+1) = 1.0;                      %replace NaN by 1.0
    X = delta*X;
    
    figure(6)
    subplot(3,1,m)
    stem(ksamples,X)
    xlabel('k')
    title(strcat('Hold number of FS Coeffs Constant - FS Coeffs for \delta = ', num2str(delta)))

    grid on
    
    m = m+1;

end

saveas(gcf,'FixNumberSamples','jpg')


%% Look at the Limiting Behavior at Duty Cycle -> 0

K = 5;

p = [0:1:K];
twopowers = 2 .^ (-p);
khigh = 20;
ksamples = [-khigh:1:khigh];

m = 1;

for d = twopowers
    
    delta = d;
    omegasamples = ksamples*(pi*delta);
    
    X = sin(omegasamples) ./ omegasamples;
    X(khigh+1) = 1.0;                      %replace NaN by 1.0
    XX = X;
    X = delta*X;
    
    figure(7)
    subplot(K+1,1,m)
    stem(ksamples,X)
    xlabel('k')
    title(strcat('FS Coeffs for \delta = ', num2str(delta)))

    grid on
    
    figure(8)
    subplot(K+1,1,m)
    stem(ksamples,XX)
    xlabel('k')
    title(strcat('FS Coeffs for \delta = ', num2str(delta)))

    grid on
    
    m = m+1;
    
end

figure(7)
set(gcf,'Position',[100 100 550 700])
saveas(gcf,'LimV1','jpg')
figure(8)
set(gcf,'Position',[100 100 550 700])
saveas(gcf,'LimV2','jpg')

