What did I do wrong with this simple filter build?

Question

I tried to put everything I have learned from people here together to code my first filter from scratch. Unfortunately, it didn't go well and I'm not getting the expected output. The math/code became quite big and messy, which makes it hard to read. But I went through it a few times and couldn't find an error. I'm not sure what I've done wrong. I feel like perhaps I have a big idea wrong about how to do this. Maybe I took a wrong turn somewhere.

Here is the circuit:

This is for a physical modeling purpose, so it is necessary to have everything in terms of the component resistor, capacitor, and inductor values as drawn.

I attempted to work through this by first solving the impedance of the parallel and series components, working out a Laplace transfer function, then substituting $s=\frac{1-z^{-1}}{T}$, expressing it in terms of per sample input/output delay, and then writing some simple C++ for it.

But the output is just an impulse of sound then it overloads. It's not filtering as expected. Any ideas what I've done wrong?

I've tried to write this out as clearly as I can in terms of the steps I've taken in the hopes it might make sense.

Parallel Component:

$\frac{1}{R_1} = \frac{1}{sL} + \frac{1}{R_d+\frac{1}{sC}}$

$R_1 = \frac{1}{\frac{1}{sL} + \frac{1}{R_d+\frac{1}{sC}}}$

$R_1 = \frac{Ls(CR_ds + 1)}{CLs^{2} + CR_ds+1}$

Series Component:

$R_2 = 2R_g$

Transfer Function:

$V_{out} = \frac{R_2}{R_2+R_1} V_{in}$

$V_{out} = \frac{2R_g}{2R_g + \frac{Ls(CR_ds + 1)}{CLs^{2} + CR_ds + 1}} V_{in}$

$V_{out}(s) = \frac{2R_g(CLs^{2} + CR_ds + 1)}{2CR_gLs^{2} + 2CR_gR_ds+CLR_ds^{2} + 2R_g + Ls} V_{in}(s)$

Substituting $s=\frac{1-z^{-1}}{T}$:

$_{Numerator} = \frac{-2 R_g (-C L z^{-2} + 2 C L z^{-1} - C L + C R_d T z^{-1} - C R_d T - T^{2})}{T^{2}}$

$_{Denominator} = \frac{2 C R_g L z^{-2} - 4 C R_g L z^{-1} + 2 C R_g L + 2 C R_g T R_d - 2 C R_g T R_dz^{-1} + C L R_d- 2C L R_dz^{-1}+C L R_dz^{-2} + 2 R_g T^2 - L T z^{-1} + L T}{T^2}$

Canceling the $1/T^{2}$:

$_{Numerator} = -2 R_g (-C L z^{-2} + 2 C L z^{-1} - C L + C R_d T z^{-1} - C R_d T - T^{2})$

$_{Denominator} = 2 C R_g L z^{-2} - 4 C R_g L z^{-1} + 2 C R_g L + 2 C R_g T R_d - 2 C R_g T R_dz^{-1} + C L R_d-2C L R_dz^{-1}+C L R_dz^{-2} + 2 R_g T^2 - L T z^{-1} + L T$

Cross Multiplying:

$_{Leftside} = 2 C R_g L V_{out}[n-2] - 4 C R_g L V_{out}[n-1] + 2 C R_g LV_{out}[n] + 2 C R_g T R_dV_{out}[n] - 2 C R_g T R_dV_{out}[n-1] + C L R_dV_{out}[n]-2C L R_dV_{out}[n-1]+C L R_dV_{out}[n-2] + 2 R_g T^2V_{out}[n] - L T V_{out}[n-1] + L TV_{out}[n]$

$_{Leftside} = V_{out}[n] (2 C R_g L + 2 C R_g T R_d +C L R_d + 2 R_g T^2 + L T) + 2 C R_g L V_{out}[n-2] - 4 C R_g L V_{out}[n-1] - 2 C R_g T R_dV_{out}[n-1] -2C L R_dV_{out}[n-1]+C L R_dV_{out}[n-2] - L T V_{out}[n-1]$

$_{Rightside} = -2 R_g (-C L V_{in}[n-2] + 2 C L V_{in}[n-1] - C LV_{in}[n] + C R_d T V_{in}[n-1] - C R_d TV_{in}[n] - T^{2}V_{in}[n])$

