Appendix

Coert Vonk

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];
    legend(strread(num2str(Rvector,3),'%s'));
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)|');
leg=strread(num2str(Rvector,4),'%s');
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);