Passive Guaranteed Simulation of Analog Audio Circuits: A Port-Hamiltonian Approach

: We present a method that generates passive-guaranteed stable simulations of analog audio circuits from electronic schematics for real-time issues. On one hand, this method is based on a continuous-time power-balanced state-space representation structured into its energy-storing parts, dissipative parts, and external sources. On the other hand, a numerical scheme is especially designed to preserve this structure and the power balance. These state-space structures deﬁne the class of port-Hamiltonian systems. The derivation of this structured system associated with the electronic circuit is achieved by an automated analysis of the interconnection network combined with a dictionary of models for each elementary component. The numerical scheme is based on the combination of ﬁnite differences applied on the state (with respect to the time variable) and on the total energy (with respect to the state). This combination provides a discrete-time version of the power balance. This set of algorithms is valid for both the linear and nonlinear case. Finally, three applications of increasing complexities are given: a diode clipper, a common-emitter bipolar-junction transistor ampliﬁer, and a wah pedal. The results are compared to ofﬂine simulations obtained from a popular circuit simulator.


Introduction
The characteristic input-to-output behavior of analog audio circuits (timbre, transitory) rests on the possibly highly nonlinear components appearing in such systems. These components make the stability of the simulations difficult to guarantee. The motivation of this work stems from the following observations: 1. Analog circuits combine energy-storing components, dissipative components, and sources. 2. Storage components do not produce energy, and dissipative components decrease it.
In this sense, analog circuits can be considered as passive systems with external power supply. We shall exploit this passivity property by transposing it to the digital domain, ensuring the stability of the simulations (see [1][2][3]).
The available approaches for the automated derivation of physical modeling and numerical simulation of audio circuits can be divided in two classes [4]: wave scattering methods (WS) and Kirchhoff's variables methods (KV). Mixed WS/KV methods have also been proposed in references [5,6]. The well-established wave-digital filter (WDF) formalism [7] belongs to the class of WS methods. For linear circuits, it provides a computationally realizable system of equations: First, by defining parametric wave variables for each elementary component and multiports (serial and parallel); The port-Hamiltonian approach is used to decompose such open physical systems in (i) a set of components that are combined according to (ii) a conservative interconnection network. These two ingredients are detailed below in the case of electronic circuits.

Components
Electronic circuit components are sorted as (or can be a combination of): n S internal components that store energy E ≥ 0 (capacitors or inductors), n D internal components that dissipate power D ≥ 0 (resistors, diodes, transistors, etc.), n P external ports that convey power S (∈ R) from sources (voltage or current generators) or any external system (active, dissipative, or mixed).
The behavior of each component is described by a relation between two power variables: the current i and the voltage v, defined in receiver convention (the received power is P = v · i).
The energy E s stored in storage component s ∈ [1 · · · n S ] is expressed as a storage function h s of an appropriate state x s : E s (t) = h s x s (t) ≥ 0. Typically, for a linear capacitor with capacitance C, the state can be the charge x = q and the positive definite function is h(q) = q 2 /(2C). Storage power variables (v s , i s ) are related to the variation of the state dx s dt and the gradient of the storage function h s (x s ), the product of which is precisely the received power: v s · i s = dE s dt = h s · dx s dt . For the capacitance, these constitutive laws are i = dq dt = dx dt and v = q/C = h . Note that these definitions apply equally for non-quadratic storage functions h(x) ≥ 0 for which h (x) is not constant.
The power D d instantaneously dissipated by the dissipative component d ∈ [1 · · · n D ] is expressed with respect to an appropriate dissipation variable w d : D d (t) ≡ D d w d (t) ≥ 0. Typically, for a linear resistance R, w can be a current w = i and D(i) = R · i 2 . As for storage components, a mapping of the dissipative power variables (v d , i d ) is provided, based on the factorization D d (w d ) = w d · z d (w d ), introducing a dissipation function z d . For the resistance, i = w and v = R · i = z(w).
The power instantaneously provided to the system through external port p ∈ [1 · · · n P ] is S p (t), and we arrange the source variables (v p , i p ) in two vectors: one is considered as an input u p , and the other as the associated output y p , so that the power received from sources on port p is S p = y p · u p = −v p · i p (receiver convention, with v p · i p the power received by the sources).