Final Equation:

$V_{out}[n] (2 C R_g L + 2 C R_g T R_d +C L R_d + 2 R_g T^2 + L T) = -2 R_g (-C L V_{in}[n-2] + 2 C L V_{in}[n-1] - C LV_{in}[n] + C R_d T V_{in}[n-1] - C R_d TV_{in}[n] - T^{2}V_{in}[n]) - (2 C R_g L V_{out}[n-2] - 4 C R_g L V_{out}[n-1] - 2 C R_g T R_dV_{out}[n-1] -2C L R_dV_{out}[n-1]+C L R_dV_{out}[n-2] - L T V_{out}[n-1])$

$V_{out}[n] = \frac{-2 R_g (-C L V_{in}[n-2] + 2 C L V_{in}[n-1] - C LV_{in}[n] + C R_d T V_{in}[n-1] - C R_d TV_{in}[n] - T^{2}V_{in}[n]) - (2 C R_g L V_{out}[n-2] - 4 C R_g L V_{out}[n-1] - 2 C R_g T R_dV_{out}[n-1] -2C L R_dV_{out}[n-1]+C L R_dV_{out}[n-2] - L T V_{out}[n-1])}{2 C R_g L + 2 C R_g T R_d +C L R_d + 2 R_g T^2 + L T }$

Code:

class PhysicalFilter{

public:

void setSampleRate(double sampleRateIn){
T = 1/sampleRateIn;
}

float filterSample(float inputSample, float C, float L, float R_d, float R_g){

input_2 = input_1;
input_1 = input;
input = inputSample;

output_2 = output_1;
output_1 = output;

float numerator = -2 * R_g * ((-C * L * input_2) + (2 * C * L * input_1) - (C * L * input) + (C * R_d * T * input_1) - (C * R_d *  T * input) - (T * T * input)) - ((2 * C * R_g * L * output_2) - (4 * C * R_g * L * output_1) - (2 * C * R_g * T * R_d * output_1) - (2 * C * L * R_d * output_1) + (C * L * R_d * output_2) - (L * T * output_1));

float denominator = (2 * C * R_g * L) + (2 * C * R_g * T * R_d) + (C * L * R_d) + (2 * R_g * T * T) + (L * T);

output = numerator/denominator; 

return output;
}

private:
float input = 0.f;
float input_1 = 0.f;
float input_2 = 0.f;
float output = 0.f;
float output_1 = 0.f;
float output_2 = 0.f;
float C = 1.f;
float L = 1.f;
float R_g = 1.f;
float R_d = 1.f;
float T = 1/44100.f;
}

Answer 1

Your analog transfer function looks OK. For the sake of clarity - and to reduce the chance of making errors - I'd just rewrite it as

$$H_a(s)=G\cdot\frac{s^2+as + b}{s^2+cs + d}\tag{1}$$

with

$$\begin{align}G&=\frac{2R_g}{R_d+2R_g}\\a&=\frac{R_d}{L}\\b&=\frac{1}{LC}\\c&=G\left(a+\frac{1}{2R_gC}\right)\\d&=G\cdot b\frac{}{}\end{align}$$

Then you can use the backward Euler transformation on the general second-order transfer function given by $(1)$, resulting in

$$H_d(z)=\frac{G}{1+cT+dT^2}\cdot\frac{1+aT+bT^2-(2+aT)z^{-1}+z^{-2}}{1-\frac{2+cT}{1+cT+dT^2}z^{-1}+\frac{1}{1+cT+dT^2}z^{-2}}\tag{2}$$

Now you can check if your discrete-time transfer function is correct.

The next step is to check your code. I would suggest to rewrite it to implement a general biquad transfer function:

$$H(z)=\frac{b_0+b_1z^{-1}+b_2z^{-2}}{1+a_1z^{-1}+a_2z^{-2}}\tag{3}$$

You can test your biquad routine by supplying some known coefficients $b_i$ and $a_i$ (e.g., for a standard low pass filter, etc.). As soon as you're convinced that the routine works properly, you can test it with the coefficients of your design. With what you have now you can't really separate the design and the implementation, so testing and debugging becomes problematic.