# Cheap non-linear zero-delay filters

DSP, Plug-in and Host development discussion.
mystran
KVRAF
6118 posts since 12 Feb, 2006 from Helsinki, Finland
Ok, so I didn't want to pollute that other thread, so I'm starting a new one.

In that other thread I mentioned you can linearize a filter around it's known state and get reasonable results for not much overhead. Here's some discussion of the topic. [edit: also, while I said "zero-delay" in the thread title, I'm really doing "topology-preservation" too; you don't need to, if you don't want to, but it kinda helps if you want to mimic analog non-linearities]

I'm not going to concern myself with BLT other than note that since this is based on trapezoidal integration, and since I'm using TDF2 BLT integrators as building blocks, we can tune filters exactly in cases where the non-linearities are negligible (which is to say that I'm not going to discuss what happens when the non-linearities force the cutoff to change and we're in a seriously warped region of the BLT frequency scale). Namely, the form of integrator building block I'm using is:

Code: Select all

``````  y[n+1] = s[n] + f*x[n]
s[n+1] = s[n] + 2*f*x[n]

where
x[n], y[n] and s[n] are input, output and state

and the tuning that maps the normalized analog frequency 1
to whatever frequency we want is
f = tan(M_PI * freq / samplerate)
``````
Few useful observations: if we think of y[n] and s[n] as continuous piece-wise linear functions y(n) and s(n), then y(n)=s(n-0.5). So the output is essentially a half-sample delayed version of the state. We're going to abuse this.

To keep things simple and the math short, I'm going to use a one-pole OTA lowpass (straight from the LM13700 datasheet) as and example and ignore all the practical details like signal levels (eg we'd really need need a ton of 1/vT factors) and auxiliary circuitry. So we have the following model for the OTA (the derivation or accuracy of which is not important) and capacitor:

Code: Select all

`````` iOta = iCtl * tanh(vIn)
dVcap/dt = 1/C * iOta
``````
And with feedback (since we want lowpass and not an integrator):

Code: Select all

`````` dVcap/dt = iCtl/C * tanh( vIn - vCap )
``````
So if we discretize, we get:

Code: Select all

`````` y[n+1] = s[n] + f * tanh( x[n] - y[n+1] )
``````
This looks like it needs a feedback solver, but here's the thing: Let's rewrite the OTA in terms of transconductance gM = iOut/vIn and we get:

Code: Select all

`````` gM(vIn) = iCtl * ( tanh( vIn ) / vIn )
or:
iOut = iCtl * T(vIn) * vIn
where:
T(x) = tanh(x) / x
``````
[footnote: tanh(x)/x=1 in the limit x -> 0, but it's easy to adapt most tanh(x) approximations to return tanh(x)/x instead and avoid this issue]
[footnote: you can do this for practically any type of memoryless non-linearity, even when "transconductance" as such wouldn't be meaningful]

Code: Select all

`````` y[n+1] = s[n] + f * T(x[n] - y[n+1]) * (x[n] - y[n+1])
``````
The point of this exercise is that we can now treat the non-linear transconductance and the actual input separately. So what we can do, is combine Euler method for the non-linearity with trapezoidal method for the linear part! In alternative interpretation we delay the transconductance by half a sample. Recall that s[n] = y[n+0.5]. For consistency, use x[n-0.5]=0.5*(x[n]+x[n-1]) for the actual input signal; everything else is available from one of the filter states:

Code: Select all

`````` y[n+1] = s[n] + f * T(x[n-0.5] - s[n]) * (x[n] - y[n+1])
``````
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. That's true at audible frequencies if we oversample to avoid aliasing. With higher gain or nastier non-linearities it deviates more, but so does aliasing increase, and once you oversample you increase both again. In practice the above works remarkably well for almost anything I've thrown at it so far (transistor ladders, diode ladders, OTA cascades... you name it).

On paper it's less accurate than fitting a linear curve directly to the tangent of the tanh() for slowly changing signals, since we are fitting a linear curve from the origin to the known operating point, but unlike the tangent fitting method, this tolerates violations of the "signal changes slowly" assumption much better; we might feed a bit too much or two little current, but most of the time the results are relatively sane (which cannot be said about tangent fitting, which can run into crazy paradoxes). You can certainly use this directly and in most cases with a bit of oversampling (eg 4 times or so usually work sensibly for reasonable input levels) it sounds quite fine (and when I say "quite fine" I mean "certainly a lot better than a traditional sequential fudge-factored filter").

Anyway, if you're NOT happy with the results (remember we're only first order as far as the non-linearities go), we can treat the calculated value as a prediction, and apply a correction step. Sensible approach would be a variation of Heun's method take the new state (and x[n+0.5]; you need one step lookahead) and recalculate the "transconductances" then redo the linear solver, then average the resulting state with the original prediction (and likewise for outputs). Since the error of the "correction" step should be opposite to the error of the "prediction" step, they should mostly cancel. As far as I can tell, this is sufficient to make it a true second-order method (don't feel like doing formal error analysis, sorry).

