Stability Characterizing Function for Electronic Circuit Design Based on Frequency-Domain Analysis With Parametric Damping

A stability characterizing function (SCF) to be used for small and large-signal stability analysis and design of single and multitransistor electronic circuits is proposed. Stability constraints based on this function can be integrated with the framework of standard computer-aided design (CAD) procedures. The method is based on the bounded-input–bounded-output (BIBO) stability criterion, where a suitable set of input and output variables is chosen at the ports of all the nonlinear intrinsic electron device (IED) models. These variables are linked by a nonlinear state perturbation (NSP) matrix, whose elements are directly computed by analyzing an associated parametrically damped (PD) (i.e., stabilized) circuit. The proposed SCF is defined in terms of the NSP matrix elements and is a scalar function of frequency and the damping parameter. Circuit stability margins are easily evaluated by visual inspection of the plot of this function. Preliminary experiment validation is carried out by applying this approach for checking the self-starting capabilities of an oscillator as well as for the small and large-signal stability analysis of a two-transistor balanced amplifier and a single-transistor PA.


Stability Characterizing Function for Electronic
Circuit Design Based on Frequency-Domain Analysis With Parametric Damping Alberto Santarelli , Member, IEEE, Leonardo Pantoli , Member, IEEE, Giorgio Leuzzi , and Fabio Filicori Abstract-A stability characterizing function (SCF) to be used for small and large-signal stability analysis and design of single and multitransistor electronic circuits is proposed.Stability constraints based on this function can be integrated with the framework of standard computer-aided design (CAD) procedures.The method is based on the bounded-input-boundedoutput (BIBO) stability criterion, where a suitable set of input and output variables is chosen at the ports of all the nonlinear intrinsic electron device (IED) models.These variables are linked by a nonlinear state perturbation (NSP) matrix, whose elements are directly computed by analyzing an associated parametrically damped (PD) (i.e., stabilized) circuit.The proposed SCF is defined in terms of the NSP matrix elements and is a scalar function of frequency and the damping parameter.Circuit stability margins are easily evaluated by visual inspection of the plot of this function.Preliminary experiment validation is carried out by applying this approach for checking the self-starting capabilities of an oscillator as well as for the small and large-signal stability analysis of a two-transistor balanced amplifier and a singletransistor PA.

I. INTRODUCTION
T HE design of analog electronic circuits consists of searching for a circuit topology and an associated set of values of the component characteristic parameters, which provide acceptable performance under a set of hard design constraints aimed at guaranteeing feasibility, reliability, and stability.
Conditions for acceptable circuit performance, as well as feasibility and reliability constraints, can be synthetically and formally defined by a set of inequalities expressed in Alberto Santarelli and Fabio Filicori are with the Electrical, Electronic, and Information Engineering Department (DEI) "Guglielmo Marconi," University of Bologna, 40136 Bologna, Italy (e-mail: alberto.santarelli@unibo.it;fabio.filicori@unibo.it).
Leonardo Pantoli is with the Department of Industrial and Information Engineering and Economics, University of L'Aquila, 67100 L'Aquila, Italy (e-mail: leonardo.pantoli@univaq.it).
Giorgio Leuzzi is with the Department of Information Engineering, Computer Science and Mathematics, University of L'Aquila, 67100 L'Aquila, Italy (e-mail: giorgio.leuzzi@univaq.it).
Color versions of one or more figures in this article are available at https://doi.org/10.1109/TMTT.2023.3305157.
Digital Object Identifier 10.1109/TMTT.2023.3305157terms of figures of merit, which are defined under periodic steady-state operating conditions [1], [2], [3].For this reason, they can be easily computed by Fourier transform (FT)-based frequency-domain circuit solution algorithms in the framework of commercial CAD tools (e.g., [2], [3]), like linearized ac analysis for small-signal operation and harmonic balance (HB) analysis for large-signals.
Stability constraint functions cannot be so directly and easily computed by using the same periodic-steady-state frequency-domain solution algorithms, since stability criteria are defined in terms of circuit/system transient responses or, equivalently, in terms of Laplace transform (LT)-based circuit/system analysis.Several methods have been proposed in the literature to solve the problem within the framework of commercial CAD tools.Linear time-invariant (LTI) stability problems, i.e., small-signal stability analysis, can be addressed in terms of frequency-domain-based open-loop function analysis [4], [5], [6], [7].For instance, the widely used method [7] analyzes the stability of multidevice power amplifiers by looking at Nyquist plots of open-loop transfer functions obtained by inserting circulators and isolators at the interface ports between passive and active subnetworks.Unconditional stability conditions [1], [8] are sometimes searched for by adding resistive feedback networks to transistors [9].
An overview of local and global stability analysis techniques can be found in [10] and [11].An analytical study of the linear time variant (LTV) stability problem, i.e., largesignal stability analysis, is addressed in [12], and further developed in [13] and [14].In the latter, an open-loop system determinant calculation procedure within the framework of frequency-analysis-based CAD tools is carried out by modifying transistor models to evaluate suitable return ratio matrices [5].Alternative nonlinear stability analysis methods are based on pole-zero identification of closed-loop transfer functions [15], [16], [17], [18], [19], [20], [21].Methods based on conversion matrix [22] and bifurcation analysis through auxiliary generators [23] have been proposed for the design of frequency dividers.An extension of the Othomo technique [7] for multitransistor circuits to the nonlinear stability analysis case has been proposed in [24].The technique described in [25] and [26], combines conversion matrix and even/odd mode theories for splitting the nonlinear stability analysis of symmetrical multitransistor circuits into a set of equivalent single-transistor problems.Time-domain circuit simulation (e.g., [27], [28], [29], [30], [31], [32], [33], [34]) represents an interesting alternative environment for the development of stability analysis methods.Some of them are based on the approximation of passive distributed circuit subsections in terms of lumped-element equivalent networks (e.g., [32], [33], [34]).
In this article, we propose a novel stability characterizing function (SCF) and a related set of stability constraints, which can be systematically computed using the same algorithms used for circuit performance evaluation, i.e., by means of FT-based periodic-steady-state analysis methods available in commercial CAD tools (e.g., [2], [3]).The positive feature of this approach is that it enables defining stability criteria valid for both single and multitransistor circuits, under both smalland large-signal operating conditions.Moreover, considering that the SCF is a real function of the frequency and an auxiliary real parameter only, the method has the advantage of performing the stability/instability check by direct inspection of a single plot of a bivariate function even in the case of multitransistor circuits.
This article is structured as follows.A proper choice of input perturbation sources and output variables to be used in the framework of the bounded-input-bounded-output (BIBO) stability criterion is first discussed in Section II.By considering that circuit instability is necessarily associated with the presence of incrementally active nonlinear devices, only the state variables of charge-conservative quasi-static intrinsic electron device (IED) models are chosen as output variables, while a set of input series equivalent voltage sources at the IED ports is adopted.
The new set of stability constraints for computer-aided design (CAD) of single and multitransistor electronic circuits is proposed in Section III.These are expressed in terms of a SCF based on the elements of a nonlinear state perturbation (NSP) matrix of transfer functions, which links the input and output variables for BIBO stability analysis.It is shown how the NSP elements along with the SCF can be directly computed through frequency-domain periodic steady-state analysis of an associated parametrically damped circuit (PDC).
The features of the PDC are introduced in Section IV, along with practical guidelines for its implementation in CAD tools for circuit analysis and design.The PDC equations are directly obtained by simply adding parametrized damping terms to all the derivative terms in the state-space linear and nonlinear differential equations, which describe the whole circuit dynamics.Practically, the implementation of the PDC can be carried out, in the framework of commercially available CAD tools for circuit analysis, by adding the same type of parametric damping modifications to all the circuit component model equations.
Preliminary validation examples of small and large-signal stability analysis are presented in Section V to confirm that the proposed stability constraints, and associated SCF, provide sufficient and necessary conditions for circuit stability of single and multitransistor circuits.The method is used for ascertaining the self-triggering of an oscillator, and for checking the small-and large-signal stability features of both a two-transistor balanced amplifier and a single-transistor PA.Finally, the possibility of integrating the new stability constraints in the framework of iterative computer-aided circuit design procedures is outlined in Section V, whereas conclusion is drawn in Section VII.Details and examples related to the parametric damping of linear and nonlinear state-space models are reported in Appendixes A and B, respectively.

