Moog Ladder Filter - Research Paper

[pscorbett]

Hello! I have been "lurking" this forum for a while, but haven't posted. I'm not sure of all the etiquette yet, so please let me know if I'm off base with anything.

I am relatively new to DSP. I have done the Max MSP / PureData thing for a couple years, and finally breaking out into other things as I'm up against the limitations of those frameworks. I am trying to make a MiniMoog Model D, and wanted to write my own implementation of the ladder filter. I was heavily inspired by Urs Heckmann's blog posts from 2016 (One-pole unlimited / monster). This really got the wheels turning for me, and I like the simplicity of the "s" terms as an accumulator for the capacitors. My degree is in Electrical Engineering, so I'm well aquatinted with the idea of the capacitor state as a function of accumulated charge. While this model may be slightly reductive, the advantage seems to be how simple the control is, particularly with setting the cutoff frequency through the "g" term. The blog articles only allude to the full implementation, and since I opted for a non-linear model with global feedback (positive feedback overdrive, not just the negative feedback resonance), I spent quite a bit of time figuring out how to set up a multi-dimensional NR solver. I was documenting this as a went, and realized I basically had half a "research paper", and that it wouldn't be too much more work to format it into something presentable. I also figured I may as well publish it somewhere in hopes that it would assist someone trying to do something similar, and answer some of the questions that I had throughout this process.

For anyone interested, I have sort of a draft copy of the paper hosted here. If you are previewing directly from GitHub, you might have to hit the "more pages" button a couples times - its 13 pages long. I'd welcome any feedback, if anyone notices anything that is not quite right. I also wanted to give Urs Heckmann and Andy Simper a heads up about this, as I cite them and their work several times. I don't have a way to reach either of them, but I know they are quite active here. Thank you both for your work, your papers and recorded talks were an invaluable resource for me!

I am also wondering if anyone has any ideas of the best place to "publish" this? My personal "portfolio" website is on the fritz, and will need to be redone before I host anything there, so I was thinking something along the likes of ResearchGate, ArXiv.org, or Academia.edu. Most of the DSP papers I've read are hosted by specific institutions like CCRMA, DAFX, and ADC, which I currently have nothing to do with (apart from enjoying their content/publications). Any advice would be appreciated!

My next steps are that I am currently working on a C++ implementation of this model (Max external first, then maybe JUCE).

[mystran]

I'm not quite sure what makes this "OTA based" model if you're still using non-linearities of the ladder. Are you talking about the idea that you can convert the differential signal path into a normal one? That's always been done in most models, since working differentially only makes sense if you bother modelling each pn-junction separately (which should still be workable in realtime, but a lot more expensive obviously).

