% 
% 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('ECGR 4124 Computer Project: z-transform ')
disp('Names: Nolan Bushnell, Federico Faggin, Gordon E. Moore')
disp(' ')

disp('PART 1: z-transform   =============================')
disp(' ')
disp('  ')


%create symbolic function
disp(' parameters and definitions:')
clear s z w W x Xz Xdtft n
syms s z w W x Xz Xdtft n
x(n)=n*(heaviside(n-1)-heaviside(n-4))
disp('***** z-transform')
Xz(z)=simplify(ztrans(x,z),'Steps',5)
nn=0:20;
xn=x(nn);
Xdtft(w)= eval(subs(Xz(exp(i*w)))); %compute DTFT
wmax=2*pi
wmin=-2*pi
% %logarithmic frequency scale
%wmin=wmax/10000
% fac=exp((log(wmax)-log(wmin))/(50));
% wn=wmin*fac.^(0:1:50);   
wn=wmin:wmax/100:wmax;  %linear freq scale
Xdtftn=double(eval(Xdtft(wn)));
%HdtftnDegrees= unwrap(atan2(imag(Hdtftn),real(Hdtftn)),5)*180/pi;
XdtftnDegrees= (atan2(imag(Xdtftn),real(Xdtftn)))*180/pi;

TFstr = char(subs(Xz)); 
z = tf('z',1)
disp('======= transfer function form: '); 
eval(['Xzsys = ',TFstr]); 
[poles,zeroes] = pzmap(Xzsys);
disp(['  poles=' num2str(poles') '  and zeroes=' num2str(zeroes') ])

 
hh=figure(1);
sp=subplot(3,1,1);stem(nn,xn,'LineWidth',3);grid on;%sp.YLim=[-10 10];
title({'Z-Transform Example' 'x[n]'});xlabel('n');  
sp=subplot(3,1,2);plot(wn,abs(Xdtftn),'LineWidth',3);grid on;
%loglog(wn,abs(Hdtftn),'LineWidth',3);grid on;
%ymax=10^(floor(log10(max(abs(Hdtftn))))+1);
sp.XLim=[wmin wmax];%sp.YLim=[ymax/1000 ymax ];
title({'|X(w)|= |X(z)| for z=e^{jw}'});xlabel('Frequency w (rad/sample)'); set(gca,'LineWidth',2); 
sp=subplot(3,1,3);plot(wn,XdtftnDegrees,'LineWidth',3);grid on;sp.XLim=[wmin wmax];
%semilogx(wn,HdtftnDegrees,'LineWidth',3);grid on;
% sp.XLim=[wmin wmax];
title('Angle of X(w) (degrees)');xlabel('Frequency w (rad/sample)');   set(gca,'LineWidth',2);
set(findall(hh,'Type','text'),'FontSize',16);
hh.Position=[360 278 550 550];
disp(' Fig. 1. Signal x[n], |X(w)|= |X(z)| for z=e^{jw}, and angle of X(w). ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q1. Fig. 1 shows discrete-time x[n], DTFT |X(w)| = |X(z)| for z=exp(jw), ')
disp('    and the angle of X(w) in degrees, where "w" denotes discrete-time')
disp('    frequency in rad/sample.   ')
disp('    Change x[n], so that x[n] = u[n] - u[n-4] to match the DTFT example')
disp('    in the lecture notes (x[n]={1,1,1,1}). ')
disp('    (Fig. 1 should change accordingly) ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q2. After making the changes in Q1, observe Fig. 1,')
disp('    and using the cursor, what is the measured magnitude of the dc value ')
disp('    of DTFT |X(w)| at w=0 rad/sample?')
disp('    Also, compute the sum of values of x[n] in Fig. 1 from n=0 to n=20 ')
disp(' ')
disp('Answer: dc value of |X(w)|= ??? ')
disp('Answer: sum of x[n] = ??? ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q3. After making the changes in Q1, observe Fig. 1,')
disp('    and using the cursor, what is the frequency in rad/sample of the  ')
disp('    first null of |X(w)| in rad/sample?  ')
disp('    Hint: first null is the positive freq. where |X(w)| first equals zero. ')
disp('    What is the corresponding continuous-time frequency in MHz,')
disp('    if the sytem has a sample rate of 10 megasample/s?')
disp(' ')
disp('Answer: first null frequency: ??? rad/sample  ')
disp('Answer: corresponding continuous-time frequency = ??? MHz  ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q4. After making the changes in Q1, observe Fig. 1,')
disp('    and using the cursor, measure |X(w)| at w=3.14 rad/sample. ')
disp('    Is this value the same as the theoretical value in the lecture notes?')
disp(' ')
disp('Answer: at w=3.14 |X(w)| = ???  ')
disp('Answer: value same as theoretical: yes/no  ')
disp(' ')
disp(' ')


