Cheap non-linear zero-delay filters

[Richard_Synapse]

I can back Mystran here, The CA3080 class of devices (eg CA3280 / LM13600 / LM13700) is most certainly very close to a Tanh funciton, to within -90 dB or so. I have not analysed the CEM3320, but if anyone wants to point me towards a schematic I'll be happy to

After doing more experiments with a PSPICE model (allowing for arbitrarily high sample rates), I stand corrected- the tanh model gives a small improvement compared to the linearized version. The error signal doesn't disappear though, so either the model needs a correction step or there's some second order effect not modelled which is relevant.

Richard

[Ichad.c]

I can back Mystran here, The CA3080 class of devices (eg CA3280 / LM13600 / LM13700) is most certainly very close to a Tanh funciton, to within -90 dB or so. I have not analysed the CEM3320, but if anyone wants to point me towards a schematic I'll be happy to

After doing more experiments with a PSPICE model (allowing for arbitrarily high sample rates), I stand corrected- the tanh model gives a small improvement compared to the linearized version. The error signal doesn't disappear though, so either the model needs a correction step or there's some second order effect not modelled which is relevant.

Richard

Hmm, just a total shot in the dark by an inept programmer - so take it with a pinch of salt. I vaguely remember reading somewhere that some chips internally also use feedback to "linearize" the circuit, if this is true - neotec's recently proposed 1pole+feedback integrators could be used to simulate it - instead of the usual intergrator - could also factor in the gain increase.

Might be wrong and or useless, just a theory

[andy-cytomic]

After doing more experiments with a PSPICE model (allowing for arbitrarily high sample rates), I stand corrected- the tanh model gives a small improvement compared to the linearized version. The error signal doesn't disappear though, so either the model needs a correction step or there's some second order effect not modelled which is relevant.

Richard

Which filter are you modelling? Can you post a picture of the PSPICE schematic?

[Richard_Synapse]

After doing more experiments with a PSPICE model (allowing for arbitrarily high sample rates), I stand corrected- the tanh model gives a small improvement compared to the linearized version. The error signal doesn't disappear though, so either the model needs a correction step or there's some second order effect not modelled which is relevant.

Richard

Which filter are you modelling? Can you post a picture of the PSPICE schematic?

It's the SH2 filter, a discrete OTA ladder. For the SPICE simulation I just used four LM13700 in series without a resonant path (since atm I'm interested only in the stage nonlinearities). The model posted here earlier was

