# Appendix

Listing of GNU/Octave code for RC low-pass filter. Used to generated the graphs in my RC low-pass filter article. Supplements the article RLC Low-pass Filter.

## Appendix

### Step response, in GNU/Octave

clc; close all; clear all; format short eng
L=47e-3; # 47mH
C=47e-9; # 47nF
#Rvector = [3.9e3]; # separate real poles
#Rvector = [2e3]; # coinciding real poles
#Rvector = [220 820 1500]; # conjugate complex poles
Rvector = [0]; # conjugate complex poles on im-axis
w=logspace(3,5,200); #*e.^(sigma*t)
f=w/(2*pi);

for R = Rvector
wn=1/sqrt(L*C);
zeta=R/2*sqrt(C/L)
t=linspace(0,2e-3,200);
if (zeta < 1) # complex conjugate poles
sigma=wn*zeta;
wd=wn*sqrt(1-zeta^2);
u=1 - sqrt(sigma^2+wd^2)/wd .* exp(-sigma*t) .* cos(wd*t-atan(sigma/wd));
hold on; h=plot(t,u); hold off;
endif
if (zeta == 1) # coinciding real poles
p=wn*(-zeta+sqrt(zeta^2-1));
u=1 + (p*t-1).*e.^(p*t);
hold on; h=plot(t,u); hold off;
endif
if (zeta > 1) # real poles
p1=wn*(-zeta+sqrt(zeta^2-1));
p2=wn*(-zeta-sqrt(zeta^2-1));
u = 1 + p2/(p1-p2)*e.^(p1*t) + p1/(p2-p1)*e.^(p1*t);
hold on; h=plot(t,u); hold off;
endif
endfor
axis([min(t) max(t) 0 2]); #1.75
xlabel('time [s]'); ylabel('|h(t)|');

t=['Step Response(t'];
t2=['), C=' num2str(C*1e9) 'nF, L=', num2str(L*1e3),'mH'];
if(length(Rvector)==1)
t=[t t2 ', R=' num2str(R/1e3) 'k\Omega']
else
t=[t ',R' t2];
endif
t = [t ];
title(t, "fontsize", 15);

### Bode magnitude, in GNU/Octave

clc; close all; clear all; format short eng
L=47e-3; # 47mH
C=47e-9; # 47nF
#Rvector = [3.9e3]; # separate real poles
#Rvector = [2e3]; # coinciding real poles
#Rvector = [220 820 1500]; # conjugate complex poles
Rvector = [0]; # conjugate complex poles on im-axis
f=logspace(1,6,200);
w=2*pi*f;

for R = Rvector
wn=1/sqrt(L*C);
zeta=R/2*sqrt(C/L)
if (zeta == 0) # complex conjugate poles on imaginary axis
u=20*log10(wn.^2) - 20*log10(wn.^2-w.^2);
hold on; h=semilogx(f,u); hold off;
hold on;
plot([wn/(2*pi) wn/(2*pi)], get(gca,'YLim'),'k--');
text(wn/(2*pi)*1.1,30,'|p|/2\pi');
hold off
poles = [-wn*zeta+sqrt(1-zeta^2)*j, -wn*zeta-sqrt(1-zeta^2)*j ];
endif
if (zeta < 1 && zeta > 0) # complex conjugate poles on left side of s-plane
u=20*log10(wn.^2) - 20*log10(sqrt((wn.^2-w.^2).^2+(2*zeta*wn*w).^2));
hold on; h=semilogx(f,u); hold off;
hold on;
plot([wn/(2*pi) wn/(2*pi)], get(gca,'YLim'),'k--');
text(wn/(2*pi),25,'\omega_n/2\pi');
hold off
poles = [-wn*zeta+sqrt(1-zeta^2)*j, -wn*zeta-sqrt(1-zeta^2)*j ];
endif
if (zeta == 1) # coinciding real poles
p=1/sqrt(L*C);
u=20*log10(p.^2) - 40*log10(sqrt(w.^2+p.^2));
hold on; h=semilogx(f,u);
plot([wn/(2*pi) wn/(2*pi)], get(gca,'YLim'),'k--');
text(wn/(2*pi),5,'|p|/2\pi');
f1=p/(2*pi);
fmax=max(f);
asymp=-40*log10((fmax-f1)/f1);
plot([min(f) f1 fmax],[0 0 asymp ],'k--');
hold off
poles=[-wn*zeta -wn*zeta];
endif
if (zeta > 1) # separate real poles
p1=wn*(-zeta+sqrt(zeta^2-1));
p2=wn*(-zeta-sqrt(zeta^2-1));
u=20*log10(p1*p2./(sqrt(w.^2+p1.^2).*sqrt(w.^2+p2.^2)));
figure(1);
hold on; h=semilogx(f,u); hold off;
hold on;
f1=-p1/(2*pi); f2=-p2/(2*pi);
plot([f1 f1], get(gca,'YLim'),'k--');
plot([f2 f2], get(gca,'YLim'),'k--');
fmax=max(f);
asymp1=0 - 20*log10((f2-f1)/f1);
asymp2=asymp1 - 40*log10((fmax-f2)/f2);
plot([min(f) f1 f2 fmax],[0 0 asymp1 asymp2],'k--');
text(f1,5,'|p1|/2\pi');
text(f2,5,'|p2|/2\pi');
hold off
poles=[p1+0j p2+0j];
endif
endfor
figure(1);
grid off;
axis([min(f) max(f) -80 40]);
xlabel('frequency [Hz]'); ylabel('20log| H(t)|');
if (zeta>=1)
leg=[leg;'asymptote'];
endif
t=['Bode Magnitude in dB(f'];
t2=['), C=' num2str(C*1e9) 'nF, L=', num2str(L*1e3),'mH'];
if(length(Rvector)==1)
t=[t t2 ', R=' num2str(R/1e3) 'k\Omega']
else
t=[t ',R' t2];
legend(leg);
endif
t = [t ];
title(t, "fontsize", 15);

This site uses Akismet to reduce spam. Learn how your comment data is processed.