Newton's Method for vectorial functions (4-pole filter)

[Urs]

Hi all,

As I'm preparing a tutorial for iteratively solving implicit filter equations, I'd thought I'd post some code here for a facts check. It's extracted from some more serious code I'm working on, stripped down to its most simple form, hacked into Maple and solved there, then copy pasted back into a header file.

The algorithm seems to work, but I'm a bit puzzled that resonance coeff has to be negative. I might have missed a - somewhere else.

Note that it's based on idealised classic 4-pole cascade filters as found in synths with CEM 3320 or SSM 2040 chips. Thus the stages are connected through inverting feedback inputs of the single lowpass filters. This is very different from ladder filters.

http://www.u-he.com/downloads/UrsBlog/f ... ectorial.h

The method will be discussed later and in depth in my blog.

Thanks,

- Urs

[hugoderwolf]

Looks good! Thanks again for putting all that out here!

I've got a similar structure working, but unexpectedly it sounds quite good when just solving the individual stages independently but in parallel. Theoretically that shouldn't work, but somehow it does. Comparing with the "correct" solution is on the list.

Regarding that sign thing: your model is a bit different to what I've come up with, as you put the minus sign inside the tanh. On my randomly cluttered sheets of paper, the sign is in the update equation (x1 = s1 - g * (tv1E) - v1E). Of course that doesn't change much, but I think that explains why you suddenly need to flip the sign on the feedback coeff. If you change everything to "+ k4pole * v4E" you should be fine.

[Urs]

unexpectedly it sounds quite good when just solving the individual stages independently but in parallel. Theoretically that shouldn't work, but somehow it does.

It's perfectly well suited for parallel processing. You can approximate the 4 tanh terms in SSE, then compute the 4 residuals in SSE, then shuffle and calculate the m-terms, then keep shuffling and calculate the y terms in SSE. Then min/max the error threshold and take the mask as a skalar for the loop, then update the state in SSE.

Because of the nature of guesses, the 4 stages are independent during the computation of each step.

[Urs]

As for the -... I'll redo this whole thing in Maple and check it out again. Or maybe I had a mental mishap when I put the first equations together.

[hugoderwolf]

What I meant was that I iterate each stage individually as they were only one pole, so 4 independent parallel running newton methods so to speak. It's as if all columns but the first in the jacobian were zeros.

Did you have a look at the assembly from that code? I could imagine that a current compiler can SSE-ize quite a bit of that code pretty well. My experience so far has been that guys like gcc or clang do a much better job at vectorizing code than I do. Using intrinsics was always slower than standard code tweaked for good compiler-optimizability.

[Ivan_C]

Thanks for sharing this Urs !

In the comments, I think there are a few wrong signs (lines 81-84), but for me in the equations everything looks fine. I'm going to run the processing class at one moment of course

I have a question related with your implementation. Is it possible to make the Newton-Raphson algorithm diverging here, for example if you have an input signal with a very high amplitude, or a cutoff frequency close to 20 kHz ?

[Urs]

You can speed convergence up by a better guess. A linear estimate works great.

I don't think it'll ever not converge. All functions are ... (Look up "stetig monoton steigend" in a german-English dictionary ).

The main problem is floating point truncation errors in the residual check. High input values or excessive resonance make it difficult. The delta method proposed by Andy might resolve this, or a dynamic adjustment of the threshold.

[Ivan_C]

I think this is "constantly monotonically increasing" in english Functions must also be Lipschitz continuous (locally) for the Newton-Raphson algorithm to work properly. The two conditions state that the absolute value derivative of the function must be strictly superior to zero AND bounded, which is the case for equations with tanh, but not with sinh for example...

What is this delta method exactly ?

[hugoderwolf]

The main problem is floating point truncation errors in the residual check. High input values or excessive resonance make it difficult. The delta method proposed by Andy might resolve this, or a dynamic adjustment of the threshold.

I think the threshold should just be relative to the absolute value, that would result in a SNR-like definition rather than a noise floor limit. For example something like (error < 0.0001 * output) would make sure the error is at least 80 dB below the signal level.

[Urs]

What is this delta method exactly ?

Check out Andy's papers about using sin() for the coefficient.

In your filter algorithm

Vout = g*(Vin-Vout)+s

substitute Vout by s + dV:

s + dV = g*(Vin-(s+dV))+s

and resolve for dV

dV = g*(Vin-s)/(g+1)

This way you don't calculate the Voltage itself, you just calculate the difference / delta between two steps.

You get Vout simply by Vout = s + dV and you update s by s += 2*dV

Andy showed that this method is numerically superior, I suppose because there's less floating point truncation.

In his paper he also demonstrates a rearrangement of terms that allows for coefficient computation with sin() rather than tan(), but I wasn't successful at applying it to my current set of equations, hence I postponed to learn it.

[Ivan_C]

Interesting idea, I'm going to have a look for that, thanks

[Urs]

I think I also found the sign error that was puzzling me.

I should rewrite each tanh term like this:

Vpositive = 0
Vnegative = Vout + Vin [- k4Pole * feedback]
tVE = tanh( VEpositive - VEnegative)
VE = g * tVE + s

And then, all inputs stay negative but feedback suddenly becomes positive

[itoa]

Guys,

What is better for solving ladder filter with NR? Iterate each stage separately (and then use new outputs for next iteration) or dealing with matrix/Jacobian approach? *better = cheaper

[Urs]

Guys,

What is better for solving ladder filter with NR? Iterate each stage separately (and then use new outputs for next iteration) or dealing with matrix/Jacobian approach? *better = cheaper

As far as I've tried, vectorial NR outperforms each stage NR and then feedback NR around it by a magnitude.

[Urs]

(Arrrgh... I still haven't posted the blog update because of a benign sign error... will try to finalize it soon)