NON-ITERATIVE SCHEMES FOR THE SIMULATION OF NONLINEAR AUDIO CIRCUITS

In this work, a number of numerical schemes are presented in the context of virtual-analog simulation. The schemes are linearly-implicit in character, and hence directly solvable without iterative methods. Schemes of increasing order of accuracy are constructed, and convergence and stability conditions are proven formally. The schemes are able to handle stiff problems very efficiently, because of their fast update, and can be run at higher sample rates to reduce aliasing. The cases of the diode clipper and ring modulator are investigated in detail, including several numerical examples


INTRODUCTION
The design of effective numerical integrators for audio rendering requires a balance between accuracy and efficiency, alongside the overriding constraint of stable operation.In virtual-analog simulation of electronic circuits, established designs, such as the trapezoid or the midpoint methods [1,2] are preeminent.These methods have the virtue of simplicity and robustness: in particular, they have the interpretation of a bilinear transform under linear conditions, in the frequency domain [3].
The trapezoid method underlies Wave Digital Filters (WDF) [4] guaranteeing passivity in the discrete case for linear systems.Wave-based methods have been extended to include many nonlinear systems, see e.g.[5,6,7,8], as well as general linear multistep integrators with variable step size [9].
Kirchhoff-domain methods include Port-Hamiltonian Systems (PHS) [10,11].These are a generalisation of Hamiltonian systems including energy storage components, dissipative components, and connection ports.They have the property of preserving passivity in the discrete setting, when discretisation is performed via the quotient method [12], which is a type of implicit numerical method.State-space models are also prominent [13,14,15].
For nonlinear systems, most passivity-preserving numerical integrators require the solution of a system of algebraic nonlinear equations at each time step.This is usually accomplished iteratively, via suitable algorithms such as Newton-Raphson.While the implementation of algebraic root-finders poses little difficulty in theory, the questions of existence and uniqueness of the update, as well as that of tolerance thresholds present problems of their own which do require some care [16].For example, while these methods are formally second-order accurate, the convergence rate could be impacted by too low a tolerance threshold.For the same reason, numerical instability can occur.
The question of accuracy may be approached in terms of orderaccuracy; higher-order accurate schemes have faster convergence rates.Standard convergence rates hold in the low-frequency limit, but for the purpose of rendering high-quality audio, one is usually interested in wideband numerical behaviour.First-order accurate schemes were shown to be less prone to spurious oscillations than the trapezoid rule, for high input voltages in the diode clipper [17].Control over numerical bandwidth expansion (aliasing) is also of prime importance for audio applications, and higher-order schemes may not have better behaviour in this respect.
In this work, the possibility of exploiting non-iterative schemes for state-space models is explored.A number of schemes is presented, such that the update can be computed explicitly, without iterative methods, though requiring the solution of a linear system.Methods of this kind are said to be linearly-(or semi-) implicit [18,19].A family of schemes of increasing order of accuracy is constructed.A scheme, formally first-order accurate, is presented.The scheme is unconditionally stable, non-iterative, and capable of reducing aliased frequencies in amplitude.A number of numerical tests, using both Matlab and C++ implementations, show that the non-iterative schemes compare favourably to trapezoid and midpoint, in that they can run at higher rates with reduced aliasing, including input/output resampling, whilst using roughly the same amount of CPU resources.Such methods have the additional advantage of sidestepping the machinery of iterative methods (including design choices such as the maximum number of iterations and tolerance).
The article is structured as follows.In Section 2, non-iterative schemes are introduced, and formal proofs of accuracy, stability and convergence are given.The properties of the trapezoid and midpoint methods are also outlined, and stability conditions given for cases of interest in virtual-analog.Section 3 presents the case of the diode clipper, described by a single stiff differential equation, and in Section 4 the more involved case of the ring modulator is presented.Numerical examples are presented throughout.

PRELIMINARIES
Consider first a basic scalar test problem of the form This equation is zero-input and nonlinear but autonomous (so that the dependence of the nonlinear function f on time is through its argument x(t) only, and not on any externally supplied control signal).It describes the time evolution of a time-dependent quantity x (such as e.g.voltage in a virtual analog application), and is initialised with x(0) = x0.Existence and uniqueness of the solution to (1) are guaranteed under the assumption of Lipschitz continuity for f [1].Furthermore, the function f is further assumed to satisfy the following conditions: sign(x) = sign(f (x)) (2a) for ϵ small, there exists The first condition is referred to as sector-boundedness, here to sector [0, ∞], the usual condition for passivity; the second condition enforces boundedness of f /x at the origin; the third condition (typical for virtual-analog systems) expresses monotonicity of f .

