Appendix A

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