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