In practice the correction step roughly doubles the CPU use. In cases where the prediction works well (and most of the time it does), it's probably better idea to double the oversampling instead, but if you hit an obviously degenerate case, then the above could solve the issue.
Last edited by mystran on Tue Jan 08, 2013 12:10 pm, edited 4 times in total.
Preferred pronouns would be "it/it" because according to this country, I'm a piece of human trash.

mystran
KVRAF
6118 posts since 12 Feb, 2006 from Helsinki, Finland
For the TL;DR folks, here's transistor ladder using the above (without correction) to play with.. note that this assume straight and linear feedback loop which is not really a very good model in practice Code: Select all

``````//// LICENSE TERMS: Copyright 2012 Teemu Voipio
//
// You can use this however you like for pretty much any purpose,
// as long as you don't claim you wrote it. There is no warranty.
//
// Distribution of substantial portions of this code in source form
// must include this copyright notice and list of conditions.
//

// input delay and state for member variables
double zi;
double s = { 0, 0, 0, 0 };

// tanh(x)/x approximation, flatline at very high inputs
// so might not be safe for very large feedback gains
// [limit is 1/15 so very large means ~15 or +23dB]
double tanhXdX(double x)
{
double a = x*x;
// IIRC I got this as Pade-approx for tanh(sqrt(x))/sqrt(x)
return ((a + 105)*a + 945) / ((15*a + 420)*a + 945);
}

// cutoff as normalized frequency (eg 0.5 = Nyquist)
// resonance from 0 to 1, self-oscillates at settings over 0.9
double cutoff, double resonance,
double * in, double * out, unsigned nsamples)
{
// tuning and feedback
double f = tan(M_PI * cutoff);
double r = (40.0/9.0) * resonance;

for(unsigned n = 0; n < nsamples; ++n)
{
// input with half delay, for non-linearities
double ih = 0.5 * (in[n] + zi); zi = in[n];

// evaluate the non-linear gains
double t0 = tanhXdX(ih - r * s);
double t1 = tanhXdX(s);
double t2 = tanhXdX(s);
double t3 = tanhXdX(s);
double t4 = tanhXdX(s);

// g# the denominators for solutions of individual stages
double g0 = 1 / (1 + f*t1), g1 = 1 / (1 + f*t2);
double g2 = 1 / (1 + f*t3), g3 = 1 / (1 + f*t4);

// f# are just factored out of the feedback solution
double f3 = f*t3*g3, f2 = f*t2*g2*f3, f1 = f*t1*g1*f2, f0 = f*t0*g0*f1;

// solve feedback
double y3 = (g3*s + f3*g2*s + f2*g1*s + f1*g0*s + f0*in[n]) / (1 + r*f0);

// then solve the remaining outputs (with the non-linear gains here)
double xx = t0*(in[n] - r*y3);
double y0 = t1*g0*(s + f*xx);
double y1 = t2*g1*(s + f*y0);
double y2 = t3*g2*(s + f*y1);

// update state
s += 2*f * (xx - y0);
s += 2*f * (y0 - y1);
s += 2*f * (y1 - y2);
s += 2*f * (y2 - t4*y3);

out[n] = y3;
}
}
``````
[/size]
Last edited by mystran on Sat Dec 13, 2014 11:34 am, edited 2 times in total.
Preferred pronouns would be "it/it" because according to this country, I'm a piece of human trash.

