Integrator filter with delayless feedback path

[antto]

neotec: THANK you
now this makes sense to me, gonna try it
(finaly this thing explained in a simple way without all the hardcore math symbols and terms)

[karrikuh]

It uses less CPU resources than I expected, and the only optimization I can think of at the moment is a nice prediction function for calculation of the initial iteration value.

You can use the technique from the KeepTopology paper for computing an initial value. This will give you perfect prediction for a linear system (no iterations necessary) and will probably still yield a good approximation for the practical nonlinear case.

[neotec]

You can use the technique from the KeepTopology paper for computing an initial value. This will give you perfect prediction for a linear system (no iterations necessary) and will probably still yield a good approximation for the practical nonlinear case.

Yep, I'm currently going through all the KeepTopology stuff again ... and discovered some neat thing. You can build a simple one-liner delayless integrator from this KeepTopology stuff (using a simple integrator, not the bilinear one and transforming the feedbacksolver a bit):

Code: Select all

buffer = (f * input + buffer) / (1.0 + f);

'buffer' equals 'output'.

[rola]

Is this the same as Newton-Raphson?

[neotec]

Na, the only thing it has in common with Newton is the iteration.

And by examining all this delayless feedback stuff ... well, I feel like *facepalm*. It's so easy that's a shame I didn't figure this out by myself

For example the delayless integrator:

Code: Select all

Starting with:
  out += f * (in - out)
Rearrange:
  out = out + f * (in - out)
The last out is the feedback, second out is our previous state (buffer):
  out = buf + f * (in - out)
Now:
  out = buf + f * in - f * out
Solving for out:
  out + f * out = buf + f * in
  (1 + f) * out = buf + f * in
            out = (buf + f * in) / (1 + f)

... ... now let's see what I can do with this stuff

Edit: So, now I have a prediction function for the filter described at the beginning of this thread that always only takes 1 or 2 iterations

[karrikuh]

Is this the same as Newton-Raphson?

No, what neotec is doing is called fixed point iteration to predict the current value of out. A fixed point of a (generally nonlinear) function f is a solution of

x = f(x), or, in our case, out = f(out).

Such a solution can often be found by iterating according to

out(i+1) = f(out(i)) with iteration i = 0,1,2,...

where you have to specify an initial value out(0). In our case evaluating f(out(i)) means just applying the normal filter equations but omitting updating the filter memory.

Note however that in general

• there may exist either zero, one or multiple fixed points for a given f
• convergence of the above iteration strongly depends on both f and the initial value out(0)

[neotec]

So, as the integrator is now delayless I started to make the feedback path delayless.

We're having 4 integrators (0 to 3) in chain:

Code: Select all

b0 = (b0 + f * input) / (1 + f)
b1 = (b1 + f * b0) / (1 + f)
b2 = (b2 + f * b1) / (1 + f)
b3 = (b3 + f * b2) / (1 + f)

For simplicity I replace (1 + f) by 't' in the following.
Our resonant ladder filter looks like:

Code: Select all

out = integrator3( integrator2( integrator1( integrator0( input - out * r ) ) ) )

Replacing the above with the formulas:

Code: Select all

out = (b3 + f * ((b2 + f * ((b1 + f * ((b0 + f * (in - out * r)) / t)) / t)) / t)) / t

And solving for out:

Code: Select all

out = (b3 * g0 + b2 * g1 + b1 * g2 + b0 * g3 + in * g4) / g5

with:

Code: Select all

g0 = t^3
g1 = f * t^2
g2 = f^2 * t
g3 = f^3
g4 = f^4
g5 = t^4 + f^4 * r

... and now we're back at having an awful lot of variables ... but reduced iteration count (1-2 using hard clipper), better stability on lower oversampling ratios and perfectly tuned when 'f = 1-exp(-2PI*fc/fs)' and 'r' still needs the gain correction factor: r = q * (1 + f * PI), where q is in [0.0, 4.0].

Another things: For full self oscillation (using only a Dirac function as input) we still need high oversampling (fs > 10*44100) ... have to check what causes this behaviour.

[mystran]

So, as the integrator is now delayless I started to make the feedback path delayless.
[...]
... and now we're back at having an awful lot of variables ... but reduced iteration count (1-2 using hard clipper), better stability on lower oversampling ratios and perfectly tuned when 'f = 1-exp(-2PI*fc/fs)' and 'r' still needs the gain correction factor: r = q * (1 + f * PI), where q is in [0.0, 4.0].

