About this zero-delay feedback again (sorry)

[christripledot]

Apologies in advance for being dense. Digital filter theory does bad things to my brain, however many papers I read. Some fantastic threads on this forum have given me a kick-start into writing a zero-delay integrator-based LPF, but I'm struggling with a few concepts. I realise this is old news now, but please show mercy!

I think I'd better introduce my confusions one by one, so here's the first:

A basic integrator LPF is given by:

Code: Select all

iin = in - buf;
buf = buf + f * iin;
// buf contains the output for this sample

Or, put another way,

Code: Select all

iin = in - buf;
out = buf + f * iin;
buf = out;

And (according to neotec's brilliantly informative post on the subject) a bilinear integrator LPF is given by:

Code: Select all

iin = in - buf;
out = buf + f * iin;
buf = f * iin + out;

In both cases, 'out' = buf + f * iin. But while the basic integrator simply stores the output in 'buf', ready for the next sample, the bilinear integrator looks like it performs another filtering step. What's going on here, in layman's terms? Is the 'out' in line 3 different to the 'out' in line 2?

[/code]

[Urs]

Hmmm... in contrast to my earlier beliefs, I'd say that a basic LPF doesn't have a feedback. It has a state.

What we're usually talking about is a more complex filter structure with at least 2 poles/integrators that takes the output (or another signal source) and feeds it back into the input (or another point).

A simple example is a 4-pole ladder filter:

output = 0
input = x - res * y4;
y1 = y1 + f * (input - y1)
y2 = y2 + f * (y1 - y2)
y3 = y3 + f * (y2 - y3)
y4 = y4 + f * (y3 - y4)
output = y4

What you can easily see is that input is composed out of the actual input sample x and the previous output sample y4. This is the lag/delay/latency.

We want the y4 at the input side to be identical to the y4 at the output side. So we

output = 0
guess = y4 // previous output sample... can be done better

while( true )
{
input = x - res * guess;
yb1 = y1 + f * (input - y1)
yb2 = y2 + f * (yb1 - y2)
yb3 = y3 + f * (yb2 - y3)
yb4 = y4 + f * (yb3 - y4)
output = y4

if ( fabs(y4 - guess) > 0.00001 ) break;

guess += (y4 - guess) * 0.1; // crap method, may work though
}

y1 = yb1
y2 = yb2
y3 = yb3
y4 = yb4

So, we calculate the filter with a guess of what y4 will be, but we do not update the state just yet. Then we compare our guess to y4. If it's close enough, then we've just calculated y4 with itself. We break out of the loop and update the state. If y4 and guess are still far apart, then we create a new (hopefully better) guess and calculate the filter again.

This method is hilariously expensive. As long as the filter is linear, you can use numerical methods to create an accurate guess without actually calculating the filter once. But as soon as non-linear elements enter the structure you can't do that anymore. You need to solve the filter iteratively, using e.g. trapezoidal integration, newton raphson or whatever other method you see fit.

Just imagine anything as simple as

yb1 = tanh( y1 + f * (input - y1) )
yb2 = tanh( y2 + f * (yb1 - y2) )
yb3 = tanh( y3 + f * (yb2 - y3) )
yb4 = tanh( y4 + f * (yb3 - y4) )

Even a Newton-Raphson of this (involves building the derivate of that whole thing) is really, really difficult and probably not fast enough to run in realtime. And then, with that amount of nonlinearity you still want to oversample 4x or so...

This is where the fun starts.

I hope this helps a bit.

Urs

[Urs]

(sorry, wife demands attention, I couldn't be bothered to write proper code in proper code tags)

[christripledot]

Please don't worry about the code tags, it was very clear. Thanks for your response!

I assume at the end of the output prediction loop you meant 'output = y4b' etc...

Yes, I understand the principle of estimating the output before performing the filtering. I've knocked up a working ZDF 4-pole already, but I wanted to post the code to check if I'm performing some unnecessary steps, but before I get there I wanted to clear up the difference between the 'regular' and the bilinear integrators. Sticking with one pole for now, and before we even consider feedback, I wondered just what is going on with the different calculation for 'buf' in the bilinear case. Why is it that the calculation for 'out' is the same for both kinds of integrator, but the calculation for 'buf' is different? Or am I just missing something obvious?

The next thing I was going to ask about has to do with iteration (or lack of it). Using the bilinear integrator, the output prediction equation can be solved to give:

out = (buf + f * input) / (1 + f)

Looking at the maths for this, it seems to be a solid solution. Does using a bilinear integrator allow us to accurately predict the output without the need to perform repeated iterations?

[Urs]

I never quite got that part either, or it must have fallen out of my head recently. I assume it has advantages for analytical methods. Might be able to clear that up shortly.

Nevertheless, iterative methods can't be avoided with a bunch of nonlinear elements. They can only be minimised by clever guesses.

[pj geerlings]

I have been meaning to bring this whole thing up again myself - so I'll join this discussion by offering the following (correct I hope) version of Urs' explanation ...

Code: Select all

guess = y4 // previous output sample... can be done better 

while( true ) { 
    input = x - res * guess; 

    yb1 = y1 + f * (input - y1) 
    yb2 = y2 + f * (yb1 - y2) 
    yb3 = y3 + f * (yb2 - y3) 
    yb4 = y4 + f * (yb3 - y4) 
    output = yb4 

    if ( fabs(yb4 - guess) < 0.00001 ) 
        break; 

    guess += (yb4 - guess) * 0.1; // crap method, may work though 
} 

y1 = yb1 
y2 = yb2 
y3 = yb3 
y4 = yb4 

my 1st question ...
How is this :

Code: Select all

out = (buf + f * input) / (1 + f) 

integrated into the code above?

peace y'all
pj

[aciddose]

there is no difference:

memory = memory + input * coefficient - memory * coefficient
or
memory = memory + (input - memory) * coefficient

http://en.wikipedia.org/wiki/Multiplication#Properties

it actually isn't really correct to think of the second version as less operations. you're just combining the operations - there are still multiple things going on in the function.

all you need to do is understand the basic properties of arithmetic and you can break down many such functions to manipulate them in more complex ways.

[christripledot]

This I know and understand. My apologies if I'm being dense, but I still think there's something screwy going on:

1. For the plain vanilla integrator:

Code: Select all

out = buf + f * iin;
buf = out;

2. For the bilinear integrator:

Code: Select all

out = buf + f * iin;
buf = f * iin + out;

The way I understood it, these two lines aren't a pair of simultaneous equations, they're algorithmic steps in pseudocode. Let's expand the 'out' in line 2 with its definition in line 1:

Code: Select all

buf = f*iin + buf + f*iin

How on earth is '2*(f*iin) + buf' equivalent to 'f*iin + buf'?!?

[aciddose]

if you don't make up weird variable names things become a lot easier to understand.

2*coefficient*delta + memory obviously just doubles the frequency. this isn't the right function.

i think what you mean is this:

Code: Select all

 temp = input + buffer * coefficient
 output = buffer - temp * coefficient
 buffer = temp

vs.

 output = buffer - input * coefficient
 buffer = output * coefficient + input

the reciprocal of the coefficient is used in the delay-less version.

(what i don't understand is why people have just 'discovered' these well-known functions recently...)

[christripledot]

All of the pseudocode I've posted has been quoted from neotec's exhaustive posts on the subject in an earlier thread (including the variable names). I thought something was up with that particular definition of the bilinear integrator; it's nice to have it confirmed.

Perhaps I've misunderstood the context in the first place. All I really wanted to know in the first place is: is there any difference in using 'bilinear' integrators (whatever they are) to build a low-pass as opposed to the usual integrator definition?

I'm guessing not; perhaps it just has to do with having an equation that can be rearranged to more easily solve the feedback when building a resonant filter?

EDIT:

Thanks, just saw your edit. Yes, either I missed a 'minus' sign or there was a typo in neotec's code. Now I see quite clearly that they are equivalent.

Re: sudden discovery of these functions:
Well, it's only in recent days that I've tried my hand at writing a resonant filter, and this forum is one of the few places with decent explanations that don't involve heavyweight mathematical notation. I'm no idiot when it comes to maths, but I'm not familiar with the notation, so I really appreciate there being so many helpful folk here willing to discuss it in plain English

[Caco]

Even a Newton-Raphson of this (involves building the derivate of that whole thing) is really, really difficult and probably not fast enough to run in realtime. And then, with that amount of nonlinearity you still want to oversample 4x or so...

This is where the fun starts.

I hope this helps a bit.

Urs

Another fun property of all this is that the more you oversample then the fewer iterations you need to do so oversampling ends up being really cheap

Generally I haven't found a huge improvement in the sound though moving to delayless filters in my code. It certainly helps but I only really notice major improvements in specific situations, such as when modulating the filter quickly, so it doesn't suddenly turn an average filter into some super analog sounding emulation. A rubbish filter made delayless still sounds rubbish, a good filter made delayless can sound better

[Urs]

Another fun property of all this is that the more you oversample then the fewer iterations you need to do so oversampling ends up being really cheap

Unless your prediction comes up with a correct initial guess most of the time

Generally I haven't found a huge improvement in the sound though moving to delayless filters in my code. It certainly helps but I only really notice major improvements in specific situations, such as when modulating the filter quickly, so it doesn't suddenly turn an average filter into some super analog sounding emulation. A rubbish filter made delayless still sounds rubbish, a good filter made delayless can sound better

Very, very true!

[cheppner]

Well, let me share some additional insights

What an Integrator does, is, it integrates. Many times we can simply assume, that in order to get the integral we just have to sum up our values.

This is called Trapezoidal Rule and it gives you an 3rd order error.

In fact, a linear interpolator looks like this:

Code: Select all

out = buf + f * in
buf = out

And, afaik a bilinear interpolator looks like this:

Code: Select all

tmp = f/2*in+buf
out = tmp+buf
buf = tmp

Now let's look at our integrators for several steps (assuming f = 1 and buffer=0 for a start):
the linear one gives

Code: Select all

f*in1, f*(in1+in2), f*(in1+in2+in3), …

the bilinear one gives

Code: Select all

f/2*in1, f*in1+f/2*in2, f*(in1+in2)+f/2*in3, …

In other words: It only integrates HALF the input sample in a step or it integrates to the middle of the sample.


Then there is a different view. What we are ultimately doing is transferring analog circuit equations to the digital domain (s- to z-plane). The question often is, how do we translate a capacitor/integrator? One answer is, with a linear interpolator or a bilinear one where the bilinear one is equivalent to using the Bilinear transform.
The latter matches amplitude and phase response of the filter, giving you an advantage in many cases.

I hope this clears up where the difference lies (in time and frequency-view).

Cheers,
- Clemens

BTW: I recommend this

[Caco]

Another fun property of all this is that the more you oversample then the fewer iterations you need to do so oversampling ends up being really cheap

Unless your prediction comes up with a correct initial guess most of the time

Cool, at the moment I average around four iterations at most filter settings but the iterations creeps up quite a bit at extreme resonance and cutoffs though.

[Urs]

Another fun property of all this is that the more you oversample then the fewer iterations you need to do so oversampling ends up being really cheap

Unless your prediction comes up with a correct initial guess most of the time

Cool, at the moment I average around four iterations at most filter settings but the iterations creeps up quite a bit at extreme resonance and cutoffs though.

Yes, high resonance becomes a cpu killer.

But take an example: In a synthesizer you actually know what's going into the filter (oscillator waveforms, white noise etc.). And once the filter goes into selfoscillation you also know what's happening. So you can take these cases into account when preparing your initial guess.

Urs