APPLICATIONS OF PORT HAMILTONIAN METHODS TO NON-ITERATIVE STABLE SIMULATIONS OF THE KORG35 AND MOOG 4-POLE VCF

This paper presents an application of the port Hamiltonian formalism to the nonlinear simulation of the OTA-based Korg35 filter circuit and the Moog 4-pole ladder filter circuit. Lyapunov analysis is used with their state-space representations to guarantee zero-input stability over the range of parameters consistent with the actual circuits. A zero-input stable non-iterative discrete-time scheme based on a discrete gradient and a change of state variables is shown along with numerical simulations. Simulations show behavior consistent with the actual operation of the circuits, e.g., self-oscillation, and are found to be stable and have lower computational cost compared to iterative methods.

One major concern is robust behaviour, and in particular the maintenance of numerical stability under nonlinear conditions; standard numerical integration strategies may not guarantee stability, and are thus prone to failure in a practical setting.Various simulation frameworks have emerged, including wave digital filters (WDFs) [12,13,14,15], a modularized simulation framework operating using wave variables and scattering operations, and newer approaches are based on discrete energy conservation and more generally passivity (also enforced using WDFs), but operate directly using physical quantities such as voltage and current.The port Hamiltonian (pH) [16] approach is one such methodology, and forms the basis for the work presented in this paper.
One drawback of the use of such methods (compared, e.g., with simpler numerical designs without a stability guarantee, such as forward Euler) is a reliance on iterative solvers such as Newton-Raphson iterations.An example of the use of Newton-Raphson iterations for modeling the Moog 4-pole ladder filter is [17].However, for a special class of pH systems as described in [18] and reviewed here in Section 2.2.1, a non-iterative scheme can be obtained by performing a change of state variables.The resulting scheme does not require an iterative method to solve for the state update.Computational cost relative to iterative methods is substantially reduced.Perhaps more importantly, many additional design considerations following from the use of iterative solvers can be sidestepped: these include convergence guarantees; the choice of tolerance and maximum number of iterations; choice of initial value; and the undesirable property (in audio) of variable load if the number of iterations changes from one time step to the next.Where applicable, this method is used in this work.
The contents of this paper are as follows.Section 2 is a review of the port Hamiltonian (pH) simulation framework, including the discretization of pH systems using a discrete gradient scheme, and an explicit scheme using a state change for a specific class of pH systems.In Section 3, the zero-input Lyapunov stability properties of the OTA based Korg35 LPF and the Moog 4-pole ladder filter along with their zero-input stable discrete-time schemes are presented.Numerical results appear in Section 4 and some concluding remarks and suggestions for future work appear in Section 5.

Mathematical Formulation
Input-State-Output Port Hamiltonian Systems [16,19] are a special case of port Hamiltonian (pH) Systems whose dynamical equations are of the form In the above, • x ∈ R N represents the state variables of the system.
• u ∈ R M represents the input sources.
• y ∈ R M represents the output flows corresponding to the input sources.
• H(x) ∈ R is the Storage Function or Hamiltonian.
• J(x) ∈ R N ×N is a skew-symmetric matrix corresponding to the power exchanges between energy storing elements.
• R(x) ∈ R N ×N is a symmetric positive semidefinite matrix corresponding to the power exchanges between energy storing and energy dissipating elements.
• G(x) ∈ R N ×M is an input matrix.
Dots above variables, as in ẋ, denote time derivatives, and ∇ represents the gradient with respect to x.Using the chain rule, the skewsymmetry of J(x), and the positive semidefiniteness of R(x) we have an energy balance If H(x) is a positive definite function, then the system is passive and is thus zero-input stable [20].
The storage function is usually the total physical energy in the system, as can be seen in the examples in [16] and in [19,21,22].However, the storage function being a true physical energy is not necessary.As long as H(x) is positive definite, passivity and thus zero-input stability can be guaranteed.
For brevity, we denote S(x) = J(x) − R(x).Also, in this work, we simply show zero-input stability directly and ignore the output flows, y, for the analyses since (3) will then hold when y is defined as in (1b).

