Adaptive Multi-rate Wavelet Method for Circuit Simulation

In this paper a new adaptive algorithm for multi-rate circuit simulation encountered in the design of RF circuits based on spline wavelets is presented. The circuit ordinary differential equations are first rewritten by a system of (multi-rate) partial differential equations (MPDEs) in order to decouple the different time scales. Second, a semi-discretization by Rothe's method of the MPDEs results in a system of differential algebraic equations (DAEs) with periodic boundary conditions. These boundary value problems are solved by a Galerkin discretization using spline functions. An adaptive spline grid is generated, using spline wavelets for non-uniform grids. Moreover the instantaneous frequency is chosen adaptively to guarantee a smooth envelope resulting in large time steps and therefore high run time efficiency. Numerical tests on circuits exhibiting multi-rate behavior including mixers and PLL conclude the paper.


Introduction
Widely separated time-scales occur in many radiofrequency (RF) circuits such as mixers, oscillators, PLLs, etc., making the analysis with standard numerical methods difficult and costly.Low frequency or baseband signals and high frequency carrier signals often occur in the same circuit, enforcing very small time-steps over a long time-period in the computation of the numerical solution.The occurrence of widely separated time-scales is also referred to as a multi-scale or multi-rate problem.Hence, classical numerical techniques, as transient analysis, result into prohibitively long run-times.
One method to circumvent the bottleneck is to reformulate the ordinary circuit DAEs in a system of partial DAEs (multi-rate PDAE).The method was first presented in [7,14], specializing to multi-rate circuits and systems in steadystate.In [22,28] the technique has been firstly applied to initial transient simulation of driven circuits with a-priori known frequencies.In [8,9,11] a generalization has been derived for circuits with a-priori unknown or time-varying frequencies.These techniques opened the door to multi-rate techniques with frequency modulated signal sources or autonomous circuits such as oscillators with a priori unknown fundamental frequency [9,11,12].
In recent time the expansion of the signal waveforms by a wavelet or spline basis instead of trigonometric basis functions became increasingly popular [5,6,10,16,17,31,32]. Trigonometric basis functions are badly suited for pulseshaped signals due to a slow convergence of the trigonometric series and the Gibbs' phenomenon, whereas spline basis functions with their compact support enable an adaptive mesh with a local refinement and hence smaller system size.It has been shown in [10] that the spectra can be easily calculated from a cubic spline basis which is very important for the acceptance of the approach from the community of circuit designers.
Essential for the success of such methods is that the spline or wavelet expansion can be adapted to the particular signal shape in an circuit simulation example.First steps have been done in [5,6], where an adaptive spline wavelet approach has been developed for transient analysis, i.e., the solution of initial value problems.Here we present an extension of this approach to multi-rate problems, where advantages of the approach can be fully explored.Sect. 2 will give a short introduction to the multi-rate circuit simulation problem.In Sect.3-5 we describe the discretization of the multi-rate circuit equations based on a spline expansion for the semi discretized problem.Sect.6 introduces the scheme for wavelet based adaptive grid generation.Numerical tests in Sect.7 show the performance of the method.

The multi-rate circuit simulation problem
We consider circuit equations in the charge/flux ori-arXiv:1604.07215v1[math.NA] 25 Apr 2016 ented modified nodal analysis (MNA) formulation, which yields a mathematical model in the form of a system of differential-algebraic equations (DAEs): Here x(t) ∈ R n is the vector of node potentials and specific branch currents and q(x) ∈ R n is the vector of charges and fluxes.The vector g(x) ∈ R n comprises static contributions, while s(t) ∈ R n contains the contributions of independent sources.The DAEs in (1) are usually solved by integration formulas for stiff systems.
To separate the different time scales the problem is reformulated as a partial differential-algebraic equation (PDAE) as where ω(τ) is an estimate of the (scaled) angular frequency.
The bivariate function x(τ,t) is related to the univariate solution x(t) of (1) as follows.For any solution x(τ,t) of (2) we get by a solution of Thus, if we choose ŝ such that then the solution of (2) provides also a solution of (1).
Although the formulation (2) is valid for any circuit, it offers a more efficient solution only for certain types of problems.This is the case if x(τ,t) is periodic in t and smooth with respect to τ.Then, a semi-discretization with respect to τ can be done resulting in a relative small number of periodic boundary problems in t for only a few discretization points τ .In the sequel we will consider (2) with periodicity conditions in t, i.e., x(τ,t) = x(τ,t + P) and suitable initial conditions x(0,t) = X 0 (t).Here P is an arbitrary but fixed period length, which results in a scaling of x(τ,t) and ω(τ) as we will see in the next section.Although P > 0 can be chosen arbitrarily only a few choices are of practical interest.These are essentially P = 1 or P = 2π if one wants to work in a fixed periodic setting or one choses P as the period of the carrier so that the scaling of x(τ,t) corresponds more to the physical setting.
As we have shown in [4,Lemma 3.2] the additional (scaled) parameter ω(τ) can be used to adapt to the instantaneous frequency f (τ) of the carrier signal.If the instantaneous frequency is known (e.g. through the sources) one can set ω(τ) = f (τ) P, to obtain optimal smoothness with respect to τ.If the instantaneous frequency is not known in advanced we consider ω(τ) as an additional unknown, which is determined by the additional smoothness condition (see [4] for details).