Finite Difference Schemes
Equation (1) will be integrated numerically using a time-stepping method with constant time step k (in seconds, with associated sample rate fs = 1/k).Here, the notation for the time step is borrowed from [1], as opposed to the h notation presented in other textbooks such as [2].Since k will never be used to denote an index, this should not generate confusion.The continuous function x(t) is approximated by a time series x n at times tn = kn for integer n.
The error E n at time step n is defined as Definitions of time difference and averaging operators are given here as Similar definitions hold when the operators are applied to the function f , rather than on x itself, so that

Trapezoid method
A standard integrator is obtained by applying trapezoidal integration of (1) on the interval tn ≤ t ≤ tn+1, yielding This is a one-step method, belonging to the more general family of implicit Adams-Moulton methods, i.e. a class of linear multistep methods [1].Detailed analysis of accuracy, convergence and stability of scheme (6) are given in many textbooks on numerical integration, see e.g.[1,2], and are briefly recalled here.In particular, the method is second-order accurate, i.e.
Since the global error is bounded, the method is, in general, zerostable (see e.g.[1] for a definition of zero-stability), and therefore convergent for sufficiently small k.Because f here also satisfies conditions (2a) and (2b), scheme (6) is in fact unconditionally stable.To see this, one expands out the operators in (6) to get Now, define α = (k/2)(f /x).From the above, one has But, because f is sector-bounded to [0, ∞], one has α ≥ 0, and hence |x n+1 | ≤ |x n |, and the solution decays monotonically in x.
Due to its stability properties and simple design, the trapezoid method is a popular choice for nonlinear systems such as those encountered in state-space models, see e.g.[13].The main drawback is its fully implicit character, requiring the solution of a nonlinear algebraic equation at each time step (in this case ( 7)).In general, a root-finding algorithm such as Newton-Raphson will be necessary.This results in various well-known practical difficulties, including the problem of choosing an appropriate threshold in the root-finding algorithm, a maximum number of iterations to preclude stalling, as well as the undesirable characteristic of variable operational cost at each time step, due to variations in the number of iterations required [16].

Midpoint method
Another popular integrator is given by the midpoint method, Like the trapezoid method, this scheme is second-order accurate, i.e. |E n | = O(k 2 ).In general, this method is zero-stable, but owing to the sector-boundedness property (2a), it is in fact unconditionally stable.This is proven easily by multiplying both sides of (10) by µ+x n , to get The midpoint rule also requires the solution of a nonlinear algebraic equation at each time-step.Under linear conditions, both the trapezoid rule and midpoint rule have the interpretation of a bilinear transformation in the frequency domain.

Non-iterative schemes
Since (2b) ensures boundedness of f near the origin, it is natural to attempt to evaluate the nonlinear function at the time step n, whilst maintaing a semi-implicit realisation via multiplication by the factor µ+x n /x n .Thus, consider the following family of schemes, approximating : Here, σ n p = σ n p (k) is a coefficient depending on the current timestep n, as well as on an order p.The coefficient σp may be chosen so that scheme (12) satisfies increasing orders of accuracy.This technique has strong links to the modified equation methods, such as the ones presented in [20].Using Taylor series arguments, one can work out expressions for σ n p as (see also Appendix A): In (13a), a ≥ 0 is a free parameter.In these schemes, the nonlinear functions are evaluated at previous time-steps, hence the update may be computed simply as Methods of this kind are sometimes referred to as semi-implicit or linearly-implicit methods [18,19], in that the implicit update is linear in x n+1 .In the vector case, such as in the ring modulator below, this translates to the solution of a linear system per update.
Using standard arguments, such as those presented in Appendix A, one may demonstrate that |E n | = O(k p ).Hence, the schemes are p th -order accurate, zero-stable and convergent in the limit k → 0. See Figure 1 for a demonstration of order of accuracy for these schemes in the case of a cubic nonlinearity.Stronger forms of stability, allowing for stable computations with a given value of k, rather than just in the limit k → 0, can be given.Stability conditions follow immediately from a cursory examination of ( 14): if βp ≥ 0 (i.e. when 1 + σp > 0), then the numerical solution is monotonically non-increasing at all time steps.By virtue of conditions (2), one has σ1 ≥ 0 and therefore the first-order accurate scheme is unconditionally stable.If f ′ ≥ f /x, then σ2 ≥ 0 and therefore the second-order accurate scheme is also unconditionally stable in this case.Similar sufficient conditions can be checked for the higher-order schemes, and this can be done on a case-by-case basis.