II. VARIABLES SELECTION FOR BIBO STABILITY ANALYSIS
Let us consider a circuit with N P1 and N P2 incrementally active one-and two-port electronic devices, respectively.Let us also assume that the IED ports are all accessible and that the extrinsic parasitic networks are all included in a single passive linear network.
According to this type of description, the circuit under test in Fig. 1 is split into two sections, the same way as done in HB-based circuit analysis and similar to [12] and [13].The first section is a nonlinear active circuit, consisting of the full set of IEDs, and the second is a linear passive network, which includes the impedance matching networks, the bias networks, the extrinsic parasitic components of the active devices, and any other linear passive component in the circuit.
All the voltages and currents at the IED ports can be grouped into two vectors consisting of N P = N P1 + 2N P2 elements, corresponding to the total number of nonlinear ports.The circuit is powered by a set of bias generators V B and signal sources a S (t).

A. Output Variables
The state variables describing each IED are related to the amount of charge in the intrinsic regions.The IED port voltages v(t) are chosen here as state variables since the functional relationship between charges and voltages at the intrinsic IED ports can be considered almost purely algebraic (quasi-static assumption) in most practical cases so that the current response can be directly expressed also in terms of voltage derivatives. 1y naming x(t) the whole set of circuit state variables, the voltage vector v(t) related to IEDs is the subset of x(t), which is relevant for circuit stability analysis.In fact, if all elements of v(t) satisfy stability criteria, then the same is necessarily true for any other circuit state variable, given that the amount of energy in the passive, linear, lossy network is attenuated over time.
According to this statement, only the perturbations of v(t) will be considered as output variables for BIBO-based stability analysis purposes, with the great advantage that the number N P of nonlinear ports is much lower than the total number N x of circuit state variables.On the other hand, it is worth noticing that the knowledge of v(t) could allow the direct computation of the entire set x(t) through standard nodal circuit analysis, once the IEDs current response i(t) is provided from nonlinear dynamic models.This choice of state variables is the same adopted in [12].

B. Input Variables and Equivalent Perturbation Sources
The multiport Thevenin theorem is applied at the IED ports, transforming the circuit in Fig. 1 into the one shown in Fig. 2, where ṽ = [ ṽ1 , . . ., ṽN P1 , ṽ1,1 , ṽ1,2 , . . ., ṽN P2 ,1 , ṽN P2 ,2 ] is the vector of voltages at the N P ports of the Thevenin-equivalent linear passive network and e(t) = E B +e S (t) are the Theveninequivalent sources, which corresponds to the turned-off voltage generators V B and a S (t).The Thevenin-equivalent linear passive network can be described by state-space linear differential equations where x LIN (t) is a vector of state variables describing the internal state of the network and A , B , C , D are matrices depending on a set of circuit design parameters.
As outlined in Section II, the nonlinear network (i.e., the set of N IEDs) is described here in terms of quasi-static nonlinear models where F and Q are vectors of memoryless nonlinear functions and The triggering of circuit instability must be sought either in internal noise generation or the possible existence of external interferences, which may affect all the circuit components.These perturbations need to be considered not only to compute the circuit performance in terms of signal-tonoise ratio but also for stability analysis purposes.By treating these noise/interference sources as if they were any other independent generator included in the circuit, it is concluded that stability analysis can be based on equivalent perturbation sources e(t) in series with the Thevenin-equivalent voltage sources e(t).
The perturbation effects on the IED port voltages, deriving from the equivalent perturbations e(t), can be described as splitting v(t) into two contributions where v(t) is the unperturbed stationary solution (IED port voltages corresponding to e(t) = 0), while v(t) describes the effects of perturbation sources e(t).Analogously, let the corresponding currents flowing into the IEDs be named i(t) = î( t) + i(t).
According to the BIBO stability criterion, instability detection can be based on the observation of the time evolution of perturbations v(t) around the stationary solution e caused by small-amplitude perturbation sources e(t).To this aim, the definition of a NSP matrix is introduced in the following, by separately considering the LTI (i.e., small signal) and LTV (i.e., large-signal) analysis cases.

