Cheap non-linear zero-delay filters

[Miles1981]

We're currently wondering what to call the method. One of the guys is finalising a set of papers that explain and compare various filter solving methods. At the moment he calls it "Cheap Method", but I'm leaning towards "Mystran's Method"…

What do you think?

If I had to give it a "proper" name, I'd probably call it "fixed-pivot one-step linearization method" or something similar, since that's what it's about: draw a line through a chosen fixed pivot point (typically origin) and your initial guess (from the previous step), and just use that for solving the system.

I was just thinking about this. It seems to be quite close to a single iteration solve of the non linear equations, doesn't it?

[Miles1981]

So, Guitar Gadgets use mystran's semi-implicit fixed-pivot trapezoidal method to simulate the EMS VCS3 filter, I said it

So I guess it is Filtrator?

[Ivan_C]

Yes it is

[cyberheater]

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 z1;
double s[4] = { 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
void transistorLadder(
    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[3]);
        double t1 = tanhXdX(s[0]);
        double t2 = tanhXdX(s[1]);
        double t3 = tanhXdX(s[2]);
        double t4 = tanhXdX(s[3]);

        // 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[3] + f3*g2*s[2] + f2*g1*s[1] + f1*g0*s[0] + 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[0] + f*xx);
        double y1 = t2*g1*(s[1] + f*y0);
        double y2 = t3*g2*(s[2] + f*y1);

        // update state
        s[0] += 2*f * (xx - y0);
        s[1] += 2*f * (y0 - y1);
        s[2] += 2*f * (y1 - y2);
        s[3] += 2*f * (y2 - t4*y3);

        out[n] = y3;
    }
}

[/size]

Hi.

I have tried to run your code but fails to compile because zi is undeclared. It that a typo. Is zi really meant to be z1?

[mystran]

I have tried to run your code but fails to compile because zi is undeclared. It that a typo. Is zi really meant to be z1?

I think it's supposed to be the other way (with "zi" for "delay of input") but .. yeah. I'll fix it to the first post.

edit: it's also worth noting that the delay/averaging of input is of questionable utility and you could just use in[n] directly (and it might or might not sound better to your ears that way)

[aciddose]

I've created a slightly tweaked version which seems to work better. Things are laid out a little differently such that it is easier to very quickly apply some SSE intrinsics on the coefficients, and if you're willing some shuffle instructions can compute the powers + estimate + state + integrate memory using SSE with very little effort.

edit: to clarify license
This code is AFAIK copyright mystran @ kvraudio.com/forum/, 2012 and I relinquish my associated and derivative rights under Creative Commons (CC0).

Please see the original code and mystran's comments about license terms: viewtopic.php?f=33&t=349859

Code: Select all

float tanhd(const float x, const float d, const float s)
{
	return 1.0f - s * (d + 1.0f) * x*x / (d + x*x);
}

float cascade_4(const float input)
{
	typedef float T;

	const T c = coefficient
	const T fb = from 0 to 4 for 100%, further "drives" feedback

	const T c2 = c * 2.0f;

	// shaper coefficients, offset, scale, shape
	// very quick approximation, close enough for me to tanh
	// yet far more flexible
	const T a = 2.0f;
	const T s = 0.1f;
	const T d = 1.0f;

	// per-sample computation
	//for (i = 0; i < steps_oversamples; i++) {
	const T t0 = tanhd(b[0] + a, d, s);
	const T t1 = tanhd(b[1] + a, d, s);
	const T t2 = tanhd(b[2] + a, d, s);
	const T t3 = tanhd(b[3] + a, d, s);

	const T g0 = 1.0f / (1.0f + c*t0);
	const T g1 = 1.0f / (1.0f + c*t1);
	const T g2 = 1.0f / (1.0f + c*t2);
	const T g3 = 1.0f / (1.0f + c*t3);
	 
	const T z0 = c*t0 / (1.0f + c*t0);
	const T z1 = c*t1 / (1.0f + c*t1);
	const T z2 = c*t2 / (1.0f + c*t2);
	const T z3 = c*t3 / (1.0f + c*t3);

	const T f3 = c       * t2*g3;
	const T f2 = c*c     * t1*g2 * t2*g3;
	const T f1 = c*c*c   * t0*g1 * t1*g2 * t2*g3;
	const T f0 = c*c*c*c *    g0 * t0*g1 * t1*g2 * t2*g3;

	const T estimate =
		     g3 * b[3] +
		f3 * g2 * b[2] + 
		f2 * g1 * b[1] + 
		f1 * g0 * b[0] + 
		f0 * input;

	// feedback gain coefficient, absolutely critical to get this correct
	// i believe in the original this is computed incorrectly?
	const T cgfbr = 1.0f / (1.0f + fb * z0*z1*z2*z3);

	// clamp can be a hard clip, a diode + highpass is better
	// if you implement a highpass do not forget to include it in the computation of the gain coefficients!
	const T xx = input - clamp(fb * estimate) * cgfbr;
	const T y0 = t0 * g0 * (b[0] + c * xx);
	const T y1 = t1 * g1 * (b[1] + c * y0);
	const T y2 = t2 * g2 * (b[2] + c * y1);
	const T y3 = t3 * g3 * (b[3] + c * y2);

	// update the stored state
	b[0] += c2 * (xx - y0);
	b[1] += c2 * (y0 - y1);
	b[2] += c2 * (y1 - y2);
	b[3] += c2 * (y2 - y3);
	//}
	
	// you must limit the compensation if feedback is clamped
	const T compensation = 1.0f + limit(fb, 0.0f, 4.0f);
	return y3 * compensation;
} 

