% 
% 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))
 
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: Closed Loop Time Response')
disp('Names: Nolan Bushnell, Federico Faggin, Gordon E. Moore')
disp(' ')
disp(' ')


disp('PART 1: Closed Loop Time Response ======')
disp(' ')
disp('  ')


clear t n T0 z s w W Dz Gp Gc gp H gprimet ghprimet gprimen ghprimen Gprime Gprimez GHprime GHprimez Gclz Gclz Golstar ZOH e E Ez Estar c cn C Cz Cstar Dz 
syms  t n T0 z s w W Dz Gp Gc gp H gprimet ghprimet gprimen ghprimen Gprime Gprimez GHprime GHprimez Gclz Gclz Golstar ZOH e E Ez Estar c cn C Cz Cstar Dz

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

disp(' ')
% disp('==== digital filter D(z):  ')
% Dz(z)=(z+1)/z 
Dz(z)=sym(1) ; %no digital filter in this project
disp(' ')
disp('==== plant function Gp(s):  ')
Gp(s)=1/(s+2) 
disp(' ')
disp('==== feedback function H(s):  ')
H(s)=sym(1/4) 
disp(' ')
disp('==== zero-order hold ZOH(s):  ')
ZOH(s)= (1-exp(-s*T0))/s
disp(' ')
disp('====  Gc(s)=Gp(s)ZOH(s):  ')
Gc(s)=Gp(s)*ZOH(s)
disp(' ')
disp('====  Gprime(s)=Gp(s)/s:  ')
Gprime(s)=(1/s)*Gp(s)
disp(' ')
disp('====  GHprime(s)=Gp(s)H(s)/s:  ')
GHprime(s)=(1/s)*Gp(s)*H(s)
 
disp(' ')
disp('===========  computed results ==================:  ')
disp(' ')
disp('====  gprimet(t):  ')
gprimet(t)=simplify(ilaplace(Gprime(s)))
gprimen(n)=gprimet(n*T0);
disp(' ')
disp('====  Gprimez(z):  ')
Gprimez(z)=simplify(ztrans(gprimen(n),n,z),'Steps',5)
disp(' ')
disp('====  ghprimet(t):  ')
ghprimet(t)=simplify(ilaplace(GHprime(s)))
ghprimen(n)=ghprimet(n*T0);
disp(' ')
disp('====  GHprimez(z):  ')
GHprimez(z)=simplify(ztrans(ghprimen(n),n,z),'Steps',5)

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q1. FIRST,  ')
disp('    change the plant function Gp(s) above to Gp(s)=30/(5*s), ')
disp('    and change the feedback above to H(s)=sym(1/10), ')
disp('    as in example 2 of the time-response lecture notes.')
disp('    Hint: "sym(1/10)" declares a symbolic fraction.')
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('    does Gprime(s) above equal Gp(s) times 1/s?')
disp('    as in example 2 of the time-response lecture notes?')
disp(' ')
disp('Answer: Gprime(s) equals Gp(s) times 1/s: yes/no  ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q4. After making the changes in Q1,')
disp('    Does GHprime(s) above equal Gp(s)H(s) times 1/s?')
disp('    as in example 2 of the time-response lecture notes?')
disp('    Hint: you may need to rearrange numerator and denominator')
disp('    to show they are equal to lecture note results')
disp(' ')
disp('Answer: GHprime(s) above equal Gp(s)H(s) times 1/s: yes/no  ')
disp(' ')
disp(' ')


disp('====  Use sylbolic solver:  ')
syms  vCs vCz vGclz vGolz vEz Rz Rstar
eqns=[  vGolz==GHprimez(z)*(z-1)/z*Dz(z), ...
        vGclz==vCz/Rz, ...
        vEz==Rz-vEz*Dz(z)*GHprimez(z)*(z-1)/z, ...
        vCs==vEz*Dz(z)*Gp(s)*ZOH(s), ...
        vCz==vEz*Gprimez(z)*(z-1)/z*Dz(z), ...
        ]

sol=solve(eqns,[vCs vCz vGclz vGolz vEz ]);
disp('=================================================== ')
disp('Closed-loop results: ')
disp(' ')
disp('==== Closed-loop pulse Transfer function Gcl(z):  ')
Gclz(z)= collect(sol.vGclz,z)
Gclstar(s)=Gclz(exp(s*T0)); %starred transform if we need it
disp(' ')
disp('==== Closed-loop output Cz(z) for input R(z) denoted Rz:  ')
Cz(z)= simplify(collect(collect(Rz*sol.vGclz,z),Rz))
disp(' ')
disp('==== Closed-loop output C(s) for input R(z) denoted Rz:  ')
C(s,z,T0,Rz)= sol.vCs
disp(' ')
disp('==== Closed-loop output Cstar(s) for input R(z) denoted Rz:  ')
Cstar(s,z,T0,Rz)= subs(subs(Cz,z,exp(s*T0)),Rz,Rstar)