A. NSP Matrix
In the case of small-signal stability analysis (i.e., LTI problem), the steady-state condition to be checked for stability is defined by the circuit solution in the absence of perturbations, i.e., when v s (t) = 0 and e(t) = 0. Let this be v(t) = V Q , where V Q is the time-invariant V B is the dependent vector of quiescent voltages across the IED ports.
When perturbations e(t) are small enough to justify linearization, the relationship between e(t) and e(t) can Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
be expressed in terms of a convolution integral where w(t) assumes the meaning of a time-invariant NSP matrix, i.e., a matrix of real impulsive responses, which describes the evolution of the state variables v(t) at the port of the nonlinear IEDs.The matrix element functions w m,n (t) should be vanishing in time (or remain amplitude-limited) as time approaches infinity in the case of a stable circuit.By naming p = σ + jω the LT-domain complex variable, state variable evolutions ( 5) can be equivalently expressed as where }, and L , L −1 are the direct and inverse-LT operators.
Coherently with the BIBO stability criterion, each complex function W m,n ( p) of the NSP-matrix W ( p) (with m, n = 1, . . ., N P ) must be evaluated.If, and only if, these functions show all poles with negative (or zero) real parts, the stationary solution v(t) = V Q will be LTI stable, i.e., the impulsive functions w m,n (t) will be vanishing or practically amplitude limited.Since the presence of a pole in the right half of the p-plane corresponds to a diverging behavior (function max), stability constraints can be put in the form |W m,n (σ with m, n = 1, . . ., N P .When large-signal sources a s (t) are also present, the unperturbed steady-state solution v(t) at the IED ports becomes time-varying (LTV problem) and the large-signal operating point (LSOP) to be checked for stability depends both on the given bias V B and signal source amplitudes a S .
In this case, the relation between the small perturbations v(t) and the state variables evolution v(t) is defined as where w(t, τ ) represents a time-variant NSP matrix of timepulse responses, depending on the actual LSOP (dependency on V B and a S is not shown explicitly for the sake of notation simplicity).By assuming a large-signal periodic regime in the t-domain (with period T s ), the matrix w(t, τ ) can be expanded in the Fourier series.By limiting the analysis to the maximum order K , we obtain where s = 2π /T s .The w k (τ ) terms represent complex NSP matrices of pulse-response-like elements, associated with each kth order harmonic component of the series expansion.By replacing ( 9) in ( 8), we obtain v k (t) being complex Fourier coefficients of v(t) given by By considering the conversion functions To the aim of stability analysis, each function element W m,n,k ( p) (with k = 0, ±1, . . ., ±K , m, n = 1, . . ., N P ) of the complex matrixes W k ( p) needs to be evaluated.If, and only if, these functions have all poles with negative (or zero) real parts, then the LSOP considered (associated with signal source amplitude a S ) will be LTV stable, leading to evanescent or amplitude-limited complex functions w m,n,k (t).Thus, a set of stability constraints can be defined at any given source amplitude a S as with m, n = 1, . . ., N P , k = 0, ±1, . . ., ±K .
It is worth noting that, the definition of the NSP matrices W k ( p) under large-signal operation could be extended in line with the principle from the strictly periodic to the more general multitone case.To this aim, the periodicity of each single elementary tone could be exploited in a way similar to how HB single-tone analysis can be extended to the multitone case (e.g., [2], [3]).A set of NSP matrices W k 1 ,...,k N T ( p), where k 1 , . . ., k N T are mixing-tone indexes, would be obtained in the presence of a N T -tones excitation.

B. FT-Based Laplace-Equivalent NSP Matrix
A method is proposed for the evaluation of the NSP matrices W ( p) (LTI case) and W k ( p) ∀k (LTV case), through purely FT-based analysis techniques.
The method is based on the introduction of a modified circuit, namely, a PDC, with the following features: 1) the PDC can be analyzed by means of standard FT-based CAD simulation tools; 2) its electrical response is dependent on a real-valued damping factor σ D ≥ 0; and 3) the frequency transfer functions between small-signal sinusoidal perturbations e and state variables deviations v from steady-state solution v are described by a matrix W( jω, σ D ) (LTI case), or a set of matrices W k ( jω, σ D ) ∀k (LTV case), which are equal, respectively, to W (σ D + jω), or W k (σ D + jω), for any choice of σ D within the Region of Convergence (RoC) of the LT-domain matrix functions.Let us consider the LTI case first.According to feature 3), linear ac-type analyses of the PDC lead to W( jω, σ D ), with Similarly, in the LTV case, provided that σ D is chosen greater than the real part σ rmp,m,n,k of the rightmost pole of any function element in any W k ( p) ∀k, i.e., σ D > σ rmp,LTV with σ rmp,LTV = σ rmp,m,n,k (RoC of the whole set of W k ( p) ∀k), Hamonic Balance-based analyses of the PDC lead, for any given LSOP with signal amplitude a S , to a set of W k ( jω, σ D ) matrices featuring with k = 0, ±1, . . ., ±K .Thanks to (16), it is concluded that PDC building rules compatible with features (14) (LTI) and ( 16) (LTV) are introduced in Section IV.Practical examples of parametrically damped (PD) linear circuit components are given in Appendix A, whereas nonlinear components are dealt with in Appendix B.

