A MATLAB script that derives a computable low pass filter from a circuit prototype. Derivation and discretisation performed with symbolic toolbox.

LowPassFilter.m

%%% Low Pass Circuit
% A computable RC circuit model using a recurrence relation derived from
% MNA, through to a transfer function which is then discretised using the
% symbolic math toolbox.

% Author: Ben Holmes
% Date: 2019/01/01
% License: GPL v3

clear;

%% MNA

% Laplace variable
syms s

% Resistor, capacitor, input, and output connections
N_R = [1 -1];
N_C = [0  1];
N_U = [1  0];
N_O = [0  1];

% Cutoff frequency
fc = 100;

% Impedances
R = 1;
C = 1/(2*pi*R*fc);

% Conductances
G_R = 1/R;
G_C = s*C;

% System matrix
S = [(N_C.'*G_C*N_C)+(N_R.'*G_R*N_R) N_U'; N_U 0];

% Input representation
X = [0; 0; 1];

%% Transfer function

% Solve for transfer functions
trfs = S\X;

% Input to output voltage relation
Vo_Vi = [N_O 0]*trfs;

% Symbolic function for fast execution (as opposed to substituting in
% values)
symTFfunc = matlabFunction(Vo_Vi);

% Number of samples and numeric sample rate
Ns = 48e3;
fs_ = 48e3;

% Frequency vector
f_vec = (0:Ns-1).*(fs_/Ns);

% Response calculation
tf_response = zeros(1,Ns);
for nn=1:Ns
    tf_response(nn) = symTFfunc(2*pi*sqrt(-1)*f_vec(nn));
end

% Plot frequency response
figure(1);
clf;
semilogx(f_vec,20*log10(abs(tf_response)));
axis tight
xlim([2 20e3])
xlabel('Frequency (Hz)');
ylabel('Amplitude (dB)');
grid on;

%% Discretise

% Discrete frequency domain variable z and symbolic sample rate
syms z fs;

% Discrete low pass transfer function
lowpassD = subs(Vo_Vi,s,2*fs*(1-z)/(1+z));

% Extract numerator and denominator from tf
[num, den] = numden(lowpassD);

% Find the coefficients of the polynomials of z
[b, bpolys] = coeffs(num,z);
[a, apolys] = coeffs(den,z);

% Normalise
a = a./b(1);
b = b./b(1);

%% Impulse response from filtering in the time domain

% Time vector
t_vec = (0:Ns-1)./fs_;

% Impulse
x = [1 zeros(1,Ns-1)];

% Numeric a and b values
b_ = eval(subs(fliplr(b),fs,48e3));
a_ = eval(subs(fliplr(a),fs,48e3));

% Run filter
y = zeros(1,Ns);
for nn=1:Ns
    if nn>1
        y(nn) = (b_(1)*x(nn) + b_(2)*x(nn-1) - a_(2)*y(nn-1))/a_(1);
    else
        y(nn) = b_(1)*x(nn)/a_(1);
    end
end

figure(1);
hold on;
semilogx(f_vec, 20*log10(abs(fft(y))));

legend('Continuous','Discrete @ 48 kHz');