% 
% copyright 2021 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
twtriangle(t)=triangularPulse(2*t) %twtriangle(t) has width=1
twpulse(t)=rectangularPulse(t)   %twpulse(t) has width=1
twsinc(t)=sinc(t/pi)   % twsinc(x)=sin(x)/x, beware of pi in maatlab sinc
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: Laplace, Fourier, and Continuous-Time Systems')
disp('Names: Nolan Bushnell, Federico Faggin, Gordon E. Moore')
disp(' ')

disp('PART 1: Laplace Transform =============================')
disp(' ')
disp('Note: maatlab Laplace transform is one-sided Laplace transform  ')
disp('  ')

disp('  ')
disp('We will use capital W as capital omega for continuous-time frequency in rad/s')
clear s t %
syms s t
disp('Check Laplace transform of dirac delta function:')
LaplaceOfDirac=laplace(dirac(t),t,s)
disp('Check Laplace transform of 2 second delayed dirac delta function:')
LaplaceOfDiracDelayedByTwo=laplace(dirac(t-2),t,s)
disp('Check Laplace transform of unit step function u(t):')
LaplaceOfStep=laplace(u(t),t,s)
disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q1. Check that the three maatlab symbolic Laplace transforms above')
disp('    are correct (LaplaceOfDirac, LaplaceOfDiracDelayedByTwo,')
disp('    and LaplaceOfStep).')
disp(' ')
disp('Answer: three symbolic Laplace transforms are correct: True/False')
disp(' ')
disp(' ')


%create symbolic function
disp(' parameters and definitions:')
clear t w W Wn x XW XWn XWnDegrees s
syms  t w W Wn x XW XWn XWnDegrees s
%x(t)=sin(3*t)/3 - (2*cos(3*t))/9 + 2/9 %+0*exp(-t/2)  % add 1e-5 to avoid divide by zero at sinc(0)
%x(t)=ilaplace((s^2+2*s+16)/(s+0.5)/(s+0.1)/(s+1));
disp('Time function:')
x(t)=(100*exp(-t))/3 - (305*exp(-t/2))/4 + (527*exp(-t/10))/12
disp('Laplace transform:')
Xs(s)=laplace(x(t),t,s)
TFstr = char(subs(Xs));         % convert symbolic Laplace into a transfer function
s = tf('s');
eval(['sys = ',TFstr]);          %sys is the transfer function
disp('maatlab transfer function form: '); 
[poles,zeroes] = pzmap(sys); 
disp(['  poles=' num2str(poles') '  and zeroes=' num2str(zeroes') ])

hh=figure(1);
Wmin=0.01; Wmax=100
sp=subplot(3,1,1);fplot(x,[0 40],'LineWidth',3);grid on;%sp.XLim=[0 max(tim)];
title({'Laplace Transform Example' 'x(t)'});xlabel('Time (s)'); set(gca,'LineWidth',2);
sp=subplot(3,1,2);bplt=bodeplot(sys,'b');
setoptions(bplt,'PhaseVisible','off');setoptions(bplt,'MagVisible','on');
setoptions(bplt,'FreqUnits','rad/s','Grid','on','YLim',[-50 50],'XLim',[Wmin  Wmax]);
title('Bode Plot Magnitude of X(s)');%title({'X(W)', ['Line2 x=' num2str(3.3)]})
figax=findall(hh,'Type','Axes');set(figax(3),'LineWidth',2);
sp=subplot(3,1,3);bplt=bodeplot(sys,'b');
setoptions(bplt,'PhaseVisible','on');setoptions(bplt,'MagVisible','off');
setoptions(bplt,'FreqUnits','rad/s','Grid','on','YLim',[-270 0],'XLim',[Wmin  Wmax]);
figax=findall(hh,'Type','Axes');set(figax(2),'LineWidth',2);
title('Bode Plot Phase of X(s)');
% sp=subplot(4,1,4);p=pzplot(sys);                   %plot poles and zeroes
% title('X(s) Pole-Zero Plot in S-Plane');%title({'X(W)', ['Line2 x=' num2str(3.3)]})
hh.Position=[360 278 550 550];
set(findall(hh,'Type','line'),'LineWidth',3);
set(findall(hh,'Type','text'),'FontSize',14);
set(findall(hh,'Type','text'),'FontWeight','bold');
disp('           Fig. 1. Laplace transform of x(t). ')
disp(' ')
disp(' ')

 disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q2. Fig. 1 shows continuous-time x(t), continuous-time frequency spectrum')