Sources
In the examples below, (1) will be modified to include a timedependent source term as follows In the discrete setting, in order to preserve the order of accuracy to at least second-order, one may approximate the source term as A stability analysis including the source term may be carried out, though it is not shown here for brevity.

CASE STUDY: DIODE CLIPPER
Outside of its relevance in the context of virtual-analog simulation, the diode clipper is particularly interesting from a numerical standpoint.Explicit numerical designs (such as e.g.Forward Euler, or the Runge-Kutta RK4 scheme) are known to fail here, unless the time step is chosen to very small values compared to the timescales of the computed solution [19].The differential equation is stiff in this case, as the exponentially-unbounded nonlinearity accounts for a fast variation in the transients of the computed solution.For these reasons, the diode clipper is used extensively as a test case for numerical designs.Yeh et.al. [19] tested several classic numerical integrators; Werner et.al. [7] presented a WDF realisation; Falaize and Hélie [10] presented a PHS discretisation; Fontana and Bozzo [16] investigated the system in terms of basins of attraction for Newton-Raphson; Holters [21] presented an antiderivative-antialiasing scheme; Parker et.al. [15] turned to neural networks.Following e.g.[19], a model, comprising input, is given in the form (15), where Here, v(t) and x(t) represent, respectively, the input and output voltages.Constants are given as: resistance R = 10 3 Ω, capacitance C = 3.3 • 10 −8 F, saturation current Is = 2.52 • 10 −9 A and thermal voltage vt = 2.6 • 10 −2 V. Here, f (x) clearly satisfies conditions (2).Furthermore, in this case f ′ ≥ f /x, and hence the second-order non-iterative scheme is also unconditionally stable.

Numerical Experiments
As a first experiment, consider Figure 2. The schemes are run using an input of the form of a sinusoid of increasing amplitude.The waveforms present in all cases comparable errors.Yeh et al. [19] observed that low-order methods are sufficient for the purpose of rendering audio signals, and this observation is confirmed here.
The bandwidth expansion due to the stiffening nonlinearity is visible in Figure 3.The figure presents the output spectrograms to an input linear sine sweep with constant peak amplitude.One observes that lower-order schemes perform somewhat better in this respect, with some evident aliasing taking place for the third and fourth-order accurate schemes.It is known that anti-aliasing can be achieved with the use of antiderivatives of the nonlinearity [22,23,21].Here, higher-order schemes work in the opposite Simulations are run at 4X audio rate, and the error is computed as the difference between each simulation, and a 100X oversampled solution.Trapezoid integration is solved using Newton-Raphson with tolerance threshold 10 −15 .Colour scheme: trapezoid (solid black), p = 1, a = 1 (blue, dash-dotted), p = 3 (solid grey), p = 4 (dashed magenta).Resampling to audio rate is achieved by means of a 12 th order Butterworth filter, with normalized cutoff frequency 0.8π.sense, through the use of higher derivatives of the nonlinear function, leading presumably to the increased effects of aliasing with order.Furthermore, higher order schemes have a higher frequency bandwidth and the lack of information of sampled discrete timestepping yields aliasing of the high-frequency spectrum (indeed, sampling the exact solution x(t) would also lead to aliasing.) Sound examples are available at [24], and Matlab sample code, illustrating the use of the non-iterative scheme p = 1 in the diode clipper case, is given in Appendix B.
The performance of the schemes depends on a number of factors.For the iterative schemes, one needs to decide on an appropriate tolerance threshold to exit the iterative loop.Generally, one wishes to compute a solution with sufficient accuracy, while avoiding an excessive number of iterations.Because the stiffening of the system due to the nonlinearity will result in a higher number of iterations, one must be careful when setting an upper bound on the number of iterations.On the other hand, the non-iterative schemes work at a fixed, predictable cost per time step.
For the new designs presented here, a higher oversampling factor is generally required, under stiff conditions, to reduce aliasing.However, given the simplicity of the update, it may still be cheaper to use the oversampled non-iterative schemes, including resampling, than the iterative schemes at audio rate, particularly for systems requiring linear system solves at each iteration, such as the ring modulator below.Leaving aside the cost of resampling, one may argue that the non-iterative schemes p = 1, p = 2 require roughly the same number of operations per time-step as the iterative schemes need per iteration-that is, two nonlinear function calls, and a similar amount of multiplies and divides.It is appropriate, then, to run the non-iterative schemes at higher rates, Sound examples available at [24].and to compare the performances of these against the trapezoid and midpoint methods.Considering now Figure 4, it is seen that the iterative schemes require an average of about 5 or 6 iterations at audio rate.It was decided then to run the non-iterative schemes at an oversampling factor of 4.
The role of the free parameter a in (13a) may be better understood by inspection of Figure 5: choosing a higher value will result in lower aliasing overall, as well as higher low-pass filtering.In general, one wishes to work at higher rates here, in order to avoid too steep a roll-off at frequencies of interest in the high range, though one may compensate by equalising the output accordingly.