I still prefer the linear version with naive saturation applied directly to the state variables, although it isn't "correct" the timbre is similar and it is a lot cheaper to compute.

This version actually works quite well though and compares well to a real analog circuit. The saturation needs to be tweaked slightly to match a particular circuit, it may still be slightly too high for most real world circuits. (see shaper "scale" coefficient)

[toasttoast]

Very nice, aciddose. Very clean and sweet.

I spent some time with your tanhd() in an online graphing program, and it didn't seem quite right unless I am completely missing the intent. Is it possible that the x*x in the numerator should be x instead? If it is, then I can get close to tanh by tweaking d and s. With the x squared up there, x below zero gives positive y rather than negative y.

Your clamp comment is intriguing. When you say diode + highpass, do you mean going through tanh and then high passing? Why do you suggest that? How would the gain coefficients need to be tweaked as a result?

Thanks for your post.

[aciddose]

tanh doesn't come close to a hard saturating diode. The diode is extremely linear up to a point, then saturates extremely fast.

You can get closer with a cosine (or parabolic) interpolation applied to a hard edge. For example 100% up to .7, then interpolate to 40% at .8.

If you really want to accurately model a diode clamp you either need to use the very expensive recursive/iterative solutions to compute the current/voltage, or you need to pre-compute this and store it in a table. I use the table, it is many times faster than any iterative solution ever could be and nearly as accurate with a simple linear interpolation.

the tanhd function is not supposed to look anything like tanh, it's supposed to look like mystran's tanh(x)/x approximation, only it looks similar while sacrificing as much accuracy as possible to minimize the expense of the computation.

The result should produce exactly the same levels and harmonics as a more accurate tanh(x)/x only with far reduced computational expense. (Do you care that the 5th harmonic is 3db less than it should be? Well I don't.)

The real benefit is that you can easily adjust the amount of saturation "scale" as well as how hard the saturation is "shape". In addition to that I add an offset to get some asymmetry.

The real advantage of doing this, however you do it is that you can precompute all the coefficient values and very quickly process the entire thing with only a few sse instructions and a couple constant values. Nothing stops you from adjusting any of the coefficents per-stage. The mm128 coefficients will store them all and process all four at once regardless.

If you put in a lot more work you can adjust the over-all scale of the transform and maintain tuning as you dynamically adjust the saturation amount, symmetry and hardness. I will work on this at a later date if motivated, don't hold your breath.

[aciddose]

Also note clamp() can just be an empty macro like:
#define clamp(X) (X)

The filter still works just fine, you'll need to adjust levels and treat it more like a minimoog filter vs. a typical rolandesque implementation.

The code I linked posted is absolutely not practical for real-world use in its current state. I've only linked it to act as a building block and tool for experimentation.

If you wanted to really target something like the original minimoog filter you should need to spend weeks bench testing and working with spice while tweaking coefficients and adding lots of highpass filters and additional non-linear sections (the output buffer in feedback path) to get the same spectra and effect that filter produces.

Likewise for any other particular filter. Keep in mind this function doesn't come anywhere even remotely close to a real world filter. It is only a very rough mold you can use to cast a new filter that you carve out the details of yourself.

[toasttoast]

Yes, I appreciate that it's not a finished filter. Just looking for some insight into your thoughts about it. Thanks again.

[aciddose]

If you're curious about the tanh(x)/x, this is because we're computing a scaling coefficient that we multiply with the delta/charge in the integrators or just as well you can also apply it as a clamp function, although it is really not like a diode.

f(x)/x * x = f(x)

So, this is just a quick way to compute the same as doing f(x). In other words clamp(x) = tanh(x) is the same as computing t[-1] = tanh(x)/x, then taking (estimate * t[-1] * feedback). So, exactly like (tanh(estimate) * feedback).

Mystran's example function is one step further toward accuracy, but this assumes you want to use exactly tanh() for some reason, while most likely you do not really want exactly that function but something approximately like it.

BTW: the "offset" I apply to get the asymmetry is something you have to really understand about typical OTA or "ladder" implementations of the circuit. In the minimoog "ladder", mismatches between transistors will make one or the other side of the ladder run at a higher current and this will approximately be like adding the offset to the input as I do, although not exactly. (You'd have to be insane to want an exact result, but it is good to note that we are not aiming for that here.)

In many OTA based designs the signal is actually inverted at each stage and so the mismatch becomes much lower. In fact any offset added initially will be matched by an inverted offset at the next stage, in theory. In addition to that however you'll often see various buffering strategies including single NPNs, darlingtons, more complex buffers, JFETs and even whole opamps.

JFETs are a common choice (see most roland OTA filters) and it is of note that the gate voltage is highly variable between fets, even across the same die.

So taking all these sorts of things into consideration you can make an attempt to more accurately model a particular circuit. The way I've naively added a bias to the input however very roughly covers them all with some minor adjustments.

[Ivan_C]

There is a typo in the first message of the topic I think. The source code Mystran has posted for the one pole OTA is the following :

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]

However, I think he wanted to write that instead :

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]-y[n+1])