disp('    magnitude |X(W)| is shown in the Bode plot of X(s), along with the phase ')
disp('    angle of X(W), where "W" denotes continuous-time frequency in rad/s.')
disp('    Use the cursor to measure the peak value of x(t)in Fig. 1 ')
disp('    and the time at which signal x(t) reaches its peak. ')
disp(' ')
disp('Answer: peak of x(t)= ???  ')
disp('Answer: time of peak of = ??? s  ')
disp(' ')
disp(' ')


 disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q3. In Fig. 1 what is the magnitude of the dc response of |X(W)| ')
disp('    in dB? ')
disp('    (approximate the dc value using the lowest frequency in the plot) ')
disp('    Convert this dB value to find |X(W)|.')
disp('    (Hint: dB=10*log10(power ratio) ')
disp('    Also,')
disp('    Approximate the area under x(t) by a triangle starting at t=0,')
disp('    reaching the peak from your previous question answer, ')
disp('    and ending at t=30.  ')
disp('    This area should approximate the dc response. ')
disp(' ')
disp('Answer: dc response of |X(W)|= ??? dB  ')
disp('Answer: converted from dB,|X(W)|= ???  ')
disp('Answer: area of triangle x(t) approximation= ??? ')
disp(' ')
disp(' ')


hh=figure(2);
opts=pzoptions;opts.XLim=[-5 5];opts.YLim=[-5 5];
p=pzplot(sys,opts);              %plot poles and zeroes
title('X(s) Pole-Zero Plot in S-Plane');set(gca,'LineWidth',2);
set(findall(hh,'Type','text'),'FontSize',14);
set(findall(hh,'Type','text'),'FontWeight','bold');
disp('           Fig. 2. Pole-zero plot of X(s). ')
disp(' ')
disp(' ')


 disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q4. Fig. 2 shows the s-plane poles of  the X(s) that was also plotted  ')
disp('    in the Bode plot of Fig. 1.  ')
disp('    Use the cursor to find the locations of the 3 poles in Fig. 2.')
disp('    Use the cursor to find the locations of the 2 zeroes in Fig. 2.')
disp(' ')
disp('Answer: three pole locations s= ???, s= ???, s= ???  ')
disp('Answer: two zero locations s= ???, s= ???  ')
disp(' ')
disp(' ')


disp('PART 2: Fourier Transform =============================')
disp(' ')
disp('We will take the approach that Fourier transform is simply the  ')
disp('Laplace transform with s=jW (since maatlab seems to do better ')
disp('symbolic computations with Laplace. ')
disp(' ')

clear t w W Wn x2 X2 X2n XWnDegrees s
syms  t w W Wn x2 X2 X2n XWnDegrees s
disp('Time-domain function: ')
x2(t)=7*exp(-t/2)*sin(5*t)
disp('Laplace transform: ')
X2s(s)=laplace(x2(t),t,s)
disp('Finding Fourier from Laplace yields: ')
X2(W)=X2s(j*W)              %Fourier transform 
fac=exp((log(Wmax)-log(Wmin))/(50));
Wn=Wmin*fac.^(0:1:50);      %logarithmic frequency scale
X2n=eval(X2(Wn));
X2nDegrees= unwrap(atan2(imag(X2n),real(X2n)))*180/pi;
disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q5. Check that the symbolic Fourier transform X2(W) of x2(t)')
disp('     above is correct.')
disp('    (you may need to search internet for a table of Fourier transforms)')
disp(' ')
disp('Answer: symbolic Fourier transform X2(W) is correct: True/False')
disp(' ')
disp(' ')