CASE STUDY: RING MODULATOR
The purpose of this section is to test the non-iterative schemes in the case of a system of nonlinear equations.A mathematical model of the system was given in [25].Similarly to the diode clipper, the ring modulator was treated in several works, see e.g.[26,27,9, where Here, x = [q ⊺ , i ⊺ ] ⊺ .The vector q = [q1, q2, q3] ⊺ contains state voltages, and the vector i = [i1, i2] ⊺ contains state currents.The output is y = q2.In (18), the linear part has been separated out from the nonlinearity, for ease of notation, but one should easily recognise that (18) is the vector equivalent of ( 15), with input um = Hmvm(t).The matrices are composed of constant coefficients from the circuit resistors, capacitors and inductors, and are as: The vector f = [f (w1), f (w2), f (w3), f (w4), 0] ⊺ includes the nonlinear current-voltage relationships depending on the voltages w = [w1, w2, w3, w4, 0] ⊺ , defined as where Finally, vc(t), vm(t) are the carrier and modulator input voltages.
Here, the nonlinear current-voltage relationships are given by the Shockley diode equation, i.e.
which clearly satisfies conditions (2).It is remarked that the vector f containing the nonlinearities is here composed of four scalar nonlinearities depending on a single scalar input, and hence conditions (2) may be checked easily componentwise.
Constants are given as:

Finite Difference Schemes
The numerical schemes used here are an extension to the vector case of the schemes presented above.In the zero-input case, the trapezoid and midpoint methods are obtained as, respectively, Only the first and second-order non-iterative methods will be considered here.They read where The first-order non-iterative method is unconditionally stable (for a ≥ 0), because the eigenvalues of the matrix σ1 are non-negative.
For σ2, a condition on k arises formally (not shown here for brevity), though practically it is met using reference time steps.Hence, both schemes can be treated as absolutely stable.For all schemes, the modulator and carrier source terms can be approximated as, respectively, Hmµ+v n m , Hcµ+v n c .Though these schemes have been written compactly as 5 × 5 systems, it is convenient to reduce the size of the update by expressing the currents in terms of the voltages.For all schemes this is accomplished as One can then work with the voltages q alone, thus reducing the size of the implicit update to 3 × 3.

Numerical Experiments
The responses to sine sweeps are visibile in Figures 6 and 7. Note that the non-iterative schemes are here run at higher rates, since the comparison here is drawn according to compute times: all simulations take roughly the same time to run in C++, as detailed below and seen in Figure 8.Clearly, the non-iterative schemes present much lower aliasing than trapezoid and midpoint.Most importantly, the run times of the non-iterative schemes are insensitive to stiffness, as visible in Figure 8.As pointed out previously, whilst both trapezoid and midpoint will require more iterations per time step as the input carrier amplitude is increased, the non-iterative schemes will always operate at the same cost.This is a particularly important aspect in view of any real-time application requiring a precise allocation of CPU resources.
Here, real-world C++ performance of two versions of the ring modulator system are compared, i.e. the trapezoid method and p = 2. first uses the Newton-Raphson iterative method and runs at an audio rate of 44.1kHz.The second runs at four times the rate, 176.4kHz.This requires both up-sampling of the input signal and down-sampling of the output, which is included in the testing.The CPU time was measured to run each system for 1 second of audio rate output, so 44100 time-steps for the iterative version and 176400 time-steps for the non-iterative.
A brief description of the algorithms in terms of their computations at each time-step is given here.The calculations are predominately algebraic, with the exception of calls to the exponential function.It is a vector system with a state size of three elements, so the calculations consist of additions and multiplications on small vectors and matrices.There is also a linear system solve on a 3 × 3 matrix, which is performed using a simple Gaussian elimination.For the iterative version, the steps are: 1. Update the current state using six vector multiplications and additions, and a matrix by vector multiplication.

Perform Newton-Raphson iterations of:
(a) Set up a 3 × 3 matrix M using a matrix by vector multiplication, and two matrix by matrix multiplications.
(b) Linear system solve on M.

DAFx.6
Proceedings of the 23 rd International Conference on Digital Audio Effects (DAFx2020), Vienna, Austria, September 2020-21 (c) A further matrix by vector multiplication, and three calls to the exponential function std::exp().
(d) Find the maximum absolute value of a vector, and check the tolerance level.
3. Update the voltage vector, and read the output.
The non-iterative version uses similar calculations, but in a single pass: 1. Update the current state using six vector multiplications and additions, and a matrix by vector multiplication.
3. Set up the 3 × 3 matrix M as above but using two extra matrix additions and matrix by vector multiplications.
4. Linear system solve on M.

Update the voltage vector, and read the output.
There is also the up-sampling and down-sampling at either end of the algorithm.For this testing a 12 th order IIR low-pass filter was used.The C++ codes were compiled using Clang with -O3 optimisation level, and it was noted that the compiler was not able to auto-vectorize any elements of the code.The executables were run on a Mac mini M1, giving the timings plotted in Figure 8.
For the non-iterative version, 22ms of the computation time was attributed to the re-sampling code.This testing shows that the noniterative scheme can be run at 4X the sample rate of the iterative version for the same level of computational cost, achieving a large reduction in aliasing artefacts.

CONCLUSIONS
In this work, non-iterative schemes of increasing order of accuracy were offered, and their performance compared to standard integrators such as the trapezoid and midpoint methods.These schemes, linearly implicit in character, require the solution of a single linear system per update, in contrast to trapezoid and midpoint for which Newton-Raphson is generally required.It was shown that, whilst the convergence properties of the schemes follow the expected trends, higher-order schemes are somewhat less suited for the purpose of rendering wideband audio.Furthermore, for higher-order schemes stability conditions can only be checked on a case-by-case basis.On the other hand, a particular form of the non-iterative schemes was given here: whilst formally first-order accurate, this scheme is always unconditionally stable and thus it preserves the passivity of the continuous system.A free parameter can be adjusted in the scheme to control the amount of low-pass filtering induced by the scheme, thus reducing aliased artefact very efficiently.Whilst this scheme must usually be run at higher rates, compared to trapezoid and midpoint, it is generally as efficient, even when the cost of resampling of the input/output is included.Furthermore, this scheme has the advantage of always running at the same cost, thus avoiding design choices such as the maximum number of iterations and tolerance.Extensions to the multivariate vector case (i.e.where the nonlinearity is composed of multivariate nonlinear functions), as well as further investigations on the role of the free parameter a, will be the subject of future work.

ACKNOWLEDGMENTS
The first author wishes to thank the anonymous reviewers for their insightful comments.This work was partly supported by the Eu-ropean Research Council (ERC), under grant 2020-StG-950084-NEMUS.

Figure 1 :
Figure 1: Global error |E n | using a cubic nonlinearity, as a function of sample rate.Here f (x) = x 3 , and x(t) = 2t + x −2 0 −1/2 sign(x0).Here, x0 = 1.3, and the error is computed at t = 0.2 s.For the trapezoid and midpoint rules, Newton-Raphson is used for the update, with a tolerance threshold set at 10 −15 .

Figure 2 :
Figure 2: Diode clipper.Response to an input sinusoid.Here, the input is v(t) = v0 sin(1000πt), with v0 as indicated in each panel.Simulations are run at 4X audio rate, and the error is computed as the difference between each simulation, and a 100X oversampled solution.Trapezoid integration is solved using Newton-Raphson with tolerance threshold 10 −15 .Colour scheme: trapezoid (solid black), p = 1, a = 1 (blue, dash-dotted), p = 3 (solid grey), p = 4 (dashed magenta).Resampling to audio rate is achieved by means of a 12 th order Butterworth filter, with normalized cutoff frequency 0.8π.

Figure 3 :
Figure 3: Diode clipper.Response spectrograms for a linear sine sweep, v(t) = sin(γ0t 2 ).Trapezoid and midpoint are run using Newton-Raphson with tolerance threshold 10 −15 , at 2X audio rate.The non-iterative schemes are run at 4X audio rate.Sound examples available at [24].

Figure 8 :
Figure 8: Ring Modulator.C++ compute times for one second of output.Here, the carrier is vc(t) = max (vc) sin(2000πt).The modulator input is vm = sin(2000πt).Colour scheme: trapezoid at audio rate (red), trapezoid at 2X audio rate (black), p = 2 at 4X audio rate (green).Linear systems are solved using Gaussian elimination, and resampling time is included.For trapezoid at 2X, 14 ms are attributed to resampling, while for p = 2, 22 ms are attributed to resampling.