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

[mystran]

Here is a basic implementation of the Cytomic SVF in C++: https://gist.github.com/hollance/2891d8 ... f0e1e55c4b

If you want, you should also be able to use the alternative bandwidth calculations for bell and shelf from the cookbook with these... though whether those are "better" are up for debate; straight BLT is too narrow at high frequencies (cramping), but the RBJ formula (in the case of the bell filter at least) is too wide, because it tries to preserve bandwidth between lower and upper band-edges which forces the lower band-edge way too low if the upper band-edge cramps against Nyquist.

One simple trick that I like that provides a sort of "middle ground" for bell filters is to compute the tan-prewarp for both the cutoff and the theoretical lower band-edge, then compute the analog bandwidth from these two prewarped frequencies to get the Q for BLT. What this does is essentially keep the lower band-edge where it should be, while ignoring the total bandwidth... so the response below cutoff is closer to what it should be, even though the response above cutoff is still heavily cramped. I'll leave the math as an exercise to reader, but it's .. actually not terribly difficult.

[JustinJ]

Here is a basic implementation of the Cytomic SVF in C++: https://gist.github.com/hollance/2891d8 ... f0e1e55c4b

Heh, I nearly posted a link to this file in your github repo a few comments back. All I'd add is that you don't necessarily need to use doubles - they work perfectly well with floats.

[Music Engineer]

Here is a basic implementation of the Cytomic SVF in C++: https://gist.github.com/hollance/2891d8 ... f0e1e55c4b

Yes - nice. Before you posted this, I was already working on my own implementation which currently looks like this:

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

Pretty much the same - which is not surprising. I might add that we don't really need to keep k and g as member variables. They can be local variables in the setup function(s). No big deal normally, but when we want to use big arrays of these filters, it may be nice to keep the object's memory footprint small.

I also added a "constant peak gain" bandpass. The original apparently corresponds to the RBJ "constant skirt gain" variant. I plotted all the responses in comparison to the RBJ cookbook responses. They are all exactly the same. This truly is the SVF-version of the RBJ filter cookbook.

What is actually the "peak" mode here (which we should not confuse with RBJs peak which is called "bell" here)? Is this a 2-pole resonator?

[Music Engineer]

Oh - and while being at it, I also did the Sallen-Key filter from here:

https://cytomic.com/files/dsp/SkfLinear ... mised2.pdf

Here it is:

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

The SKF also produces exactly the same output as the SVF in lowpass mode when one uses Q = 1 / (2 - 2*resonance) to convert between Q and resonance parameters.

[Music Engineer]

A slight variation is possible if you take the highpass directly from the solver as the "voltage" after the feedback sum, but before the first integrator, which then gives you HP, BP, LP directly as the basis (ie. without summing input). If one solves for the HP term, then BP and LP are simple to compute simply as feed-forward integrators.

The cool thing about this variation is that if you have an analog prototype (such as those given in the cookbook) of the form (a*s^2 + b*s + c)/(s^2 + s/Q + 1) then a,b,c give you the mixing gains for HP,BP,LP (respectively) directly. [...] once there's non-linearities (eg. Q limiting by additional BP feedback) the "inside the loop" HP can be a bit "better" [...]

That sounds indeed desirable. How would the coefficient-computation and sample-processing have to be modified to achieve that? I really don't know enough about circuit modeling to understand how I would "take the highpass directly from the solver" etc.

Is that variation also what is described in that paper?:

https://cytomic.com/files/dsp/SvfLinear ... dalSin.pdf

on page 9, the Mathematica output at the top for the highpass does indeed produce "m0 = 1, m1 = 0, m2 = 0" which suggests that the mixing coeffs assume HP, BP, LP as the signals to be mixed. I have also tried to implement the conversion formulas from the DF biquad coeffs - but something seems to be wrong with them. If the paper does indeed refer to a variation of the filter structure, that could be a plausible explanation for that.

[mystran]

That sounds indeed desirable. How would the coefficient-computation and sample-processing have to be modified to achieve that? I really don't know enough about circuit modeling to understand how I would "take the highpass directly from the solver" etc.

There's actually no "circuit modelling" even required here, we can do it all purely in terms of abstract signals and integrators ... and even though Andy writes in terms of circuit quantities, I think that mostly just makes things more confusing for linear filters... but that's a matter of opinion I guess. I'll post the "purely in terms of signals" version here so you can compare equivalence to Andy's version to perhaps gain more insight.

Anyway, regarding HP the idea is that you take the ODE with 2-dimensions and add a third algebraic dimension for the HP node (ie. feedback summing opamp's output) turning the thing into DAE (which makes little difference in terms of solving; it's still just Ax=b, except now you only integrate some of the values; I'm not going to write that part out right now, but perhaps another day) ... and then if you solve for this new variable first, this thing (which I've posted many times) falls out:

