% 
% copyright 2021 by Thomas P. Weldon
% all rights reserved
%

clear all
close all hidden
clc
sympref('HeavisideAtOrigin',1); %set equal to 1 as in unitstep


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: Digital Filter Design')
disp('Names: Nolan Bushnell, Federico Faggin, Gordon E. Moore')
disp(' ')

disp('PART 1: Designing a "digital RC filter"  =============================')
disp(' ')
disp('  ')

disp('Theory:  ')
disp('   In the later sections we will use impulse invariance and bilinear')
disp('   transforms to convert the analog filter into a digital filter.')
disp('   So, we first design an analog filter.')
disp('   ')

disp('   Consider an RC lowpass filter where H(s)=1/(1+sRC)')
disp('   ')

clear R C Ts fs s n t w W z h hn H Hbilin himpinvn Himpinv
syms  R C Ts fs s n t w W z h hn H Hbilin himpinvn Himpinv

disp(' Continuous-time RC filter Laplace domaintransform of impulse response:')
H(s)=1/(1+s*R*C)
disp(' Continuous-time RC filter impulse response:')
h(t)=ilaplace(H,s,t)

disp(' RC filter constants:')
R=sym(200)
C=sym(1/2000)
disp('these can be converted to maatlab transfer functions:')
syms s;TFstr = char(eval(vpa(H))); s=tf('s'); eval(['tfH = ',TFstr]); 
tfH=series(tfH,tf([1],[1])); syms s; % fix a bug for when transfer function equals a number
 

disp(' ')
disp('====Digital filter designs ==========:')
disp('Sample period and sample rate:')
Ts=sym(1/10);
fs=1/Ts
nn=0:50;
tim=vpa(nn.*Ts);