C. SCF and BIBO Constraints
According to the BIBO criterion, inequalities ( 7) and ( 13) define the set of small-and large-signal stability constraints for electronic circuit design.Thanks to features (14), (16) of the associated PDC, the LT-based stability constraints can be evaluated in terms of FT-based functions obtained by analyzing the σ D -dependent PDC with frequency-domain ac or HB circuit analysis.In particular, the LT-domain W m,n (σ + jω), or W m,n,k (σ + jω) functions in (7), or (13), can be replaced with the associated equivalent FT-domain functions Stability constraints can be verified across finite ranges of both perturbation frequencies f = ω/2π and damping factors σ D .Frequencies can be swept from a minimum value f min , above the bias network cutoff region, up to a maximum value f max near the highest maximum oscillating frequency of all the electron devices, i.e., f max ≃ f max osc .On the other hand, the potentially unstable region of the LT p-plane can be practically limited to the range of damping factors 0 ≤ σ D ≤ σ D,max , σ D,max being a value corresponding to the incrementally passive behavior of any IED PD model over the whole frequency axis (thus leading to a necessarily stable PDC).It is worth noting that, from a practical point of view, σ D,max can be searched for as a value large enough to make any IED PD model capable of providing available power gain less than one for any choice of source and load terminations at any frequency.In cases where the S-matrix of the IED does exist (i.e., if the device is stable on 50 ), σ D,max could also be searched for as a σ D value large enough to make the PD IED model unconditionally stable for any source and load termination at any frequency.To this aim, known necessary and sufficient conditions for unconditional stability could be practically used (e.g., [5], [8]).
All these boundary values can be sought at the beginning of the stability analysis process, by exploiting well-known formulas based on small-signal S-parameters of electron devices (well-known equations to determine f max osc and/or to check the unconditional stability of transistors can be found for instance in [1]).
Since stability constraints must be necessarily satisfied by also considering nonnegligible technological uncertainties in circuit implementation and inaccuracies in component models, an adequate margin is introduced.In addition, a normalization term is conveniently adopted to choose a unique stability margin ϵ lim > 0 shared by all the stability constraints (e.g., ϵ lim = 0.01 in this work).On this basis, ( 7) and ( 13) can be practically replaced by with m, n = 1, . . ., N P and k = 0, ±1, . . ., ±K .
For a simpler and more compact circuit stability analysis, an overall SCF can be introduced.In the LTI case, this can be defined as where the min (•) operator selects the minimum value of its argument across the m, n index space. 2According to the set of constraints (18) and the SCF definition (20), it is concluded that the quiescent condition V Q of the circuit under test is LTI stable if, and only if Similarly, in the LTV case and for any given source amplitude a S , the SCF can be defined as where the min (•) operator is here applied across the m, n, k index space.According to the set of constraints (19) and the SCF definition (22), it is concluded that the circuit is LTV stable under large-signal operation around the quiescent condition V Q if, and only if, S C ( f, σ D ) > ϵ lim for any σ D > 0 and f > 0, which is a condition formally identical to (21).Practically, (22) will be checked for a discrete set of source available powers associated with a S .
Constraint (21) can be considered a necessary and sufficient condition for BIBO stability.In fact, if (21) is violated then at least one element of the NSP matrix will be unbounded for some positive values of f , σ D and the circuit will be BIBO unstable (necessary condition).Conversely, if (21) is satisfied then all the elements of the NSP matrix will be bounded over the entire f > 0, σ D > 0 domain and the circuit will be BIBO stable (sufficient condition).
On this basis, circuit stability/instability can be simply verified by directly examining the SCF plots, drawn over an adequately dense and wide sampling grid of the variables f and σ D .In both LTI and LTV cases, when stability constraint (21) is not met, the location of the rightmost pole in the LT-plane causing unstable behavior (i.e., violating the stability margin) is identified by the values of f and σ D , which minimize the SCF.This can also be directly seen by visual inspection of SCF plots shown for the application examples provided in Section V.
Depending on the actual application context, algorithms for function minimization could be alternatively exploited for the verification of the SCF-based constraint (21), by avoiding the more computationally expensive function evaluation over dense grids in the f , σ D domain.