mystran
KVRAF
6118 posts since 12 Feb, 2006 from Helsinki, Finland
Oh and AQ: if you use this you own me a Poly-Ana. Preferred pronouns would be "it/it" because according to this country, I'm a piece of human trash.

mystran
KVRAF
6118 posts since 12 Feb, 2006 from Helsinki, Finland
• if you want to filter stereo signals, you can treat the signals and states as vectors, and evaluate "transconductances" with the Pythagorean norm sqrt(L*L+R*R) which ofcourse is specifically why I was approximating tanh(sqrt(x))/sqrt(x) in the first place... this isn't really circuit modelling anymore (good luck with analog version) but unlike a dual-mono implementation, distortion tends to localize at the input signals rather than the individual speakers.. I personally like the effect and I've used it for traditional filters before I even touched any zero-delay stuff (eg Sweep uses it for example; my dev-builds for new version already use a better "OTA SVF" model with the same approach to stereo)
• if you still feel you need to iterate (which personally think is quite futile if you apply the correction step, since trapezoidal is only second-order anyway), then my educated guess is that iteration on the transconductance is likely to converge faster than iteration on the signal itself... but like I said, I don't see the point without a better integrator
• note that without the correction step, the cost here is almost exactly the cost of solving a linear zero-delay filter + the cost of evaluating the non-linearities once.
Preferred pronouns would be "it/it" because according to this country, I'm a piece of human trash.

mr.bungle
KVRist
95 posts since 25 Jul, 2007 from Finland
I tried the filter in my polysynth. Hot damn. I think it's the best sounding thing in a VST so far.

Here's a build: http://dl.dropbox.com/u/13976481/MauSynth_BETA_RC1.dll KVRian
1081 posts since 8 Feb, 2012 from South - Africa
Thanks for all the detailed explenations! The way in which you make it sound so easy though - is kinda scary My one brain-cell is learning to do simple SSE and the other one is taking a nap Your 4L2 plug is pretty dang useful as well.

Regards
Andrew

camsr
KVRAF
7055 posts since 17 Feb, 2005
Wow this filter tunes spot on. I was able to use full resonance while tracking the fundamental frequency with the MauSynth and it stayed sharp from the bass to the treble. Nice work guys Z1202
KVRian
1271 posts since 12 Apr, 2002
mystran wrote:On paper it's less accurate than fitting a linear curve directly to the tangent of the tanh() for slowly changing signals, since we are fitting a linear curve from the origin to the known operating point, but unlike the tangent fitting method, this tolerates violations of the "signal changes slowly" assumption much better; we might feed a bit too much or two little current, but most of the time the results are relatively sane (which cannot be said about tangent fitting, which can run into crazy paradoxes).
Nice! (well, I didn't verify or try this myself, but from what you're saying seems really nice).

Regards,
{Z}

mystran
KVRAF
6118 posts since 12 Feb, 2006 from Helsinki, Finland

edit: also added a remark in the original post about "topology-preserving" since that's what it is.. when I wrote the thing originally, I didn't realize anyone would want to do "zero-delay" without doing topology-preserving, but apparently that's also possible Preferred pronouns would be "it/it" because according to this country, I'm a piece of human trash.

Big Tick
KVRAF
3340 posts since 29 May, 2001 from New York, NY
Thanks ! There are a few x[x] which should be x[n].

Also, something doesn't make sense to me: if

Code: Select all

``x[n-0.5]=0.5*(x[n]+x[n-1])``
and

Code: Select all

``y[n+1] = s[n] + f * T(x[n-0.5] - s[n]) * (x[n] - y[n+1])``
then shouldn't you have

Code: Select all

`` t = T(0.5*(x[n] + x[n-1]) - s[n] )``

Code: Select all

