% 
% 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 tuening it in.') 
disp(' ======================================================.')
disp(' ')
disp('Digital Control Computer Project: Open Loop Digital Control Systems')
disp('Names: Nolan Bushnell, Federico Faggin, Gordon E. Moore')
disp(' ')
disp(' ')

disp('PART 1: Analysis of Open Loop Digital Control Systems =============================')
disp(' ')
disp('  ')


clear t n T0 z s w W Gp Gc gp gprimet gprimen Gprimez Golz Golstar ZOH e E Ez Estar c cn C Cz Cstar Dz 
syms  t n T0 z s w W Gp Gc gp gprimet gprimen Gprimez Golz Golstar ZOH e E Ez Estar c cn C Cz Cstar Dz

disp('===========  system parameters ==================:  ')


disp('==== plant function Gp(s):  ')
Gp(s)=1/(s+5)
disp('==== zero-order hold ZOH(s):  ')
ZOH(s)= (1-exp(-s*T0))/s
disp('====  Gc(s)=Gp(s)(1-exp(-sT0)/s:  ')
Gc(s)=((1-exp(-s*T0))/s)*Gp(s)
disp('====  Gprime(s)=Gp(s)/s:  ')
Gprimes(s)=(1/s)*Gp(s)
disp('==== Sample period and rate:  ')
T0=1/7         %sample period
fs=1/T0        %sample rate
% disp('==== D(z):  ')
Dz(z)=sym(1);  %we do not use digital filters in this project
disp('==== input e(t):  ') 
e(t)=heaviside(t)-heaviside(t-5*T0) 
disp('==== input E(s):  ')
E(s)=laplace(e)
en(n)=subs(e(n*T0));
disp('==== z-transform of e(t) is denoted as Ez(z):  ')
Ez(z)= simplifyFraction(ztrans(en))
 
 
disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q1. FIRST,  ')
disp('    change the plant function Gp(s) above to Gp(s)=1/(s+3), ')
disp('    as in the open-loop example 1 from the lecture notes.')
disp('    ALSO, ')
disp('    change the input signal e(t) above to be the unit step, ')
disp('    as in the open-loop example 1 from the lecture notes.')
disp('    Hint: see maatlab heaviside() function')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q2. After making the changes in Q1,')
disp('    does Gc(s) above equal the zero-order hold times the plant function')
disp('    Gp(s)?')
disp(' ')
disp('Answer: Gc(s) equals zero-order hold times Gp(s): yes/no  ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q3. After making the changes in Q1,')
disp('    note that the Laplace transform of e(t) is denoted as E(s).')
disp('    Does E(s) above equal Laplace transform of the unit step?')
disp(' ')
disp('Answer: E(s) equal Laplace transform of unit step: yes/no  ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q4. After making the changes in Q1,')
disp('    note that the z-transform of e(t) is denoted as Ez(z).')
disp('    Does Ez(z) above equal  the z-transform of the unit step?')
disp(' ')
disp('Answer: Ez(z) equals z-transform of a unit step: yes/no  ')
disp(' ')
disp(' ')

 
disp('===========  computed results ==================:  ')
T0save=T0;          %save value of T0
clear T0; syms T0; %we need a symbolic result in T0
disp('====  gprimet(t):  ')
gprimet(t)=simplify(ilaplace(Gprimes(s)))
gprimen(n)=gprimet(n*T0);
disp('====  Gprimez(z):  ')
%Gprimez(z)=simplify(ztrans(gprimen(n),n,z),'Steps',5)
Gprimez(z)=simplifyFraction(ztrans(gprimen(n),n,z))
T0=T0save; %restore T0


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q5. After making the changes in Q1,')
disp('    does Gprime(s) above equal Gp(s) times 1/s?')
disp('    Does Gprimez(z) above equal G`(z) in example 1 of ')
disp('     the open-loop lecture notes?')
disp('    Hint: you may need to multiply numerator and denominator')
disp('    by exp(-3T0) or some similar thing to verify the equality')
disp(' ')
disp('Answer: Gprime(s) equals Gp(s) times 1/s: yes/no  ')
disp('Answer: Gprimez(z) same as lecture notes: yes/no  ')
disp(' ')
disp(' ')



disp('====  Use sylbolic solver:  ')
T0save=T0;          %save value of T0
clear T0; syms T0; %we need a symbolic result in T0
syms  vCs  vGolz EstarIn 
eqns=[  vGolz==Gprimez(z)*(z-1)/z*Dz(z), ...
        vCs==EstarIn*Gc(s), ...
        ]