A. PD Circuit
The original circuit under test is modified in order to satisfy FT-domain to LT-domain equivalence conditions, leading to verification (14) (LTI case) and ( 16) for any choice of k indexes and signal amplitude a S (LTV case).
The PDC of the Thevenin-equivalent linear network in Fig. 2 is considered first.Original port voltages ṽ(t), currents i ( t), and state variables x LIN (t) are replaced by corresponding modified PDC variables ṽ(t, σ D ), i(t, σ D ), and x LIN (t, σ D ), which are dependent on a damping parameter σ D .The model of the PDC in the FT domain must equal the Laplace-transformed set of equations ( 1), (2) for any σ = σ D within the RoC of the LT functions.To this aim, it is sufficient to modify the original circuit so that a term σ D x LIN (t, σ D ) is added to all the time derivatives in (2), leading to the following PDC linear model, i.e., Examples of how linear lumped and distributed circuit components can be practically modified to achieve this result and proof of FT-domain to LT-domain equivalence are given in Appendix A.
By proceeding in a similar way, the model of the PD nonlinear IED network is obtained in terms of the intrinsic transistor port voltages v(t, σ D ), currents i(t, σ D ), and charges Q(v(t, σ D )).To achieve FT-domain to LT-domain correspondence a term σ D Q needs to be added in (3) to all the charge derivatives, i.e., with dQ(v(t, σ D )) = C(v(t, σ D ))v(t, σ D ).See Appendix B for major details on the PD nonlinear IED network model IV-B and the proof of FT-domain to LT-domain equivalence under linearized conditions.It is worth observing that the quantities σ D x LIN (t, σ D ) in (24) and σ D Q(t, σ D ) in IV-B are damping terms since, as observed in (15) and (17), they shorten the duration of the circuit pulse responses.For very large positive values of σ D , the time derivatives in (24) and IV-B become negligible with respect to the algebraic damping terms, leading to a stable, almost memoryless PDC due to the very short duration of pulse responses.Parametrical damping with increasing positive values of σ D modifies pole locations so that unstable poles are moved from the right to the left half-plane of the LT domain and any unstable circuit is transformed into a stable one.For this reason, the damping parameter σ D can be also considered as a circuit stabilizing parameter.Thanks to this behavior, the PDC can be used for stability analysis in the FT domain, by detecting the unstable pole crossing the imaginary axis of the LT plane for some value of σ D .
For practical implementation, in the framework of existing CAD tools, the PDC can be directly obtained by replacing all the circuit component models with the corresponding PD ones, with unmodified circuit topology, as shown in Appendix A.
The full equivalence of the nonlinear IED network PDC model IV-B to the aim of LTI and LTV stability analysis is mathematically proved in Appendix B.

B. NSP Matrix Computation With FT-Based Analysis of PDCs
The circuit schematic for the computation of the NSP matrix W( jω, σ D ) (LTI case), or set of matrices W k ( jω, σ D ) (LTV case), by using FT-domain-based periodic steady-state analysis algorithms (ac in the LTI case or HB in the LTV one), is shown in Fig. 3.This circuit is made up of three parts, to be conveniently placed on the same schematic page of the commercial CAD tool adopted.The circuit for performance analysis is represented in Fig. 3(c), whereas the associated PDC for stability analysis is shown in Fig. 3(a).Smallsignal ac generators e at angular frequency ω are used as perturbation sources of the current steady-state solution in the PDC [Fig.3(a)].This circuit is connected to Fig. 3(c) and the associated auxiliary PDC is represented by Fig. 3(b).In fact, circuits [Fig.3(b) and (c)] are both needed for the calculation of the Thevenin-equivalent sources e appearing in the PDC [Fig.3(a)].The functions of the three circuit parts are described in more detail in the following.
The steady-state voltages v at the IED ports are computed by means of the unperturbed original circuit shown in Fig. 3(c).When an ac analysis is executed (with signal generators a S off), the steady-state voltages will be the time-invariant quiescent values v = V Q (LTI case), whereas, when a power-swept multitone HB analysis is carried out (for instance, two-tone mixer-mode HB with one large-signal tone at an angular frequency S and a small-signal tone at ω), the steady-state voltages v(t) will represent time varying LSOPs depending on the current levels a S of large-signal generators (LTV case).Thevenin-equivalent voltage sources e are applied in the PDC [Fig.3(a)] by means of voltage-controlled voltage sources (VCVSs).According to the Thevenin-equivalent circuit representation shown in Fig. 2 and by imposing the conservation of the steady-state voltages v at the ports of the PD nonlinear IED network, the appropriate e values are: e = v + ṽ.The v contribution is copied from circuit [Fig.3(c)], whereas ṽ is evaluated by means of the auxiliary PDC [Fig.3(b)] and copied to Fig. 3(a).In circuit [Fig.3(b)], the steady-state solution v is considered as an input.By using VCVSs, the voltages v in Fig. 3(c) are copied to the PD-IED ports in Fig. 3(b) leading to the computation of the corresponding σ D -dependent currents i flowing into the PD-IED network.These are in turn copied with current-controlled current sources (CCCSs) and forced to flow into the Thevenin-equivalent PD linear passive network in Fig. 3(b), leading to the evaluation of the corresponding port voltages ṽ.Finally, the Thevenin-equivalent sources e = v + ṽ are inserted into Fig.3(a) by using VCVSs, which copy voltages v from (c) and ṽ from Fig. 3(b).
It is worth noting that, apart from the preliminary creation of a PD component model library as described in Section IV-A, the implementation and analysis of schematics shown in Fig. 3 only require standard components and ac/HB simulation methods commonly available in commercial CAD tools (e.g., [2], [3]).

C. Procedure for Stability Analysis
Stability analysis is carried out by simulation of the circuit in Fig. 3 with parametric sweeps over the damping parameter σ D , the perturbation frequency f and, in the case of LTV analysis, over the large-signal levels a S .The N P ac sources e provide small-signal sinusoidal perturbations that can be switched on one at a time, by repeating the parametric sweeps for the evaluation of different columns of the NSP matrix (14) (LTI case) or set of matrices (16) (LTV case).Simple ac simulations are sufficient in the case of small-signal stability analysis (LTI case), whereas multitone HB analysis is used for the large-signal case.Standard mixer-like HB simulation engines provided by commercial CAD tools (e.g., [2], [3]) can Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
be used for this aim and HB power sources adopted for both the large-signal tone(s) a s (t) and the small-signal perturbations e (for instance, with amplitudes as low as −30 dBm).
After completing the NSP matrix(es) computation, necessary and sufficient BIBO stability conditions (21) are easily evaluated in terms of the SCFs (20), or (22).
It is worth noting that, although the schematics related to the application examples provided in this work were created manually, dedicated software-based procedures could be implemented for an automatic netlist generation of the auxiliary PDC and PDC in Fig. 3.

