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);