hh=figure(3);
sp=subplot(3,1,1);fplot(x2,[0 10],'LineWidth',3);grid on;%sp.XLim=[0 20];
title({'Fourier Transform Example' 'x2(t)'});xlabel('Time (s)');set(gca,'LineWidth',2);
sp=subplot(3,1,2);loglog(Wn,abs(X2n),'LineWidth',3);grid on;sp.XLim=[Wmin Wmax];
title({'Magnitude of Fourier Transform' '|X2(W)|'});xlabel('Frequency W (rad/s)'); set(gca,'LineWidth',2); 
sp=subplot(3,1,3);semilogx(Wn,X2nDegrees,'LineWidth',3);grid on;sp.XLim=[Wmin Wmax];
title('Angle of X2(W) (degrees)');xlabel('Frequency W (rad/s)');   set(gca,'LineWidth',2);
set(findall(hh,'Type','text'),'FontSize',16);
hh.Position=[360 278 550 550];
disp('           Fig. 3. Fourier transform of x2(t). ')
disp(' ')
disp(' ')


 disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q6. Fig. 3 shows continuous-time x2(t), continuous-time frequency')
disp('    spectrum magnitude |X2(W)| is shown  along with the phase ')
disp('    angle of X(W), where "W" denotes continuous-time frequency in rad/s.')
disp('    Use the cursor to measure the frequency of damped sinusoid of x2(t). ')
disp('    Use the cursor to measure the frequency of the peak of |X2(W)| . ')
disp(' ')
disp('Answer: frequency of x2(t) in Hz = ??? Hz  ')
disp('Answer: frequency of x2(t) in rad/s = ??? rad/s  ')
disp('Answer: frequency of peak of |X2(W)| in rad/s = ??? rad/s  ')
disp(' ')
disp(' ')



disp('PART 3: Continuous-Time Systems  =============================')
disp(' ')
disp(' ')



%create symbolic function
disp(' parameters and definitions:')
clear t w W Wn x XW XWn XWnDegrees s S Gps H gp h
syms  t w W Wn x XW XWn XWnDegrees s S Gps H gp h
disp('======= symbolic function form: '); 
Gps(s)=10/(s*(s+1))  
Hs(s)=sym(4/10) 
gp(t)=ilaplace(Gps(s),s,t)*u(t)
h(t)=ilaplace(Hs(s),s,t)%*u(t)
TFstrGp = char(subs(Gps)); 
TFstrH = char(subs(Hs)); 
s = tf('s');
%disp('======= transfer function forms: '); 
eval(['Gpsys = ',TFstrGp]);
eval(['Hsys = ',TFstrH]) ;
disp('======= System plant transfer function Gp(s): '); 
Gpsys=Gpsys
disp('======= System feedback transfer function H(s): '); 
Hsys =Hsys 
disp('======= System open loop gain transfer function Gol(s)=Gp(s)H(s): '); 
Golsys=series(Gpsys,Hsys)
disp('======= System closed loop gain transfer function Gcl(s): '); 
Gclsys = feedback(Gpsys, Hsys)

disp('Data check:  ')
abs(freqresp(Golsys,0.3))
abs(freqresp(Gclsys,0.3))

disp('System block diagram:  ')

%