V. EXPERIMENT VALIDATION
All the stability test cases provided are based on BFR92A medium-power silicon bipolar transistors from Infineon.The nonlinear Gummel-Poon model provided by the foundry was modified into a quasi-static version with parametrical damping.
It is worth noting that the PDC, based on components described in Section IV-A and Appendix A might involve some approximations.Although the PD models might not be fully adequate for accurate performance prediction, especially at high operating frequencies, 3 they can be conveniently used for stability analysis, where accuracy requirements are weaker when considering adequately large stability margins, especially when the corresponding inaccuracies are of a pessimistic type (e.g., losses in passive components act as additional damping and quasi-static IED models tend to overestimate gain).
In the following validation examples, stability analysis was carried out by using a commercial CAD tool [3], by adopting an adequately dense grid of damping factors σ D (e.g., 10 points/ • ) and frequencies (e.g., f = 2 MHz) to detect the minima of the SCFs (20), or (22).As far as the stability margin is concerned, the chosen value of ϵ lim should be neither too small nor too large.An excessively large stability margin might involve practical limitations on the possibility of achieving near-optimal circuit performance in the circuit design phase.On the other hand, too small values of ϵ lim could lead to insufficient margins in consideration of technology process dispersion and possible inaccuracies in component models.Too small values of ϵ lim could as well lead to using unnecessarily dense grids in the f , σ D domain for detecting instabilities.Since ϵ lim is based on a normalized definition [see (20) or (22)] of SCF, which reaches near unity in stable situations, values in the range 10 −2 ÷ 10 −1 can be considered reasonable.A stability margin ϵ lim = 10 −2 was chosen in (21) for the following examples.This relatively small value, chosen for stability analysis method validation purposes, is coherent with a hypothesis of amplitude variations and uncertainties of the SCF (due to pseudorandom variations of technological parameters or other error sources) not greater than 1%, 3 The experiment results presented in this section were obtained by leaving line transitions (e.g., bends and T-junctions) unmodified in the PDC and Auxiliary PDC, whereas lines were damped in terms of attenuation blocks (see Fig. 8, bottom-right case, in Appendix A) after estimation of the corresponding delays.Accuracy in the estimation of starting frequencies of oscillations in all cases was presented to prove that the performance predictions were negligibly affected by these approximations.corresponding to very accurate implementation technology and component models.This choice allows the achieving of reliable prediction of both stable and unstable circuit behavior when using an adequately dense discretization grid in the f , σ D domain.Instead, larger values of ϵ lim (e.g., 0.05 or 0.1) could be considered as more appropriate when stability must be guaranteed for design purposes, by assuming larger margins for tolerances of implementation technologies and component models.
Three validation examples with increasing order of complexity in the required instability detection capabilities are presented.Despite the nominal frequency of operation of the examples being in the RF range, significant reactive effects interact in all cases with nonlinearity to determine the overall circuit behavior.

A. Single-Transistor Oscillator
The self-starting capability (i.e., instability) of a single-transistor hybrid oscillator designed at 430 MHz was verified by means of an LTI stability analysis.See the measured output spectrum and a photograph of the circuit implementation in Fig. 4(a) [34].
The SCF (20) was computed over a dense grid of values of the damping factor σ D and frequency f , with σ D,max = 1ns −1 .Fig. 4(b) shows the S C -plot in the relevant frequency range of 420 MHz ÷ 440 MHz, while a contour plot displaying iso-level lines of S C in the same region is shown in Fig. 4(c).As expected, the presence of a local minimum, whose value (3.11 • 10 −3 ) is less than ϵ lim = 10 −2 (corresponding to the smallest iso-level line shown in the plot), points out the presence of an unstable pole at about σ = σ rmp,LTI = 0.5 ns −1 and f = f rmp,LTI = 428 MHz in the grid of values considered.The measured output spectrum of the circuit confirms the steady-state oscillation at nearly the same frequency.

B. Single-Transistor PA
The stability analysis of a single-transistor hybrid largesignal amplifier, designed at the nominal frequency of 560 MHz, was considered [33].The output spectrum measurement (P in = 10 dBm) and a photograph of the circuit implementation is shown in Fig. 5(a).This circuit is LTI stable, but it shows LTV instability arising at half of the nominal frequency for larger input power levels [22], [35], [36].The three SCF graphs in Fig. 5(b)-(d) correspond to the three different input power levels: −35(small-signal), −4, and 7 dBm, respectively.Fig. 5(d) shows that at 7 dBm input power, the SCF has a local minimum, whose value (0.12 • 10 −3 ) is smaller than ϵ lim at the zero gradient point σ D = σ rmp,LTI = 0.1 ns −1 and f = f rmp,LTI = 288 MHz, i.e., nearly half of the operational frequency.This is coherent with the experiment results in Fig. 5(a).On the other hand, Fig. 5(c) shows a constrained minimum greater than ϵ lim at σ D nearly equal to zero, suggesting that the LSOP corresponding to the input power level −4 dBm is stable.
According to Fig. 5(b) and (c) no local minimum can be found in the same region of the LT plane at reduced input power (−4 dBm) and under small-signal operation (−35 dBm), which implies that, by lowering input power, the potentially unstable pole is shifted into the negative σ region.This suggests stable circuit operation, coherently with the experiment results.