dV0/dt = f * tanh( 1.1 * in - r' * tanh(1.96 * V3) - v0)
dV1/dt = f * tanh( 1.1 * v0 - v1 )
dV2/dt = f * tanh( 1.1 * v1 - v2 )
dV3/dt = f * tanh( 1.1 * v2 - v3 )

Meanwhile I performed some more tests, applying a "proper" solver to the OTA equations above. They match to -80 db or so, hence the solver as such does not seem to be at fault.

Richard

[andy-cytomic]

Mystran, I think there is another correction for your first post. It looks like you used a pure integrator to update the "s" variable instead of a one pole low pass:

...Now the feedback dependence is linear, so we can implement this as::

Code: Select all

 t = T(0.5*(x[n] + x[n-1]) - s[n]) 
 y[n+1] = (s[n] + f*t*x[n]) / (1 + f*t) 
 s[n+1] = s[n] + 2*f*t*x[n] 

Note that technically it this only reasonable when the signal changes slowly....

correcting both the use of 'n' as previously discussed and using a low pass yields:

Code: Select all

 t = T(0.5*(x[n] + x[n-1]) - s[n-1]) 
 y[n] = (s[n-1] + f*t*x[n]) / (1 + f*t) 
 s[n] = s[n-1] + 2*f*t*(x[n] - y[n])

A quick note also on using 0.5*(x[n] + x[n-1]). If you plot the frequency response of s[n]/x[n] you will see a large high frequency boost as f gets higher. For a frequency response function like (x[n]-s[n])/x[n] where you are actually shaping the high pass output like the above you can normalise the amplitude of this boost at nyquist by using a factor like this:

Code: Select all

 m = min (1-f/2, 1/2) 
 t = T(m*x[n] + (1-m)*x[n-1] - s[n-1])

Here are some plots to show what I mean, the blue is the ideal response (x[n] - y[n])/x[n], green the (0.5*(x[n] + x[n-1]) - s[n-1])/x[n], the mustard/yellow one is (x[n] - s[n-1])/x[n], and the red is the corrected one (m*x[n] + (1-m)*x[n-1] - s[n-1])/x[n] and on the bottom plot please note the green and red are on top of each other:




When the input to the non-linearities is only an s[n-1] term on its own then there isn't much you can do, but from my tests it's seems to give slightly better results putting up with this boost in frequency to get a 1/2 sample delay than using the y[n-1] value with the correct frequency response and being a whole sample behind, but it probably depends on the particular case. As always with DSP you never get something for nothing, there are always tradeoffs.

[andy-cytomic]

For completeness here are the low pass and band bass responses as well. The low pass is s[n-1]/x[n] and the band pass (s1[n-1]-s2[n-1])/x[n], and the cutoff frequencies are 2^-4, 2^-2.5 and 2^-1.5, which at a sample rate of 44.1khz come out to around 2.7k, 7.8k, and 15.6k

[andy-cytomic]

....
4. Fire up your favourite Computer Algebra System, and solve the whole thing. For Maxima as an example:

Code: Select all

  declare([s0,s1,s3,s4,in], mainvar);
  solve([
     y0 = (s0 + f * t0 * (in - r*y3)) * g0,
     y1 = (s1 + f * t1 * y0) * g1,
     y2 = (s2 + f * t2 * y1) * g2,
     y3 = (s3 + f * t3 * y2) * g3
    ], [y0, y1, y2, y3]);

Take one of the solutions (I'm using y3) and write it down in code. In most cases it doesn't really matter which one you pick, so whichever looks the least evil is usually a good candidate, as long as it appears in the right-hand side of at least one other equation, so we can hope to simplify the remaining set of equations. Naturally you can pull out more than one if they are easy to calculate. Whatever seems to result in the smallest amount of total work.

Hey Mystran, since you already have everything grouped with as few terms as possible you don't need to involve a maths package and then try and undo the mess it delivers and pull back out the common terms manually, you can substitute directly as follows, which is manual but not too error prone since since it's mostly copy and paste. You could even get the maths package to do these steps for you:

Code: Select all

y3 = (s3 + f*t3*y2)*g3

now bring the g3 out the front as it's a little neater:
y3 = g3*(s3 + f*t3*y2)

insert what y2 is directly in the above:
y3 = g3*(s3 + f*t3*g2*(s2 + f*t2*y1))

then insert y1:
y3 = g3*(s3 + f*t3*g2*(s2 + f*t2*g1*(s1 + f*t1*y0)))

then insert y0:
y3 = g3*(s3 + f*t3*g2*(s2 + f*t2*g1*(s1 + f*t1*g0*(s0 + f*t0*(in - r*y3)))))

now you have y3 on the right hand side and you can eliminate it in the usual way, just delete all the 's' terms and only keep the scaling terms of y3 on the right hand side:
y3 = g3*(s3 + f*t3*g2*(s2 + f*t2*g1*(s1 + f*t1*g0*(s0 + f*t0*in)))) / 
(1 + g3*(f*t3*g2*(f*t2*g1*(f*t1*g0*(f*t0*r)))))

now mystran did a slightly different grouping but the one I'm showing makes more sense due to the placement of brackets by the result of the direct substitution above:

f3 = f*t3*g2*g3;
f2 = f*t2*g1*f3;
f1 = f*t1*g0*f2;
f0 = f*t0*f1;

y3 = (g3*s3 + f3*s2 + f2*s1 + f1*s0 + f0*in) / (1 + r*f0)

here is the solution mystran posted previous which gives the same result:

f3 = f*t3*g3
f2 = f*t2*g2*f3
f1 = f*t1*g1*f2
f0 = f*t0*g0*f1

y3 = (g3*s3 + f3*g2*s2 + f2*g1*s1 + f1*g0*s0 + f0*in) / (1 + r*f0) 

[andy-cytomic]

From what I can tell from a real OTA ladder filter, the nonlinearity is not tanh. For the diode clipper it's most likely the wrong model to choose as well. Of course it may sound perfectly fine, just saying unlikely to be 100% authentic if that's what you're after.

Like I tried to point out (perhaps not explicitly enough) it depends on the OTA.

For something like CA3080 (or even LM13700 as long as you leave the diode linearization unconnected) tanh() is quite reasonably model, since the whole this is just another long-tailed pair plus a few current mirrors. So if you built an OTA ladder with 4x CA3080 for the stages and another for the resonance control, then you should get roughtly tanh() all the way.

Ofcourse, most OTA cascades in the wild are not built out of CA3080s, but something like CEM3320 so the question of "what are the correct non-linearities" then becomes "how are the CEM3320 gaincells implemented".

I can back Mystran here, The CA3080 class of devices (eg CA3280 / LM13600 / LM13700) is most certainly very close to a Tanh funciton, to within -90 dB or so. I have not analysed the CEM3320, but if anyone wants to point me towards a schematic I'll be happy to

I've now done further investigations into the shape of the CA3080 and it's actually not quite tanh as I previously posted. Antti put me on to this great (and free download as a pdf) book which goes into more detail on how much slight variations in manufacturing effect performance: http://www.designinganalogchips.com/ I thought the tolerances in chips were better than they are. Even with a four transistor current mirror you still have problems in the real word due to very slight mismatching, so depending on the chip it's within only about -50 dB or so of a tanh.

[andy-cytomic]

In emails with Mystran he has said that his method of using f(x)/x for non-linearities is already contained in the source code of certain circuit simulation packages, possibly as fallback methods for stability. Maths is a bit tricky like that, there are lots of really smart people that have already covered most of this ground long ago. I think i've only ever come up with one cool idea that I haven't seem yet in all my time doing DSP, and this a smoothing algorithm, most of the other stuff I do is just fast approximations - and the maths behind those is already discovered, it's more the choice it picking which approximation works best for which case. These choices matter and contribute to the tone of your work.

So in summary, everything presented in this thread is old news. Circuit simulation packages have been doing this stuff better than us since computers were first invented, and most of the maths was sorted out long before computers even existed.

Now there are still some people doing odd things I've noticed, mostly due to the fact they are uncomfortable with currents. The integration of the sallen key filter is the most striking example where people try and remove integrating the current over the capacitor by solving for a different set of equations and so introduce many more terms. When I get a little more time I'll write a paper for everyone that runs through exactly how to apply arbitrary numerical integration schemes to basic (filter) circuits. Hopefully such a summary will be easy to follow and show the equivalence between standard trapezoidal integration of circuits and the so called "TPT", and so dispel that this is anything new as either.

[mystran]

In emails with Mystran he has said that his method of using f(x)/x for non-linearities is already contained in the source code of certain circuit simulation packages, possibly as fallback methods for stability.

Well I haven't really found any source code utilizing the code, but in a book (and I regret forgetting which one it was, couldn't find it later) about circuit simulation methods there was a comparison between various iterative methods, including fixed-point, secant, this method and Newton. The author referred to the method mentioned here as the "BLAS-3" method and basically said that it's not been studied too carefully, but suspected the convergence (with iteration) should be somewhere between linear and that of secant method (which is about sqrt(2)).

Personally, from empirical observations, I suspect the convergence is damn near linear. As far as I'm concerned, the only advantage over Newton (or secant) is that a truncated iteration (most notably single evaluation only) tends to result in something that isn't totally useless (as could happen with Newton).

[mystran]

....
4. Fire up your favourite Computer Algebra System, and solve the whole thing. For Maxima as an example:

Code: Select all

  declare([s0,s1,s3,s4,in], mainvar);
  solve([
     y0 = (s0 + f * t0 * (in - r*y3)) * g0,
     y1 = (s1 + f * t1 * y0) * g1,
     y2 = (s2 + f * t2 * y1) * g2,
     y3 = (s3 + f * t3 * y2) * g3
    ], [y0, y1, y2, y3]);

Take one of the solutions (I'm using y3) and write it down in code. In most cases it doesn't really matter which one you pick, so whichever looks the least evil is usually a good candidate, as long as it appears in the right-hand side of at least one other equation, so we can hope to simplify the remaining set of equations. Naturally you can pull out more than one if they are easy to calculate. Whatever seems to result in the smallest amount of total work.

Hey Mystran, since you already have everything grouped with as few terms as possible you don't need to involve a maths package and then try and undo the mess it delivers and pull back out the common terms manually, you can substitute directly as follows, which is manual but not too error prone since since it's mostly copy and paste. You could even get the maths package to do these steps for you:

Yeah.. you're right. For whatever reason I didn't realize at a time that there's a much better way to let the math package (say Maxima) do the work for you:

Code: Select all

get_lu_factors(lu_factor(matrix(
[ 1+f*t1,      0,      0, t0*f*r ],
[  -f*t1, 1+f*t2,      0,      0 ],
[   0,     -f*t2, 1+f*t3,      0 ],
[   0,         0,  -f*t3, 1-f*t4 ])));

I'll save you the output since it doesn't format nicely for the forum (it does in the wxMaxima window though), but it's rather trivial to translate to an efficient implementation (lu_factor on "generalring" doesn't do any simplification so the structure of the reduction remains visible).

That said, in this particular case it's probably a bit of an overkill.

Anyway, should probably update the original algorithm for a more efficient reduction, in case someone still cares..

[Z1202]

Hopefully such a summary will be easy to follow and show the equivalence between standard trapezoidal integration of circuits and the so called "TPT", and so dispel that this is anything new as either.

Actually, in the VAFilterDesign book the bilinear integrators are introduced via the trapezoidal integration. It's absolutely clear that the two methods are absolutely equivalent, except for the "instantaneously unstable" cases.

The clear benefit of the block-diagram approach (for me) is that, being visual, it's more intuitive and thus sometimes can help avoiding certain kinds of math errors (YMMV).

Edit: one example showing the intutivity, is the replacement of the DF1 integration with TDF2 integration, thus halving the state variable count. Of course, this is also possible analytically (and IIRC I'm doing this in the book), but the idea is way more obvious in the graphical form.

PS. It continues to amaze me, how much emotional effort you are putting into trying to prove I didn't come up with anything new. If I ever claimed to have created anything new with the TPT, it was the methodology not the results.

[andy-cytomic]

In emails with Mystran he has said that his method of using f(x)/x for non-linearities is already contained in the source code of certain circuit simulation packages, possibly as fallback methods for stability.

Well I haven't really found any source code utilizing the code, but in a book (and I regret forgetting which one it was, couldn't find it later) about circuit simulation methods there was a comparison between various iterative methods,

Hey Mystran, sorry for mis-paraphrasing you. It was a while since you sent me that email, and what you actually said was pretty vague not mentioning source code or even a manual just "found some indications that apparently a pre-historic (is pre-SPICE) circuit simulator called BIAS...", so sorry for that and thanks for clearing it up!

[mystran]

...

Block diagrams are a tool of the devil. Get out of this thread, you heretic!

The serious part then is that the only way I personally can read a block-diagram is to convert it into equations (or something similar), and it usually takes several attempts before I'm confident that I actually got it right (which doesn't necessarily mean I actually got it right). I'd rather describe them as frightening rather than intuitive. But I also can't do geometry without converting it to algebra first, so maybe I'm just stupid. YMMV.

[mystran]

In emails with Mystran he has said that his method of using f(x)/x for non-linearities is already contained in the source code of certain circuit simulation packages, possibly as fallback methods for stability.

Well I haven't really found any source code utilizing the code, but in a book (and I regret forgetting which one it was, couldn't find it later) about circuit simulation methods there was a comparison between various iterative methods,

Hey Mystran, sorry for mis-paraphrasing you. It was a while since you sent me that email, and what you actually said was pretty vague not mentioning source code or even a manual just "found some indications that apparently a pre-historic (is pre-SPICE) circuit simulator called BIAS...", so sorry for that and thanks for clearing it up!

No big deal. Was it BIAS? See, I've even forgot the name. Anyway, I was intentionally vague back then, because I'd just lost track of the source (not source code, but the source of information) and I kinda hate it when I can't give a proper reference.