% File: "Lab4from3_01.m" 2-18-2014avo
% program derived from Couch example3_01.m matlab for use
% with EELE445 lab 4 on Natural Sampling
% function lab4

clear;
clc

fm =2.5;  %sine wave input frequency KHz
k=2;    %factor to oversample above nyquist rate
fs = fm*1.5*k;    %sample frequency KHz  choose Fs/Fm to be an integer
order=4;   %harmonics of Fs to consider in spectrum 
% f = -2*fm/fs:fm/fs:5*fs;
f = -2*fs/(k*fm):fm/(2*k*fs):order*(fs+fm);
d = 1/5;  %duty cycle of sampling waveform

%***************input 1+cos(2pi*2.5e3*t) input spectrum ******************

W = zeros(length(f),1);
for (i = 1:1:length(f))
   if (f(i) == 0.0)   %DC
    W(i) = 1;
   end;
   if (abs(f(i)) == fm)  %sinwave
   W(i) = .5;
   end;   
 end;


% d is the sample pulse width divided by the sampling period.
n = 0:1:order;
x = pi*n*d;
a = SA(x);

% Generating spectrum of Natural Pulse PAM signal.
Ws = zeros(length(f),1);
nfs=find(f==fs);  %find f=0 frequency
nfdc=find(f==0);    %find f=fs frequency
% Generating Shifted versions of W due to sampling.
for (i = 1:1:length(n))
  Wtemp = zeros(length(f),1);
  for j = n(i)*(nfs-nfdc)+1:1:length(W)
    Wtemp(j) = W(j - n(i)*(nfs-nfdc));
  end;
  Ws = Ws + d*a(i)*Wtemp;
end;

%******************** format screen output *******************
out=[f(:),Ws(:)];
disp(sprintf(['input frequency = ', num2str(fm), ' KHz \r' ]))
disp(sprintf(['sample frequency = ', num2str(fs), ' KHz \r']))
disp('   FREQ KHz    |CN|')
disp( out(find(abs(out(:,2))>0),:) ) %list only nonzero spectrum


subplot(211);
h=stem(f,abs(W));
set(h(1),'MarkerSize',3,'MarkerFaceColor','red','Marker','square')
xlabel('     f (KHz) -->');
ylabel('|W(f)|');
axis([-2*fm 2*fm 0 max(abs(W))]);
title('Cn Magnitude Spectrum for the AnalogSignal');

subplot(212);
h=stem(f,Ws);
set(h(1),'MarkerSize',3, 'MarkerFaceColor','red','Marker','square')
xlabel('f (KHz) -->');
ylabel('|Ws(f)|');
axis([-.25 max(f)  0 1.1*max(abs(Ws))]);
title(['Cn Magnitude Spectrum of Naturally Sampled PAM Signal for f>0, fs=',num2str(fs)]);
subplot(111);


% function y = SA(x)
% % File: SA.M
% % CALL: y = SA(x)
% % This function computes sin(x)/x.
% % x is assumed to be a "scalar" or a vector;
% 
% x = x(:);
% y = zeros(length(x),1);
% 
% for i = 1:length(x)
%   if (x(i) == 0)
%     y(i) = 1;
%   else
%     y(i) = sin(x(i))/x(i);
%   end;
% end;