Conservative Interconnection
The interconnection of the components is achieved by relating all the voltages and currents through the application of the Kirchhoff's laws to the interconnection network (schematic). This defines a conservative interconnection, according to Tellegen's theorem recalled below (see also [26] and [27] §9.4).
Theorem 1 (Tellegen). Consider an electronic circuit made of N edges defined in same convention (here receiver), with individual voltages v = (v 1 , · · · , v N ) and currents i n = (i 1 , · · · , i N ) which comply with the Kirchhoff's laws. Then v · i = 0.
A direct consequence of (2) is that no power is created nor lost in the structure: v · i = ∑ N i=1 P n = 0, with P n = v n · i n the power received by edge n, thus defining a conservative interconnection (Tellegen's theorem is a special case of a more general interconnection structure, namely, the Dirac structure (see [23] §2.1.2 for details)). Now, denote (v s , i s ), (v d , i d ), and (v p , i p ) the sets of all the power variables associated with storage components, dissipative components, and sources (respectively), and v = (v s , v d , v p ) , i = (i s , i d , i p ) the vectors of all the power variables. Then, Tellegen's theorem restores the power balance (1) with where ∇H : R n S → R n S denotes the gradient of the total energy E = H(x) = ∑ n S s=1 h s (x s ) with respect to (w.r.t.) the vector of the states [x] s = x s , and function z : R n D → R n D denotes the collection of The above description of storage components, dissipative components, and source, along with the conservative interconnection stated by the Kirchhoff's laws, constitute the minimal definition of a port-Hamiltonian system (PHS) (see [23] §2.2). In this work, we focus on circuits that admit an explicit realization of PHS, for which the quantities b = (b 1 , · · · , b N ) = ( dx dt , w, −y) (with b n = v n or b n = i n ) can be expressed as linear combinations of the remaining N powers variables organized in the dual vector a = (a 1 , · · · , a N ) = (∇H(x), z(w), u) (with a n = i n if b n = v n or a n = v n if b n = i n ): Then, a · b = a · J · a = 0 from Tellegen's theorem, so that the matrix J is necessarily skew-symmetric (J = −J). More precisely, we consider the following algebraic-differential system of equations  where matrices J x , J w , J y are skew-symmetric. The significance of the structure matrices is the following: J x ∈ R n S ×n S expresses the conservative power exchanges between storage components (this corresponds to the so-called J matrix in classical Hamiltonian systems); J w ∈ R n D ×n D expresses the conservative power exchanges between dissipative components; J y ∈ R n P ×n P expresses the conservative power exchanges between ports (direct connections of inputs to outputs); K ∈ R n S ×n D expresses the conservative power exchanges between the storage components and the dissipative components; G x ∈ R n S ×n P expresses the conservative power exchanges between ports and storage components (input gain matrix); G w ∈ R n D ×n P expresses the conservative power exchanges between ports and dissipative components (input gain matrix).
Property 1 (Power Balance). The variation of the total energy E = H x of a system governed by (5) is given by (1), with D = z(w) · w ≥ 0 the total dissipated power, and S = u · y the total power incoming on external ports.

Proof.
We have a · b = dE dt + D − S. Now a · b = a · J · a = 0 since J is skew-symmetric.

Remark 1 (Power variables).
This work is devoted to the treatment of electronic circuits for which power variables are chosen as current and voltage. However, all the aforementioned definitions apply equally to multiphysical systems, provided an adapted set of power variables, generically denoted by flux (currents, velocities, magnetic flux variations) and efforts (voltages, forces, magnetomotive force), the product of which is a power (see [23] Table 1.1). This follows the bond-graph modeling approach [28,29], on which the PHS formalism is built (see [23] §1.6 and 2.1). The treatment of multiphysical audio systems in the PHS formalism can be found in [30] (electromechanical piano that includes mechanical, electrical, and magnetic phenomena) and [31] ( §III.B) (modulated air flow for musical acoustics applications that includes mechanical and acoustical phenomena).

Example
Consider the resistor-inductor-capacitor (RLC) circuit in Figure 1, with n S = 2, n D = 1, and n P = 2, described as follows. For the linear inductance L, the state and the positive definite function can be the magnetic flux x 1 = φ and h 1 (φ) = φ 2 /(2L), so that v L = dh 1 /dx 1    Applying Kirchhoff's laws to this simple serial circuit yields From the constitutive laws of components, this equation restores the form (5) exactly, block by block. It provides the algebraic-differential equations that govern the system with input u and output y.
This work aims at simulating such passive systems by firstly generating Equation (5) associated to a given circuit, and secondly by deriving its numerical version so that a discrete power balance is satisfied.
Remark 2 (Reduction). The system (5) can be reduced by decomposing function z into its linear and nonlinear parts. See Appendix A for details.

Generation of Equations
This section provides a method to translate the description of a circuit (components and interconnections) from a netlist in a Spice-style [32] to the Formulation (5). Compared to standard methods that express all the currents as a function of all the voltages (see [18,27,32]), Formulation (5) expresses vector b of selected power variables (voltage or current) as a function of the vector a of complementary power variables (if [b] n is a voltage of a branch, [a] n is the associated current with receiver convention). To derive the matrix J that relates the voltages and the currents arranged in vectors a and b according to Kirchhoff's laws (as in example §2.2), we propose a two-step method: Step 1 : from a netlist (L) to a graph (G) that represents the Kirchhoff's laws for a chosen orientation (convention); Step 2 : from (G) to the skew-symmetric matrix J in (5).
Step 1 is standard. The presentation focuses on convention choices and details our procedure. In step 2, we propose an algorithm that analyzes if Formulation (2) is available (that is, the circuit is realizable into the PHS formalism) and delivers the matrix J in this case. Otherwise, the circuit corresponds to an implicit formulation that is not addressed in this paper. In practice, such cases appear for serial(/parallel) connection of voltage(/current)-controlled components. In this case, port-Hamiltonian Formulation (5) requires extension (see [22,23]).

Netlists
Each line of a netlist describes an element of the corresponding schematic, with: identification label, list of connection nodes, type of element, and list of parameters. We divide netlists into two blocks: internal components (dissipative and storage) and external ports (supplies and ground). In the first block (components), each line includes a reference to the appropriate entry in the dictionary and a list of the parameters for the corresponding model. Each line of the second block (external ports) provides the label of the externalized node, the type of supply (voltage or current), and the symbol∼if the supply is modulated (typically, the input signal), or a value if constant (typically, a battery).
As an example, the netlist corresponding to the circuit in Figure 2 is given in Table 1. Here, the components are given lines 1 to 3 (gray). The first two lines describe dipoles: a linear capacitor between N 1 and N 2 with label C1 and capacitance value 20e −9 F; and a resistor between N 3 and N 4 with label R1 and resistance value 1.5e 3 Ω. The third line describes a npn bipolar-junction transistor. From the dictionary (Appendix B), the base terminal appears to be connected to the circuit's node N 2 , the emitter terminal to N 3 , and the collector terminal to N 5 . For this component, the list of parameters is: forward and reverse common emitter current gain, reverse saturation current, and thermal voltage. External ports are given lines 4 to 7 . Line 4 describes a constant 9 V voltage supply (labeled Vcc) on the circuit's node N 4 . 5 describes a modulated voltage supply (here considered as the input signal) on N 4 . 6 describes a constant 0 A current supply on node N 3 ; this permits the recovery of the voltage on that node N 3 as an output to the circuit. 7 describes the connection of the circuit's node N 5 to the ground.  Table 1). Table 1. Example of a netlist corresponding to the circuit in Figure 2. The grey part corresponds to the components, and the other elements correspond to the external ports, or sources (as in figure 2).

. Graph
A graph G = {N, B} is defined by two lists of nodes N (also called vertices) and branches B (also called edges), with B ⊂ N 2 (each element of B is an object defined on two elements of N, see [33] for details). The dictionary (Appendix B) encodes the graph of each elementary component. The branches of such an elementary graph contain the constitutive laws of the corresponding component: • Dipoles are made of two nodes and a single branch, defining a single couple of state x and storage function h(x) (storage component), or dissipative variable w and scalar relation z(w) (dissipative components). • More generally, n-ports multipole are made of n nodes and at least n − 1 branches, defining n − 1 couples of variables and functions. Typically, the graph for the bipolar junction is made of two branches (base-emitter and base-collector).
The graph corresponding to a given circuit is derived from its netlist description in two steps: 1. build the internal graph by connecting the elementary graph of the components from the first block of the netlist, 2. introduce a reference node N 0 (or datum, see [27] §10) to define the external branches from the second block.
Typically, N 0 corresponds to the ground or any local electrostatic potential which does not impact the currents nor voltages. Then, N is built from the list of nodes appearing at least once in the netlist, plus the reference node N [N 0 , N 1 , · · · , N n N ]. According to Section 2, the set of branches is organized as B = {B S , B D , B P }, with B S the n S energy storage branches, B D the n D dissipative branches, and B P the n P sources.
As an example, the construction of the graph in Figure 2 from its netlist 1 is as follows. Firstly, the internal graph is built. It is made of n N = 5 nodes {N 1 , · · · , N 5 } and four branches B S = {C 1 } and B D = {R 1 , Q 1,bc , Q 1,be }. Secondly, we introduce the (virtual) reference node N 0 to define the four branches corresponding to the external ports B P = {IN, Vcc, OUT, GRD}.

Kirchhoff's Laws on Graphs
We assign to each branch b both a voltage v b and a current i b in receiver convention, the direction of the branch indicating the direction of the current. Note that the power supplied to the system on port p is the power emitted by the port branch S p = u p · y p = −v p · i p . For a circuit made of n N + 1 nodes and n B = n S + n D + n P branches, we define: the set of electrostatic potentials on the nodes e = (e 1 · · · e n N ) , the set of voltages v = (v 1 · · · v n B ) , and the set of currents i = (i 1 · · · i n B ) . The orientation of an entire graph is encoded in its incidence matrix Γ ∈ R (n N +1)×n B , defined below [27] ( §9).
As an example, the incidence matrix for the circuit described in Table 1 is given equation below (0 are replaced by dots). Notice the grey columns correspond to the components, and the other columns correspond to the external ports, or sources (as in Table 1 and Figure 2).
Since the reference potential e 0 does not influence the voltages nor the currents, it is not taken into account in the Kirchhoff's laws, and we define the reduced incidence matrix Γ ∈ R n N ×n B obtained by deleting the row corresponding to the datum N 0 in Γ. This leads to the following matrix formulation of Kirchhoff's Voltage Law (KVL) and Kirchhoff's Current Law (KCL) [27] ( §10), from which the structure (5) is derived.

Realizability Analysis
The PHS structure (5) relies on (i) an arrangement of currents i and voltages v in two vectors a and b and (ii) a set of linear relations encoded in the skew-symmetric matrix J that corresponds to the conservation laws (8) applied on (i, v). For storage and sources components, step (i) is straightforward with the constraints given in Table 2. For dissipative components, this step is achieved by selecting each component as voltage-controlled or current-controlled in order to satisfy a criterion on the matrix description of the interconnection scheme. This realizability criterion is given in Section 3.2.1, assuming the control type of every edge is known. A method of choosing the control type of dissipative edges so as to satisfy the realizability criterion is addressed in Section 3.2.2. This leads to Algorithm 1, which solves (i) and (ii).

A Criterion for Realizability
In Section 3.1, the set of edges B has been partitioned based on the differentiation between internal edges (or component edges {B S , B D }, grey) and external edges (or ports edges B P , white), in order to build the complete graph from the netlist. In this section, we are interested in the PHS Formulation (5) associated to a given complete graph. To that end, we leave out from here the differentiation between internal and external edges in order to focus on the differentiation between voltage-controlled and current-controlled edges.
Suppose the control type of every edge is known, and the set of edges is split according to B = {B 1 , B 2 }, with B 1 the set of n 1 voltage-controlled edges and B 2 the set of n 2 current-controlled edges (see Table 2). Correspondingly, the sets of power variables are split as v = (v 1 , v 2 ) and i = (i 1 , i 2 ) , and we defineã = (v 1 , i 2 ) andb = (i 1 , v 2 ) . Since the reference potential e 0 defined on node N 0 does not influence the voltages nor the currents, it is not considered in the sequel, and the incidence matrix splits as follows: This leads to a rewrite of the Kirchhoff laws (8) as: Proposition 1 (Realizability). If γ 2 is invertible, then the port-Hamiltonian structure (5) provides a realization of the graph Proof. From the relation on the voltages in (9), we get v 1 = γ 1 e and v 2 = γ 2 e. From the relation on the currents (10), we get The PHS (5) is obtained by rearranging the edges according to their role with respect to the power balance, according to the permutation of vector elements Π(ã) = dx dt , w, y = a (and correspondingly Π(b) = (∇H, z, u) = b), which is also applied on rows and columns ofJ to yield a = J b.
From the invertibility condition on γ 2 in Proposition 1, we state the following remark, which is used in the sequel to derive the realizability analysis algorithm.

Remark 3 (Necessary condition for realizability).
A necessary condition for the graph G to be realizable as a PHS (5) is that it includes as many current-controlled edges as nodes n N , with γ 2 ∈ R n N ×n N .

Algorithm
This section introduces an algorithm that selects the appropriate control type for each dissipative edge so that the partition B = {B 1 , B 2 } satisfies Proposition 1. From Remark 3, the total number n 2 of current-controlled edges should be exactly equal to the number of nodes n N . From the special structure of the incidence matrix Γ, this in turn ensures that the potential on each node is uniquely defined by a linear combination of the voltages v 2 (elements in a associated to current-controlled edges); i.e., γ 2 is invertible.
Consider the current-controlled edge b from node i to node j in Figure 3. If the potential on node j is known, the remaining potential is obtained from e i = v b ± e j , where the sign depends on the orientation. In this case, we say edge b imposes the potential on node i. Now, the objective is to perform this analysis globally so that e is a linear combination of v 2 with e = (γ 2 ) −1 · v 2 . To that end, we introduce the realizability matrix Λ defined element-wise as follows Then, a given graph is realizable if the type of each resistor can be selected so that the following set of constraints is fulfilled.
(C1) The potential on each node n ∈ [1, · · · , n N ] is uniquely defined so that ∑ Constraints (C1-C2) ensure that γ 2 is invertible, so that e = (γ 2 ) −1 · v 2 . Constraint (C3) ensures that the reference potential on datum does not contribute to the system's dynamics. Constraint (C4) ensures that inputs of voltages-controlled edges are explicitly given by a linear combination of the nodes potentials so that v 1 = γ 1 e = (γ −1 2 · γ 1 ) · v 2 . To build and analyze the matrix Λ, we start from the adjacency matrix A of the graph, defined as follows: Then the non-zero elements in A are analyzed to cope with the realizability constraints (C1-C4). This yields Algorithm 1. The PHS (5) is finally recovered as discussed in the proof of Property 1.
Algorithm 1: Analysis of realizability. If successfully complete, the resulting PHS structure is given by the procedure in the proof of Proposition 1 Input: A graph G = {N, (B x , B w , B y )} corresponding to the interconnection of storage, dissipative, and source edges. Output: Sets of voltage-controlled edges B 1 and current-controlled edges B 2 .

Example
As an example, the realizability analysis for the system in Figure 2 with the choice of inputs/outputs in Table 1 is as follows. In step 1, the realizability matrix Λ is initialized with the adjacency matrix A, which is built by taking the absolute value of incidence matrix [A] n,b = abs [Γ] n,b : In The realizability matrix after step 11 in Algorithm 1 is After step 19, the algorithm concludes that the potential on node N 1 is imposed by edge B IN so that the potential on node N 2 is imposed by the capacitor B C1 . After step 28, the algorithm concludes that the potential on node N 4 is imposed by the edge B Vcc so that the resistor is current-controlled (so as to impose the potential on node N 3 ).
This concludes the realizability analysis. To recover the associated port-Hamiltonian structure, we return to the incidence matrix Γ. With the new edges ordering B = {B 1 , B 2 } prescribed by the above analysis, it is rewritten as Finally, the structure (5) is recovered by computing the matrix γ = γ −1 2 · γ 1 in (11) with

Guaranteed-Passive Simulation
This section is devoted to the discrete-time simulation of the algebraic-differential system (5); that is, the computation of x(k) ≡ x(k · T) from u(k) ≡ u(k · T), with k ∈ N, for the constant sampling frequency f s = 1/T. First, we present the design of a numerical scheme that properly transposes the power balance (1) to the discrete time domain: this choice makes the passivity property preserved, from which stability issues stem. Second, a numerical method is used to solve the implicit equations due to the numerical scheme (on x) and the algebraic equations (on w).

Numerical Scheme
To ensure the stable simulation of stable dynamical system dx dt = f(x), many numerical schemes focus on the approximation quality of the time derivative (or integration), combined with operation of the vector field f. Here, we adopt an alternate point of view, by transposing the power balance (1) into the discrete time-domain to preserve passivity. This is achieved by numerical schemes that provide a discrete version of the chain rule for computing the derivative of the composite function E = H(x). This is the case of the forward difference scheme, for which first order approximation of the differential applications dx(t, dt) = dx dt (t) · dt and dH(x, dx) = ∇H(x) · dx on the sample grid t ≡ kT, k ∈ Z are given by where, for mono-variate energy storing components (H(x) = ∑ n S n=1 h n (x n )), the n-th coordinate is given by A discrete chain rule is indeed recovered so that the following substitution in (5) dx Remark 4 (Multi-variate components). The case of mono-variate energy storing components covers most of the applications in electronics. Additionally, a generalization of the discrete gradient for multi-variate Hamiltonians such that Equations (20) and (21) are satisfied is given in Appendix C.
In this paper, we consider the class of the PHS composed of a collection of linear energy storing components, with quadratic Hamiltonian h n (x n ) = x 2 n 2C n (C n is a capacitance or an inductance and we define Q = diag(C 1 · · · C n S ) −1 ). Then the discrete gradient (22) reads which restores the midpoint rule that coincides in this case with the trapezoidal rule. For nonlinear cases, (22) leads to another numerical scheme depending on the nonlinearity, still preserving passivity (see (25) and §4.3).

Solving the Implicit Equations
Injecting the numerical scheme (26) in (5) and solving for the quantity δx(k) = x(k + 1) − x(k) leads to the following energy-preserving numerical system: where matrices are related to J in (5) as follows.
Given u(k), the solution of (27) is obtained from the solution of the static nonlinear implicit function f(w(k)) = p x(k), u(k) , with Remark 5 (Explicit mapping). From the global inverse function theorem (see [34]), there exists an explicit mapping w(k) = f −1 p x(k), u(k) provided the Jacobian matrix J f (w(k)) = I d − (J w − 1 2 K · Q · D · K) · J z (w(k)) is invertible for all w(k), connecting the proposed method to the K method [18,19]. This is true since Q, D and the Jacobian of z (for the components of the dictionary in Table B1) prove positive definite, and J w is skew-symmetric.
In this paper, we use the Newton-Raphson algorithm, which iteratively approximates the nearest root of function r : w(k) ∈ R n D → r w(k) ∈ R n D with the following update rule: w n+1 (k) = w n (k) − J r (w n (k)) −1 . r w n (k) , where J r (w) is the Jacobian matrix of r w(k) = f w(k) − p(k). Once a solution w(k) to the implicit equation is available, the output and state updates are given by: Finally, denoting by n t the number of time-steps and n NR the number of Newton-Raphson iterations per time-step, the simulation is performed according to Algorithm 2.
Algorithm 2: Simulation, with n t the number of time-steps and n NR the (fixed) number of Newton-Raphson iterations.

Comparison with Standard Methods
In this section, the proposed approach (PHS structure combined with the discrete gradient method) is compared with two standard methods: the trapezoidal rule (average of the vector field at x(k) and x(k + 1), used in the WDF approach [7]) and the midpoint rule (evaluation of the vector field at , suitable for any differential-algebraic system of equations). Both are known to preserve the passivity of linear undamped systems (see [35] for a detailed analysis). The updates associated with these three methods are given in Table 3. These methods are applied on the same conservative system dx dt = J x · ∇H(x), the power balance of which is given by dE dt = 0 (with D = S = 0). The comparison measure is then the relative error on energy ε(k) = for k ≥ 0. Table 3. Updates for the three methods considered in §4.3. PHS stands for port-Hamiltonian system.

Method Update
Trapezoidal rule Midpoint rule First, notice that for quadratic Hamiltonian H(x) = x ·Q·x 2 with linear gradient ∇H(x) = Q · x, the three methods yield the same update: As a consequence, these three methods induce the same frequency warping (see [36] for the analysis of the bilinear transform derived from the trapezoidal rule).
We focus on the nonlinear case. For comparison, we choose a simple nonlinear conservative system with state x = (x 1 , x 2 ) , non-quadratic Hamiltonian and canonical skew-symmetric matrix J x = 0 −1 1 0 .     Table 3, for a nonlinear conservative system dx dt = J x · ∇H(x) with H(x) given in (30) . We see from Figure 4d that the error associated with the proposed method (PHS approach combined with the discrete gradient method) is low compared to the two other methods (with machine precision 10 −16 ). The accumulation of these errors is responsible for the apparently unstable behavior of the trapezoidal rule.
In each case, the resulting implicit equations are solved by Python iterative solver (see [37]). The Python code is available at the url given in [38]. In order to exhibit the behavior of each method in the worst case, simulations are performed with an especially low sample rate of f s = 10 Hz. The results for each method are given in Figure 4, with comparison in Figure 4d. We see that the error of the proposed method is low (close to machine precision 10 −16 ) compared to standard methods.

Applications
This section is devoted to the simulation of three analog audio circuits by the application of Algorithms 1 and 2. Those circuits are a diode clipper, a common-emitter BJT audio amplifier, and a wah-pedal as a full device. Results obtained with (i) the method in Section 4 and (ii) with the offline circuit simulator LT-Spice [32] are compared.

Diode Clipper
Diode clipper circuits can be found in several audio-distortion devices. They are made of one resistor and two diodes (n S = 0, n D = 3) connected to the ground in reversed bias (see Figure 5a). The external ports are the input/output and the ground (n P = 3). The resistor is current-controlled and the ground is removed. The vectors (a, b) and the structure J returned by Algorithm 1 are: The simulation is performed according to Algorithm 2 at the sample rate f s = 96 kHz, with three Newton-Raphson iterations (shown to be enough to converge in practice). We apply a linearly increasing 1 kHz sinusoidal excitation u IN = −v IN during 10 ms with maximum amplitude 2 V (i OUT = 0 A, v GRD = 0 V). The output y OUT = −v OUT is given in Figure 5b. We see the signal is clamped between ±0.6 V, in accordance with LT-Spice results.

Common-Emitter BJT Audio Amplifier
Common-emitter bipolar-junction transistor (BJT) amplifiers are widely used as amplification stages in analog audio processing. They are made of two capacitors (n S = 2), two resistors, and one NPN transistor which is made of two nonlinear dissipative branches (n D = 4, see Figure 6a and the dictionary in Table B1). The external ports are the input/output signals, the 9 V supply, and the ground (n P = 4). Note that the ground is removed. The resistor Rc is current-controlled, and the resistor Rf is voltage-controlled. The vectors (a, b) and the structure J returned by Algorithm 1 are given in Equations (31) and (32).
The system is reduced according to Appendix A, and the simulation is performed according to Algorithm 2 at the sample rate f s = 384 kHz, with 10 Newton-Raphson iterations. The reason for increasing the sample-rate and the number of Newton-Raphson iterations is twofold. Firstly, it attenuates the effect of aliasing for input signals limited to the audio range (see [39] for details). Secondly, it ensures that the iterative solver converges, which is difficult due to the numerical stiffness of the problem; that is, the Lipschitz constant associated to the inital Cauchy problem is very high, see [40]. At first, we turn the supply v VCC = −9 V on, and we wait 0.3 s for the system to reach its steady state. Then, we apply a 10 ms sinusoidal excitation u IN = −v IN at 1 kHz with linearly increasing amplitude between 0 V and 0.2 V (i OUT = 0 A). The resulting output y OUT = −v OUT is given in Figure 6b. We see that the signal is amplified between 0 V and 9 V, with a strong asymmetrical saturation, in accordance with LT-Spice results. Additionally, spectrograms obtained for an exponential chirp on the audio range are given in Figure 7 (see Figure 16 in [39] for comparison).

Wah Pedal
This section addresses the simulation of a full device (namely the Dunlop Cry-Baby wah pedal) to be used in real time. The circuit is given in Figure 9. It provides a continuously varying characteristic wah filtering of the input signal. This circuit has been treated with the nodal discrete K-method in [41] and with the PHS framework in [42]. It is composed of n S = 7 storage branches (6 capacitors and 1 inductor), n D = 18 dissipative branches (11 resistors, 1 PN diode, 2 NPN transistors and a potentiometer), and n P = 3 ports (input/output signals and battery, discarding the 5 grounds). The wah parameter is the potentiometer's coefficient α. Notice this circuit includes several edges that do not contribute to the device input-to-output behavior, as analyzed in [41]. In this work, we consider the complete original schematic. From Algorithm 1, the resistors R 1 , R 6 · · · R 9 and R 11 are considered as conductances, and the others as resistances. The structure J is not shown here. The sets of PHS variables are:ẋ where w R is the set of dissipative states and z R the set of characteristics according to each resistor's type, and The system is reduced according to Appendix A (with potentiometer's time varying resitors kept in w, z). Firstly, we realize an offline simulation (in Python) with Algorithm 2 for the sampling rate f s = 96 kHz, and three Newton-Raphson iterations. We apply a white noise normalised to 1 V on the input u IN = −v IN (i OUT = 0 A). The magnitudes of transfer functions obtained from fast Fourier transform are given in Figure 8 for the two extreme positions of the pedal. These results are in accordance with LT-Spice.  Secondly, a VST plugin [43] to simulate the Cry-Baby in real-time is made from Algorithm 2. First, a C++ code is automatically generated; Second, this code is encapsulated in a Juce template to compile the audio plugin (see [44]). The sample rate f s is imposed by the host digital audio workstation (here Ableton Live!), and we force five Newton-Raphson iterations. The simulation performed well (audio examples are available at the url [45]). The CPU load on a laptop (Macbook 2.9 GHz Intel Core i7 with 8Go RAM) is 37% for f s = 96 kHz, and 20% for f s = 48 kHz.
Remark 6 (Time-varying stability). The use of the Newton-Raphson method can hamper the stability of the numerical solution for time-varying systems, especially in the case of fast variations (here, of the potentiometer) for which the Jacobian matrix of the implicit Function (28) can be ill-conditioned. For linear storage components, a solution is to use the K-method instead (see Remark 5).

Conclusions
We have established a method to automatically recast an analog audio circuit into PHS formalism, which guarantees passivity of the continuous time model. The generation of the PHS from a given schematic lies on two points: 1. the graph theory to describe the interconnection network of a given circuit's schematic, 2. a dictionary of elementary components which are conformable with PHS formalism.
Then, we transposed this physical principle to the digital domain by properly defining the discrete gradient of the Hamiltonian, such that a discrete time version of the power balance is satisfied. The resulting stable numerical scheme is of second order (restoring the midpoint rule for linear systems). It has been shown that the K-method is always applicable to PHS (providing efficient implementations of the implicit relation due to the proposed numerical scheme).
Offline simulations are consistent with LT-Spice results. The whole method allows the automatic generation of C++ simulation code to be used in the core of a real-time VST audio plug-insimulating the Dunlop Cry-Baby wah pedal.
A first perspective on this work is to consider higher-order numerical schemes (namely, the class of Runge-Kutta schemes). Moreover, it would be possible to symmetrize the roles of the voltages and the currents at the interconnection by applying the Cayley transform to the PHS structure, thus adopting wave variables, with possible connection with the WDF formalism. Additionally, an automated analysis of the original schematic could be developed, so as to identify the unimportant or degenerate states and to reduce the dimensionality of the system. Finally, it could be possible to exploit the compatibility of the proposed method with the K-method to alleviate the numerical cost due to Newton-Raphson iterations.
with z L a diagonal matrix whose elements are the resistance or conductance of linear dissipative components, and z N a collection of the nonlinear dissipative relations. Correspondingly, the structure matrices are decomposed as Defining L = (K L , −K LN , −G L ), and M = z L −1 − J LL −1 a positive definite matrix, the system (5) withJ a skew-symmetric matrix given bỹ and R is a symmetric positive definite matrix given by Indeed, matrixJ (respectively R) corresponds to the conservative (respectively resistive) interconnection of dynamical storage components, nonlinear dissipative components, and sources.
(z D (w D ) .w D ≥ 0). As in LT-Spice simulators, a minimal conductance G min is added in Table B1. This helps convergence in the simulation process. NPN junctions are passive 3-ports, with dissipated power D Q = v B · i B + v C · i C + v E · i E ≥ 0. Here we use the Ebers-Moll model, which preserves this passivity property. I S is the saturation current, β R and β F are respectively the reverse and forward common emitter current gains, and v t is the thermal voltage. The corresponding voltage-controlled dissipative characteristic z Q (w Q ) = [i BC , i BE ] is given in Table B1, denoting α R = β R +1 β R , α F = β F +1 β F , and including minimal conductances. Note that D Q = z Q (w Q ) .w Q ≥ 0. In our final simulations, some resistors are added to model the resistance of contacts in the nonlinear components, choosing the same values as in LT-Spice models.
Appendix B.4. Incidence Matrices Γ Incidence matrices for 2-ports, potentiometer P, and transistor Q with conventions of Table B1.
For mono-variate components, (22) and (C1) coincide and yield (discrete) constitutive laws that are insensible to the ordering of the state variables.