disp('======= Bilinear transform design:')
disp('====  Hbilin(z):  ')
Hbilin(z)=collect(simplify(subs(H,s,(2/Ts)*(z-1)/(z+1)),'Steps',5),z)
hbilin(n)=vpa(simplify(iztrans(Hbilin(z),z,n),'Steps',5))
olddigits=digits(); digits(5);vpa(eval(Hbilin(z)))
digits(olddigits);
disp('this can be converted to a maatlab transfer function:')
syms z;TFstr = char(eval(vpa(Hbilin))); z=tf('z',eval(Ts)); eval(['tfHbilin = ',TFstr]); 
tfHbilin=series(tfHbilin,tf([1],[1])); syms z; % fix a bug for when transfer function equals a number
[poles,zeroes] = pzmap(tfHbilin);
disp(['tfHbilin poles=' num2str(poles') '  and zeroes=' num2str(zeroes') ])
disp(['  magnitude of poles=' num2str(abs(poles')) ])


disp('======= Impulse invariance  design:')
himpinv(n)=Ts*h(n*Ts);
disp('====  Himpinv(z):  ')
Himpinv(z)=collect(simplify(ztrans(himpinv(n),n,z),'Steps',5),z)
olddigits=digits(); digits(5);vpa(eval(Himpinv(z)))
digits(olddigits);
disp('this can be converted to a maatlab transfer function:')
syms z;TFstr = char(eval(vpa(Himpinv))); z=tf('z',eval(Ts)); eval(['tfHimpinv = ',TFstr]); 
tfHimpinv=series(tfHimpinv,tf([1],[1])); syms z; % fix a bug for when transfer function equals a number
[poles,zeroes] = pzmap(tfHimpinv);
disp(['tfHimpinv poles=' num2str(poles') '  and zeroes=' num2str(zeroes') ])
disp(['  magnitude of poles=' num2str(abs(poles')) ])



disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q1. FIRST, ')
disp('    change the RC filter parameters above, ')
disp('    to R=sym(500) and C=sym(1/1000).')
disp('    For the system above, what is the sample period Ts?')
disp('    For the system above, what is the sample rate fs?')
disp(' ')
disp('Answer: sample period = ??? s ')
disp('Answer: sample rate = ??? sample/s ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q2. After making the changes in Q1, ')
disp('    for the continuous-time RC filter parameters above, ')
disp('       what is the RC time constant in seconds?')
disp('       What is the theoretical 3 dB bandwidth in rad/s?')
disp('       Does the theoretical impulse invariance equal "tfHimpinv" above?')
disp(' ')
disp('Answer: RC time constant = ??? s ')
disp('Answer: theoretical 3 dB bandwidth = ??? rad/s ')
disp('Answer: theoretical impulse invariance equal "tfHimpinv": yes/no ')
disp(' ')
disp(' ')



 
hh=figure(1);
sp=subplot(3,1,1);plot(tim,eval(h(tim)),'LineWidth',3);grid on;%sp.YLim=[-10 10];
title({['RC Filter Impulse Responses,  R=' num2str(eval(R)) '  C=' num2str(eval(C))  '  fs=' num2str(eval(1/Ts)) ' samples/s' ], 'h(t)' });xlabel('Time (s)');  
sp=subplot(3,1,2);stem(nn,eval(hbilin(nn)),'LineWidth',3);grid on;%sp.YLim=[-10 10];
title({'Bilinear transform design hbilin[n]'});xlabel('n');  
sp=subplot(3,1,3);stem(nn,eval(himpinv(nn)),'LineWidth',3);grid on;%sp.YLim=[-10 10];
title({'Impulse invariance design himpinv[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);
set(findall(hh,'Type','text'),'FontSize',14);
hh.Position=[360 278 550 550];
disp(' Fig. 1. Impulse responses.  Continuous-time h(t), bilinear transform hbilin[n], ')
disp('         and impulse invariance himpinv[n]. ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q3. After making the changes in Q1, ')
disp('    Fig. 1 shows the impulse responses, ')
disp('    where the top plot is the continuous-time impulse response h(t), ')
disp('    the middle plot is bilinear transform impulse response hbilin[n],  ')
disp('    and the bottom plot is the impulse response for the impulse  ')
disp('    invariance design himpinv[n]. ')
disp('      Using the cursor, measure:')
disp('        the time for h(t) to reach 37% of its initial value?')
disp('        the value of n for himpinv[n] to reach 37% of its initial value?')
disp('    Why is the signal level of himpinv[n] 1/10 as large as h(t)? ')
disp(' ')
disp('Answer: time for h(t) to reach 37% = ??? s ')
disp('Answer: value of n for himpinv[n] to reach 37%: n= ???  ')
disp('Answer: why himpinv[n] is 1/10 as large as h(t): short explanation  ')
disp(' ')
disp(' ')

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/250:wmax;  %linear freq scale
 
hh=figure(2);
sp=subplot(3,1,1);plot(wn/Ts,abs(eval(subs(H(i*wn/Ts)))),'LineWidth',3);grid on;
sp.XLim=eval([wmin/Ts wmax/Ts]);sp.YLim=[0 1.25 ];
title({['RC Filter Frequency Responses,  R=' num2str(eval(R)) '  C=' num2str(eval(C)) '  fs=' num2str(eval(1/Ts)) ' samples/s' ], 'Continuous-time |H(W)|' });
xlabel('Frequency W (rad/s)');  
sp=subplot(3,1,2);plot(wn,abs(eval(subs(Hbilin(exp(i*wn))))),'LineWidth',3);grid on;
sp.XLim=[wmin wmax];sp.YLim=[0 1.25 ];
title({'Bilinear transform |Hbilin(z)| for z=e^{jw}'});xlabel('Frequency w (rad/sample)'); set(gca,'LineWidth',2); 
sp=subplot(3,1,3);plot(wn,abs(eval(subs(Himpinv(exp(i*wn))))),'LineWidth',3);grid on;
sp.XLim=[wmin wmax];sp.YLim=[0 1.25 ];
title('Impulse invariance design |Himpinv(z)| for z=e^{jw}');xlabel('Frequency w (rad/sample)');   set(gca,'LineWidth',2);
set(findall(hh,'Type','text'),'FontSize',14);
hh.Position=[360 278 550 550];
disp(' Fig. 2. Frequency responses.  Continuous time |H(W)|, bilinear transform  ')
disp('         design = |Hbilin(z)| for z=e^{jw}, impulse invariance design ')
disp('         |Himpinv(z)| for z=e^{jw}. ')
disp(' ')
disp(' ')
 


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q4. After making the changes in Q1, ')
disp('    Fig. 2 shows the frequency responses of continuous-time |H(W)|, bilinear  ')
disp('         transform |Hbilin(z)| for z=e^{jw}, and impulse invariance design ')
disp('         |Himpinv(z)| for z=e^{jw}, where "w" denotes discrete-time')
disp('         frequency in rad/sample, and "W" denotes continuous-time')
disp('         frequency in rad/s.   ')
disp('    Using the cursor, measure:')
disp('        the 3 dB bandwidth of continuous-time |H(W)| rad/s?')
disp('        the 3 dB bandwidth of bilinear transform |Hbilin(z)| in rad/sample?')
disp('        the 3 dB bandwidth of impulse invariance |Himpinv(z)| in rad/sample?')
disp(' ')
disp('Answer: 3 dB bandwidth of continuous-time |H(W)|= ??? rad/s ')
disp('Answer: 3 dB bandwidth of bilinear transform |Hbilin(z)| = ??? rad/sample ')
disp('Answer: 3 dB bandwidth of impulse invariance |Himpinv(z)|= ??? rad/sample')
disp(' ')
disp(' ')

 
disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q5. After making the changes in Q1, ')
disp('    in Fig. 2,  ')
disp('    Is aliasing visible in the bilinear transform design |Hbilin(z)| ?')
disp('    Is aliasing visible in the impulse invariance design |Himpinv(z)| ?')
disp(' ')
disp('Answer: aliasing visible in |Hbilin(z)|: yes/no ')
disp('Answer: aliasing visible in |Himpinv(z)|: yes/no ')
disp(' ')
disp(' ')

disp(' =============== Next, change the sampling period =============== ')

disp('Sample period and sample rate:')
Ts=sym(1/40);
fs=1/Ts
nn=0:50;
tim=vpa(nn.*Ts);

 
disp('======= Bilinear transform design:')

disp('====  Hbilin(z):  ')
Hbilin(z)=collect(simplify(subs(H,s,(2/Ts)*(z-1)/(z+1)),'Steps',5),z)
hbilin(n)=vpa(simplify(iztrans(Hbilin(z),z,n),'Steps',5))
olddigits=digits(); digits(5);vpa(eval(Hbilin(z)))
digits(olddigits);
disp('this can be converted to a maatlab transfer function:')
syms z;TFstr = char(eval(vpa(Hbilin))); z=tf('z',eval(Ts)); eval(['tfHbilin = ',TFstr]); 
tfHbilin=series(tfHbilin,tf([1],[1])); syms z; % fix a bug for when transfer function equals a number
[poles,zeroes] = pzmap(tfHbilin);
disp(['tfHbilin poles=' num2str(poles') '  and zeroes=' num2str(zeroes') ])
disp(['  magnitude of poles=' num2str(abs(poles')) ])


disp('======= Impulse invariance  design:')
himpinv(n)=Ts*h(n*Ts);
disp('====  Himpinv(z):  ')
Himpinv(z)=collect(simplify(ztrans(himpinv(n),n,z),'Steps',5),z)
olddigits=digits(); digits(5);vpa(eval(Himpinv(z)))
digits(olddigits);
disp('this can be converted to a maatlab transfer function:')
syms z;TFstr = char(eval(vpa(Himpinv))); z=tf('z',eval(Ts)); eval(['tfHimpinv = ',TFstr]); 
tfHimpinv=series(tfHimpinv,tf([1],[1])); syms z; % fix a bug for when transfer function equals a number
[poles,zeroes] = pzmap(tfHimpinv);
disp(['tfHimpinv poles=' num2str(poles') '  and zeroes=' num2str(zeroes') ])
disp(['  magnitude of poles=' num2str(abs(poles')) ])

hh=figure(3);
sp=subplot(3,1,1);plot(wn/Ts,abs(eval(subs(H(i*wn/Ts)))),'LineWidth',3);grid on;
sp.XLim=eval([wmin/Ts wmax/Ts]);sp.YLim=[0 1.25 ];
title({['RC Filter Frequency Responses,  R=' num2str(eval(R)) '  C=' num2str(eval(C)) '  fs=' num2str(eval(1/Ts)) ' samples/s' ], 'Continuous-time |H(W)|' });
xlabel('Frequency W (rad/s)'); 
sp=subplot(3,1,2);plot(wn,abs(eval(subs(Hbilin(exp(i*wn))))),'LineWidth',3);grid on;
sp.XLim=[wmin wmax];sp.YLim=[0 1.25 ];
title({'Bilinear transform  |Hbilin(z)| for z=e^{jw}'});xlabel('Frequency w (rad/sample)'); set(gca,'LineWidth',2); 
sp=subplot(3,1,3);plot(wn,abs(eval(subs(Himpinv(exp(i*wn))))),'LineWidth',3);grid on;
sp.XLim=[wmin wmax];sp.YLim=[0 1.25 ];
title('Impulse invariance design |Himpinv(z)| for z=e^{jw}');xlabel('Frequency w (rad/sample)');   set(gca,'LineWidth',2);
set(findall(hh,'Type','text'),'FontSize',14);
hh.Position=[360 278 550 550];
disp(' Fig. 3. Frequency responses, with new sample rate.  Continuous time |H(W)|,  ')
disp('          bilinear transform design = |Hbilin(z)| for z=e^{jw}, ')
disp('          impulse invariance design |Himpinv(z)| for z=e^{jw}. ')
disp('          ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q6. After making the changes in Q1, ')
disp('    and with the higher sample rate used in Fig. 3,  ')
disp('       What is the higher sample rate used in Fig. 3?')
disp('       What is the old sample rate in Fig. 2?')
disp('    Has |Himpinv(z)| aliasing been reduced in Fig. 3 as compared to Fig.2?')
disp(' ')
disp('Answer: higher sample rate in Fig. 3: ???samples/s ')
disp('Answer: old sample rate in Fig. 2: ??? samples/s ')
disp('Answer: |Himpinv(z)| aliasing been reduced: yes/no ')
disp(' ')
disp(' ')



disp('PART 2: Designing a "digital Butterworth filter"  =============================')
disp(' ')
disp('  ')
format SHORTE

disp('Sample period and sample rate:')
Ts=sym(1/44100);
fs=1/Ts
nn=0:50;
tim=vpa(nn.*Ts);

[z,p,k]=butter(2,2*pi*9000,'s'); %design butterworth

[n,d]=zp2tf(z,p,k);
disp('Continuous-time Butterworth:')
tfHbutter=tf(n,d)

%freqresp(Hbutter,ww)
% hh=figure(4)
% bode(n,d)

disp('======= Bilinear transform  design:')
[nd,dd]=bilinear(n,d,eval(1/Ts));
% freqz(nd,dd,256) 
tfHbutterBilin=tf(nd,dd,eval(Ts))
% pzmap(tfHbilin)
% zgrid;

disp('======= Impulse invariance  design:')
[nd,dd]=impinvar(n,d,eval(fs));
% freqz(nd,dd,256)
tfHbutterImpinv=tf(nd,dd,eval(Ts))
[top,bot]=residuez(nd,dd);
% pzmap(tfHimpinv)
% zgrid;


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/250:wmax;  %linear freq scale



fr=freqresp(tfHbutter,eval(wn/Ts)); nf=1:length(fr); freqrespButterH(nf)=fr(1,1,nf);
fr=freqresp(tfHbutterBilin,eval(wn/Ts));nf=1:length(fr); freqrespButterHbilin(nf)=fr(1,1,nf); 
fr=freqresp(tfHbutterImpinv,eval(wn/Ts)); nf=1:length(fr); freqrespButterHimpinv(nf)=fr(1,1,nf);


hh=figure(4);
sp=subplot(3,1,1);plot(wn/Ts,abs(freqrespButterH),'LineWidth',3);grid on;
sp.XLim=eval([wmin/Ts wmax/Ts]);sp.YLim=[0 1.25 ];
title({['Butterworth Filter Frequency Responses  '  ' fs=' num2str(eval(1/Ts)) ' samples/s' ], 'Continuous-time |tfHbutter(W)|' });
xlabel('Frequency W (rad/s)');  
sp=subplot(3,1,2);plot(wn,abs(freqrespButterHbilin),'LineWidth',3);grid on;
sp.XLim=[wmin wmax];sp.YLim=[0 1.25 ];
title({'Bilinear transform |tfHbutterBilin(z)| for z=e^{jw}'});xlabel('Frequency w (rad/sample)'); set(gca,'LineWidth',2); 
sp=subplot(3,1,3);plot(wn,abs(freqrespButterHimpinv),'LineWidth',3);grid on;
sp.XLim=[wmin wmax];sp.YLim=[0 1.25 ];
title('Impulse invariance design |tfHbutterImpinv(z)| for z=e^{jw}');xlabel('Frequency w (rad/sample)');   set(gca,'LineWidth',2);
set(findall(hh,'Type','text'),'FontSize',14);
hh.Position=[360 278 550 550];
disp(' Fig. 4. Butterworth frequency responses.  Continuous time |tfHbutter(W)|,  ')
disp('         bilinear  transform design = |tfHbutterBilin(z)| for z=e^{jw}, ')
disp('         impulse invariance design |tfHbutterImpinv(z)| for z=e^{jw}. ')
disp(' ')
disp(' ')
 

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q7. Fig. 4 shows  Butterworth filter design frequency responses.    ')
disp('    The top plot is continuous time frequency response |tfHbutter(W)|,  ')
disp('    middle is bilinear transform design = |tfHbutterBilin(z)| for z=e^{jw}, ')
disp('    the bottom plot is impulse invariance design |tfHbutterImpinv(z)|  ')
disp('    for z=e^{jw}. ')
disp('      ')
disp('    Before making changes,  ')
disp('    observe the aliasing apparent in the impulse invariance design.  ')
disp('    Then,  ')
disp('    change the arguments in the "[z,p,k]=butter(...)" command line above to')
disp('    design a 4th-order lowpass Butterworth with 3 dB bandwidth of 2200 Hz.')
disp(' ')
disp(' ')

 
disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q8. After making the changes in Q7, ')
disp('    in Fig. 4, using the cursor, measure:')
disp('        the 3 dB bandwidth of continuous-time |tfHbutter(W)| in rad/s?')
disp('        the 3 dB bandwidth of bilinear transform |tfHbutterBilin(z)| in rad/sample?')
disp('        the 3 dB bandwidth of impulse invariance |tfHbutterImpinv(z)|in rad/sample?')
disp(' ')
disp('Answer: 3 dB bandwidth of |tfHbutterImpinv(z)|= ??? rad/s ')
disp('Answer: 3 dB bandwidth of |tfHbutterImpinv(z)| = ??? rad/sample ')
disp('Answer: 3 dB bandwidth of |tfHbutterImpinv(z)| = ??? rad/sample ')
disp(' ')
disp(' ')


hh=figure(5);
%pzm=pzmap(Hzsys)
pzp=pzplot(tfHbutterBilin);
setoptions(pzp,'FreqUnits','auto','Grid','off','YLim',[-1.5 1.5],'XLim',[-1.5  1.5]);
xlabel('Real(z)');ylabel('Imag(z)');
title({ 'Bilinear Transform Butterworth Filter' ' Pole-Zero Plot ' });
%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. 5. Poles and zeroes of bilinear transform Butterworth in z-plane. ')
disp(' ')
disp(' ')

 
disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q9. Fig. 5 shows the poles and zeroes of the bilinear transform  ')
disp('    Butterworth filter design, tfHbutterBilin. ')
disp('    Based on the poles shown in Fig. 5,  ')
disp('    is the filter stable? ')
disp(' ')
disp('Answer: filter in Fig. 5 is stable: yes/no ')
disp(' ')
disp(' ')

hh=figure(6);
sp=subplot(3,1,1);stepplot(tfHbutter,eval(nn.*Ts));grid on;
set(gca,'LineWidth',2);set(findall(gcf,'type','line'),'linewidth',3);sp.YLim=[0 1.25];
title({['Butterworth Filter Step Responses '  '  fs=' num2str(eval(1/Ts)) ' samples/s' ], 'step(t)' });
xlabel('Time (s)');  
sp=subplot(3,1,2);stepplot(tfHbutter,eval(nn.*Ts));grid on;
set(gca,'LineWidth',2);set(findall(gcf,'type','line'),'linewidth',3);sp.YLim=[0 1.25];
title({'Bilinear transform design stepbilin[n]'});xlabel('Time (s)');  
sp=subplot(3,1,3);stepplot(tfHbutter,eval(nn.*Ts));grid on;
set(gca,'LineWidth',2);set(findall(gcf,'type','line'),'linewidth',3);sp.YLim=[0 1.25];
title({'Impulse invariance design stepimpinv[n]'});xlabel('Time (s)');  
set(findall(hh,'Type','text'),'FontSize',14);
hh.Position=[360 278 550 550];
disp(' Fig. 6. Butterworth step responses.  Continuous-time step(t), bilinear  ')
disp('         transform stepbilin[n], and impulse invariance stepimpinv[n]. ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q10. Fig. 6 shows the step responses of the Butterworth filter.')
disp('     The top is the continuous-time step(t), middle is bilinear  ')
disp('      transform design step response stepbilin[n], and the bottom plot ')
disp('      is impulse invariance design step response stepimpinv[n]. ') 
disp('    using the cursor, measure:')
disp('        the percent overshoot of continuous-time step(t)?')
disp('         the percent overshoot of bilinear transform stepbilin[n]?')
disp('         the percent overshoot of impulse invariance stepimpinv[n]?')
disp(' ')
disp('Answer: percent overshoot of continuous-time step(t)= ??? % ')
disp('Answer: percent overshoot of bilinear transform stepbilin[n] = ??? % ')
disp('Answer: percent overshoot of impulse invariance stepimpinv[n] = ??? % ')
disp(' ')
disp(' ')