sol=solve(eqns,[ vGolz  vCs])
T0=T0save; %restore T0
disp('====  End sylbolic solver:  ')

disp(' ')
disp(' ')

disp('=================================================== ')
disp('Open-loop pulse transfer functions ')

disp('==== Open-loop pulse Transfer function Gol(z):  ')
Golz(z)= simplify(sol.vGolz,'Steps',5)
Golstar(s)=Golz(exp(s*T0)); 
disp('==== Open-loop output C(s) for generic input E(s) with star transform Estar(s):  ')
subs(simplify(sol.vCs,'Steps',5))
% disp('==== open-loop Cstar(s) = Estar*Gol(z) with z=exp(s*T0):  ')
% Cstar(s,Estar)=subs(subs(Estar*sol.vGolz,z,exp(s*T0)))  
% disp('==== Open-loop pulseTransfer Cstar(s)/Estar(s)=Gol(z) with z=exp(s*T0):  ')
% subs(Golz(z),z,exp(s*T0))
disp('==== Estar(s) for above given e(t):  ')
Estar(s)=Ez(exp(s*sym(T0)))
disp('==== Open-loop output C(s) for given Estar(s):  ')
C(s)= simplify(subs(simplify(sol.vCs,'Steps',5),EstarIn,Estar))
disp('==== Open-loop output Cz(z) for given input Ez(z):  ')
Cz(z)= simplify(Ez(z)*Golz(z),'Steps',5)
disp('==== Open-loop output c[n]  given input Ez(z):  ')
cn(n)= iztrans(Cz)
disp(' ')
disp(' ')

disp(['==== Unit step response c[n] at n=100 is ' num2str(eval(cn(100))) ])
disp(['==== Openloop dc gain (z=1) from Gol(z)= ' num2str(eval(Golz(1))) ])
disp(['==== Openloop dc gain (s=0) from Golstar(s)= ' num2str(eval(Golstar(0))) ])
disp(' ')
disp(' ')



disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q6. After making the changes in Q1,')
disp('    does the open-loop pulse transfer function Golz(z) equal')
disp('    the open-loop transfer function in example 1 of the ')
disp('    open-loop lectuer notes?')
disp('    What is the dc (z=1) response of the open-loop pulse-transfer ')
disp('    function Golz(z) above?')
disp(' ')
disp('Answer: Golz(z) same as lecture notes: yes/no  ')
disp('Answer: at z=1 Golz(z)=:  ???  ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q7. After making the changes in Q1,')
disp('    what is the dc (s=0) response of the  plant function Gp(s) above?')
disp('    Does the dc response of the open-loop pulse transfer function Golz(z)  ')
disp('    above equal the dc response of the plant Gp(s) to within 10%?')
disp(' ')
disp('Answer: at s=0 Gp(s)=: ???  ')
disp('Answer:  dc value of Golz(z)= dc value of Gp(s) +/-10% =: yes/no  ')
disp(' ')
disp(' ')
 



nn=0:1:floor(2/T0);
tim=nn.*T0;
en=double(eval(e(tim)));
Wn=-4*pi*fs:pi*fs/25:4*pi*fs;
Wn=Wn+1e-6; %suppress divide-by zero errors
wn=Wn*T0;
GpWn=double(eval(Gp(i*Wn)));
%Estarn=double(eval(Estar(i*Wn)));
ZOHn=double(eval(ZOH(i*Wn))); 
GolzWn=eval(Golz(exp(i*wn)));
Golstarn=eval(Golstar(i*Wn));
Cn=double(eval(C(i*Wn)));
c(t)=ilaplace( subs(C));
n2=0:1:floor(5*2/T0);
tim2=n2.*T0/5;
ctn=real(eval(c(tim2)));
cnn=real(eval(cn(nn)));
gp(t)=ilaplace(Gp);
gpn=gp(tim2);