Code: Select all

            float hp = (in - (g+r)*z1 - z2) / (1 + g*(g+r));
            float bp = z1 + g*hp; 
            float lp = z2 + g*bp;
            // state variable update
            z1 = 2*bp - z1;   // equivalent to: z1 += 2*g*hp
            z2 = 2*lp - z2;    // equivalent to: z2 += 2*g*bp
            
            // this numbering will make sense below
            float out = a2*hp + a1*bp + a0*lp;

The coefficients are g=tan(pi*f/fs) and r=1/Q. The 1/(1+g*(g+r)) term (and optionally (g+r) term too) should be precomputed for static and perhaps even "modulation rate" filters (it's safe to interpolate, even though that results in a bit of extra damping); I left it here "as-is" just so it's easier to read and because for an audio-rate modulated filter you can just as well write it like this directly.

For emphasis, this isn't really different from Andy's filter, you can obtain one from the other simply by rearranging the terms, but I like this form because it's sort of clear in terms of what is going on: we solve the feedback sum (hp node) and then integrate twice.

The "analog prototype response" for this thing is (a2*s^2+a1*s+a0)/(s^2+s/Q+1), so if we have any analog prototype (such as those listed in RBJ Cookbook) where the denominator is already in the form 1/(s^2+s/Q+1) then basically all you have to do is take the scalar multiplies from numerator as a2,a1,a0 and you're set.

If we have a peaking filter, then RBJ gives (s^2+A/Q*s+1)/(s^2+s/(A*Q)+1) so deminator is wrong, so what do we do?!? PANIC. Fortunately we can fix this by substituting Q'=AQ and we get (s^2+A^2/Q'*s+1)/(s^2+s/Q'+1) and now it's in the correct form again when we use r=1/Q'=1/(A*Q).

What about shelves though?
We could panic again, but given RBJ lowshelf A*(s^2+sqrt(A)/Q*s+A)/(As^2+sqrt(A)/Q*s+1) we can also do a similar trick, but this with "s" where we substitute s' = s/sqrt(A) so that the denominator multiplies into the desired form and this time we then set g=tan(pi*f/fs)/sqrtA() to scale the cutoff (doing the scaling after tan() pre-warp keeps the pre-warp frequency as-is, which is what RBJ also does).

Symbolically the results are exactly the same whether you do it this way or go through the direct form coefficients, but if you're doing things numerically then the direct form coefficient solve is a bit less stable, so it's helpful to start directly from the analog prototype if possible.

[stoopicus]

This won't be any help to the OP but really simple question about expressing these as VA or SVF vs Direct Form: With VA and SVF the emphasis seems to be on delay-free feedback, but don't they just move the delay terms in to the integrators instead of the direct forward and feedback path?

i.e. aren't they equivalent? Well, I guess they have to be, as they are modeling more or less the same response. Hmm.

[mystran]

This won't be any help to the OP but really simple question about expressing these as VA or SVF vs Direct Form: With VA and SVF the emphasis seems to be on delay-free feedback, but don't they just move the delay terms in to the integrators instead of the direct forward and feedback path?

i.e. aren't they equivalent? Well, I guess they have to be, as they are modeling more or less the same response. Hmm.

The whole "zero delay" thing is bit of a misnomer, though what it really means is that "delay-free loops" are allowed in the signal-flow graph and the significance of this is that rather than building filters out of unit delays, we can build filters out of "implicit" numerical integrators.

Implicit integration methods are useful, because unlike explicit methods, they can be A-stable, which essentially means that every stable continuous-time system maps to a stable discrete-time system without any additional precautions. Trapezoidal method is particular nice, because besides being the most accurate A stable linear multi-step method (ie. the "best" simple method), it also maps unstable system to unstable system and in fact preserves the response exactly up to frequency warping (by arctan, so we can compensate); it is equivalent to BLT.

Ordinarily implicit methods present a problem, because given that they contain a "delay-free path" from input to output, if you put them into feedback loops you have a "delay-free loop" where you cannot pick any order of operations that allows you to actually implement the thing directly.

If you look at classical VA filters like Chamberlin or the various transistor-ladder models in literature before this whole "ZDF" business, then what you see is that they are still built out of integrators, but they have to add additional unit delays into the loop (or in case of Chamblin, make at least one of the integrators explicit) which then runs into all the stability problems of explicit methods... and a lot has been published about how to minimize these issues under the assumption that every loop must have a delay.

So the two insights of the whole ZDF business are that (1) we can solve the simultaneous equations instead and (2) this often isn't nearly as expensive as you might think. The first one is fairly obvious ("to someone skilled in the art"), but the second one isn't and it's the second insight that's most likely the reason why we've only been using these recently. Classically the whole approach was dismissed as "requires solving an implicit system" with the assumption that this would necessarily be too expensive.

Solving a matrix equation Ax=b is expensive in the general case, sure... but when the system is small and very sparse (ie. most of the matrix is just zeroes) and especially if we precompute the divisions as coefficients (or have modern hardware where division isn't crazy expensive) ... then it's actually not bad at all. In fact, with something like an SVF using trapezoidal method, it turns out that it's roughly as much computation as a direct form biquad...

... but that's very much not obvious if you start from the classical midset of "implicit methods are expensive because they require solving the matrix equation Ax=b" and because the insights weren't quite obvious and this suddenly allowed us to completely ignore the old "every loop must contain a delay" rule ... what is essentially just trapezoidal integration with a few optimisations (eg. unrolling the solver and using TDF2-form integrators taking advantage of a fixed time-step) got a bunch of names, including "zero-delay filters."

So "zero-delay" here refers to loops in the signal-flow graph, rather than whether or not there are any delays involved at all.

[stoopicus]

Thank you, that's helpful. I was coming at this from the other direction in that I was looking at how (in this case) Pirkle expresses the VA/ZDF topology in his book and comparing it to TDF2 and thinking that while they were clearly different, they weren't *that* different (as in, considering it is a different technique, I was expecting more differences.)

I still need to read Zavalishin's book, which I am guessing would clear this up for me more completely.

[mystran]

Thank you, that's helpful. I was coming at this from the other direction in that I was looking at how (in this case) Pirkle expresses the VA/ZDF topology in his book and comparing it to TDF2 and thinking that while they were clearly different, they weren't *that* different (as in, considering it is a different technique, I was expecting more differences.)

I want to add that when we talk about "TDF2" in the context of ZDF, that's different from "TDF2" in the direct form sense. The trapezoidal integrator itself can be written in a form that's essentially a one-pole DF1 or (in the case of fixed time-step) TDF2.. but that's just a way to express the integrator itself as part of the signal-flow (which is solved as simultaneous equations) and there's no actual "TDF2 filters" as such in a ZDF scheme.

[stoopicus]

Yes - this is what I suspect Pirkle is doing (ch12 of his effects book if you have it). He's expressing the integrator itself that way, if I am reading it right.

[mystran]

I should probably add that the standard trapezoidal rule leads to DF1 form.
We compute h[n]/2*(x[n]+x[n-1]) where h[n] is the time-step which can vary, then add this to the previous output y[n-1] to obtain the new output y[n].

We can't just take any random direct form here, because in the time-varying case the values we actually store as state variables are important and we actually need trapezoidal rule in the time-domain, not just the correct response.

However, if we assume that h[n] is constant (ie. this is valid if the time-step is fixed), then what we can do instead is store as state:
u[n] = y[n] + h[n]/2*x[n]

So for the next time-step, we can then compute (well, really solve from the simultaneous equations) output (before updating state as above separately) as
y[n] = u[n-1] + h[n]/2*x[n] == ( y[n-1] + h[n-1]/2*x[n-1] ) + h[n]/2*x[n-1]

This is the TDF2 form. I left the indexes for h[n] there in order to emphasize that this is equivalent to the standard trapezoidal rule if and only if h[n] == h[n-1], that is the time-step is fixed... but for audio purposes we generally want fixed time-step anyway, because non-uniform sampling is such a huge pain, so TDF2 is generally a win, because now we don't need to store previous input and previous output separately.

The other two direct forms (DF2 and TDF1) are not terribly useful here, because they only give us correct trapezoidal rule in the LTI case... which kinda spoils all the fun.

[Music Engineer]

The "analog prototype response" for this thing is (a2*s^2+a1*s+a0)/(s^2+s/Q+1), so if we have any analog prototype (such as those listed in RBJ Cookbook) where the denominator is already in the form 1/(s^2+s/Q+1) then basically all you have to do is take the scalar multiplies from numerator as a2,a1,a0 and you're set.

If we have a peaking filter, then RBJ gives (s^2+A/Q*s+1)/(s^2+s/(A*Q)+1) so deminator is wrong, so what do we do?!? PANIC. Fortunately we can fix this by substituting Q'=AQ and we get (s^2+A^2/Q'*s+1)/(s^2+s/Q'+1) and now it's in the correct form again when we use r=1/Q'=1/(A*Q).

What about shelves though?
We could panic again, but given RBJ lowshelf A*(s^2+sqrt(A)/Q*s+A)/(As^2+sqrt(A)/Q*s+1) we can also do a similar trick, but this with "s" where we substitute s' = s/sqrt(A) so that the denominator multiplies into the desired form and this time we then set g=tan(pi*f/fs)/sqrtA() to scale the cutoff (doing the scaling after tan() pre-warp keeps the pre-warp frequency as-is, which is what RBJ also does).

Thanks for elaborating! I implemented it today. Here it is:

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

works perfectly. Much leaner and cleaner than my older first implementation that I did back in the day when Vadim's book came out. I posted this here:

viewtopic.php?p=5243733#p5243733

The per-sample update algorithm is the same but my coefficient calculations were quite a mess back then (at least for some of the modes) because I parameterized it in terms of bandwidth (in octaves) rather than Q. And the API was rubbish, too. The new code is much better. I also agree that it's nicer to mix highpass, bandpass and lowpass rather than input, bandpass and lowpass.

[stoopicus]

I should probably add that the standard trapezoidal rule leads to DF1 form.

yep I think I see that now, thanks again.

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.

Ok, time for me to just actually read Vadim's book. Thanks again for the help above.

[mystran]

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.