How can I learn how to implement filters from scratch? (C++)

[stoopicus]

I should clarify that the form Pirkle is using for the integrator is exactly the same as Vadim uses on the cover of his book.

That's the TDF2 form that is valid under "fixed time-step" assumption. That's what you normally want.

Yep, coming full circle, that's what I was getting at with my original question above around how it looked a lot to me like it was actually equivalent to TDF2.

[Music Engineer]

I'd really like to add some sort of "getTransferFunctionAt(complex z)" to the SVF that somehow computes the value of the transfer function for a given z but I'm somehow stumped. I got as far as writing down the following:

Code: Select all

  h[n] = s*x[n] - s*c*u[n-1] - s*v[n-1]   <-->   H(z) = s*X(z) - s*c*z^-1 * U(z) - s*z^-1 * V(z)
  b[n] = g*h[n] + u[n-1]                  <-->   B(z) = g * H(z) + z^-1 * U(z)
  l[n] = g*b[n] + v[n-1]                  <-->   L(z) = g * B(z) + z^-1 * V(z)
  u[n] = 2*b[n] - u[n-1]                  <-->   U(z) = 2 * B(z) - z^-1 * U(z)
  v[n] = 2*l[n] - v[n-1]                  <-->   V(z) = 2 * L(z) - z^-1 * V(z)

where h[n], b[n], l[n] are the time domain highpass, bandpass and lowpass signals and u[n], v[n] are the integrator states and H(z), etc. are the corresponding z-transforms of these signals (Warning: in this notation, H(z) is *not* a transfer function but the z-trafo of the highpass signal h[n]). I renamed the variables a bit compared to the source code to make them shorter (because expressions will get long and messy anyway). So far so correct? But what now? I'd like to find expressions for H(z)/X(z), B(z)/X(z) and L(z)/X(z). These would then be the transfer functions of the highpass, bandpass and lowpass. Right? To obtain the overall transfer function, I just form a weighted sum with the a-coefficients. But so far, I had no success finding these expressions. It gets messy really quick and I'm not really sure, how I could explain this task to a computer algebra system (like Sage - which is what I typically use). How would one go about this? I tried the following snippets:

Code: Select all

var("d g s c X H B L U V")              # d = z^-1
e1 = H == s*X - s*c*U*d - s*V*d
e2 = B == g*H + U*d
e3 = L == g*B + V*d
e4 = U == 2*B - U*d
e5 = V == 2*L - V*d
solve([e1,e2,e3,e4,e5],[H/X,B/X, L/X])

  ->  TypeError: H/X is not a valid variable.


var("d g s c X H B L U V") 
e1 = H == s*X - s*c*U*d - s*V*d
e2 = B == g*H + U*d
e3 = L == g*B + V*d
e4 = U == 2*B - U*d
e5 = V == 2*L - V*d
HX = H/X
BX = B/X
LX = L/X
solve([e1,e2,e3,e4,e5],[HX,BX, LX])

  -> TypeError: H/X is not a valid variable.


var("d g s c X H B L U V HX BX LX") 
e1 = H  == s*X - s*c*U*d - s*V*d
e2 = B  == g*H + U*d
e3 = L  == g*B + V*d
e4 = U  == 2*B - U*d
e5 = V  == 2*L - V*d
e6 = HX == H/X
e7 = BX == B/X
e8 = LX == L/X
solve([e1,e2,e3,e4,e5,e6,e7,e8],[HX,BX,LX])

  -> []

here:

https://sagecell.sagemath.org/ The stuff after the arrow -> are the results I get - which are no results. I'm using d = z^-1, btw. I'm doing something wrong. Am I even on the right track or is this all pure nonsense?

[stoopicus]

I'd really like to add some sort of "getTransferFunctionAt(complex z)" to the SVF that somehow computes the value of the transfer function for a given z but I'm somehow stumped.

Do you mean as a function of ω or of z? In terms of something that is not immediately determinable without additional work, I would want this as a function of ω (i.e. return the complex magnitude and phase response at any given ω). Isn't the transfer function as a function of z simply the output of the filter?

(Not asking rhetorically; I am assuming I am the one missing something obvious here)

[Music Engineer]

Do you mean as a function of ω or of z? In terms of something that is not immediately determinable without additional work, I would want this as a function of ω (i.e. return the complex magnitude and phase response at any given ω).

In general, as function z, i.e. for an arbitrary point in the z-plane. This can then be used to compute the frequency response in terms of w by just passing z = exp(i*w).

Isn't the transfer function as a function of z simply the output of the filter?