Semi-discretization
We discretize (2) with respect to τ (Rothe's method) using Gear's BDF1 method of order s and obtain Here we have to determine an approximation X k (t) of x(τ k ,t) from known, approximative solutions X k−i (t) at previous time steps τ k−i .i = 1, . . ., s.As it is well known, the BDF coefficients α k i are chosen such that for any polynomial p of degree up to s the derivative is computed exactly, i.e., Assuming that the τ k−i , ω k−i and X k−i (t) are known and fixed for i > 0 we define Then X k is the solution of the periodic boundary value problem The new problem ( 9) is closely related to the original periodic steady state problem of the circuit, only modified by the additional 'source term' of ( 6).

Multi-rate source term
Except for autonomous systems we have to choose a multi-rate source term ŝ(τ,t), which satisfies for the given source s(t) the relation (5), even if the frequency term ω(τ) is not known in advance.Therefore, we consider first a reference source term for constant ω(τ).We assume that there is a reasonable multi-rate term s(τ,t) which satisfies the relation s(t) = s(t, ωt) i.e., (5) for a constant ω(τ) = ω.For any non-constant ω(τ) we can than choose ŝ(τ,t) = s(τ,t + σ(τ)) where If ω(τ) is fixed, then σ(τ k ) needed for the computation of f k can be computed exactly or at least with arbitrary accuracy.If ω(τ) is not known, then an approximation σ k ≈ σ(τ k ) has to computed together with ω k .We will use the quadrature formula In particular, with weights W = 0, W = 1 and W = 1  2 we obtain the left and right rectangular rule and the trapezoidal rule, respectively. 2If W > 0 we have to take into account σ k and thus f k (x,t) depends on the unknown ω = ω k and with ( 8) and ( 12) the derivative can be calculated as

A spline Galerkin method for the periodic boundary value problems
In [5,6], an adaptive spline wavelet method for the initial value problem on the circuit equations (1) has been developed.For the periodic problem we will modify this approach by using a periodic basis.Then the periodic boundary conditions (10) are fulfilled automatically and has not to be enforced by additional equations.
We want to approximate the solution of ( 9) by a periodic spline function of order m.For given grid points t with 0 = t 0 < t 1 < . . .< t N = P we consider all m − 2-times differentiable, P-periodic functions, which are polynomials of degree less than m at each subinterval (t ,t +1 ).The break points t are called spline knots, which are periodically extended by t kN+ = t + k P. Note, that the periodicity condition implies that the periodic spline is a piecewise polynomial also with respect to the extended grid.
A stable and computational efficient basis for the linear space of spline functions is constituted by the B-splines N m (t), which are uniquely determined by their minimal support [t ,t +m ] and the partition of unity ∑ i N m i (t) = 1 (for normalization).For more information on spline functions and B-splines as well as for efficient computational methods we refer the reader to the detailed description in [18,30].
To expand the periodic solution of (9), we need periodic basis functions, which we obtain by the periodized Bsplines which form a basis for the space of periodic spline functions.
Now we have to find a spline function which approximates the solution of ( 9).Since c ∈ R n (n is the number of equations and unknowns in the circuit equations (1)) we have to determine n × N coefficients c k,i .Thus, we have to derive from ( 9) n × N conditions, which ensure that x indeed approximates the solution of ( 9).Taking the integral over N subintervals we obtain = ω q k x(t ) − q k x(t −1 ) + t t −1 f k x(t),t dt for = 1, . . ., N. This N vector valued equations determine the vector coefficients c k .Note, that the splitting points tk do not coincide with the spline knots t k , but they have to be chosen in relation to the spline grid.In particular, we need due to periodicity that tN − t0 = P.
The nonlinear system (15) can be solved by Newton's method.If ω = ω(τ) is known in advance we have to solve the linear system where j is the Newton count.Starting from a sufficiently good initial guess, e.g.c c c (k,0) = c c c (k−1) , one obtains usually after sufficiently many steps a good approximation c c c (k) = c (k,J) of the solution of (15).In order to set up the linear system we have to compute the right hand side 3 where the Jacobians C = D x q(x) and G = D x g(x,t) are available in any circuit simulator, which is based on the MNA formulation (1).The integrals can be computed by a suitable quadrature rule.In practical computations where the spline order is chosen usually as 3 or 4 one can use the Simpson rule or the two point Gauß quadrature.
If ω k = ω(τ k ) has to be determined during the computation, linearization of (15) results in the underdetermined system under condition (11).Here z z z = D ω F (c c c, ω) =1,...,N , where with ( 13) For computational reasons we replace ( 11) by a similar condition on the spline coefficients, namely Following [4] we obtain the solution of ( 17) which satisfies ( 18 where in addition to ( 16) ω k is determined by the iteration

Adaptive grid generation
It is still open, how the spline knots t i are chosen.In particular, for functions with sharp transients it is important that we put more grid points at the location of this refinement.Since those locations are often not known in advance, it is important that the grid is generated automatically by an adaptive refinement.For this refinement we use spline wavelets for non uniform grids introduced in [1].
Starting on some coarse initial grid we solve (15) by Newton's method and obtain a first approximation, which we denote as X(t).Now we apply a fast wavelet transform introduced in [3] and obtain The spline X(t) is an approximation of X(t) using only the even spline knots t 2k .The ψ k are wavelets which carry detail information.In the setting of [1] the wavelets are compactly supported splines, with a chosen number of vanishing moments, which are localized near t 2k+1 .Thus, the coefficients d k are a measure of the local approximation error, which describes also the smoothness of X(t).Therefore we insert additional knots in the neighborhood of t 2k+1 if the coefficient d k exceeds a given threshold.The spline representation of X(t) for the new grid can then be computed efficiently by the Oslo algorithm [15,20].These refined spline expansion is used as a new initial guess for another Newton iteration.
This approach is repeated several times, with Newton tolerances and wavelet threshold reduced for each refinement.This leads to better and better approximations, which provide more and more information for the grid refinements leading to an almost optimal grid.Furthermore, this approach is computational efficient, since the first Newton iterations are computed on a relatively small grid.For the final refinements usually one Newton update is sufficient, since we have already an excellent initial guess.The whole process is stopped if the Newton iteration after the last refinement gets below a given error bound.These error bound is indeed a very good estimate of the achieved error.For details on the refinement algorithm we refer to [3].
In Sect. 5 we suggested to use the solution X k−1 of the previous time step τ k−1 as initial guess (predictor) for Newton's method.In fact any reasonable predictor should be based on solutions from previous time steps such that the following approach should be used for other predictors, too.Due to the smoothness in τ the difference between X k−1 and X k should be small, and we can expect fast convergence of Newton's method.However, the refinements in the adaptive method described above will increase the number of spline knots and thus the size of the problem.Therefore, we will use an approximation of X k−1 on a coarser grid as initial guess.An efficient approximation can be achieved by wavelet thresholding as follows.From the wavelet expansion (19) of X k−1 we remove terms with small wavelet coefficients d k , i.e., an approximation is obtained as If the threshold ε is chosen in the order of magnitude of the error tolerance of the method, than we can still expect the initial guess Y (t) to be as good as X (t) itself, while the computational cost is reduced.Using the approach from [3] removing ψ k (t) is equivalent to remove the knot t 2k+1 .Thus, the method will remove spline knots, which became dispensable due to a change in signal shape over several time steps τ k , or which where not needed from the beginning due to the rough error estimate in the refinement described above.This leads to an almost optimal grid for any given error tolerance.
For more details on this grid coarsening algorithm we refer to [3].

Numerical tests
The above algorithm is implemented in C++, using a C++ library for circuit description and device models.

Gilbert mixer
Our first example is the simulation of a downconverting Gilbert mixer.The inputs are a sinusoidal frequency of f rf = 99.9MHz and a digital reference signal of fundamental frequency f lo = 100MHz, which switches between −5V and 5V.Since the driving frequencies are known, we have chosen ω(τ) = 1 and carrier period to P = 10 −8 s, which corresponds to the frequency of 100MHz.In Fig. 1 shows the resulting output signal.In τ-direction one can see the smooth sinusoidal envelope of frequency f out = f lo − f rf = 0.1MHz, while the t direction shows the influence of the carrier signal.In particular, at the jump discontinuities of the input signal one can observe an overshoot of the output signal, due to capacitive effects.The simulation with 28 envelope time steps for 20000 periods needed 10s, while a transient simulation was done in 8min.In Fig. 2 the corresponding adaptive spline grids for the envelope time steps τ k are shown.As expected, it can be observed that the adaptive wavelet refinement leads to a fine grid at the jump locations, while only few spline knots are needed to represent the smooth parts in between.Furthermore we can see that roughly 30 envelope time steps are needed for an interval which contains 2200 oscillations of the reference signal.Obviously the spline grid does not change much with since the locations of sharp transients do not change with τ.However, we will consider examples in the sequel, where an optimal grid is not known in advance, but may even change with τ.

Colpitts Oscillator
In our next example we consider the start up phase a 3MHz Colpitts-Quartz oscillator.To avoids unwanted effects from numerical damping, we have used trigonometric splines [30, Sect.10.8] instead of the usual polynomial splines (see [2] for details).The transistor used as feedback amplifier introduces some nonlinear effects into the output (see Fig. 3), which results in a sharp edge near t = 0.24µs.The simulation with 202 envelope time steps for 7500 periods needed 10s, while a transient simulation was done in 165min.The spline grids in Fig. 4 are adapted to this edge, when it becomes prominent at τ ≈ 5ms, and the grids are even adapted to a slight change of location in the following envelope time steps, which shows the perfect functioning of our approach.

A Phase Locked Loop
Our next example is a Phase Locked Loop (PLL), which generates an output signal whose phase is related to the phase of the input "reference" signal.As depicted in Fig. 5, this is achieved, by comparing the phases of the "reference" signal and the output or "feedback" signal in a phase detector.The result is than filtered, in order to stabilize the behavior of the PLL, and fed to a Voltage Controlled Oscillator whose frequency depends on the (filtered) phase difference.The output is fed back (possibly to a frequency divider) to the phase detector.If the PLL is locked, the phases of reference and feedback signal are synchronized.Thus both signal will have (almost) the same local frequency.We have tested a PLL (containing 205 MOSFETs) with a frequency modulated sinusoidal signal with center frequency 25kHz, which corresponds to 800kHz oscillator frequency due to a frequency of factor 32.The baseband signal is also sinusoidal with frequency 10Hz and frequency deviation 100Hz.The estimated instantaneous frequency in Fig. 6 (determined based on ω(τ)) shows after a short startup phase a good agreement with the baseband signal.This is also reflected by the control signal of the VCO in Fig. 7, which is obtained by filtering the output of the phase detector (Fig. 8).Note, that the envelope corresponds to the baseband signal, while the carrier signal is due to the filtering almost constant.This allows to do the shown envelope simulation using only 160 envelope time steps, while a corresponding transient analysis would contain 5000 oscillations.The wavelet envelope simulation of this circuit was done in 7 min.However, a comparable transient simulation on the same simulator needed 33 hours.We will have a closer look to the start-up phase before the PLL is locked, which gives a good impression how the adaptive wavelet method works.For this we have a closer look at the phase detector output for t near 33ms, where sharp edges are present as one can see in Fig. 9.The change in the signal shape is due to changes in the phase difference (with τ) during transient response.In Fig. 10 one can see the corresponding spline grid.Obviously, the grid is refined around the sharp edges in the signal, as expected b b b = F(c c c, ω) = F (c c c, ω) =1,...,N as well as the Jacobian (with respect to c c c) A = D c F(c c c, ω) with the matrix block entries ) by d d d c = b b b − d ω z z z, where b b b = A A A −1 b b b and z z z = A A A −1 z z z are computed by solving the corresponding linear systems and