disp(' ')
disp('==== Sample period and rate:  ')
T0=1/10         %sample period
fs=1/T0        %sample rate

disp(['==== Closed-loop dc gain (z=1) from Gcl(z)= ' num2str(eval(Gclz(1))) ])
disp(['==== Closed-loop dc gain (s=0) from Gclstar(s)= ' num2str(eval(Gclstar(0))) ])

 disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q5. After making the changes in Q1,')
disp('    does the closed-loop pulse transfer function Gclz(z) above equal')
disp('    the closed-loop transfer function Gcl(z)=C(z)/R(z) in example 2  ')
disp('    of the time-response lecture notes?')
disp('    What is the dc (z=1) response of the closed-loop pulse-transfer ')
disp('    function Gclz(z) above?')
disp(' ')
disp('Answer: Gclz(z) same as lecture notes C(z)/R(z): yes/no  ')
disp('Answer: at z=1 Gclz(z)=: ???  ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q6. After making the changes in Q1,')
disp('    does Cz(z) equal C(z) in example 2  ')
disp('    of the time-response lecture notes?')
disp('    Note that Rz in the above result denotes R(z) from lecture notes')
disp(' ')
disp('Answer: Cz(z) same as lecture notes C(z): yes/no  ')
disp(' ')
disp(' ')
  
nn=0:1:floor(5/T0);
tim=nn.*T0;
Wn=-4*pi*fs:pi*fs/25:4*pi*fs;
Wn=Wn+1e-6; %suppress divide-by zero errors
wn=Wn*T0;

wmax=4*pi; wmin=-wmax;
wmax=pi;
wmin=wmax/1000;
fac=exp((log(wmax)-log(wmin))/(100));
wn=wmin*fac.^(0:1:100);      %logarithmic frequency scale
Wn=wn/T0;

GpWn=double(eval(Gp(i*Wn)));
%Estarn=double(eval(Estar(i*Wn)));
ZOHn=double(eval(ZOH(i*Wn))); 
Gclzn=eval(Gclz(exp(i*wn)));
GclznDegrees= unwrap(atan2(imag(Gclzn),real(Gclzn)))*180/pi;

Gclstarn=eval(Gclstar(i*Wn));
n2=0:1:4*max(nn);
tim2=n2.*T0/4;
% gp(t)=ilaplace(Gp);
% gpn=gp(tim2);

hh=figure(1);
sp=subplot(3,1,1);loglog(wn,abs(Gclzn),'LineWidth',3);grid on;
ymax=10^(floor(log10(max(abs(Gclzn))))+1);
sp.XLim=[wmin wmax];sp.YLim=[ymax/100 ymax ];
title({'|Gcl(w)|= |Gcl(z)| for z=e^{jw}'});xlabel('Frequency w (rad/sample)'); set(gca,'LineWidth',2); 
sp=subplot(3,1,2);semilogx(wn,GclznDegrees,'LineWidth',3);grid on;sp.XLim=[wmin wmax];sp.YLim=[-270 90 ];
title('Angle of Gcl(w) (degrees)');xlabel('Frequency w (rad/sample)');   set(gca,'LineWidth',2);
sp=subplot(3,1,3);loglog(Wn,abs(Gclstarn),'LineWidth',3);grid on;
ymax=10^(floor(log10(max(abs(Gclstarn))))+1);
sp.XLim=[wmin/T0 wmax/T0];sp.YLim=[ymax/100 ymax ];
title({'Starred transform |Gcl^*(W)|'});xlabel('Frequency W (rad/s)'); set(gca,'LineWidth',2); 
set(findall(hh,'Type','text'),'FontSize',16);
hh.Position=[360 278 550 650];
disp(' Fig. 1. Closed-loop pulse transfer function |Gcl(w)|, angle of Gcl(w), ')
disp('         and starred transform. ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q7. After making the changes in Q1,')
disp('    in Fig. 1 above, use the cursor tool to measure the following: ')
disp('    What is the 3 dB bandwidth of pulse transfer function |Gcl(w)|?')
disp('    What is the corresponding continuous-time bandwidth in rad/s?')
disp('    Hint: lowpass bandwidth is from zero rad/sample to the 3 dB point')
disp('    Hint: Gcl(w) is the ratio of ac-volts-out/ac-volts-in')
disp('    Hint: use sample rate to find corresponding continuous-time bandwidth')
disp(' ')
disp('Answer:  3 dB bandwidth of |Gcl(w)| = ??? rad/sample ')
disp('Answer:  corresponding continuous-time bandwidth = ??? rad/s ')
disp(' ')
disp(' ')


