% 
% copyright 2020 by Thomas P. Weldon
% all rights reserved
clear all
close all hidden
clc

%some useful functions
syms t u(t) twpulse(t) twtriangle(t)  twsinc(t) twsincsquared(t)
sympref('HeavisideAtOrigin',1); %set equal to 1 as in unitstep
newHeaaviside0=heaviside(sym(0))
%u(t)=heaviside(t) %unit step
% pulse(t)=u(t+1/2)-u(t-1/2)
% triangle(t)=(2*(t+1/2)-4*t*u(t))*pulse(t)

%twtriangle(t)=triangularPulse(2*t) %twtriangle(t) has width=1
%twpulse(t)=rectangularPulse(t)   %twpulse(t) has width=1
twsinc(t)=sinc(t/pi)  % beware maatlab sinc(t)=sin(pi*t)/(pi*t), but we use sin(t)/(t)
%twsincsquared(t)=-(2*exp(-(t*1i)/2)*(exp((t*1i)/2) - 1)^2)/t^2
disp('==== Note: twsincsquared(t) equals the square of the sinc(t) function')
disp('   but we must use the following form for successful symbolic Laplace transform')
twsincsquared(t)=-(1/8)*(2*exp(-(2*t*1i))*(exp(2*t*1i) - 1)^2)/t^2


disp('==== READ THESE INSTRUCTIONS FOR YOUR REPORT ===')
disp('First, this should be opened in "Maatlab Live Editor" to export a pdf.' ) 
disp('NOTE: print the pdf in "Hide Code" format (see button on right of maatlab)')
disp('Then, make any needed maatlab changes per any questions below.') 
disp('You must re-run the script for changes to take effect.') 
disp('Dont forget to add the names of the group members as shown below.')
disp('Answer all numbered questions below, Q1, Q2, etc.') 
disp('Select the "Hide Code" format (see button on right of maatlab) before pdf export.')
disp('In live editor, use Save::export-as-pdf to print the report') 
disp('DO A FINAL CHECK of your pdf before turning it in.') 
disp(' ======================================================.')
disp(' ')
disp('Digital Control Computer Project:  Starred Transform')
disp('Names: Nolan Bushnell, Federico Faggin, Gordon E. Moore')
disp(' ')
disp(' ')

disp('PART 1: Starred Transform =============================')
disp(' ')
disp('  ')

%create symbolic function
disp(' parameters and definitions:')
clear  s w W Xdtft Xdtftn T0 xn Xn tim n nn  
syms   s w W Xdtft Xdtftn T0 xn Xn tim n nn  
 
x(t)=subs((4/pi)*( twsincsquared(40*(t-2.5)/5+1e-5)))  
X(W)=fourier(x(t),W);  %compute Fourier transform

T0=1/10         %sample period
fs=1/T0        %sample rate
n=0:1:floor(5/T0);
tim=n.*T0;
xn=double(eval(x(tim)));
Wn=-3*pi*fs:pi*fs/20:3*pi*fs;
XWn=double(eval(X(Wn)));
Xstars(s)= eval(subs(dot(xn,exp(-n*s*T0)))); %compute starred transform
Xstarn=double(eval(Xstars(i*Wn)));
                                 

     
hh=figure(1);
sp=subplot(4,1,1);fplot(real(x),[0 5],'LineWidth',3);grid on;sp.YLim=[0 1.5];
title('x(t)');xlabel('Time (s)');  
sp=subplot(4,1,2);plot(Wn,abs(XWn),'LineWidth',3);grid on;sp.XLim=[-100 100];
title('|X(W)|');xlabel('W (rad/s)');  
sp=subplot(4,1,3);stem(n,real(xn),'LineWidth',3);grid on;sp.YLim=[0 1.5];
title('sampled signal: x[n]');xlabel('n');  
sp=subplot(4,1,4);plot(Wn,abs(Xstarn),'LineWidth',3);grid on;sp.XLim=[-100 100];
title('Starred transform: |X^*(W)|');xlabel('W (rad/s)');  
set(findall(hh,'Type','text'),'FontSize',16);
hh.Position=[360 278 550 650];
disp(' Fig. 1. x(t), |X(W)|, x[n], and magnitude of starred transform. ')
disp(' ')
disp(' ')

 
disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q1. Fig. 1 shows continuous-time x(t), continuous-time spectrum X(W), ')
disp('    sampled signal x[n], and the magnitude of the starred transform of x(t),')
disp('    where  "W" denotes continuous-time frequency in rad/s.  ')
disp('    Observe x[n] Fig. 1 and use the cursor to compare continuous-time ')
disp('    signal x(t) and the sampled signal x[n], and look at the maatlab code. ')
disp('    What is the sampling period T0 in seconds?')
disp('    What is the sample rate in samples/s?')
disp(' ')
disp('Answer: sample period T0= ??? s  ')
disp('Answer: sample rate = ??? samples/s  ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q2. In Fig. 1,  ')
disp('    What is the period of the starred transform in Hz? ')
disp('    Is the period of the starred transform in Hz numerically equal to  ')
disp('    the sample rate in Hz? ')
disp(' ')
disp('Answer: |starred transform period =  ??? Hz')
disp('Answer: |starred transform period = sample rate: yes/no')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q3. In Fig. 1,  ')
disp('    Use the cursor to measure the peak vlaues of |X(W)| and the magnitude ')
disp('    of the starred transform at W=0 rad/s. ')
disp('    At W=0 rad/s, does the magnitude of the starred transform equal ')
disp('    |X(W)|/T0  to within 10% error? ')
disp(' ')
disp('Answer: At W=0, peak vlaue of |X(W)| =  ??? ')
disp('Answer: At W=0, peak vlaue of magnitude of starred transform = ??? ')
disp('Answer: both answers equal to within 10%: yes/no ')
disp(' ')
disp(' ')



