3
$\begingroup$

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:

enter image description here

enter image description here

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;
}
$\endgroup$
  • 1
    $\begingroup$ Your confusing time and frequency domain. You should NEVER have to divide. This is a simple difference equation, the "denominator" should be subtracted, not divided. Try your code first with a very simple first order filter $\endgroup$ – Hilmar Jan 5 at 22:45
  • 1
    $\begingroup$ @Hilmar: At a first glance I thought the same, but the denominator is just a constant, so it's correct. $\endgroup$ – Matt L. Jan 6 at 6:30
3
$\begingroup$

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.

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Thanks Matt! I followed that approach and also the "how to apply a biquad" here: dsprelated.com/freebooks/filters/BiQuad_Section.html and it seems to have worked nicely. I haven't hooked it up to a spectrogram but it's functioning and outputting audio that sounds like one would expect. Appreciate the guidance. It's much easier when you do it with standardized coefficients like that and it also makes it more efficient because then I only have to calculate the sets of coefficients when changing parameters. Obviously this makes it much cheaper per sample for processing. Thanks again. $\endgroup$ – mike Jan 5 at 23:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.