You can approximate tanh() efficiently in various ways. Probably the fan favorite is to approximate sinh() by some low order polynomial (cubic if you're really cheap, but something like x+x^3/6+x^5/120 is quite workable) and then compute tanh=sinh/sqrt(1+sinh^2) using said sinh approximation. One can even use hardware rsqrt to compute 1/sqrt if speed is more important than accuracy... but one thing to watch out is that if one approximates tanh() then Newton should really use the proper derivative of the approximation rather an approximation (perhaps derived) of the derivative of the original function, otherwise you'll hurt convergence.... which can be an issue in general when it comes to sigmoids.

With these filters, a lot of the performance comes down to how fast you can solve Ax=b (and also how to tweak Newton to converge in the least number of iterations, but that's somewhat more "black magic"). In most cases LU is your method of choice, but don't assume the first pivoting order you'll pick is the fastest. Since you do definitely have to unroll the LU into straight-line code and then optimize that, I'd highly suggest writing a code-generator or you'll quickly go nuts.

edit: The problem of finding the optimal pivots that minimize LU fill-in (which is basically what "fastest possible order" comes down to) is NP-complete problem (and if you want to further optimize for minimum arithmetics taking into account the fact that not all fill-in is created equal, that definitely does not help) ... but as with many NP-complete problems some heuristics (or just a couple of rounds of brute-force trial and error) can still be worthwhile to get a "reasonably good" solver rather than something approaching the worst-case.

[pscorbett]

Wow, great comment! And thank you for taking the time to read the paper! Admittedly, me calling it OTA based is a bit lazy, but I'm an electrical engineer first and foremost, and I was using the schematic symbol for an OTA, and I think my brain would have broken if I'd tried to call it anything else. For what its worth, I plan to do the Juno filter next, which is actually OTAs

I've seen the "x+x^3/6+x^5/120" tanh() approximation before, and I haven't tried it yet. I've always heard tanh() was taxing, but also assumed that the Cmath library was using a LUT or something. Is it actually more expensive that this method, with three exponentials and a square root?

The thought hadn't even crossed my mind to roll my own Ax=b solver, or code generator. I am not done my C++ implementation yet, but I was planning on using a standard linear algebra library like Armadillo or Eigen. I have no idea the ramifications of using them in a real-time application though. Is this a bad idea?

[Z1202]

I have an impression that "OTA" is often used to label the type of nonlinearity occurring in the discrete-component model of the SSM-based Prophet 5 filter, which is still a ladder, but the nonlinearity is qualitatively different. No idea, how correct both are, but the model is effectively using a differential transistor pair to drive the capacitor loading (compared to the Moog ladder, where the capacitor is driven by the difference of two differential transistor pairs). I believe I used the term OTA when talking about this model earlier, possibly some other people picked this up, or used the same term independently, or it could mean smth else.
NB. didn't read the OP paper (yet).

Edit: just checked my book, I was also using the term OTA in that sense in the book, possibly the term was picked from there

[pscorbett]

The term OTA stands for Operational Transconductance Amplifier, which is basically just differential amplifier with an additional input that controls the output current. This makes them very useful for VCAs and filters. The commonly used OTAs in audio tend to be LM17300s and CA3080s, but both the SSM2040 and IR3109 "filter" chips appear to be basically and OTAs followed by buffers, in the configuration that I borrowed from Urs and Andy in my paper (you have to "bring your own capacitors"). The Roland Juno VCF "chips" such as the 80017A were just a daughter board with a similar OTA/buffer structure.

A transistor ladder is definitely NOT an OTA of course, that's why the title of the paper says it's an "OTA Abstraction". Urs switched his nomenclature from calling it an OTA to a Voltage Controlled Current Source (VCCS), but I like just referring to it as an OTA* (caveat).

The one pole formulas laid out in the paper (and in Urs's blog) are as follows:
I_linear = g*(V+ - V-)
I_OTA = g*tanh(V+ - V-)
I_ladder = g*(tanh(V+) - tanh(V-))

So the formula of the nonlinearity is reflective of the correct type of filter, but the structure is abstracted into something easier to analyze and manipulate

Sorry, Z1202, I don't know from your user name (new here)... which book? Are you Vadim by chance??
(If so, I cited the Art of VA filter Design as well )

[Z1202]

The one pole formulas laid out in the paper (and in Urs's blog) are as follows:
I_linear = g*(V+ - V-)
I_OTA = g*tanh(V+ - V-)
I_ladder = g*(tanh(V+) - tanh(V-))

That's what I was referring to, in the formulas above I_OTA (approximated as a single differential transistor pair, with a buffered output - which I guess is probably correct at least for the first OTAs) has only one tanh, while the transistor ladder has two, corresponding to the two differential transistor pairs before and after the capacitor.

Sorry, Z1202, I don't know from your user name (new here)... which book? Are you Vadim by chance??
(If so, I cited the Art of VA filter Design as well )

Happy that the book continues to serve the community

[pscorbett]

Very cool! It's been a fantastic resource for me over the last few months. I keep referring back to it Its interesting to hear hear some context on it. Was it originally you that devised this structure of the one-pole "OTA block" for all three filter types? I can't say that I traced the full history back in terms of dates. I assumed that Urs borrowed if from Andy Simper, and I know you briefly mention OTAs in the start of the "Ladder Filter" chapter. Section 6.10 was very helpful when I was trying to figure out the multidimensional NR Thank you for your contributions!!

[mystran]

I've seen the "x+x^3/6+x^5/120" tanh() approximation before, and I haven't tried it yet. I've always heard tanh() was taxing, but also assumed that the Cmath library was using a LUT or something. Is it actually more expensive that this method, with three exponentials and a square root?

Oh, there's no exponentials here. As a rule of thumb, square root on a modern desktop CPU costs about the same as division. Since we're doing both here, we might be able to save a few cycles by using rsqrt (=lower precision 1/sqrt) and a Newton round or two to refine, but even with div+sqrt it's not too bad.

The suggested polynomial "sinh" is x+(x^3)/6+(x^5)/120 (or similar). That's a totally terrible approximation for sinh as such (with sinh growing exponentially, really any polynomial is going to be terrible), but it doesn't matter too much for tanh, because sinh/cosh approaches 1/1 quick enough that we really just need to get the shape of the knee right. So code would look something like (with divisions written in such a way that the compiler can optimize to reprocidural multiplication even without fast-math):

Code: Select all

float x2 = x*x;
float sh = x*(1+x2*((1/6.f) + x2*(1/120.f)));
float th = sh / sqrt(1+sh*sh);

[Z1202]

Was it originally you that devised this structure of the one-pole "OTA block" for all three filter types?

No, that definitely wasn't me, I prefer thinking in terms of continuous-time block diagrams. I can assume that someone who prefers to think in terms of circuit diagrams came up with the idea to represent an inverting summator + tanh + integrator as an OTA + capacitor + buffer instead.

@mystran: I think there is a +/* typo in your sh formula.

[mystran]

@mystran: I think there is a +/* typo in your sh formula.

Thanks.. fixed.. I actually thought I'd fixed that already (after proofreading), but must have cancelled edit or something.

[mystran]

No, that definitely wasn't me, I prefer thinking in terms of continuous-time block diagrams. I can assume that someone who prefers to think in terms of circuit diagrams came up with the idea to represent an inverting summator + tanh + integrator as an OTA + capacitor + buffer instead.

The only place I remember seeing this is Urs' article. For what it's worth, I prefer to think in terms of differential equations, so definitely wasn't me either.

edit: Now that I think about it... I might have used such equivalence a long time ago... but I don't remember writing about it publicly so I'd still imagine this is either from Urs or from someone else.. idk.. in any case, it is arguably kinda easier to just work with some mathematical formalism (eg. integrators) directly, whether in terms of equations or in terms of block diagrams (which can then be converted to equations in a straight-forward fashion)...

[Z1202]

edit: Now that I think about it... I might have used such equivalence a long time ago... but I don't remember writing about it publicly so I'd still imagine this is either from Urs or from someone else.. idk.. in any case, it is arguably kinda easier to just work with some mathematical formalism (eg. integrators) directly, whether in terms of equations or in terms of block diagrams (which can then be converted to equations in a straight-forward fashion)...

I can imagine that it can be instructive to compare such circuit diagrams to the circuit diagrams of the real filters. But beyond that comparison, I also agree that it doesn't seem to give any special insights and rather unnecessarily obfuscates things. YMMV.

[pscorbett]

That's a good! I did notice the book explicitly stated integrators in the system diagrams. My gripe with the structure I used was that the development of the "s" term ("state accumulator"?) was a bit handwavey. I'm used to s being used as a Laplace variable, and so it wasn't such a stretch to consider it as such a charge accumulator or integrator here (although in the Laplace domain the integrator would be 1/s). The nice thing about it was that it could be completely isolated from the non-linearity, which is why I stuck with it.

I also didn't derive the update step for "s" myself: s[n+1] = 2y[n] - s[n]
It looks like it is just the result of one of the Euler's methods, although I could be wrong. To me, its very interesting to see the different approaches people take, and preferences, for developing their filter models!

@mystran Thanks for the explanation of this polynomial approximation! I will have to look into incorporating this. Right now, I am just getting my head around CMake hell and trying to get my compilers to work with the SDK and libraries... there will be further optimization opportunities once that is figured out

[Z1202]

s as a state has little to do with Laplace transform, it's a very different variable. The integrators in the book are almost exclusively trapezoidal, not Euler.

[mystran]

That's a good! I did notice the book explicitly stated integrators in the system diagrams. My gripe with the structure I used was that the development of the "s" term ("state accumulator"?) was a bit handwavey. I'm used to s being used as a Laplace variable, and so it wasn't such a stretch to consider it as such a charge accumulator or integrator here (although in the Laplace domain the integrator would be 1/s). The nice thing about it was that it could be completely isolated from the non-linearity, which is why I stuck with it.

Not sure which book you're referring to (Vadim's?), but usually we do use "s" for the Laplace variable and "1/s" for an integrator. When discretizing continuous-time systems (most commonly with trapezoidal) some people then use unrelated (well, sort of) state variables such as s1, s2 for the integrators, though these days I personally tend to follow the convention of calling them z1, z2 instead (but really that's totally arbitrary). These variables are the "implementation details" for the "abstract" integrators.

I also didn't derive the update step for "s" myself: s[n+1] = 2y[n] - s[n]
It looks like it is just the result of one of the Euler's methods, although I could be wrong.

It depends on the surrounding context. Often this kind of operation (specifically with the "2" term) is found in code for trapezoidal when using the common "trick" of implementing the integrator in "TDF2" form (valid in fixed sample-rate case), where the state variable actually stores the sum of the previous input AND the contribution of the previous input to the next time-step (which allows for just one state variable per integrator).