% 
% 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: FIR Filters and Digital Non-Foster Circuits')
disp('Names: Nolan Bushnell, Federico Faggin, Gordon E. Moore')
disp(' ')

disp('PART 1: Designing FIR filters using windowing method  =============================')
disp(' ')
disp('  ')



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

disp('====Parameters forfilter design ==========:')
disp('Sample period, sample rate, and filter order and 3 dB bandwidth:')
Ts=sym(1/100)
fs=1/Ts
butterOrder=15
omega3dB=30 % change this
disp('window length = "windowLength"+1, since window id from n=0 to n=windowLength :')
disp('            Note: variable "windowLength" must be an even number ')
windowLength=50;% change this
windowLength=floor(windowLength/2)*2 %force to be an even number
nn=0:150;
tim=double(nn.*eval(Ts));
nn2=0:5*(length(nn)-1);
tim2=double(nn2/5.*eval(Ts));
wmax=2*pi;wmin=-2*pi;
wmax=pi;wmin=-pi;
% %logarithmic frequency scale
%wmin=wmax/10000
% fac=exp((log(wmax)-log(wmin))/(50));
% wn=wmin*fac.^(0:1:50);   
wn=wmin:wmax/500:wmax;  %linear freq scale
 

[z,p,k]=butter(butterOrder,omega3dB,'s'); %sol design butterworth
[n,d]=zp2tf(z,p,k);
disp('Continuous-time Butterworth:')
tfHbutter=tf(n,d)
hButtert= impulse(tfHbutter,tim2);
hDesiredButtern= eval(Ts)*impulse(tfHbutter,tim);
[M,I] = max(hDesiredButtern);
count=0;
for n3=1:length(nn)
    rectWindow(n3)=1;
    hdes(n3)=0;
    triangWindow(n3)=2*(n3-1)/windowLength; %triang window=0 at n=0, and n=n3-1, pwak=2
    if n3>windowLength/2+1 triangWindow(n3)=2*(windowLength/2-1-count)/windowLength; count=count+1; end;
    if n3>windowLength+1 rectWindow(n3)=0; triangWindow(n3)=0; end;
    if n3>windowLength rectWindow(n3)=0; triangWindow(n3)=0; end;
end

%hDesiredButtern=hDesiredButtern.*0+1
hRectn=rectWindow.*hDesiredButtern';
hTriangn=triangWindow.*hDesiredButtern';
denz=rectWindow.*0; denz(1)=1; %denominator for conputing freqz
tfHRectz=tf(hRectn,denz,double(eval(Ts)));
tfHTriangz=tf(hTriangn,denz,double(eval(Ts)));

 
%impulse(tfHbutter,tim)

 
% fr=freqresp(tfHbutter,eval(wn/Ts)); nf=1:length(fr); freqrespButterH(nf)=fr(1,1,nf);
% fr=freqresp(tfHRectz,eval(wn/Ts));nf=1:length(fr); freqrespRectWin(nf)=fr(1,1,nf); 
% fr=freqresp(tfHTriangz,eval(wn/Ts)); nf=1:length(fr); freqrespTriangWin(nf)=fr(1,1,nf);
%dB
fr=freqresp(tfHbutter,eval(wn/Ts)); nf=1:length(fr); freqrespButterH(nf)=10*log10(abs(fr(1,1,nf)).^2+1e-8);
fr=freqresp(tfHRectz,eval(wn/Ts));nf=1:length(fr); freqrespRectWin(nf)=10*log10(abs(fr(1,1,nf)).^2+1e-8); 
fr=freqresp(tfHTriangz,eval(wn/Ts)); nf=1:length(fr); freqrespTriangWin(nf)=10*log10(abs(fr(1,1,nf)).^2+1e-8);


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q1. FIRST, ')
disp('    change the Butterworth 3 dB bandwidth "omega3dB" to 20 rad/s above, ')
disp('    and the windowlength "windowLength" to 140.')
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('    what is the order of the continuous-time Butterworth filter? ')
disp(' ')
disp('Answer: order = ??? ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q3. After making the changes in Q1, ')
disp('    the un-windowed "desired filter impulse response" is defined above  ')
disp('    as "hDesiredButtern" and is equal to Ts*h(n*Ts). ')
disp('     Does Ts*h(n*Ts) equal the impulse-invariance design? ')
disp('     Note: windows are later applied to "hDesiredButtern."')
disp(' ')
disp('Answer: Ts*h(n*Ts) equal the impulse-invariance design: yes/no ')
disp(' ')
disp(' ')