Discrete Gradient Scheme
The pH system may be discretized while preserving passivity by following the discrete gradient method outlined in [19] and [18].The key is to ensure a discrete version of the chain rule is satisfied by the discrete gradient, i.e., For a quadratic storage function, H(x) = 1 2 x T Px, with P positive definite, the choice of can be verified to satisfy the discrete chain rule.For a separable storage function H(x) = i Hi(xi), the choice of ∇ d H(x, δx) such that satisfies a discrete chain rule.For general storage functions, a definition of the discrete gradient can be found in [18].The resulting discretization of the scheme with time step k and time index n is This ensures that an energy balance is satisfied under zero-input conditions, so that the discretization is zero-input stable.This scheme is in general implicit and will require an iterative method such as, e.g., Newton-Raphson in order to solve for δx n at each time step.If the storage function is quadratic, H(x) = 1 2 x T Px, with positive definite P then the discrete gradient is (5) and the discrete gradient rule yields a scheme that requires a single inversion of an n × n matrix at each time step to calculate the update, δx t .The scheme above is first-order accurate.[18] describes a two-stage scheme that is second-order accurate.For the sake of brevity, this work will only show firstorder accurate schemes, but they can be easily made second order accurate by using the scheme in [18].

Explicit Scheme
To avoid the need to employ an iterative method such as Newton-Raphson to compute the state update when the storage function is not quadratic, we perform a change of variables so that the storage function is quadratic in this transformed state space.In [18], a variable change is proposed for pH systems with storage function satisfying For separable storage functions, H(x) = i Hi(xi), where H(x) is, as always, positive definite so that Hi(0) = 0, this change of state is given by For general storage functions, [18] provides a suitable change of variables as well.This transformation yields a pH system in x, where, with • denotes function composition and J f (x) is the Jacobian of f in the following convention, This finally yields the discretization The conversion back to the original state variables can be done using x = g(x).This scheme is only first-order accurate but the two-stage form in [18] can achieve second-order accuracy.

Circuit Analysis
A reduced form of the OTA-based MS20 filter circuit is shown in Figure 1.RC stages followed by buffers are realized using OTAs.
A control circuit is able to change R in order to vary the cutoff frequency, ω = 1/RC.The buffer of the second stage of filter is combined with a voltage divider of gain K, and used to change the resonance parameter of the filter.Applying Kirchhoff's Current Law (KCL) and the ideal op-amp laws to the input of the op-amp gives where ID pair is the current through the pair of cascaded diodes with the direction shown in Figure 1.Also, due to the non-inverting amplifier used, we must have |vf| ≥ |Kv2|, with equality when the diodes conduct enough to be approximated by a short circuit.Similar to the circuit analyses in [23] and [24], the current due to the forward-biased diode can be considered to dominate the saturation current of the reverse-biased diode.This, along with |vf| ≥ |Kv2| gives where λ = sgn(v2) and sgn(•) is the sign function.Solving explicitly for vf in (17) using the Lambert W function with G = 1 + R2/R1 ≈ 4 as the non-inverting amplifier's gain, applying KCL at the other nodes of the circuit, and scaling the capacitor voltages to obtain a dimensionless system except for time, we get where u is the scaled input voltage, u = vin/Vref, x1 and x2 are the scaled capacitor voltages, xi = vi/Vref with Vref = 3nVT .VT = 0.02585V is the thermal voltage of a diode and n is the emission coefficient of the diode assumed to be approximately 1. ω = 1/RC is the cutoff frequency in radians per second.α = KG ≥ 0 is the resonance parameter of the circuit.η(x2) is the nonlinearity due to the clipping non-inverting amplifier in the feedback path of the cicuit and is given by where λ = sgn(x2) and β = IsatR2/Vref.It can also be shown that x2η(x2) ≥ 0.
The output of the system is the voltage across the second capacitor, vout = v2 = x2Vref

Lyapunov Stability Analysis
It will be convenient to write the dynamical equations of the system in a different form.The system is where γ(x2) = α − η(x2)/x2.To ensure this is well defined for x2 = 0, we let In practice, the division by x2 is ill-conditioned for small x2.A limiting approximation for small x2 is 3αβ/(4 + 4β).It grows asymptotically to 3α/4 as |x2| → ∞.A plot of this function can be seen in Figure 2 A candidate Lyapunov function for this system is which is positive definite.Testing for Lyapunov stability of the zero-input system with this, we have This is non-positive whenever For this to hold unconditionally, we must have which is the condition for the zero-input system to be stable.Stability analysis of the linearization of ( 19) about x = 0 gives the same stability condition.This is also consistent with the small signal analysis in [25] when assuming the clipping diodes do not conduct any current at all (β = 0).If ω and α are not time-varying, the zero-input system is autonomous.Then LaSalle's invariance principle [20] implies the origin of the zero-input system is globally asymptotically stable for the condition (27).Further, if ω and α are time-varying such that ( 27) is satisfied at all times, the system is also uniformly stable [20], since the Lyapunov function H(x) is independent of time.

