Listing of GNU/Octave code for RC low-pass filter. Used to generated the graphs in my RC low-pass filter article. This supplements the article RC Low-pass filter.
Appendix A
Unit Step Response in GNU/Octave
GNU/Octave code for RC low-pass filter
clc; close all; clear all; format short eng R=100; # 100Ohm C=470e-9; # 470nF w=logspace(3,5,200); f=w/(2*pi); t=linspace(0,2e-3,200); p=-1/(R*C); u=1-e.^(p*t); h=plot(t,u); axis([min(t) max(t) 0 2]); #1.75 xlabel('time [s]'); ylabel('|h(t)|'); t=['Step Response(t), C=' num2str(C*1e9) 'nF, R=' num2str(R) '\Omega'] title(t, "fontsize", 15);
Frequency Response in GNU/Octave
GNU/Octave code for RC low-pass filter
clc; close all; clear all; format short eng R=100; # 100Ohm C=470e-9; # 470nF f=logspace(1,6,200); w=2*pi*f; p=-1/(R*C); u=-20*log10(sqrt(1+(w*R*C).^2)); h=semilogx(f,u); hold on; wn=-p; 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=-20*log10((fmax-f1)/f1); plot([min(f) f1 fmax],[0 0 asymp ],'k--'); hold off poles=[-p -p]; figure(1); grid off; axis([min(f) max(f) -80 40]); xlabel('frequency [Hz]'); ylabel('20log| H(t)|'); leg=[strread(num2str(R,1),'%s');'asymptote']; t=['Bode Magnitude in dB(f), C=' num2str(C*1e9) 'nF, R=' num2str(R) '\Omega']; title(t, "fontsize", 15); hold off;
Nyquist Diagram in GNU/Octave
GNU/Octave code for RC low-pass filter
clc; close all; clear all; format short eng pkg load control R=100; # 100Ohm C=470e-9; # 470nF p=-1/(R*C); H = tf([-p], [ 1 -p ]); [mag, phi, w] = bode(H); nyquist(H); h=gcf; axis ([-0.2, 1.2, -.7, .7], "square");
Square Wave in GNU/Octave
clf; t=linspace(0,2e-3,1e4); # t from 0 to 2 msec, with 10,000 steps f=2e3; # input frequency [Hz] R=100; # 100 Ohms C=470e-9; # 470 nF M=1e5; # number of harmonics ut=0; # input signal (square wave) yt=0; # output signal w=2*pi*f; # omega fc=1/(2*pi*R*C) # cutoff (-3dB) frequency for n=1:2:M, nwt = n*w*t; nwRC = n*w*R*C; ut = ut + 4/pi * sin(nwt) / n; argH = 1 / sqrt( 1 + nwRC^2 ); angH = -atan(nwRC); yt = yt + 4/pi * argH * sin(nwt + angH) / n; end plot(t*1e3,ut, 'b-',t*1e3,yt) title('Square Wave input to RC filter') xlabel('t [msec]') ylabel('[Volt]') grid on; legend('u(t)','y(t)') saveas(1,"square.svg")