disp('  ')
disp('  ')
disp('==== Step response =================  ')
T0save=T0;          %save value of T0
clear T0; syms T0; %we need a symbolic result in T0

clear stepCs stepct stepctall
syms stepCs stepct stepctall

disp('==== Closed-loop output step response Cs(s) for r(t)=u(t) :  ')
stepCs(n,s,T0)=iztrans(partfrac(C(s,z,T0,z/(z-1)),z))
disp('==== Closed-loop output step response c(t) for r(t)=u(t) :  ')
stepct(n,t,T0)=ilaplace(stepCs(n,s,T0))
%stepct(t,T0)=ilaplace(stepCs(1,s,T0))
disp('==== Closed-loop output step response c[n] for r(t)=u(t) :  ')
stepcn(n,T0)=iztrans(Gclz(z)*z/(z-1))

T0=T0save; %restore T0

%the step response is the sum of each sample value times the delayed laplace
%sample values all equal one for a unit step input  
Digits = 32
stepctall(t)=sym(0)
for n2=0:max(nn)
   stepctall(t)=eval(stepctall(t)+heaviside(t-n2*T0)*subs(ilaplace(stepCs(n2,s,T0)),t,t-n2*T0));
end

%Digits = 32
stepctn=vpa(stepctall(tim2+T0/1e6));
stepcnn=eval(stepcn(nn,T0));
 
hh=figure(2);
sp=subplot(2,1,1);stem(nn,stepcnn,'LineWidth',3);grid on;%sp.XLim=[0 max(tim)];
title({'Second System Step Response  c[n]' });xlabel('n');  
sp=subplot(2,1,2);plot(tim2,stepctn,'LineWidth',3);grid on;%sp.XLim=[0 max(tim)];
title({'Step Response c(t)' });xlabel('n');  
% sp=subplot(4,1,3);loglog(Wn,abs(closedLoopGainStarWn),'LineWidth',3);grid on;%sp.XLim=[max(closedLoopGainStarWn)/100 max(closedLoopGainStarWn)];
% sp.YLim=[maxG/1000 maxG];
% title({'Magnitude of Closed Loop Gain Starred |C*(W)/R*(W)|'});xlabel('Frequency (rad/s)'); set(gca,'LineWidth',2); 
% sp=subplot(4,1,4);semilogx(Wn,closedLoopGainStarWnDegrees,'LineWidth',3);grid on;%sp.XLim=[Wmin Wmax];
% title('Angle of C*(W)/R*(W) (degrees)');xlabel('Time (s)');   set(gca,'LineWidth',2);
%hh.Position=[360 278 550 650];
set(findall(hh,'Type','text'),'FontSize',16);
disp('           Fig. 2. Closed-loop step response c[n] and c(t). ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q8. Fig. 2 above shows the closed-loop step response of the system ')
disp('    when the continuous-time input r(t) is a unit step,')
disp('    and c(t) is the continuous-time output (inverse Laplace of C(s)),')
disp('    and c[n] is the output at the sample instants (inverse z-transform')
disp('    of C(z)) ')
disp('    In Fig. 2 above, use the cursor tool to measure the following: ')
disp('    Does c[n] equal c(t) to within 5% at the sample instants? ')
disp(' ')
disp('Answer: c[n] equal c(t) to within 5%: yes/no ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q9. In Fig. 2 above, use the cursor tool to measure the rise time of')
disp('    c(t) from 10% to 90% of final value (from 1 volt to 9 volts)? ')
disp('    Does the final value of c(t) approximately equal the theoretical ?')
disp('    dc gain (Gclz(z) at z=1) ?')
disp(' ')
disp('Answer: rise time = ??? s ')
disp('Answer: final value approx. equal to dc gain: yes/no ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q10. Recall that the rise time of a simple lowpass RC filter is ')
disp('     approximately 2.2 divided by the 3 dB bandwidth in rad/s.')
disp('     Is the rise time measured in Q9 within 15% of 2.2 divided by  ')
disp('     the measured 3 dB bandwidth in Q7?')
disp(' ')
disp('Answer: yes/no ')
disp(' ')
disp(' ')