Another things: For full self oscillation (using only a Dirac function as input) we still need high oversampling (fs > 10*44100) ... have to check what causes this behaviour.

If you want you can replace all the integrators with the bilinear variant. The feedback is solvable for any linear system and since bilinear transform maps an analog filter exactly (except warping the frequency axis) to a digital filter you can have not just perfect tuning (using tan(w/2) prewarp as usual) but also exact feedback gain.

[neotec]

If you want you can replace all the integrators with the bilinear variant. The feedback is solvable for any linear system and since bilinear transform maps an analog filter exactly (except warping the frequency axis) to a digital filter you can have not just perfect tuning (using tan(w/2) prewarp as usual) but also exact feedback gain.

Thanks, I will try the bilinear integrators.

[neotec]

So, I've tried the bilinear integrators. Using bilinear transformed H(s)=1/s in transposed DF II:

Code: Select all

out = buf + f * in
buf = f * in + out

transformed into a delayless one pole lowpass:

Code: Select all

pin = in - (f * in + buf) / (1 + f)
out = buf + f * pin
buf = f * pin + out

and 'f' is now:

Code: Select all

f = tan(PI * fc / fs)

Works great. Perfect tuning, constant resonance gain (no more corrections), thanks mystran.

The interresting thing is that you can put all kinds of non-linearities on the buffers and still have a totally predictable filter (using e.g. the formula I provided)

Code: Select all

buf = saturate(f * pin + out)

ossom^^ thanks

Edit: *gnah* double posting *gnah* sry

[Urs]

Yep.

[mystran]

The interresting thing is that you can put all kinds of non-linearities on the buffers and still have a totally predictable filter (using e.g. the formula I provided)

Wait a minute. This would work if you have a nice analytic formula for the non-linearity, and just throw it in with the rest of the stuff before you solve the feedback, but doesn't it fall apart if you need stuff like absolute values or hard-clip guards with the non-linearity?

[neotec]

Wait a minute. This would work if you have a nice analytic formula for the non-linearity, and just throw it in with the rest of the stuff before you solve the feedback, but doesn't it fall apart if you need stuff like absolute values or hard-clip guards with the non-linearity?

The feedback solver only works with buffer values, input and resonance factor. As long as you don't mess around with the filter input (input - r * feedback) or with the integrator input you can do everything you want with the buffers.

The feedback solver always predicts perfectly, no matter what I do with the buffers. Currently I removed all clipping from the filter input, set a hard-clipper [-2,2] on the first buffer and a soft-clipper on the last buffer. Works perfectly and my iterating test filter outputs an average iteration count of exactly 1.0.

The filter also doesn't explode, no matter what I choose for 'r' (tried 1000.0^^). Just put the hard clippers on the first buffers always, so the other buffers smooth out the 'corners'.

So my filter looks like this at the moment (using the stuff from my previous posts):

Code: Select all

double out = b3 * g0 + b2 * g1 + b1 * g2 + b0 * g3 + input * g4;
            
final double tin = (input - out * r);

final double in0 = tin - (f * tin + b0) * f2;
final double c0 = b0 + f * in0;
b0 = hardClip(f * in0 + c0);

final double in1 = c0 - (f * c0 + b1) * f2;
final double c1 = b1 + f * in1;
b1 = f * in1 + c1;

final double in2 = c1 - (f * c1 + b2) * f2;
final double c2 = b2 + f * in2;
b2 = f * in2 + c2;

final double in3 = c2 - (f * c2 + b3) * f2;
final double c3 = b3 + f * in3;
b3 = softClip(f * in3 + c3);

// f2 = 1.0 / (1.0 + f)
// g0..4 were pre-divided by g5

[antto]

tan(pi*fc/fs)

tan() grows to +inf then goes from -inf back to +inf .. is that correct?!
or is this meant to work only _with_ oversampling?

[mystran]

Oh I see. I thought you ment "buffer" as in "buffer amplifier" as you'd have between stages in an analog ladder. I realize now that you really ment the state variables. So as long as you don't pass the signal through any non-linearities (during a single sample) you can distort the state variables when you store them for the next sample all you want. Sure, that works.