C. Balanced Amplifier
As a final validation example, stability analysis was carried out for a hybrid paralleled-transistor balanced amplifier, theoretically designed for large-signal operation at the nominal frequency of 590 MHz.See the measured output spectrum at P in = 10 dBm and a photograph of the circuit implementation in Fig. 6(a) [26], which shows large-signal (LTV) instability at nearly half the nominal frequency [22], [35], [36].The stability analysis approach in [26] revealed the presence of an odd-mode kind of large-signal instability, involving signals flowing into the internal circuit loop.To demonstrate the ability of the SCF-based stability analysis to detect this sort of event, the SCF (22) (K = 8) was computed by PDC analysis based on two-tone mixer-mode HB simulation in the presence of a large-signal input tone at 590 MHz with the Fig. 5. PDC-based stability analysis of a single-transistor hybrid PA ( f 0 = 560 MHz).(a) Measurement of the output spectrum (P in = 10 dBm) and inset photograph from [33].SCF (b) under small-signal excitation (stable LTI-case); (c) at input power level −4 dBm (stable LTV case), and (d) at input power level 7 dBm (unstable LTV case) with detail of the constant-level contour plot (e).
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.PDC-based stability analysis of a paralleled-transistor hybrid balanced amplifier ( f 0 = 590 MHz).(a) Measurement of the output spectrum (P in = 10 dBm) and inset photograph from [26].SCF S of the balanced amplifier at input power levels b) −35 dBm and (c) 5 dBm with detail of the iso-level contour plot (d).Analysis confirms that the circuit is LTI stable but shows LTV instability at large-signal input power levels and 288 MHz.power level of 5 dBm.Fig. 6(c) and (d) shows the presence of a local minimum, whose value (0.46 • 10 −3 ) is smaller than ϵ lim at the zero gradients at σ D = σ rmp,LTI = 0.16 ns −1 and f = f rmp,LTI = 288 MHz, coherently with the experiment evidence.This corresponds to a pole with a positive real part, thus leading to LTV unstable circuit behavior.
A new set of SCFs ( 22) was evaluated at the input power level of −35 dBm (small-signal operation).The total absence of local minima in the σ D ≥ 0 plane [see Fig. 6(b)] confirms that, by lowering the signal amplitude, the unstable pole has Fig. 7. Flowchart describing the proposed analysis method as part of a general circuit design flow.According to the outlined procedure, the stability check could be treated run-time and together with other performance goals.If the method is used for a one-shot stability check of a given circuit, as in the present work, the feedback loop must be neglected.
been shifted into the left-hand side of the LT plane (LTI stable circuit).Similar conclusions were drawn in [26] for the same amplifier both experimentally and by means of an alternative stability analysis method, whose validity is however limited to symmetrical circuit topologies.
Coherently with the theoretical outcome outlined in Section III-C, the three experiment validation examples presented confirm that ( 21) is a necessary and sufficient constraint for circuit stability.In fact, in all the situations where ( 21

VI. INTEGRATION WITH CIRCUIT DESIGN PROCEDURES
The proposed method paves the way for the direct integration of stability constraints in the framework of iterative circuit design procedures carried out with commercially available CAD tools (e.g., [2], [3]).In fact, the evaluation of stability constraints (21), either based on SCFs (20) or (22), can be executed in the framework of commercial CAD tools together with the estimation of circuit performance goals and the verification of other constraints like reliability and physical feasibility.
According to the flowchart shown in Fig. 7, the library of PD component models for the chosen technology process is initially built and the design procedure gets started with an initial circuit design using unmodified component models.
As shown in Fig. 3, besides the circuit with unmodified component models, the perturbed PDC and auxiliary PDC need to be built for SCF-based stability checking.As outlined above, the netlist generation of these PD circuits could be automatically executed by a special-purpose software routine, based on the original layout with unmodified component models.Performance figures of merit and SCF evaluation over a dense grid in the f , σ D domain can now be evaluated through standard ac-and HB-analysis-based methods, as outlined in Section IV.In case performance goals and stability constraints were not satisfied, circuit variants could be checked iteratively after any manual topology change and/or optimization-driven automatic parameter value update.
For the validation examples presented in this work, PDC and Auxiliary PDC were manually built, and the flowchart was followed all the way down to the stability check without implementing the feedback loop, i.e., without the objective of further proceeding to any circuit modification.

VII. CONCLUSION
BIBO stability constraints for small and large-signal circuit analysis and design have been proposed.The method is compatible with circuits described by nonlinear time-domain differential equations, possibly including linear distributed components modeled by delay or convolution integral operators.The approach is compatible with commercial CAD tools used for circuit performance analysis.
Provided that a library of PD component models is preliminarily built for the adopted process technology, any steady-state solution of circuits, either single-or multitransistor-based, can be checked for stability by visual inspection of a single 2-D graph of the proposed scalar SCF S C ( f, σ D ).The latter can be obtained with the same analysis techniques adopted for circuit performance evaluation, i.e., either through frequency-domain ac (small-signal) or multitone quasi-periodic HB (large-signal) simulations with parametric sweeps of the perturbation source frequency f and of the damping parameter σ D .
The SCF expression consists of a normalized equation based on the elements of an NSP matrix, which links small-signal perturbations connected in series to all the intrinsic transistor ports with the corresponding port voltages.The BIBO stability constraints are evaluated in a systematic way and stability is checked by verifying that the function minimum is greater than a specified stability margin.This makes the method potentially compatible with automatic iterative circuit design procedures, a feature not easily shared by other methods in the literature.In fact, the creation of PD versions of the original circuit, needed for the method execution, could be made automatic through dedicated software routines embedded into CAD tools.
Single and multitransistor experiment validation examples confirm the validity of the approach under both small-and large-signal conditions.

APPENDIX A PDC BUILDING RULES AND PD MODELS OF LINEAR COMPONENTS
The PDC is built according to the following guidelines to guarantee coherence with the original problem.
1) The PDC has the same topology as the original circuit (in such a way that the Kirchhoff voltage and current laws are the same).2) energy-storage (i.e., frequency-dependent) circuit elements are replaced by associated PD components, whose behavior depends on the damping factor σ D .In particular, the PDC response is damped by acting on all the components that store energy either in the form of an electric field (like capacitors), magnetic field (inductors), or both (TLIN-based components, active devices).3) Memory-less component models remain unmodified.4) The steady-state solution v at the N P ports of the IEDs keep unmodified for each damping factor σ D value.For each component type, an associated PD component is now defined along with the corresponding model.