Or maybe something like that :

Code: Select all

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

[Ichad.c]

Code: Select all

// initialization (the resonance factor is between 0 and 8 according to the article)
f = tan(Pi * cutoff/sampleRate);
r = (7.f * resonance + 0.5f);

// the input x[n+1] is given by 'in', and x[n] by zi

// input with half delay
float ih = 0.5f * (in + zi);

// evaluate the non-linear factors
float t0 = f*tanhXdX(ih - r * s[3]);
float t1 = f*tanhXdX(s[1]-s[0]);
float t2 = f*tanhXdX(s[2]-s[1]);
float t3 = f*tanhXdX(s[3]-s[2]);
float t4 = f*tanhXdX(s[3]);

// This formula gives the result for y3 thanks to MATLAB
float y3 = (s[2] + s[3] + t2*(s[1] + s[2] + s[3] + t1*(s[0] + s[1] + s[2] + s[3] + t0*zi)) + t1*(2*s[2] + 2*s[3]))*t3 + s[3] + 2*s[3]*t1 + t2*(2*s[3] + 3*s[3]*t1);

y3 /= (t4 + t1*(2*t4 + 4) + t2*(t4 + t1*(t4 + r*t0 + 4) + 3) + 2)*t3 + t4 + t1*(2*t4 + 2) + t2*(2*t4 + t1*(3*t4 + 3) + 2) + 1;

// Other outputs
float y2 = (s[3] - (1+t4+t3)*y3) / (-t3);
float y1 = (s[2] - (1+t3+t2)*y2 + t3*y3) / (-t2);
float y0 = (s[1] - (1+t2+t1)*y1 + t2*y2) / (-t1);
float xx = (zi - r*y3);

// update state
s[0] += 2 * (t0*xx + t1*(y1-y0));
s[1] += 2 * (t2*(y2-y1) - t1*(y1-y0));
s[2] += 2 * (t3*(y3-y2) - t2*(y2-y1));
s[3] += 2 * (-t4*(y3) - t3*(y3-y2));

zi = in;
out = y3*r;

Shouldn't xx = (zi - r*y3); actually be xx = (in - r*y3); and ditto for your y3 calculation "t0*zi" should be "t0*in", and r = (7.f * resonance + 0.5f); should be r = (17.f * resonance + 0.5f); with a loop gain of 8 it barely gives a squelch, or am I missing something here?

Regards
Andrew

[Ivan_C]

About zi/in you are right, it doesn't make any sense to take x[n+1] as the reference somewhere and x[n] in the other place. I think I followed a little too much mystran's code last summer, and I didn't checked the change he made there a few months ago.

About the resonance, the real problem comes from the consideration of the terms inside the tanh in the original equations, which were omitted for simplification purposes. I saw that if I use the exact dividers from the equations, I got something which saturates too easily, so I did a few tests with other values, and finally I dismissed these dividers. I imagine that you have to match the input amplitude values of the original circuits, and a few things in the modeling if you want to keep them as they are supposed to be, and having the resonance control working with the range from the Fontana's article.

However, something I decided to keep later and which is not included in this source code is the ratio between the dividers in each tanh. If you look into the original equations, you will see that sometimes the divider is 2*Vt, sometimes 2*n*Vt, sometimes 6*n*Vt or something like that.

Sorry if I didn't update the code sooner

PS. I have updated the original source code I posted.

[Ichad.c]

However, something I decided to keep later and which is not included in this source code is the ratio between the dividers in each tanh. If you look into the original equations, you will see that sometimes the divider is 2*Vt, sometimes 2*n*Vt, sometimes 6*n*Vt or something like that.

Okay, that makes a bit more sense, should have read the Fontana paper with more care Btw, in your original equations:

Code: Select all

dvc1/dt = f.(tanh((vin-vout)/(2*Vt)) + tanh((vc2-vc1) / (2*gamma)))

2*Vt = 1, so that can be omitted. Think I need to throw this into Maxima and see what pops out. Always actually wondered how different in output Maxima vs. Matlab is. Finding an efficient form is quite tedious. Once Maxima's 'simplifier' substituted 4 adds for 3 divides

Cheers
Andrew