disp('PART 2: z-transform poles and zeroes  =============================')
disp(' ')
disp('  ')


%create symbolic function
disp(' parameters and definitions:')
clear s z w W x Hz Hdtft n
syms s z w W x Hz Hdtft n
%h(n)=1*(heaviside(n-1)-heaviside(n-4))
h(n)=1*(1*((4/5)^n) + 2*((1/3)^n))  %sol
disp('***** z-transform')
Hz(z)=simplify(ztrans(h,z),'Steps',5)
nn=0:20;
hn=h(nn);
Hdtft(w)= eval(subs(Hz(exp(i*w)))); %compute DTFT
wmax=2*pi
wmin=-2*pi
% %logarithmic frequency scale
%wmin=wmax/10000
% fac=exp((log(wmax)-log(wmin))/(50));
% wn=wmin*fac.^(0:1:50);   
wn=wmin:wmax/100:wmax;  %linear freq scale
Hdtftn=double(eval(Hdtft(wn)));
%HdtftnDegrees= unwrap(atan2(imag(Hdtftn),real(Hdtftn)),5)*180/pi;
HdtftnDegrees= (atan2(imag(Hdtftn),real(Hdtftn)))*180/pi;

TFstr = char(subs(Hz)); 
z = tf('z',1)
disp('======= transfer function form: '); 
eval(['Hzsys = ',TFstr]); 
[poles,zeroes] = pzmap(Hzsys);
disp(['  poles=' num2str(poles') '  and zeroes=' num2str(zeroes') ])

 
hh=figure(2);
sp=subplot(3,1,1);stem(nn,hn,'LineWidth',3);grid on;%sp.YLim=[-10 10];
title({'Z-Transform Example' 'h[n]'});xlabel('n');  
sp=subplot(3,1,2);plot(wn,abs(Hdtftn),'LineWidth',3);grid on;
%loglog(wn,abs(Hdtftn),'LineWidth',3);grid on;
%ymax=10^(floor(log10(max(abs(Hdtftn))))+1);
sp.XLim=[wmin wmax];%sp.YLim=[ymax/1000 ymax ];
title({'|H(w)|= |H(z)| for z=e^{jw}'});xlabel('Frequency w (rad/sample)'); set(gca,'LineWidth',2); 
sp=subplot(3,1,3);plot(wn,HdtftnDegrees,'LineWidth',3);grid on;sp.XLim=[wmin wmax];
%semilogx(wn,HdtftnDegrees,'LineWidth',3);grid on;
% sp.XLim=[wmin wmax];
title('Angle of H(w) (degrees)');xlabel('Frequency w (rad/sample)');   set(gca,'LineWidth',2);
set(findall(hh,'Type','text'),'FontSize',16);
hh.Position=[360 278 550 550];
disp(' Fig. 2. Signal h[n], |H(w)|= |H(z)| for z=e^{jw}, and angle of H(w). ')
disp(' ')
disp(' ')



disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q5. Fig. 2 shows discrete-time h[n], DTFT |H(w)| = |H(z)| for z=exp(jw), ')
disp('    and the angle of H(w) in degrees, where "w" denotes discrete-time')
disp('    frequency in rad/sample.   ')
disp('    Observe Fig. 2,')
disp('    and using the cursor, what is the measured magnitude of the dc value ')
disp('    of DTFT |H(w)| at w=0 rad/sample?')
disp('    Also, compute the sum of values of h[n] in Fig. 2 from n=0 to n=20 ')
disp(' ')
disp('Answer: dc value of |H(w)|= ???  ')
disp('Answer: sum of h[n] = ???   ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q6. Hz(z) above is the symbolic z-transform of symbolic function h[n]. ')
disp('    Hzsys above is the system transfer function form of H(z).')
disp('    Comparing Hz(z) to Hzsys, does  Hz(z) equal Hzsys?   ')
disp('    Hint: you might have to rearrange them to verify if they are equal')
disp(' ')
disp('Answer: does  Hz(z) equal Hzsys: yes/no  ')
disp(' ')
disp(' ')


hh=figure(3);
%pzm=pzmap(Hzsys)
pzp=pzplot(Hzsys);
setoptions(pzp,'FreqUnits','auto','Grid','off','YLim',[-1.5 1.5],'XLim',[-1.5  1.5]);
title('Z-plane Pole-Zero Plot of H(z)| ');%title({'X(W)', ['Line2 x=' num2str(3.3)]})
% set(findobj(gca, 'type', 'line'),'markersize', 12)
% set(findobj(gca, 'type', 'line'),'linewidth', 2)
% figax=findall(hh,'Type','Axes');set(figax(2),'LineWidth',2);
set(findall(hh,'Type','text'),'FontSize',14);
set(findall(hh,'Type','text'),'FontWeight','bold');
set(findobj(gca, 'type', 'line'),'markersize', 10)
set(findobj(gca, 'type', 'line'),'linewidth', 2)
%figax=findall(hh,'Type','Axes');set(figax(2),'LineWidth',2);
grid off
%hh.Position=[360 278 550 550];
disp(' Fig. 3. Poles and zeroes of H(z) in z-plane. ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q7. Fig. 3 shows the pole-zero plot of H(z) for the h[n] in Fig. 2.')
disp('    Note the unit circle, and locations of poles and zeroes.')
disp('    Use the cursor to find the locations of the 2 zeroes in Fig. 3.')
disp(' ')
disp('Answer: two zero locations z= ???, z= ??? ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q8. Fig. 3 shows the pole-zero plot of H(z) for the h[n] in Fig. 2.')
disp('    Note the unit circle, and locations of poles and zeroes.')
disp('    Use the cursor to find the locations of the 2 poles in Fig. 3.')
disp(' ')
disp('Answer: two zero locations z= ???, z= ???  ')
disp(' ')
disp(' ')


%create symbolic function
disp(' parameters and definitions:')
clear s z w W x2 X2z X2dtft n
syms s z w W x2 X2z X2dtft n
%h(n)=1*(heaviside(n-1)-heaviside(n-4))
x2(n)=1*(heaviside(n)-heaviside(n-5))  %sol
disp('***** z-transform')
X2z(z)=simplify(ztrans(x2,z),'Steps',5)
nn=0:20;
x2n=x2(nn);
X2dtft(w)= eval(subs(X2z(exp(i*w)))); %compute DTFT
wmax=2*pi
wmin=-2*pi
% %logarithmic frequency scale
%wmin=wmax/10000
% fac=ex2p((log(wmax)-log(wmin))/(50));
% wn=wmin*fac.^(0:1:50);   
wn=wmin:wmax/100:wmax;  %linear freq scale
X2dtftn=double(eval(X2dtft(wn)));
%HdtftnDegrees= unwrap(atan2(imag(Hdtftn),real(Hdtftn)),5)*180/pi;
X2dtftnDegrees= (atan2(imag(X2dtftn),real(X2dtftn)))*180/pi;