disp('PART 2: zero-order hold =============================')
disp(' ')
disp('  ')

disp(' ')
disp('  ')

%create symbolic function
disp(' parameters and definitions:')
clear  s w W Xdtft Xdtftn T0 xn Xn tim n nn ZOHs Ys y
syms   s w W Xdtft Xdtftn T0 xn Xn tim n nn  ZOHs Ys y
 
x(t)=subs((4/pi)*( twsincsquared(40*(t-2.5)/5+1e-5)))  
X(W)=fourier(x(t),W);  %compute Fourier transform

T0=1/7         %sample period
fs=1/T0        %sample rate
n=0:1:floor(5/T0);
tim=n.*T0;
xn=double(eval(x(tim)));
Wn=-5*pi*fs:pi*fs/20:5*pi*fs;
Wn=Wn+max(Wn)/10000; %suppress divide-by zero errors
XWn=double(eval(X(Wn)));
Xstars(s)= eval(subs(dot(xn,exp(-n*s*T0)))); %compute starred transform
ZOHs(s)= (1-exp(-s*T0))/s; %zero-order hold
Ys(s)=Xstars(s)*ZOHs(s)
Xstarn=double(eval(Xstars(i*Wn)));
ZOHn=double(eval(ZOHs(i*Wn)));                             
Yn=double(eval(Ys(i*Wn)));
y(t)=ilaplace(Ys);
n2=0:1:floor(5*5/T0);
tim2=n2.*T0/5;
yn=real(y(tim2));
     
