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


disp('PART 1: Analysis of Closed Loop Digital Control Systems ======')
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   
disp(' ')
disp('==== plant function Gp(s):  ')
Gp(s)=1/(s+2) 
disp(' ')
disp('==== feedback function H(s):  ')
H(s)=sym(1/10)
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('==== input e(t):  ')
% e(t)=heaviside(t)-heaviside(t-5*T0) 
% e(t)=heaviside(t) %sol
% en(n)=subs(e(n*T0));
% disp('==== z-transform of e(t) is denoted as Ez(z):  ')
% Ez(z)= simplifyFraction(ztrans(en))
% disp('==== input E(s) denoted as E(s):  ')
% E(s)=laplace(e)


 
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)=1/s, ')
disp('    and change the digital filter Dz(z) above to Dz(z)=1+1/z, ')
disp('    as in the closed-loop example 1 from the lecture notes.')
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('    does Gprime(s) above equal Gp(s) times 1/s?')
disp('    as in the closed-loop lecture notes?')
disp('    Does GHprime(s) above equal Gp(s) times 1/s times H(s)?')
disp(' ')
disp('Answer: Gprime(s) equals Gp(s) times 1/s: yes/no  ')
disp('Answer: GHprime(s) equals Gp(s) times 1/s times H(s): yes/no  ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q4. After making the changes in Q1,')
disp('    does Gprimez(z) above equal G`(z) in example 1 of ')
disp('       the closed-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: Gprimez(z) same as lecture notes: yes/no  ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q5. After making the changes in Q1,')
disp('    does GHprimez(z) above equal GH`(z) in example 1 of ')
disp('       the closed-loop lecture notes? ')
disp('       Here, GH(z) denotes the GcH with an overbar in 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: GHprimez(z) same as lecture notes: 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)*Gc(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 given input R(z) denoted Rz:  ')
%Cz(z)= collect(Rz*sol.vGclz,'Steps',5)
%Cz(z)= simplify(collect(collect(Rz*sol.vGclz,z),Rz))
Cz(z)= simplify(collect(collect(Rz*sol.vGclz,z),Rz))


disp(' ')
disp('==== Closed-loop output C(s) for given input R(z) denoted Rz:  ')
C(s)= simplify(collect(collect(sol.vCs,z),Rz))
C(s)= subs(simplify(collect(collect(sol.vCs,Rz),z)),z,exp(s*T0))

disp(' ')
disp('==== Closed-loop output Cstar(s) for given input R(z) denoted Rz:  ')
Cstar(s)= 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('Q6. 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 1  ')
disp('    of the closed-loop 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('Q7. After making the changes in Q1,')
disp('    what is the dc (s=0) response of the  plant function Gp(s) above?')
disp('    Does this equal the dc responce in Q6?')
disp(' ')
disp('Answer: at s=0 Gp(s)= ???  ')
disp('Answer: this is equal to dc response of Q6: yes/no  ')
disp(' ')
disp(' ')
 

nn=0:1:floor(2/T0);
tim=nn.*T0;
nn=0:1:floor(2/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/10000;
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:floor(5*2/T0);
tim2=n2.*T0/5;
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('Q8. 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 does |Gcl(w)| equal at the 3 dB point?')
disp('    What is the phase of Gcl(w) at the 3 dB point?')
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(' ')
disp('Answer:  3 dB bandwidth of |Gcl(w)| = ??? rad/sample ')
disp('Answer:  at 3 dB point |Gcl(w)| = ??? ')
disp('Answer:  at 3 dB point, phase of Gcl(w) = ??? degrees ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q9. 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 starred transform?')
disp('    Does this bandwidth in rad/s equal the answer in Q8 divided by ')
disp('    sample period T0?')
disp(' ')
disp('Answer: 3 dB bandwidth of starred transform = ??? rad/s ')
disp('Answer: equals Q8 answer divided by T0: yes/no ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q10. What is the theoretical phase at the 3 dB point of a')
disp('    bode plot of a function with a single pole, such as 1/(s+12)?  ')
disp('    Hint: do not forget to include the +/- sign ')
disp(' ')
disp('Answer: ??? degrees ')
disp(' ')
disp(' ')