hh=figure(4);
Wmin=0.1; Wmax=1000
sp=subplot(3,1,1);splt=stepplot(Gclsys,'b');grid on;%sp.XLim=[0 max(tim)];
title({'Closed-loop step response'});xlabel('Time (s)'); set(gca,'LineWidth',2);
sp=subplot(3,1,2);bplt=bodeplot(Golsys,'r:',Gclsys,'b',Gpsys,'c--');
setoptions(bplt,'PhaseVisible','off');setoptions(bplt,'MagVisible','on');
setoptions(bplt,'FreqUnits','rad/s','Grid','on','YLim',[-60 60],'XLim',[Wmin  Wmax]);
title({'Magnitude of open-loop (red-dotted),'  'closed-loop (blue), and plant (cyan dash)'});%title({'X(W)', ['Line2 x=' num2str(3.3)]})
figax=findall(hh,'Type','Axes');set(figax(3),'LineWidth',2);
sp=subplot(3,1,3);bplt=bodeplot(Golsys,'r:',Gclsys,'b',Gpsys,'c--');
setoptions(bplt,'PhaseVisible','on');setoptions(bplt,'MagVisible','off');
setoptions(bplt,'FreqUnits','rad/s','Grid','on','YLim',[-270 0],'XLim',[Wmin  Wmax]);
figax=findall(hh,'Type','Axes');set(figax(2),'LineWidth',2);
title({'Phase of open-loop (red-dotted),'  'closed-loop (blue), and plant (cyan dash)'});
% sp=subplot(4,1,4);p=pzplot(sys);                   %plot poles and zeroes
% title('X(s) Pole-Zero Plot in S-Plane');%title({'X(W)', ['Line2 x=' num2str(3.3)]})
hh.Position=[360 278 550 640];
set(findall(hh,'Type','line'),'LineWidth',3);
set(findall(hh,'Type','text'),'FontSize',14);
set(findall(hh,'Type','text'),'FontWeight','bold');
disp('           Fig. 4. System analysis. ')
disp(' ')
disp(' ')


 disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q7. Fig. 4 shows continuous-time step response, and the Bode plots of')
disp('    open-loop gain Gol(s)=Gp(s)H(s), closed-loop gain Gcl(s), and ')
disp('    plant function Gp(s). Feedback H(s) is not plotted.')
disp('    Change the parameters to Gps(s)= 100/(s*(s+11)) and H(s)=sym(1/10) to')
disp('    match the example in the lecture notes. ')
disp('    (The results in Fig. 4 should match the lecture notes)')
disp(' ')
disp(' ')


 disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q8. After making the changes in Q7, ')
disp('    use the cursor to measure the magnitude of the open-loop gain Gol(s)')
disp('    in dB at a frequency of 0.3 rad/s.  ')
disp('    Convert this dB value to find |Gol(s)|.')
disp('    Use the matlab command window to find . ')
disp('    (Hint: dB=10*log10(power ratio) ')
disp(' ')
disp('Answer: |Gol(s)| at 0.3 rad/s in dB = ??? dB  ')
disp('Answer: converted from dB, |Gol(s)| = ???   ')
disp(' ')
disp(' ')

 disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q9. After making the changes in Q7, ')
disp('    use the cursor to measure the magnitude of the closed-loop gain Gcl(s)')
disp('    in dB at a frequency of 0.3 rad/s.  ')
disp('    Convert this dB value to find |Gcl(s)|.')
disp('    (Hint: dB=10*log10(power ratio) ')
disp(' ')
disp('Answer: |Gcl(s)| at 0.3 rad/s in dB = ??? dB  ')
disp('Answer: converted from dB, |Gcl(s)|  = ???  ')
disp(' ')
disp(' ')

 disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q10. After making the changes in Q7, ')
disp('    having changed the feedback to H(s)=1/10,')
disp('    is the magnitude of the closed-loop gain |Gcl(s)| at at 0.3 rad/s ')
disp('    approximately equal to 1/|H(s)|?')
disp('    (Hint: dB=10*log10(power ratio) ')
disp(' ')
disp('Answer: at 0.3 rad/s |Gcl(s)| approx. equal to 1/|H(s)|: True/False')
disp(' ')
disp(' ')