TFstr = char(subs(X2z)); 
z = tf('z',1)
disp('======= transfer function form: '); 
eval(['X2zsys = ',TFstr]); 
[poles,zeroes] = pzmap(X2zsys);
disp(['  poles=' num2str(poles') '  and zeroes=' num2str(zeroes') ])

 
hh=figure(4);
sp=subplot(3,1,1);stem(nn,x2n,'LineWidth',3);grid on;%sp.YLim=[-10 10];
title({'Z-Transform Example' 'x2[n]'});xlabel('n');  
sp=subplot(3,1,2);plot(wn,abs(X2dtftn),'LineWidth',3);grid on;
%loglog(wn,abs(Hdtftn),'LineWidth',3);grid on;
%ymax=10^(floor(log10(max(abs(Hdtftn))))+1);
sp.XLim=[wmin wmax];%sp.YLim=[ymax/1000 ymax ];
title({'|X2(w)|= |X2(z)| for z=e^{jw}'});xlabel('Frequency w (rad/sample)'); set(gca,'LineWidth',2); 
sp=subplot(3,1,3);plot(wn,X2dtftnDegrees,'LineWidth',3);grid on;sp.XLim=[wmin wmax];
%semilogx2(wn,HdtftnDegrees,'LineWidth',3);grid on;
% sp.X2Lim=[wmin wmax];
title('Angle of X2(w) (degrees)');xlabel('Frequency w (rad/sample)');   set(gca,'LineWidth',2);
set(findall(hh,'Type','text'),'FontSize',16);
hh.Position=[360 278 550 550];
disp(' Fig. 4. Signal x2[n], |X2(w)|= |X2(z)| for z=e^{jw}, and angle of X2(w). ')
disp(' ')
disp(' ')



hh=figure(5);
%pzm=pzmap(Hzsys)
pzp=pzplot(X2zsys);
setoptions(pzp,'FreqUnits','auto','Grid','off','YLim',[-1.5 1.5],'XLim',[-1.5  1.5]);
title('Z-plane Pole-Zero Plot of X2(z)| ');%title({'X(W)', ['Line2 x=' num2str(3.3)]})
% set(findobj(gca, 'type', 'line'),'markersize', 12)
% set(findobj(gca, 'type', 'line'),'linewidth', 2)
% figax=findall(hh,'Type','Axes');set(figax(2),'LineWidth',2);
set(findall(hh,'Type','text'),'FontSize',14);
set(findall(hh,'Type','text'),'FontWeight','bold');
set(findobj(gca, 'type', 'line'),'markersize', 10)
set(findobj(gca, 'type', 'line'),'linewidth', 2)
%figax=findall(hh,'Type','Axes');set(figax(2),'LineWidth',2);
grid off
%hh.Position=[360 278 550 550];
disp(' Fig. 5. Poles and zeroes of X2(z) in z-plane. ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q9. Fig. 5 shows the pole-zero plot of X2(z) for the x2[n] in Fig. 4.')
disp('    Note the unit circle, and locations of poles and zeroes.')
disp('    Use the cursor to find the locations of the 4 zeroes in Fig. 3.')
disp(' ')
disp('Answer: four zero locations z= ??? +/- j???, z= ??? +/- j???  ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q10. Recall that the unit circle in Fig. 5 corresponds to the frequency')
disp('    response in Fig. 4 with DTFT X2(w) = X2(z) for z=exp(jw).')
disp('    Calculate the 4 angles in radians corresponding to the 4 zeroes ')
disp('    locations in Fig. 5. (two angles are positive, two angles are negative)')
disp('    Using the cursor, measure the 4 frequencies in rad/sample  ')
disp('    between -pi and pi where |X2(w)| in Fig. 4 equals zero. ')
disp('    Are the 4 angles numerically equal (+/-10%) to the 4 frequencies? (ignore units)')
disp(' ')
disp('Answer: angles of 4 zero locations: +/- ??? radians, and +/- ??? radians ')
disp('Answer: frequencies where|X2(w)|=0 : +/- ??? rad/sample, and +/- ??? rad/sample ')
disp('Answer: angles of zeroes equals frequencies to within 10%: yes/no ')
disp(' ')
disp(' ')