The transfer function is the ratio of the z-transforms of the output signal and the input signal of the filter, i.e. H(z) = Y(z)/X(z). For a simple direct form first order filter, we can derive it as follows:

Code: Select all

  y[n] = b0*x[n] + b1*x[n-1] - a1*y[n-1]                Translate to z-domain
  
  Y(z) = b0*X(z) + b1*z^-1*X(z) - a1*z^-1*Y(z)          Move Y(z) terms to left and X(z) to right
  Y(z) + a1*z^-1*Y(z) = b0*X(z) + b1*z^-1*X(z)          Factor out Y(z) and X(z)
  Y(z) (1 + a1*z^-1) = X(z)  (b0 + b1*z^-1)             Solve for X(z)/Y(z) and call that H(z)
  
  H(z) = Y(z)/X(z) = (b0 + b1*z^-1) / (1 + a1*z^-1)     Done

And this generalizes straightforwardly to higher order direct form filters because you can just as easily factor out Y(z) and X(z) there - just the polynomials in z^-1 that constitute the other factor will get longer. But for a more complicated filter like the SVF, it is not so clear anymore how one would do this in general. There are so many intermediate signals and multiple feedback paths that makes it a total mess to try to bring it into this form. It's all intermingled such that you cannot just simply factor out X(z) and Y(z) anymore. I suppose there must be some way to let a computer algebra system solve this in general - I just haven't figured out how, i.e. what commands I would have throw at the CAS to get what I want.

[mystran]

But for a more complicated filter like the SVF, it is not so clear anymore how one would do this in general.

The most straight-forward method that works with whatever structure is to write it in state-space form, from which you'll basically get it as C((zI-A)^-1)B+D (see https://en.wikipedia.org/wiki/State-spa ... esentation for which matrix is which).

Wikipedia explains the whole thing with more readable math than what I can do here. It's written in terms of continuous time filters, but it works exactly the same for discrete time too. If you write the state-space in terms of delays, then just swap z in place of s.

For trapezoidal filters though, you could just as well take the continuous-time state-space directly and do the BLT substitution s<-2/T*(z-1)/(z+1). However.. in the case of the SVF, we already know that the continuous-time response is (a2*s^2+a1*s+a0)/(s^2+s/Q+1) so you can just as well use that for the BLT substitution.

[Music Engineer]

OK - thanks. I will look into this.

[S0lo]

But for a more complicated filter like the SVF, it is not so clear anymore how one would do this in general. There are so many intermediate signals and multiple feedback paths that makes it a total mess to try to bring it into this form. It's all intermingled such that you cannot just simply factor out X(z) and Y(z) anymore. I suppose there must be some way to let a computer algebra system solve this in general - I just haven't figured out how, i.e. what commands I would have throw at the CAS to get what I want.

You can derive 3 frequency domain equations directly from the SVF topology. Which then you can use to solve for the transfer function(s). Replace the integrators with a factor of wc/s (in case you want the prewarp built in, use: g/s, I think)

[Music Engineer]

OK - I finally figured it out by a mix of manual derivations and letting the CAS do certain steps for me:

https://github.com/RobinSchmidt/RS-MET/ ... ctions.txt

...but I really think, there should be a way to let the CAS do things like that automatically. The thing is that I do not want to solve for any (subset) of the variables but rather for an expression in some of the variables - namely, for a quotient of two of them. I tried to define this quotient as an additional variable and then set up an additional equation for that (like H = Y/X) but somehow, that also didn't work.

[S0lo]

OK - I finally figured it out by a mix of manual derivations and letting the CAS do certain steps for me:

https://github.com/RobinSchmidt/RS-MET/ ... ctions.txt

...but I really think, there should be a way to let the CAS do things like that automatically. The thing is that I do not want to solve for any (subset) of the variables but rather for an expression in some of the variables - namely, for a quotient of two of them. I tried to define this quotient as an additional variable and then set up an additional equation for that (like H = Y/X) but somehow, that also didn't work.

I believe the final low pass form should be something like this:

ddd.png

Give or take. Thats assuming your we're talking about the exact same SVF mentioned in Vadim's book.

[Music Engineer]

But that's the s-domain transfer function and I'm talking about the z-domain transfer function here and all that I have available to compute it are the member variables of this class:

https://github.com/RobinSchmidt/RS-MET/ ... leFilter.h