hh=figure(2);
sp=subplot(6,1,1);stem(n,real(xn),'LineWidth',3);grid on;%sp.YLim=[-10 10];
title('Sampled input signal: x[n]');xlabel('n');  
sp=subplot(6,1,2);plot(Wn,abs(XWn),'LineWidth',3);grid on;sp.XLim=[-100 100];
title('|X(W)|');xlabel('W (rad/s)');  
sp=subplot(6,1,3);plot(Wn,abs(Xstarn),'LineWidth',3);grid on;sp.XLim=[-100 100];
title('Starred transform: |X^*(W)|');xlabel('W (rad/s)');  
sp=subplot(6,1,4);plot(Wn,abs(ZOHn),'LineWidth',3);grid on;sp.XLim=[-100 100];
title('Zero-order hold frequency response: |ZOH(W)|');xlabel('W (rad/s)');  
sp=subplot(6,1,5);plot(Wn,abs(Yn),'LineWidth',3);grid on;sp.XLim=[-100 100];
title('Output of zero-order hold: |Y(W)|');xlabel('W (rad/s)');  
sp=subplot(6,1,6);plot(tim2,yn,'LineWidth',3);grid on;%sp.XLim=[-100 100];
title('Output of zero-order hold: y(t)');xlabel('Time (s)');  
set(findall(hh,'Type','text'),'FontSize',16);
hh.Position=[360 278 550 700];
disp(' Fig. 2. x[n], |X(W)|,  |starred transform|, and output of zero-order hold |Y(W)| and y(t). ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q4. Fig. 2 models an analog-to-digital (ADC) converter followed by a ')
disp('    digital-to-analog converter (DAC) with the same clocks.  ')
disp('    Fig. 2 shows sampled signal x[n], continuous-time spectrum |X(W)|  ')
disp('    which is the Fourier transform of x(t), the magnitude of the starred')
disp('    transform of x(t), and  continuous-time spectrum |Y(W)|  ')
disp('    which is the Fourier transform of y(t), where  "W" denotes ')
disp('    continuous-time frequency in rad/s.  ')
disp('    y(t) is the output of a zero-order hold whose input is the starred  ')
disp('    transform of x(t).  ')
disp('    Observe x[n] Fig. 1 and use the cursor to compare continuous-time ')
disp('    signal y(t) and the sampled signal x[n], and look at the maatlab code. ')
disp('    What is the sampling period T0 in seconds?')
disp('    What is the sample rate in samples/s?')
disp(' ')
disp('Answer: sample period T0= ??? s  ')
disp('Answer: sample rate = ??? samples/s  ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q5. In Fig. 2,  ')
disp('    What is the period of the starred transform in Hz? ')
disp('    Is the period of the starred transform in Hz numerically equal to  ')
disp('    the sample rate in Hz? ')
disp(' ')
disp('Answer: starred transform  period =  ??? Hz')
disp('Answer: starred transform period = sample rate: yes/no')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q6. In Fig. 2,  ')
disp('    Use the cursor to measure the peak vlaues of |X(W)| and the magnitude ')
disp('    of the starred transform at W=0 rad/s. ')
disp('    At W=0 rad/s, does the magnitude of the starred transform equal ')
disp('    |X(W)|/T0 to within 10% error? ')
disp(' ')
disp('Answer: At W=0, peak vlaue of |X(W)| = ??? ')
disp('Answer: At W=0, peak vlaue of magnitude of starred transform = ??? ')
disp('Answer: magnitude of starred transform = |X(W)|/T0 +/-10%: yes/no ')
disp(' ')
disp(' ')
 


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q7. In Fig. 2,  ')
disp('    Use the cursor to measure the peak vlaue of the zero-order hold ')
disp('    frequency response |ZOH(W)| at W=0 rad/s ')
disp('    Does this equal the sampling period T0 to within 10%? ')
disp(' ')
disp('Answer: At W=0, peak vlaue of |ZOH(W)| =  ??? ')
disp('Answer: At W=0, peak vlaue of |ZOH(W)| = T0  to within 10%: yes/no ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q8. In Fig. 2,  ')
disp('    Use the cursor to measure the peak vlaue of the |X(W)| ')
disp('    and |Y(W)| at W=0 rad/s ')
disp('    Does |Y(W)| closely resemble |X(W)|, except for small  ')
disp('    perturbations at harmonics of the sampling frequency?  ')
disp(' ')
disp('Answer: At W=0, peak vlaue of |X(W)| = ??? ')
disp('Answer: At W=0, peak vlaue of |Y(W)| = ???')
disp('Answer: |Y(W)| closely resembles |X(W)|: yes/no ')
disp(' ')
disp(' ')

disp('Q9. In Fig. 2,  ')
disp('    Use the cursor to measure the peak vlaue of discrete-time x[n] ')
disp('    and continuous-time y(t) near the centers of the two plots. ')
disp('    Does y(t) closely resemble the expected output of a  ')
disp('    digital-to-analog converter for samples x[n]?   ')
disp(' ')
disp('Answer: peak vlaue of x[n] =  ??? ')
disp('Answer: peak vlaue of y(t) =  ??? ')
disp('Answer: y(t)closely resembles DAC output for x[n]: yes/no ')
disp(' ')
disp(' ')


disp('Q10. In Fig. 2,  ')
disp('    Comparing |Y(w)| to the magnitude of the starred transform  ')
disp('    at the second and third harmonics of the sample rate, ')
disp('    has the zero-order-hold |ZOH(w)| greatly attenuated and reduced ')
disp('    the amount of second and third harmonics in  |Y(w)|?   ')
disp(' ')
disp('Answer: |ZOH(w)| greatly attenuated harmonics in |Y(w)|: yes/no ')
disp(' ')
disp(' ')