DAFx.3
Proceedings of the 24 th International Conference on Digital Audio Effects (DAFx20in21), Vienna, Austria, September 8-10, 2021 For the system is dissipative for only a part of the phase space as shown in Figure 3. Simulation for α in this range shows that the system self oscillates, i.e., the state trajectories fall into a limit cycle.This behavior is consistent with that of the actual circuit.Rigorously proving the existence of a limit cycle will require the use of the Poincaré-Bendixson criterion [20].Practically, this would involve finding a closed bounded subset, M ⊂ R 2 , of the phase plane so that the only equilibrium points within it are unstable and all trajectories starting in M stay in M for all future time.The equilibrium point of the system at the origin is unstable for α satisfying (28).The set M is normally chosen as M = {x : H(x) ≤ c}, where c is large enough so that Ḣ(x) ≤ 0 for all x such that H(x) = c.However due to the nature of the regions that the phase plane is split into, as can be seen in Figure 3, such a c does not exist.This is because no matter how large a c is chosen, the curve H(x) = c always goes through the strip of the phase space in which Ḣ(x) ≥ 0 so that Poincaré-Bendixson's criterion fails to show the existence of a limit cycle.A more sophisticated way of searching for a trapping region, M , will be required to be able to use Poincaré-Bendixson's criterion to rigorously prove the existence of a limit cycle.For α ≥ 8, Ḣ(x) ≥ 0 unconditionally and the Lyapunov analysis fails.Simulation shows unbounded growth of the state for these values of α.

Port Hamiltonian Form
The previous conclusions can alternatively be arrived at by rewriting (21) in the form of (1a), x 1 Lyapunov stability analysis would proceed by requiring be non-positive.This happens when γ(x2) − 2 ≤ 0 which is identical to the conditions obtained from the earlier analysis.

Discrete Time Scheme
Discretization of ( 29) can be done directly using the discrete gradient scheme.Since the storage function function is already quadratic, no state change is required.One of the forms of the (first-order accurate) discretization is with A(x) and G as in (21).

Moog 4-pole Ladder Filter
From the circuit analyses done in [5,17,26,27,28], we have the dimensionless version of the system dynamics of the Moog 4-Pole Ladder Filter as Similarly to the identity on the feedback loop in [17], we write where This allows us to represent the system in (32) as Stability analysis of the linearization of the system [29] shows it is stable for 0 ≤ α ≤ √ 2, i.e., 0 ≤ r ≤ 1.

Lyapunov Stability Analysis
A more detailed stability analysis can be found in [30].
The state is first transformed by a diagonal scaling with d unspecified for now.The zero-input system then becomes where Choosing a candidate Lyapunov function [30] where the system can be written in the form in (1a) with zero-input.
For this to be stable, the symmetric part of S(w) must be negative semidefinite.This is guaranteed when Notice that this choice of Lyapunov function recovers the full range of values allowed for the resonance parameter α, in contrast with previous choices of Lyapunov function [17].

Discrete-Time Scheme
The storage function H(w) satisfies the condition ( 9), the procedure in Section 2.2.1 can be used to obtain an explicit scheme that does not require the use of Newton's iterations.The required change of state is where, with g = f −1 , The matrices in these expressions are as follows: where and where and where γ is defined in (34).The limiting values of each of the functions are needed to avoid ill-conditioning due to division by small numbers.This happens when the state variables are small.From these we finally have the discretization

NUMERICAL RESULTS
All simulations were performed in MATLAB at a sampling rate of 44.1 kHz on a Windows machine with an Intel Core i7-1065G7 CPU.Audio examples can be found on this repository 1 .

Korg35 Rev. 2
The results of simulating the circuit for different resonance parameters with zero input are shown in Figures 4, 5, and 6.In all three figures, the energy balance is satisfied with the total energy varying from its constant value by only floating-point precision.
Figure 4 shows the dissipative behavior of the system for α < 2.
The fixed point at the original is an asymptotically stable spiral.Figures 5 and 6 show the self-oscillating behavior when α is increased beyond 2. The system tends towards the limit cycle regardless of whether the initial state is inside or outside the periodic orbit.Although it is expected that unbounded growth would occur for α ≥ 8, ill-conditioning effects appear from α ≥ αmin ≈ 7.7 even though double-precision floats are used for calculations.The Lambert-W function is computed using Fritsch's iterations as in [23] and [31].A MATLAB implementation of Fritsch's iterations can be found on the website accompanying [23].Processing a 12.8 second input audio sample takes about 2.0 seconds.