A. PD Linear Capacitor Model
A linear capacitor C is described in the Laplace domain by with Q( p) = C V ( p).According to the parametric damping criteria defined in Section IV-A, the associated PD capacitor model is defined as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where V( jω, σ D ) and I( jω, σ D ) are σ D -dependent voltage and current phasors and Q( jω, σ D ) = CV( jω, σ D ).Since (27) corresponds in the time domain to with q = Cv, the PD capacitor is described by the model shown in Fig. 8, where an additional resistor with conductance σ D C is connected in parallel to the standard capacitor.The current flowing through the resistor can be considered a damping current i DMP (t, σ D ) given by The FT-domain PD capacitor model ( 27) is clearly equivalent to the LT-domain model (26)

B. PD Linear Inductor Model
A linear inductor L is described in the Laplace domain by with ( p) = L I ( p).According to the parametric damping criteria defined in Section IV-A, the associated PD inductor model is defined as or in the time domain According to (32), the PD inductor model is obtained by connecting a resistor with a value σ D L in series to the standard inductor, as shown in Fig. 8.The voltage across the resistor can be considered as a damping voltage v DMP (t, σ D ) given by LT/FT-domain equivalence conditions are clearly satisfied by (30) and (31) for σ = σ D .

C. PD Transmission Line (TLIN) and Other Linear Distributed Component Models
Based on the PD models associated with linear capacitors and inductors, the σ D -dependent PD TLIN model can be easily obtained by adding series and parallel damping resistors to the inductor and capacitor per unit length in the standard LC model of the lossy line, as shown in Fig. 8.The same procedure can be followed for any lumped-component-based equivalent circuit of common lossy line transitions (e.g., bends and T-junctions).
More in general, model parametric damping can be also applied to linear component models whose dynamic behavior is described not only by state-space derivative operators, but also by delay operators, like pulse response functions with associated time-domain convolution integrals and, coherently, also PD transfer functions.
A general-purpose procedure is outlined here for the σ Ddependent PD-model construction of any linear distributed component described by the S-parameter matrix.S β ( jω).The frequency-domain S-matrix model may derive from direct measurements or, alternatively, by S-parameter simulation of black-box linear models provided by the component foundry.
By using numerically efficient inverse-FFT algorithms, or equivalently time-domain simulations of the single component, the scattering-wave pulse responses s β,i j (t) = F −1 {S β,i j ( jω)} for each i, j = 1, 2 can be calculated.As is well-known, these time-domain functions relate reflected and incident waves (typically defined by adopting 50 normalization) through causal convolution integrals, i.e., The PD elements of the lossy distributed component Sparameter matrix S β ( jω, σ D ) are then obtained, through DFT-based numerical methods, according to It is straightforwardly verified that S β ( p) = S β ( jω, σ D ) when σ = σ D .In fact, for each element of the matrix, it holds Here, the attenuation has been equivalently split into two equal factors (e −σ D τ/2 ) at the beginning and at the end of the TLIN, so preserving the line symmetry.

D. PD Model of the Linear Passive Network
According to the rules given, the associated PD passive linear network is created by replacing each elementary energy-storage component in the original linear passive network with its PD counterpart.Since the original network and the PD one share the same topology, the respective state space models (1) and ( 2) and ( 23) and ( 24) necessarily entail the same matrices A , B , C , D .Thanks to this feature, the Laplace-transformed model of (1) and ( 2 (47)

B. LTV Case
In where sums are limited to the maximum order K , were ĝk and Ĉk are complex Fourier coefficient matrices.By replacing (49) in (48) and L-transforming it is obtained Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Fig. 1 .
Fig. 1. Circuit to be checked for stability.A number N of electron devices, either one-or two-port are considered.Dashed lines indicate the second port of IEDs if this is present.

Fig. 3 .
Fig. 3. Circuit used for both performance and stability analysis using frequency-domain CAD tools.This is made of three parts on the same schematic page, namely, (a) PDC with perturbation sources, (b) auxiliary PDC, and (c) unperturbed original circuit.(c) Computes the unperturbed steady-state solution v. (b) evaluates voltages ṽ.(a) Allows for computation of the NSP matrix (14) (LTI case) or set of matrices (16) (LTV case).

Fig. 6 .
Fig. 6.PDC-based stability analysis of a paralleled-transistor hybrid balanced amplifier ( f 0 = 590 MHz).(a) Measurement of the output spectrum (P in = 10 dBm) and inset photograph from[26].SCF S of the balanced amplifier at input power levels b) −35 dBm and (c) 5 dBm with detail of the iso-level contour plot (d).Analysis confirms that the circuit is LTI stable but shows LTV instability at large-signal input power levels and 288 MHz.
D )} of the PDC are damped (or, equivalently, parametrically stabilized) versions of the original w(t) functions, i.e., w(t, σ D ) provided that σ D is chosen greater than the real part σ rmp,m,n of the rightmost pole of any function element in W ( p), i.e., σ D > σ rmp,LTI with σ rmp,LTI = σ rmp,m,n [RoC of W ( p)].In other words, according to (14) the FT-based NSP matrix W( jω, σ D ) of the PDC coincides with the LT-based W ( p) of the original circuit for each σ = σ D > σ rmp,LTI .It is also worth noticing that, thanks to (14), the pulse response functions w(t, σ D ) = F −1 {W( jω, σ the LTV case, linearization of the IED model (3) around any generic LSOP v(t) (stationary solution) leads to the time-variant linear relationship are small-amplitude deviations of instantaneous voltages with respect to v(t), with ĝ(t) = (∂ F(v)/∂v) v(t)= v(t) and Ĉ(t) = C( v(t)).By assuming a periodic regime with period T s ( s = 2π/T s ) the time-variant conductance ĝ(t) and capacitance Ĉ(t) can be expanded in the Fourier series as