For the derivations, I renamed scl to s and gpr to c. I don't assume any knowledge about the original design parameters. It's actually possible to reconstruct the design parameters from the filter's coefficients - but it's a bit tricky because the mapping is different for different filter types, so one must first figure out, which type it is (which is indeed possible by looking at the pattern of the a0, a1, a2 coeffs) - I dabbled with the idea of reconstructing the s-domain design parameters from the filter coeffs - with some good success actually - and if you have them back, you can evaluate H(s) with an appropriate bilinear transform formula to compute s from z - but I consider that cheating. I made it almost work, though.

[mystran]

But that's the s-domain transfer function and I'm talking about the z-domain transfer function here

Substitute BLT s<-2/T*(z-1)/(z+1) and ask your CAS to simplify.

[S0lo]

But that's the s-domain transfer function and I'm talking about the z-domain transfer function here and all that I have available to compute it are the member variables of this class:

https://github.com/RobinSchmidt/RS-MET/ ... rMystran.h

Ahh My bad. Here is the one I got a few years back:

zdom.png

Edit: got it using BLT.

[stoopicus]

In general, as function z, i.e. for an arbitrary point in the z-plane. This can then be used to compute the frequency response in terms of w by just passing z = exp(i*w).
...
The transfer function is the ratio of the z-transforms of the output signal and the input signal of the filter, i.e. H(z) = Y(z)/X(z).

Ahh, derp. Yes, obviously, my mistake. Sorry about that.

[JustinJ]

But that's the s-domain transfer function and I'm talking about the z-domain transfer function here and all that I have available to compute it are the member variables of this class:

https://github.com/RobinSchmidt/RS-MET/ ... rMystran.h

For the z-transform, here's what I'm using. My variables are as named the same as in the SvfLinearTrapOptimised2 paper.

Code: Select all

float SVFFilter::magnitude(double frequency, double sampleRate) const
{
    auto z = std::exp(std::complex<double>(0.0, -2.0 * M_PI) * frequency / sampleRate);

    const double gsq = g * g;
    const auto zsq = z * z;
    const auto twoz = z * 2.0;
    const double gm1 = g * m1;
    const double gk = g * k;
    const double twogsq = gsq * 2.0;

    auto numerator = gsq * m2 * (zsq + twoz + 1.0) - gm1 * (zsq - 1.0);
    auto denominator = gsq + gk + zsq * (gsq - gk + 1.0) + z * (twogsq - 2.0) + 1.0;

    return std::abs(static_cast<double>(m0) + (numerator / denominator));
}

Shouldn't be hard to alter it to work with what you've got.

[Music Engineer]

OK - thanks! Meanwhile, I have finally figured out how to let Sage do the whole messy derivation. This code works:

Code: Select all

var("d  s c g  X H B L U V  TF_H TF_B TF_L")
e1 = H == s*X - s*c*U*d - s*V*d
e2 = B == g*H + U*d
e3 = L == g*B + V*d
e4 = U == 2*B - U*d
e5 = V == 2*L - V*d
e6 = TF_H == H/X
e7 = TF_B == B/X
e8 = TF_L == L/X
sH = solve([e1,e2,e3,e4,e5, e6],[TF_H, H,B,L,U,V])
sB = solve([e1,e2,e3,e4,e5, e7],[TF_B, H,B,L,U,V])
sL = solve([e1,e2,e3,e4,e5, e8],[TF_L, H,B,L,U,V])
sH[0][0], sB[0][0], sL[0][0]

and it gives the result:

Code: Select all

(TF_H ==  (d^2     - 2*d       + 1)*s  /  (d^2 - 2*((c*g - g^2)*d^2 - (c*g + g^2)*d)*s - 2*d + 1),
 TF_B == -(d^2*g   -             g)*s  /  (d^2 - 2*((c*g - g^2)*d^2 - (c*g + g^2)*d)*s - 2*d + 1),
 TF_L ==  (d^2*g^2 + 2*d*g^2 + g^2)*s  /  (d^2 - 2*((c*g - g^2)*d^2 - (c*g + g^2)*d)*s - 2*d + 1))

My mistake was that I tried to let Sage only to solve for the variable that I was interested in. But with n equations, one must solve for n variables and then throw away the results for the variables that one is not interested in. When thinking about it, it's pretty obvious - trying to solve for less variables then there are equations gives an overdetermined system, i.e. more equations than degrees of freedom. That will in general give an empty solution set (except in exceptional situations) - which is what I've got (after fixing the more stupid type errors). I made the mistake of assuming that solving for just one variable is somehow "easier" than solving for multiple variables. I have learned a new thing about using a CAS effectively (and learned it the hard way). I updated my textfile with the derivation accordingly:

https://github.com/RobinSchmidt/RS-MET/ ... ctions.txt