Moog 4-Pole Ladder Filter
The results of simulating the circuit for different resonance parameters with zero input are shown in Figures 7 and 8.This energy balance can be seen to be satisfied with the total energy varying from its constant value by only floating-point precision.Figure 7 shows the dissipative behavior of the system for r < 1.When r is increased beyond 1, the circuit appears to self-oscillate in a periodic manner as in Figure 8. Poincaré-Bendixson's criterion does not apply to systems of dimension other than 2.
Processing a 12.8 second input audio sample takes about 2.8 seconds.Discretizing (41) without a state change and using Newton-Raphson iterations at each time step yields a processing time of about 7.8 seconds for the same 12.8 second input.Iterations are stopped after either 10 iterations are complete, or the relative error is 10 −8 .

CONCLUDING REMARKS
This paper is a study of the application of using a Port Hamiltonian framework, alongside a discrete gradient scheme and a state change to the modeling of the classic Korg35 (Rev.2) and Moog 4-Pole Ladder filters.This allows guarantees on zero-input stability of the discretization.The simulations display behavior, such as self-oscillation, consistent with that of the actual circuits.The use of the state change allows for a discrete-time scheme which does not require an iterative method to calculate the state updates.Instead it requires the solution of a linear system of size equal to the dimension of the system only once at each time step.This results in a considerable savings compared with schemes requiring iterative solvers, and ensures existence and uniqueness of computed solutions by construction.
Future work on the Moog 4-pole filter would involve exten-Figure 6: Output, phase portrait, energy variation and relative error in total energy for the Korg35 for a cutoff frequency of 10 Hz and resonance, α = 2.5.Initial state is the circular point in the phase portrait.
sions to the more realistic case of time-varying parameters.Due to the dependence of the storage function of the Moog filter on the resonance parameter, the arguments used to prove time-varying stability for the Korg35 filter can not be extended to this case.For uniform stability [20], positive definite functions W1(w) and W2(w), both independent of time, t, and therefore independent of the resonance parameter, must be found such that W1(w) ≤ H(w, α(t)) ≤ W2(w) (54a) ∂H ∂t (w, α(t)) + ∇H(w, α(t)) T F(w, α(t)) ≤ 0, where ∇H is the gradient with respect to w, and ẇ = F(w; α(t)) is the state dynamics of the system.Bounded-Input-Bounded-Output (BIBO) stability in the presence of input u must be examined.This may be approached using the concept of input-to-state stability [20].Aliasing behavior of these simulations has not been addressed here, and methods to suppress them also needs to be examined.We refer the reader to recent work on antialiasing methods in virtual analog [32,33,34,35], including in the context of port-Hamiltonian methods [36].
Copyright: © 2021 Mohammed Danish et al.This is an open-access article distributed under the terms of the Creative Commons Attribution 3.0 Unported License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Figure 3 :
Figure 3: Regions of phase space where the Korg35 system necessarily dissipates.Also shown is the velocity field of the system for x1 large compared to x2. ±x * 2 are the values of x2 for which (26) is an equality.

DAFx. 4
Proceedings of the 24 th International Conference on Digital Audio Effects (DAFx20in21), Vienna, Austria, September 8-10, 2021 where u = vin/Vref is the scaled input voltage, xi = vi/Vref, i = 1, 2, 3, 4 are the scaled capacitor voltages, Vref = 2VT , with VT = 0.02585V being the thermal voltage of a diode as before.ω is the cutoff frequency of the filter in radians per second.α 4 = 4r where r is the resonance parameter of the filter.The output of the filter is the voltage, v4 = x4Vref, across the fourth capacitor in the circuit.

Figure 4 :
Figure 4: Output, phase portrait, energy variation and relative error in total energy for the Korg35 for a cutoff frequency of 10 Hz and resonance, α = 1.9.Initial state is the circular point in the phase portrait.

Figure 5 :
Figure 5: Output, phase portrait, energy variation and relative error in total energy for the Korg35 for a cutoff frequency of 10 Hz and resonance, α = 2.5.Initial state is the circular point in the phase portrait.