Contents

Initial stuff

clear,warning off all
s = tf('s');

Make prototype low-pass filter

[z p k]=buttap(2);%Get info for a 2nd order butterworth filter LPF
proto=tf(k*poly(z),poly(p));%Make tf

% Get frequency response info
ww=.01:.1:1e3;
[mag phase]=bode(proto,ww);

disp('Low-pass prototype:'),proto %#ok<NOPTS>

% Display frequency response
figure(1);clf;subplot(211),semilogx(ww,20*log10(mag(:,:)))
title('Bode Diagram'),ylabel('Magnitude (dB)')
axis([0 1e2 -30 10])
subplot(212),semilogx(ww,phase(:,:)/180)
ylabel('Phase (pi radians)'),xlabel('Frequency (rps)')
axis([1e-2 1e2 -1.2 .2])

Low-pass prototype:
 
Transfer function:
        1
-----------------
s^2 + 1.414 s + 1
 

Transmission coefficient and reflection coefficient

t_squared=1/(s^4+1); %Transmission coefficient (2nd order butterworth LPF)
rho_squared=1-t_squared;

pees=pole(rho_squared);
rho=s^2/((s-pees(1))*(s-pees(2)));

disp('rho(s):'),rho %#ok<NOPTS>

% Display reflection coefficient stuff
figure(1);clf;subplot(211)
pzmap(rho_squared)
title('Pole-Zero Map for |\rho(j\omega)|^2')
ax=axis;subplot(212),pzmap(rho)
title('Pole-Zero Map for \rho(s)'),axis(ax)

rho(s):
 
Transfer function:
       s^2
-----------------
s^2 + 1.414 s + 1
 

Getting z11(s) and applying Foster's

R1=5e4;
% z11=(R1*(1+rho))/(1-rho);
z11=R1*(2*s^2+s*sqrt(2)+1)/(s*sqrt(2)+1); %Fixed Matlab glitch

% za=z11-s*(sqrt(2))*R1;
za=R1/(s*sqrt(2)+1); %Fixed Matlab glitch

ya=1/za;

disp('z11(s):'),z11 %#ok<NOPTS>
L=sqrt(2)*R1;
fprintf('Remove pole at infinity with %0.4g H inductor\n',L)
disp('za(s):'),za %#ok<NOPTS>
disp('Invert')
disp('ya(s):'),ya %#ok<NOPTS>
C=sqrt(2)/R1;R2=R1;
fprintf('Implement with %0.3d F capacitor and %0.1g Ohm terminating resistor\n',[C R2])

% Show circuit
bg=[1 1 1];figure(1);clf;
a=imread('Figures/prototype.png','BackgroundColor',bg);
image(a)
axis off;

z11(s):
 
Transfer function:
100000 s^2 + 7.071e004 s + 50000
--------------------------------
          1.414 s + 1
 
Remove pole at infinity with 7.071e+004 H inductor
za(s):
 
Transfer function:
   50000
-----------
1.414 s + 1
 
Invert
ya(s):
 
Transfer function:
1.414 s + 1
-----------
   50000
 
Implement with 2.828e-005 F capacitor and 5e+004 Ohm terminating resistor

Transform into 60Hz notch filter

wc=2*pi*60;
bw=2*pi*3;
w1=1;

% The inductor becomes an inductor and a cap in parallel
L1=w1*bw*L/wc^2;C1=1/(w1*bw*L);

% The capacitor becomes a cpacitor in series with an inductor
C2=w1*bw*C/wc^2;L2=1/(w1*bw*C);

fprintf('We now have R1=%0.3d Ohm, R2=%0.3d Ohm, C1=%0.3d F,\nL1=%0.3d H, C2=%0.3d F, L2=%0.3d H\n',[R1 R2 C1 L1 C2 L2])
% format this line better

% Show circuit
bg=[1 1 1];figure(1);clf;
a=imread('Figures/notch.png','BackgroundColor',bg);
image(a)
axis off;

We now have R1=50000 Ohm, R2=50000 Ohm, C1=7.503e-007 F,
L1=9.378e+000 H, C2=3.751e-009 F, L2=1.876e+003 H

Simulate inductors

% Floating inductor (L1) L1=C_L1*R_L1^2
C_L1=1e-8;R_L1=sqrt(L1/C_L1);

% Grounded inductor (L2) L2=C_L2*R_L2^2
C_L2=1e-7;R_L2=sqrt(L2/C_L2);

fprintf('Finally, we get R1=%.0g Ohm, R2=%.0g Ohm, C1=%0.3g F, C2=%0.3g F,\n R_L1=%.0g Ohm, R_L2=%.0g Ohm, C_L1=%0.3g F, C_L2=%0.3g F\n',...
    [R1 R2 C1 C2 R_L1 R_L2 C_L1 C_L2])

% Show circuit
bg=[1 1 1];figure(1);clf;
a=imread('Figures/simulated.png','BackgroundColor',bg);
image(a)
axis off;

% Show Ks:1
bg=[1 1 1];figure(2);clf;
a=imread('Figures/Ks1.png','BackgroundColor',bg);
image(a)
axis off;

Finally, we get R1=5e+004 Ohm, R2=5e+004 Ohm, C1=7.5e-007 F, C2=3.75e-009 F,
 R_L1=3e+004 Ohm, R_L2=1e+005 Ohm, C_L1=1e-008 F, C_L2=1e-007 F

Theoretical Results Compared to Data

[num den]=lp2bs(k*poly(z),poly(p),wc,bw);
theory=tf(num,den);

disp('Theoretical Notch Filter:'),theory %#ok<NOPTS>

ww=[1e1:1:1e4];
[mag phase]=bode(theory,ww);

fData=[2*pi*10 2*pi*70 2*pi*80 2*pi*82 2*pi*84 2*pi*100 2*pi*1e3];
magData=[-9 -9 -9.5 -26 -9.3 -9 -9];
% The flat part was not actually measured. It just seemed flat, so I'm
% making it so in the graph.


% Display frequency response
figure(1);clf;semilogx(ww/(2*pi),20*log10(mag(:,:)),fData/(2*pi),magData)
title('Bode Diagram'),ylabel('Magnitude (dB)'),xlabel('Frequency (Hz)')
axis([1e1 1e3 -30 10]),legend('Theoretical','Measured','Location','SouthEast')

Theoretical Notch Filter:
 
Transfer function:
s^4 - 5.4e-012 s^3 + 2.842e005 s^2 - 7.702e-007 s + 2.02e010
------------------------------------------------------------
  s^4 + 26.66 s^3 + 2.846e005 s^2 + 3.789e006 s + 2.02e010
 

Clean up

warning on all
clear

Published with MATLAB® 7.7