hh=figure(1);
sp=subplot(4,1,1);plot(Wn,abs(GpWn),'LineWidth',3);grid on;sp.XLim=[-max(Wn) max(Wn)];
title('Plant function |Gp(s)| at s=jW');xlabel('W (rad/s)');  grid on;
maxy=floor(10*max(abs(GpWn)))/10+0.1;sp.YLim=[0 maxy];
sp=subplot(4,1,2);plot(Wn,abs(ZOHn),'LineWidth',3);grid on;sp.XLim=[-max(Wn) max(Wn)];
title('Zero-order hold |ZOH(s)| at s=jW');xlabel('W (rad/s)');  
maxy=floor(10*max(abs(ZOHn)))/10+0.1;sp.YLim=[0 maxy];
sp=subplot(4,1,3);plot(wn,abs(GolzWn),'LineWidth',3);grid on;sp.XLim=[-max(wn) max(wn)];
maxy=floor(10*max(abs(GolzWn)))/10+0.1;sp.YLim=[0 maxy];
title('Open-loop pulse transfer function |Golz(z)| at z=e^j^w');xlabel('w (rad/sample)');  
sp=subplot(4,1,4);plot(Wn,abs(Golstarn),'LineWidth',3);grid on;sp.XLim=[-max(Wn) max(Wn)];
maxy=floor(10*max(abs(Golstarn)))/10+0.1;sp.YLim=[0 maxy];
title('Open-loop starred transform |Gol^*(s)|=|Golz(e^s^T^o)| at s=jW');xlabel('W (rad/s)');  
set(findall(hh,'Type','text'),'FontSize',16);
hh.Position=[360 278 550 550];
disp(' Fig. 1. Plant function |Gp(s)|, Zero-order hold |ZOH(s)|, open-loop  ')
disp('         pulse transfer function |Golz(z)|, and starred transform |Gol^*(s)|')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q8. After making the changes in Q1,')
disp('    in Fig. 1 above use the cursor to measure the following:  ')
disp('    What is the 6 dB bandwidth of the plant function |Gp(s)|?')
disp('    What is the 6 dB bandwidth of the zero-order hold |ZOH(s)|?')
disp('    What is the 6 dB bandwidth of starred transform |Golstar(s)|?')
disp('    Hint: lowpass bandwidth is from zero rad/s to the 6 dB point')
disp('    Hint: Gp(s) is the ratio of ac-volts-out/ac-volts-in')
disp(' ')
disp('Answer:  6 dB bandwidth of |Gp(s)|= approx ??? rad/s ')
disp('Answer:  6 dB bandwidth of |ZOH(s)|= approx ??? rad/s ')
disp('Answer:  6 dB bandwidth of |Golstar(s)|= approx ??? rad/s ')
disp(' ')
disp(' ')
 

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q9. After making the changes in Q1,')
disp('    in Fig. 1, ')
disp('    is the bandwidth of|Gp(s)| or the bandwidth of |ZOH(s)|')
disp('    the dominant factor in the overall open-loop bandwidth')
disp('    observed in the starred transform |Golstar(s)|?')
disp(' ')
disp('Answer: bandwidth of|Gp(s)| or |ZOH(s)| dominates: |Gp(s)| or |ZOH(s)|    ')
disp(' ')
disp(' ')

     
hh=figure(2);
sp=subplot(4,1,1);stem(nn,real(en),'LineWidth',3);grid on;sp.XLim=[-max(nn) max(nn)];
title('Sampled input signal: e[n]');xlabel('n');  
sp=subplot(4,1,2);plot(tim2,abs(gpn),'LineWidth',3);grid on;sp.XLim=[-max(tim) max(tim)];
title('Plant impulse response gp(t)');xlabel('W (rad/s)');  
sp=subplot(4,1,3);plot(tim2,ctn,'LineWidth',3);grid on;sp.XLim=[-max(tim) max(tim)];
title('Output response c(t)');xlabel('Time (s)');  
sp=subplot(4,1,4);stem(nn,real(cnn),'LineWidth',3);grid on;sp.XLim=[-max(nn) max(nn)];
title('Sampled output response: c[n]');xlabel('n');  
set(findall(hh,'Type','text'),'FontSize',16);
hh.Position=[360 278 550 700];
disp(' Fig. 1. Sampled input signal e[n], plant impulse response gp(t),. ')
disp('        continuous-time output response c(t), and sampled output c[n]. ')
disp(' ')
disp(' ')



disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q10. After making the changes in Q1,')
disp('    in Fig. 1, use the cursor to measure: ')
disp('    What is the decay time (90% to 10%) in gp(t)?')
disp('    What is the rise time (10% to 90%) of c(t)?')
disp('    Is c[n] equal to the samples of c(t), such that c[n]=c(n*T0)? ')
disp(' ')
disp('Answer: decay time (90% to 10%) in gp(t): ??? s   ')
disp('Answer: rise time (10% to 90%) of c(t): ??? s   ')
disp('Answer: does  c[n]=c(n*T0) +/-10%: yes/no  ')
disp(' ')
disp(' ')