hh=figure(1);
sp=subplot(3,1,1);plot(tim2,hButtert,'b',tim,triangWindow,'r:',tim,rectWindow,'c:','LineWidth',3);grid on;
sp.XLim=[0 max(tim2)];
title({['Butterworth Filter Impulse Responses,  omega3dB=' num2str(omega3dB) '  windowLength' num2str(windowLength)  '  fs=' num2str(eval(1/Ts)) ' samples/s' ], 'h(t) in solid blue, and windows (dotted curves) ' });xlabel('Time (s)');  
sp=subplot(3,1,2);stem(nn,hRectn,'LineWidth',3);grid on;sp.YLim=[-0.05 0.1];
title({'Rectangular windowed design hrect[n]'});xlabel('n');  
sp=subplot(3,1,3);stem(nn,hTriangn,'LineWidth',3);grid on;sp.YLim=[-0.05 0.1];
title({'Triangular windowed design design htriang[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), rectangular windowed hrect[n], ')
disp('         and triangular windowed htriang[n]. ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q4. After making the changes in Q1, ')
disp('    Fig.1 shows the impulse response of the continuous-time filter h(t), ')
disp('    rectangular windowed hrect[n], and triangular windowed htriang[n].')
disp('      Using the cursor, measure:')
disp('        what is the peak value of h(t)?')
disp('        what is the peak value of hrect[n]?')
disp('    Does Ts=(peak of hrect[n])/(peak of h(t) )? ')
disp(' ')
disp('Answer: the peak value of h(t) = ???  ')
disp('Answer: peak value of hrect[n] = ???  ')
disp('Answer: Ts=(peak of hrect[n])/(peak of h(t) ) : yes/no  ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q5. After making the changes in Q1, ')
disp('    in Fig.1 , using the cursor, measure:')
disp('        what is the peak value of the triangular window?')
disp('    Does the triangular window have the same width as the rectangular window?')
disp(' ')
disp('Answer: peak value of the triangular window = ??? ')
disp('Answer: triang. window same width as rect. window +/-5%: yes/no  ')
disp(' ')
disp(' ')

 
hh=figure(2);
sp=subplot(3,1,1);plot(wn/Ts,(freqrespButterH),'LineWidth',3);grid on;
sp.XLim=eval([wmin/Ts wmax/Ts]);sp.YLim=[-70 0];
title({['Butterworth FilterFrequency Responses, ', ' omega3dB=' num2str(omega3dB) '  windowLength' num2str(windowLength)  '  fs=' num2str(eval(1/Ts)) ' samples/s' ], 'h(t) in solid blue, and windows (dotted curves) ' });xlabel('Time (s)');  
xlabel('Frequency W (rad/s)');  ylabel('dB');
sp=subplot(3,1,2);plot(wn,(freqrespRectWin),'LineWidth',3);grid on;
sp.XLim=[wmin wmax];sp.YLim=[-70 0];
title({'Rectangular window design |Hrect(z)| for z=e^{jw}'});
xlabel('Frequency w (rad/sample)');  ylabel('dB'); set(gca,'LineWidth',2); 
sp=subplot(3,1,3);plot(wn,(freqrespTriangWin),'LineWidth',3);grid on;
sp.XLim=[wmin wmax];sp.YLim=[-70 0];
title('Triangular window design  |Htriang(z)| for z=e^{jw}');
xlabel('Frequency w (rad/sample)');  ylabel('dB');  set(gca,'LineWidth',2);
set(findall(hh,'Type','text'),'FontSize',14);
hh.Position=[360 278 550 550];
disp(' Fig. 2. Frequency responses in dB.  Continuous time |H(W)|,  ')
disp('          rectangular window design |Hrect(z)| for z=e^{jw}, ')
disp('          triangular window design  |Htriang(z)|  for z=e^{jw}. ')
disp(' ')
disp(' ')



disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q6. After making the changes in Q1, ')
disp('    Fig.2 shows the frequency response of the continuous-time filter |H(W)|, ')
disp('    of the rectangular window design  |Hrect(z)|, and triangular ')
disp('     window design  |Htriang(z)|  for z=e^{jw}. ')
disp('     Using the cursor, measure:')
disp('       what is the 3 dB bandwidth of |H(W)| in rad/s?')
disp('          Hint: 3 dB bandwidth is frequency where the plot is -3 dB ')
disp('       what is the 3 dB bandwidth of |Hrect(z)| in rad/sample?')
disp('    Does the second answer equal the first answer times Ts (+/-5%)?')
disp(' ')
disp('Answer: 3 dB bandwidth of |H(W)| = ??? rad/s  ')
disp('Answer: 3 dB bandwidth of |Hrect(z)| = ??? rad/sample  ')
disp('Answer: second answer equals first answer times Ts, (+/-5%): yes/no  ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q7. After making the changes in Q1, ')
disp('    in Fig.2,  using the cursor, measure:')
disp('       what is the 40 dB bandwidth of |H(W)| in rad/s?')
disp('          Hint: 40 dB bandwidth is frequency where the plot is -40 dB ')
disp('       what is the 40 dB bandwidth of |Hrect(z)| in rad/sample?')
disp('       what is the 40 dB bandwidth of |Htriang(z)| in rad/sample?')
disp(' ')
disp('Answer: 40 dB bandwidth of |H(W)| = ??? rad/s  ')
disp('Answer: 40 dB bandwidth of |Hrect(z)| = ??? rad/sample  ')
disp('Answer: 40 dB bandwidth of |Htriang(z)| = ??? rad/sample  ')
disp(' ')
disp(' ')

disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q8. After making the changes in Q1, ')
disp('    observe Fig.2,  ')
disp('       which design gives more attenuation at high frequencies,')
disp('       the triangular windowed design or the rectangular windowed design? ')
disp('      (At high frequencies near 3.1 rad/sample)')
disp(' ')
disp('Answer: which design gives more attenuation: triangular/rectangular  ')
disp(' ')
disp(' ')



disp('PART 2: Digital non-Foster impedances  =============================')
disp(' ')
disp('  ')

 
clear R C Ts fs s n t w W z h hn hd hrect htriang win H Hbilin himpinvn Himpinv
syms  z C Ts s Zc Hc

disp('====Parameters Digital non-Foster design ==========:')
disp('Sample period, sample rate, and capacitance:')
fs=200*10^(6)
Ts=1/fs
C=-80*10^(-12)
disp('Find Hc(z) for the digital capacitor:')
Hc(z)=C*(1-1/z)/Ts;
disp('Hc(z)=')
olddigits=digits(); digits(5);vpa(eval(Hc(z)))
digits(olddigits);
disp('Find Zc(s):')
Zc(s)=subs( (s*Ts/((1-1/z)*Hc(z))), z, exp(s*Ts));
olddigits=digits(); digits(5);vpa(eval(Zc(s)))
digits(olddigits);

wmax=pi;wmin=-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
 

 
hh=figure(3);
%sp=subplot(3,1,1);
plot(wn/Ts,real(vpa(Zc(i*wn/Ts+pi/10^6))),'r',wn/Ts,imag(vpa(Zc(i*wn/Ts+pi/10^6))),'b',wn/Ts,imag(vpa(1./((-C)*(i*wn/Ts+pi/10^6)))),'b:', 'LineWidth',3);grid on;
title({['Digital negative capacitor,  C=' num2str(C)  '  fs=' num2str(fs) ' samples/s' ], 'Digital Re\{Zc(s)\} in solid red, Im\{Zc(s)\} solid blue, ', 'Im\{positive capacitance\} in dotted red' });xlabel('Time (s)');  
xlim([0 double(wmax/Ts)]);ylim([-1000 1000]);
xlabel('Frequency W (rad/s)');  ylabel('ohms'); set(gca,'LineWidth',2); 
set(findall(hh,'Type','text'),'FontSize',14);

%hh.Position=[360 278 550 550];
disp(' Fig. 3. Impedance Zc(s) of a digital negative capacitor.  Solid red  ')
disp('          curve is real part of Zc(s), solid blue curve is imaginary ')
disp('          part of Zc(s).  Dotted blue curve shows the corresponding ')
disp('          reactance of a positive capacitance of the same value.')
disp(' ')
disp(' ')




disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q9. The input impedance of a digital negative capacitor is ')
disp('    shown in Fig.3, where the solid red curve is real part of Zc(s),  ')
disp('    the solid blue curve is the imaginary part of Zc(s), and the  ')
disp('    dotted blue curve shows the corresponding  reactance of a ')
disp('    positive capacitance of the same value.')
disp('    In the plotted example, what is the design value for C in pF?')
disp('        Hint: C should be negative')
disp(' ')
disp('Answer: design value for C = -??? pF  ')
disp(' ')
disp(' ')


disp('QQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQQ ')
disp('Q10. In Fig.3, ')
disp('    What is theoretical impedance of the negative capacitance at 100 megaradian/s? ')
disp('        Hint: impedance should be purely imaginary, ... +j times something')
disp('    Using the cursor in Fig. 3, ')
disp('    what is the measured reactance (solid blue curve) at 100 megarad/s? ')
disp('        Hint: should be fairly close to theory, +/-20% ')
disp(' ')
disp('Answer: theoretical impedance of C at 100 megarad/s= approx +j??? ohms ')
disp('Answer: measured impedance of C at 100 megarad/s= approx +j??? ohms ')
disp(' ')
disp(' ')