``t = T(s[n] + 0.5*(x[n] + x[n-1]) ``

mystran
KVRAF
6118 posts since 12 Feb, 2006 from Helsinki, Finland
Big Tick wrote:Thanks ! There are a few x[x] which should be x[n].

Also, something doesn't make sense to me: if

Code: Select all

``x[n-0.5]=0.5*(x[n]+x[n-1])``
and

Code: Select all

``y[n+1] = s[n] + f * T(x[n-0.5] - s[n]) * (x[n] - y[n+1])``
then shouldn't you have

Code: Select all

`` t = T(0.5*(x[n] + x[n-1]) - s[n] )``

Code: Select all

``t = T(s[n] + 0.5*(x[n] + x[n-1]) ``
Yes the sign of s[n] was obviously wrong. Corrected in OP, thanks for noticing.

That's what you get for proof-reading your own text. Regading x[x] I found one and correct it.
Preferred pronouns would be "it/it" because according to this country, I'm a piece of human trash.

Big Tick
KVRAF
3340 posts since 29 May, 2001 from New York, NY
Can you 2x oversample it, and then get the actual half-sample delay from the previous (oversmpled) sample ? Or is this useless ? KVRian
1081 posts since 8 Feb, 2012 from South - Africa
Big Tick wrote:Can you 2x oversample it, and then get the actual half-sample delay from the previous (oversmpled) sample ? Or is this useless ?
Good Q -> sounds sneaky, but wouldn't that be dependant on your oversampling technique?

Random add-on: Since there is the word "cheap" in the heading, you could use 1 of these 2 approximations for tan(x):
A = -15*x+x^3;
B = 3*(-5+2*x^2);
tan_out = A/B;

A = 5*(-21*x+2*x^3);
B = 105-45*x^2+x^4;
tan_out = A/B;
I generally use nr2 - works well - even without oversampling. Haven't checked it out with mystran filter above though. nr1 might track equally as well if you oversample.

Just my 2cents.

Regards
Andrew

mystran
KVRAF
6118 posts since 12 Feb, 2006 from Helsinki, Finland
Big Tick wrote:Can you 2x oversample it, and then get the actual half-sample delay from the previous (oversmpled) sample ? Or is this useless ?
You mean like 2x oversample input, instead of taking average?

I've not tried it but it sounds rather useless. In the worst-case the high-freq attenuation from averaging might reduce aliasing slightly. In the best case (when you're oversampling enough that aliasing isn't an issue) the signal is smooth enough that you could forget the delay and just use the input sample (I've tried both ways and audible difference is generally speaking none).

Mostly I do it for theoretical consistency (ie looks good on paper). It's certainly not a critical component of the method, so I wouldn't waste too many cycles on it (you can always find something else to improve instead).

Oh and remember that aliasing in a non-linear recursive filter isn't really just an audible problem: once you start getting significant aliasing, your filter state essentially becomes bogus, and you enter a garbage-in-garbage-out feedback loop. In my book it doesn't really matter how accurately we recycle garbage (not in DSP anyway).
Preferred pronouns would be "it/it" because according to this country, I'm a piece of human trash.

mystran
KVRAF
6118 posts since 12 Feb, 2006 from Helsinki, Finland
Ichad.c wrote: Random add-on: Since there is the word "cheap" in the heading, you could use 1 of these 2 approximations for tan(x):
A = -15*x+x^3;
B = 3*(-5+2*x^2);
tan_out = A/B;

On top of that, you probably have an exponentiation inside the tangent, so what you really want is tan(pow(2,octaves+log2_basetune)) or something sufficiently similar. If you run this once per sample (and I certainly do), you almost certainly want to optimize the whole thing in one way or another, but that's kinda different topic. edit: I like your approach though.. Pade approximations are generally what I always try first whenever I need to approximate something nasty. edit2: I can't resist commenting on the "different topic" just enough to mention that in general the accuracy requirements for tuning should be optimized in terms of octaves (ie on the logarithmic scale) because the perception of pitch is such that 10 cents deviation at 100Hz sounds roughly as bad as 10 cents deviation at 2kHz and the absolute errors are quite different. 