Nonlinear ZDF SVF part 2

Previously we derived a nonlinear zero delay feedback filter where the nonlinear tanh() elements were placed at the inputs of the integrators. This is not ideal, since it will lead to instabilities at high input amplitudes even with implicit integration.

Another way to model nonlinearities present in active analog integrators is to place the tanh() clipping at the “output” of the integrator, where it will limit the system to the $[-1, 1]$ interval.

Nonlinear “derivative”

Expressing tanh()-limited integration as a differential equation is not really possible with traditional methods since it breaks the fundamental theorem of calculus and the definition of the derivative as a linear difference quotient.

However, it is easy to express as an equation of discrete integration where for example

taken as the linear slope approximation (Euler integration)

becomes

which is more than enough for the needs of dealing with approximations of differential equations as difference equations.

Nonlinear SVF difference equation

Let’s now use this nonlinear discrete integrator as a basis for building the discrete-time version of the state-variable filter.

Taking the state-variable filter differential equation

and applying the trapezoidal rule to discretize it as a difference equation yields the following system.

The latter difference equation can be then nonlinearized as the following system.

If we want to use this as a practical time-stepping integration scheme, we need a way to solve for $bp_{t+1}$ to calculate the filter state from. It is possible to solve it from this system numerically with Newton-Raphson iteration as an implicit integration method.

By substitution $a=\frac{1}{2}h$ and gathering terms the equation for $bp_{t+1}$ becomes

Further

and substituting $C_{t} = bp_{t} - alp_{t}-a\sqrt{2}bp_t+au_{t}+au_{t+1}$ as a constant

Using tanh summation identity

and gathering terms to LHS

Further substitutions $D_t=lp_{t}+\frac{1}{2}hbp_{t}$, $E_t=tanh(C_{t})$ and $x=bp_{t+1}$ reduce the equation to

This still looks more complicated than it has to be so substituting the functions

and notating it as a function of x

Newton-Raphson needs a derivative of this function

and of the substituted functions

The iteration for solving $x=bp_{t+1}$ can now be carried over as

from where the full filter state can be integrated as

The rough pseudocode for achieving the full implicit time-stepping is as follows.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 a = 0.5*dt; D_t = lp + a*bp; E_t = tanh(bp - a*lp - a*fb*bp + a*u_t1 + a*input); x_k = tanh(bp + dt*(-lp - fb*bp + input)); // newton-raphson for(ii=0; ii < 32; ii++) { T2 = tanh(a*x_k + D_t); T3 = tanh(-a*fb*x_k - a*T2); DT2 = a - a*T2*T2; DT3 = (1.0 - T3*T3)*(-a*fb - a*DT2); fx = x_k + E_t*x_k*T3 - T3 - E_t; Dfx = 1.0 + E_t*(T3 + x_k*DT3) - DT3; x_k2 = x_k - fx/Dfx; // breaking limit if(abs(x_k2 - x_k) < 1.0e-15) { x_k = x_k2; break; } x_k = x_k2; } lp = tanh(lp + 0.5*dt*(bp + x_k)); bp = x_k; hp = -lp - fb*bp + input; // set input at t-1 u_t1 = input;

Conclusions

This method is implemented in my state-variable filter module for VCV Rack together with other methods from this blog. The performance is excelent with a constant Q factor across the frequency spectrum and practically unconditional stability, although the SVF filter structure together with nonlinear tanh() integration is stable even with the most basic explicit integration methods.

The mathematical derivation can be carried over to other nonlinear filter structures where stability is not so great (such as the Moog ladder filter) with the same basic steps so hopefully this serves the reader mainly as a good example and a learning stepping stone for further experimentation.

Happy hacks.