# Circuit Distortion Analysis Based on the Simplified Newton's Method

Authors: M. M. Gourary; S. G. Rusakov; S. L. Ulyanov; M. M. Zharov; B. J. Mulvaney
Journal: Journal of Electrical and Computer Engineering
 (2011)

---

## Abstract

A new computational technique for distortion analysis of nonlinear circuits is presented. The new technique is applicable to the same class of circuits, namely, weakly nonlinear and time-varying circuits, as the periodic Volterra series. However, unlike the Volterra series, it does not require the computation of the second and third derivatives of device models. The new method is computationally efficient compared with a complete multitone nonlinear steady-state analysis such as harmonic balance. Moreover, the new technique naturally allows computing and characterizing the contributions of individual circuit components to the overall circuit distortion. This paper presents the theory of the new technique, a discussion of the numerical aspects, and numerical results.


---

## Body

## 1. Introduction

RF circuits are generally designed to be linear with respect to the signal path. However, the desired signal may be weakly distorted due to nonlinearities of the circuit components. Analyzing this nonlinear distortion is an important problem in the design of RF circuits [1, 2]. Another important practical task is the computation of distortion contributions due to individual nonlinearities of the circuit. This information is essential to determine critical circuit elements and improve the design.The traditional approach to measuring distortion is to apply one or more pure test tones to the circuit’s input and determine harmonics or intermodulation products at the output [3, 4]. This approach is applied to both time-invariant circuits (such as amplifiers) and circuits with periodic excitation (such as mixers).Consequently, distortion characteristics can be obtained using quasiperiodic steady-state nonlinear analysis such as harmonic balance technique [5–11] or mixed frequency time method [5, 12]. Such methods can be applied to circuits with strong nonlinear behavior but require essential computational efforts. Also these methods do not provide the computation of individual distortion contributions.For distortion analysis of time-invariant weakly nonlinear circuits, the approach based on Volterra series is more efficient [3, 4, 13–15]. According to this approach, each nonlinearity in the circuit is presented by a power series expansion up to the third order, and each term of expansion is associated with a source. The analysis leads to three linear systems, one for each order of nonlinear circuit response, that are solved successively beginning with the first order. This analysis implemented in Spice simulators (e.g., the DISTO analysis) provides the distortion characterization after the computation of DC operating point.The essential disadvantage of this approach is the necessity to compute power series coefficients for each nonlinearity in the circuit. This requires that the simulator has explicitly coded second- and third-order derivatives of the device models. With modern device models, which include nonlinear dependencies in several variables as well as different model behavior in different regions of device operation, computation of high-order derivatives is extremely difficult, if not impossible. Also this requirement limits the introduction of new models that may prevent the distortion analysis with a wide class of behavioral models.In comparison with conventional nonlinear distortion analysis based on computing distortion about the DC operating point, distortion analysis of communication circuits often requires determining the distortion about a periodically time-varying operating point. In this case, the input signal is considered as a small excursion about the periodically varying operating point, and the distortion characteristics can be obtained from the simulation of weakly nonlinear behavior with respect to the input signal [16]. Thus, these circuits must be considered as periodically time-varying systems. Conventional distortion analysis can be extended to periodically time-varying nonlinear circuits using time varying Volterra series [17]. However, the same difficulty with device models persists.Thus, the introduction of an efficient special-purpose distortion mode into circuit simulators is a current challenge in RF and microwave CAD engineering.We focus on the following practical problems for distortion analysis. First, in order to avoid the analytic computation of second and third derivatives of device models, we need an alternative to the Volterra series technique that retains the same overall accuracy. Second, we need to handle the problem of distortion of periodic time-varying circuits. The third practical problem we wish to address is the computation of the individual contributions of each component of the circuits to the overall distortion. This of course has important applications in guiding circuit design. The Volterra series technique supports such computation, so we need to retain this capability.This paper discusses a new approach for the distortion analysis [18] based on the simplified Newton’s method [19]. The computational efforts of the new approach are practically the same as for distortion analysis based on Volterra series, and radically lower in comparison with complete multitone nonlinear steady-state analysis. The new approach supports the computation of individual distortion contributions. A matrix form for the characterization of individual contributions is proposed that is well suited for circuit simulators, and the computational technique for extracting individual contributions is presented.It is important to note that the proposed approach does not require computing high-order derivatives of device model nonlinearities, and thus there is no need to code the second and third derivatives for all the device models. This approach provides basically the same order of accuracy as the Volterra series, due to properties of the simplified Newton’s method. The technique can be readily implemented in a general-purpose circuit simulator.A similar approach to distortion analysis has been proposed [20] which is based on linear-centric models and successive chord method. This likewise avoids the computation of second and third derivatives of device models but differs from the presented technique in computation of third-order and higher-distortion components, essential for determining the IM3 metric. The approach does not take into account the interaction of nonlinearities in the computation of distortion contributions.The paper is organized as follows. Section2 describes the formulation of the new approach for distortion analysis and presents the basic computational procedure. Section 3 discusses the application of this approach to the periodic distortion analysis of nonlinear circuits in the frequency domain in the framework of the harmonic balance method. Section 4 is devoted to the problem of computing individual contributions. Examples of computing distortion metrics and individual contributions for some typical circuits are given in Section 5.


## 2. Formulation of the Approach

### 2.1. Basic Model for Distortion Analysis

The distortion analysis is intended to provide a measure of the distortion products when one or more pure sinusoids are applied to the input of the circuit. In RF circuits, there may be an extra periodic excitation that determines the periodically varying operating point.In this case, the nonlinear circuit can be described by a system of differential-algebraic equations [1, 5](1)ddtq(x(t))+i(x(t))+ũ(t)+u(t)=0,
where q is the vector of charges, x is the vectors of nodal voltages, ũ(t) is a large periodic excitation with frequency ω̃, and u(t)=u1(t)+u2(t)+⋯+uM(t) is a vector of small periodic excitations corresponding to the set of independent sources with incommensurate frequencies ωi,i=1,…,M.Equation (1) corresponds to the problem of distortion analysis for periodically time-varying circuits. The basic model for conventional distortion analysis about DC operating point can be obtained as the particular case of (1) with ũ(t)=0, that is,(2)ddtq(x(t))+i(x(t))+u(t)=0.


### 2.2. Features of Volterra Series-Based Technique

The purpose of the nonlinear distortion analysis is the evaluation of the deviation of the circuit behavior from the desired linear behavior. There are effective methods to provide this evaluation, such as Volterra series [13, 15] or alternative method of direct calculation of nonlinear responses [3].These methods present the approximate solution of (1) or (2) in the form(3)x(t)≈x(n)(t)=x̃(t)+Δx(0)(t)+⋯+Δx(n-1)(t),
where x̃(t) is the periodic solution without small excitations, Δx(0)(t) is the linear response of the system, that is, the result of periodic AC analysis [1], and Δx(k)(t) (k>0) is the deviation of (k+1)th order with respect to the norm of small excitations ∥Δx(k)(t)∥=O(∥u(t)∥k+1). Here, the Landau symbol O means the asymptotic order.Each deviation is computed by solving the linear system in the frequency domain. The matrices of linear systems are obtained using linear AC analysis. The rhs vectors (equivalent current sources) are obtained by the computation of polynomials of previous-order deviations. The polynomials coefficients are defined by Taylor expansion of voltage-current and voltage-charge characteristics of circuit devices. The accuracy of the approximate solution (3) is characterized by the expression(4)‖x(n)(t)-x(t)‖=O(‖u(t)‖n+1).Usually it is assumed that nonlinearities are sufficiently small to neglect terms higher than third order. Thus, the analysis leads to successively solving the set of linear systems for each order beginning from the first order. The first-order solutionx(1)(t) is then used to calculate the equivalent current sources to determine the second-order solution x(2)(t). Both the first- and second-order solutions are used to compute the third-order solution x(3)(t). Thus, only three steps are needed to obtain third-order solution.The main advantage of Volterra series in the comparison with complete nonlinear analysis is the reduction of computational efforts. This results in a much smaller dimension of the small signal linear systems in the comparison with the dimension of the linear system solved at the iteration step of nonlinear analysis.The essential disadvantage of the Volterra series is the necessity of computing high-order power series coefficients for each nonlinearity in component models. This requirement limits the introduction of new device models because it leads to the very complex practical problem of analytically deriving high-order derivatives of the models.


### 2.3. Foundation of Computational Procedure Based on Simplified Newton’s Method

Below, we discuss the new approach that provides the similar accuracy (4) by the solving of the same number of linear systems. But unlike the Volterra series, this approach does not require complicated procedures to evaluate residuals in terms of high-order coefficients of Taylor expansions of circuit nonlinearities.Firstly, we can point out that any computational method to determine steady-state solution of (1) reduces the problem to the system of nonlinear algebraic equations(5)f(z)=w,
where w, z are N-dimensional vectors of input and internal signals in frequency or time domain. Let the vector z(0) be a solution of the problem obtained from (5) by neglecting small excitations, that is,(6)f(z(0))=0.
The nonlinear algebraic equation (5) are typically solved by the Newton’s method [19](7)J(z(k))Δz(k)=w-f(z(k)),z(k+1)=z(k)+Δz(k),k=0,1,2,…,
where(8)J(z(k))=[∂f1(z(k))∂z1⋯∂f1(z(k))∂zN⋮⋱⋮∂fN(z(k))∂z1⋯∂fN(z(k))∂zN]
is the Jacobian matrix of the system (5), and Δz(k) is the vector of corrections.Here, we apply the simplified Newton’s method [19] that is similar to the Newton’s method (7) except that the Jacobian matrix is kept fixed and equal to J(z(0)). Thus, the method does not require any update of the Jacobian matrix.The numerical scheme contains one Newton step and (k-1) steps of the simplified Newton’s method. For our consideration, it is important that the reduction of error can be estimated as follows [19, 21]: (9)‖z(k)-z*‖=O(‖z(0)-z*‖k+1),k=1,2,…,
where z* is the exact solution.Using the estimate∥z(0)-z*∥=O(∥w∥), we have from (9) for k=n that(10)‖z(n)-z*‖=O(‖w‖n+1).Thus, the simplified Newton’s method provides the same accuracy order as the Volterra series method (4) after the same number of steps and under the same limitations on the level of input amplitude. Hence, three steps of the method (7) are sufficient to provide the third-order solution.To obtain correctionsΔz(k) up to the third order, the computational procedure after the determination of the initial guess by (6) includes the following three steps:(1)
(11)J(z(0))Δz(0)=w,(2)
(12)J(z(0))Δz(1)=w-f(z(1)),(3)
(13)J(z(0))Δz(2)=w-f(z(2)).
Step 1 of this technique coincides with the Volterra series method while steps 2 and 3 differ by the evaluation of rhs vectors. Rather than expansion of nonlinear functions in power series (forming equivalent current sources), the rhs vectors for systems (12), (13) are computed by the same technique as nonlinear analysis, described in Section 3. As a result, the new approach avoids the computation of Taylor expansions of circuit nonlinearities.


### 2.4. Numerical Aspects

One shortcoming of the approach is connected with the possible growth of numerical noise in (12), (13) for cases of small deviations Δz(0), Δz(1). This is particularly critical for step 3. To decrease this effect, we propose another form for the rhs vector (13).The first-order Taylor expansion off(z(2)) yields(14)f(z(2))=f(z(1)+Δz(1))=f(z(1))+J(z(1))Δz(1)+O(‖Δz(1)‖2).Taking into account that∥Δz(1)∥=O(∥w∥2), we can evaluate O(∥Δz(1)∥2)=O(∥w∥4). So neglecting this term does not change the order of error at step 3 (see (10) for n=3). The substitution (14) into (13) leads to the modified expression of step 3(15)J(z(0))Δz(2)=w-f(z(1))-J(z(1))Δz(1).Subtracting (12) from (15), a more convenient form of this expression is obtained (16)J(z(0))Δz(2)=[J(z(0))-J(z(1))]Δz(1).
In many cases, it is efficient to perform step 2 by the equation (17)J(z(0))Δz(1)=w-φ,
where vector φ approximates f(z(1)) with the accuracy(18)‖φ-f(z(1))‖=O(‖w‖3).If step 2 is performed by (17) instead of (12), then previous considerations lead to the following expression instead of (16):(19)J(z(0))Δz(2)=[J(z(0))-J(z(1))]Δz(1)+φ-f(z(1)).Expression (16) is useful for the analysis of contributions. Expressions (17), (19) will be used in Section 3 to reduce the number of linear systems to be solved.


## 2.1. Basic Model for Distortion Analysis

The distortion analysis is intended to provide a measure of the distortion products when one or more pure sinusoids are applied to the input of the circuit. In RF circuits, there may be an extra periodic excitation that determines the periodically varying operating point.In this case, the nonlinear circuit can be described by a system of differential-algebraic equations [1, 5](1)ddtq(x(t))+i(x(t))+ũ(t)+u(t)=0,
where q is the vector of charges, x is the vectors of nodal voltages, ũ(t) is a large periodic excitation with frequency ω̃, and u(t)=u1(t)+u2(t)+⋯+uM(t) is a vector of small periodic excitations corresponding to the set of independent sources with incommensurate frequencies ωi,i=1,…,M.Equation (1) corresponds to the problem of distortion analysis for periodically time-varying circuits. The basic model for conventional distortion analysis about DC operating point can be obtained as the particular case of (1) with ũ(t)=0, that is,(2)ddtq(x(t))+i(x(t))+u(t)=0.


## 2.2. Features of Volterra Series-Based Technique

The purpose of the nonlinear distortion analysis is the evaluation of the deviation of the circuit behavior from the desired linear behavior. There are effective methods to provide this evaluation, such as Volterra series [13, 15] or alternative method of direct calculation of nonlinear responses [3].These methods present the approximate solution of (1) or (2) in the form(3)x(t)≈x(n)(t)=x̃(t)+Δx(0)(t)+⋯+Δx(n-1)(t),
where x̃(t) is the periodic solution without small excitations, Δx(0)(t) is the linear response of the system, that is, the result of periodic AC analysis [1], and Δx(k)(t) (k>0) is the deviation of (k+1)th order with respect to the norm of small excitations ∥Δx(k)(t)∥=O(∥u(t)∥k+1). Here, the Landau symbol O means the asymptotic order.Each deviation is computed by solving the linear system in the frequency domain. The matrices of linear systems are obtained using linear AC analysis. The rhs vectors (equivalent current sources) are obtained by the computation of polynomials of previous-order deviations. The polynomials coefficients are defined by Taylor expansion of voltage-current and voltage-charge characteristics of circuit devices. The accuracy of the approximate solution (3) is characterized by the expression(4)‖x(n)(t)-x(t)‖=O(‖u(t)‖n+1).Usually it is assumed that nonlinearities are sufficiently small to neglect terms higher than third order. Thus, the analysis leads to successively solving the set of linear systems for each order beginning from the first order. The first-order solutionx(1)(t) is then used to calculate the equivalent current sources to determine the second-order solution x(2)(t). Both the first- and second-order solutions are used to compute the third-order solution x(3)(t). Thus, only three steps are needed to obtain third-order solution.The main advantage of Volterra series in the comparison with complete nonlinear analysis is the reduction of computational efforts. This results in a much smaller dimension of the small signal linear systems in the comparison with the dimension of the linear system solved at the iteration step of nonlinear analysis.The essential disadvantage of the Volterra series is the necessity of computing high-order power series coefficients for each nonlinearity in component models. This requirement limits the introduction of new device models because it leads to the very complex practical problem of analytically deriving high-order derivatives of the models.


## 2.3. Foundation of Computational Procedure Based on Simplified Newton’s Method

Below, we discuss the new approach that provides the similar accuracy (4) by the solving of the same number of linear systems. But unlike the Volterra series, this approach does not require complicated procedures to evaluate residuals in terms of high-order coefficients of Taylor expansions of circuit nonlinearities.Firstly, we can point out that any computational method to determine steady-state solution of (1) reduces the problem to the system of nonlinear algebraic equations(5)f(z)=w,
where w, z are N-dimensional vectors of input and internal signals in frequency or time domain. Let the vector z(0) be a solution of the problem obtained from (5) by neglecting small excitations, that is,(6)f(z(0))=0.
The nonlinear algebraic equation (5) are typically solved by the Newton’s method [19](7)J(z(k))Δz(k)=w-f(z(k)),z(k+1)=z(k)+Δz(k),k=0,1,2,…,
where(8)J(z(k))=[∂f1(z(k))∂z1⋯∂f1(z(k))∂zN⋮⋱⋮∂fN(z(k))∂z1⋯∂fN(z(k))∂zN]
is the Jacobian matrix of the system (5), and Δz(k) is the vector of corrections.Here, we apply the simplified Newton’s method [19] that is similar to the Newton’s method (7) except that the Jacobian matrix is kept fixed and equal to J(z(0)). Thus, the method does not require any update of the Jacobian matrix.The numerical scheme contains one Newton step and (k-1) steps of the simplified Newton’s method. For our consideration, it is important that the reduction of error can be estimated as follows [19, 21]: (9)‖z(k)-z*‖=O(‖z(0)-z*‖k+1),k=1,2,…,
where z* is the exact solution.Using the estimate∥z(0)-z*∥=O(∥w∥), we have from (9) for k=n that(10)‖z(n)-z*‖=O(‖w‖n+1).Thus, the simplified Newton’s method provides the same accuracy order as the Volterra series method (4) after the same number of steps and under the same limitations on the level of input amplitude. Hence, three steps of the method (7) are sufficient to provide the third-order solution.To obtain correctionsΔz(k) up to the third order, the computational procedure after the determination of the initial guess by (6) includes the following three steps:(1)
(11)J(z(0))Δz(0)=w,(2)
(12)J(z(0))Δz(1)=w-f(z(1)),(3)
(13)J(z(0))Δz(2)=w-f(z(2)).
Step 1 of this technique coincides with the Volterra series method while steps 2 and 3 differ by the evaluation of rhs vectors. Rather than expansion of nonlinear functions in power series (forming equivalent current sources), the rhs vectors for systems (12), (13) are computed by the same technique as nonlinear analysis, described in Section 3. As a result, the new approach avoids the computation of Taylor expansions of circuit nonlinearities.


## 2.4. Numerical Aspects

One shortcoming of the approach is connected with the possible growth of numerical noise in (12), (13) for cases of small deviations Δz(0), Δz(1). This is particularly critical for step 3. To decrease this effect, we propose another form for the rhs vector (13).The first-order Taylor expansion off(z(2)) yields(14)f(z(2))=f(z(1)+Δz(1))=f(z(1))+J(z(1))Δz(1)+O(‖Δz(1)‖2).Taking into account that∥Δz(1)∥=O(∥w∥2), we can evaluate O(∥Δz(1)∥2)=O(∥w∥4). So neglecting this term does not change the order of error at step 3 (see (10) for n=3). The substitution (14) into (13) leads to the modified expression of step 3(15)J(z(0))Δz(2)=w-f(z(1))-J(z(1))Δz(1).Subtracting (12) from (15), a more convenient form of this expression is obtained (16)J(z(0))Δz(2)=[J(z(0))-J(z(1))]Δz(1).
In many cases, it is efficient to perform step 2 by the equation (17)J(z(0))Δz(1)=w-φ,
where vector φ approximates f(z(1)) with the accuracy(18)‖φ-f(z(1))‖=O(‖w‖3).If step 2 is performed by (17) instead of (12), then previous considerations lead to the following expression instead of (16):(19)J(z(0))Δz(2)=[J(z(0))-J(z(1))]Δz(1)+φ-f(z(1)).Expression (16) is useful for the analysis of contributions. Expressions (17), (19) will be used in Section 3 to reduce the number of linear systems to be solved.


## 3. Distortion Analysis in Frequency Domain

### 3.1. Distortion Analysis of Periodic Circuits

We now apply the simplified Newton’s method described in Section2 to compute distortion of a periodic circuit simulated with the harmonic balance method. This can be used in place of a multitone harmonic balance simulation in situations where some of the driving signals are small and subject to mild nonlinearity. The quasiperiodic harmonic balance (HB) method [5] transforms (1) to the following nonlinear algebraic system in the frequency domain(20)H(X)=I(X)+Ω⋅Q(X)+Ũ+∑i=1MUi=0,
where vectors X, I(X), Q(X), Ũ, and U are frequency domain representation of vectors x(t), i(x), q(x), ũ(t), and u(t), and Ω is a block-diagonal matrix of combinational frequencies of fundamentals ω̃, ωi.Solving the nonlinear system (20) by the Newton’s method results in the linear system for corrections ΔX(k),(21)JHB(X(k))ΔX(k)=-H(X(k)),
where JHB is the multitone harmonic Jacobian matrix [5].The vector of unknowns contains all harmonics of the linear combinations with integer coefficients of all fundamentalsω=k̃ω̃+ν, where ν denotes the small signal combinational frequency of the form(22)ν=k1⋅ω1+k2⋅ω2+⋯+kM⋅ωM.The number of complex unknowns in (20) or (21) is equal to N·(2K̃·K+K̃+K+1), where N is the circuit dimension, K̃ is the number of large signal harmonics, and K is the number of positive small signal combinational frequencies ν.Now, we consider the approximate solution of (20) by the simplified Newton’s procedure (7) with fixed Jacobian matrix.The initial step (6) is solving the reduced problem that is in our case a single tone HB problem(23)H̃(X̃)=Ĩ(X̃)+Ω̃⋅Q̃(X̃)+Ũ=0.
The nonlinear system (23) is obtained from (20) by applying only the large signal excitation.After determining the initial guessX(0)=X̃, one can apply the iterative process (7) that for system (20) is presented in the form(24)JHB(X̃)ΔX(k)=-H(X(k)).
The Jacobian matrix JHB(X̃) corresponds to (20) linearized about the single tone steady-state solution, and the system (24) presents the problem of the periodic small signal analysis. This problem can be decomposed into the set of low dimensional problems [22] for each small signal combinational frequency(22)(25)JSS(νj)ΔX(k)(νj)=-H(k)(νj).The number of unknowns in each system (25) is equal to N·(2K̃+1), and there are K systems with different νj.The entries of the matrixJSS(νj) can be defined by the Fourier transform of periodically varying conductance G(t) and capacitance C(t) matrices. But there is no need to obtain the small signal matrix in the frequency domain because systems (25) can be solved by efficient Krylov subspace methods with matrix-vector multiplication in the time domain [6].Thus, each step of the simplified Newton’s iterative process includes two numerical procedures:(i)
evaluation of the full rhs vectorH(X(k)) by the same algorithm that is applied in the solving of multitone HB problem,(ii)
solving the linear systems (25) with rhs vectors defined by the frequency decomposition of the full rhs vector for the required small signal frequency set.The evaluation of the full rhs vector for the modified step 3 (16) is performed by the corresponding HB expression(26)H(2)=[JHB(X(1))-JHB(X̃)]ΔX(1),
where JHB(X(1)) is harmonic Jacobian matrix (21) that is computed at the step 2 simultaneously with the computation of H(X(1)).


### 3.2. The Numerical Scheme

Before proceeding with the distortion analysis, it is first necessary to determine the harmonics that must be considered in the computations.The number of large signal harmonics is determined during the solving of the single-tone HB problem (23). To minimize the number of small signal harmonics that provides the accuracy order (10), we take into account the sinusoidal form of small input signals ui=Uiexp(±jωit). In such case, the minimal sets of small signal combinational frequencies at the nth step ν(n)={ν1(n),ν2(n),…} are defined by the following recurrent expressions:(27)ν(1)={…,ωi,…},ν(n+1)={…,|νj(n)±ωi|,…}.In the case of two small excitations, frequency sets (27) of second and third orders result in the forms(28)ν(2)={0,2ω1,2ω2,ω1+ω2,|ω1-ω2|},(29)ν(3)={ω1,ω2,2ω1+ω2,ω1+2ω2,|2ω1-ω2|,|2ω2-ω1|,3ω1,3ω2}.Note that at the last step, it is sufficient to determine only those harmonics of (29) that are required by the user. For example, if only intermodulation distortion is needed, then frequency set (29) is reduced to only one frequency ν(3)={|2ω1-ω2|}.Small signal systems at step 2 are solved only for second-order harmonics (28). It is equivalent to the performing of step 2 by (17), where vector φ represents the rhs vector with zero third-order harmonics (29). So in this case, the third step must be performed by (19), that is, vector (26) must be added by the magnitudes of third-order harmonics obtained from the rhs vector of step 2, H(1)(νj(3)).The computational procedure includes the following steps.Step 1 (Periodic large signal solution).
Solve the system (23) by the single-tone HB method.
Save in memory the solution plus the conductance and capacitance matrices obtained at the last Newton stepx̃(t),G̃(t), C̃(t).Step 2.
This step coincides with the standard periodic small signal analysis [23]. Here, the linear systems corresponding to the fundamentals ω1, ω2 are solved,
(30)JSS(ω1)ΔX(0)(ω1)=-U1,JSS(ω2)ΔX(0)(ω2)=-U2.Step 3.
(a) Determine the first-order solution in time domain by applying inverse Fourier transformΓ-1,
(31)x(1)(t)=x̃(t)+Γ-1ΔX(0).
(b) Compute time domain vectors of circuit charges, currents, and admittance matrices,(32)q(x(1)(t)),i(x(1)(t)),G(x(1)(t)),C(x(1)(t)).
(c) Compute the residual vector in the frequency domain by applying the Fourier transformΓ and using (20),(33)H(X(1))=Γi(x(1)(t))+ΩΓq(x(1)(t))+Ũ+U1+U2.
(d) Decompose the vectorH(X(1)) into H(1)(νj(2)) and H(1)(νj(3)).
(e) Determine corrections of solution by solving the systems(34)JSS(νj(2))ΔX(1)(νj(2))=-H(1)(νj(2)),
for all of 5 small signal combinational frequencies νj(2) from the set (28).Step 4.
(a) Transform the second-order correction to time domain by applying inverse Fourier transform(35)Δx(1)(t)=Γ-1ΔX(1).
(b) Compute vector (26) by matrix-vector multiplication in time domain and Fourier transform to frequency domain(36)H(2)=Γ[(G̃(t)-G(x(1)(t)))Δx(1)(t)]+jΩΓ[(C̃(t)-C(x(1)(t)))Δx(1)(t)].
(c) Decompose vectorH(2) into H(2)(νj(3)).
(d) Determine corrections by solving systems(37)JSS(νj(3))ΔX(2)(νj(3))=-H(2)(νj(3))-H(1)(νj(3)),
for all required third-order harmonics from (29).The conversion of signals from time to frequency domain and vise versa is performed using the multidimensional fast Fourier transform, which results in some aliasing error due to the finite number of sample points. Defining this number in accordance with the Nyquist frequency (≥2n+1), the aliasing error depends on harmonics higher than n. Thus, the required accuracy order (10) is provided for n=3 if the number of sample points is greater or equal to 7 for each small-signal fundamental.


### 3.3. DC Operating Point

The method for distortion analysis of weakly nonlinear circuits described by (2) can be obtained as a special case of the above method for periodically time-varying circuits.In this case, the initial step of the analysis is determining the vector of DC solution. Thus, the vectorX(0) contains N nonzero components corresponding to the DC solution vector, and the N×N Jacobian matrix corresponds to the AC matrix JSS(ν)=G+jνC, where matrices G and C are computed at the DC operating point. The computational procedure then follows the one presented above.


### 3.4. Comparison with Linear-Centric Approach

The linear-centric approach for distortion analysis of time-varying and weakly nonlinear circuits has been described in [20]. The main distinction in the numerical strategy of the approach presented above in comparison with the linear-centric approach is the following: our numerical procedure includes the third step of the recursive process (11)–(13) in the form given by expression (19).The linear-centric technique [20] also avoids the computation of high-order derivatives of circuit nonlinearities. The linear-centric model uses the successive chord iterative method, exploiting the constant Jacobian matrix for solving the nonlinear equations of harmonic balance, and requiring only one linear system solution to estimate distortion components.It can be mentioned that this computational procedure is equivalent in practice to two first steps of our approach. However, as it follows from (10), two steps do not provide the desired accuracy of computations. In the general case, solving the linear problem (13) is required to compute third-order distortion component with required accuracy.The following simple illustrative example demonstrates that the linear-centric model can lead to incorrect results due to limitations of number of solved linear problems.Let a nonlinear resistor with voltage-current relationI=V2 be excited by unit DC current source and a small sinusoidal signal. The corresponding equation has the form(38)v2-1-a⋅cos⁡(ωt)=0.
The DC solution of (38) is v(0)=1.The first-order correction is determined from the linear system obtained by linearization (38) (39)2Δv(0)=a⋅cos⁡(ωt).
Therefore, the first-order solution is expressed as(40)v(1)=v(0)+Δv(0)=1+a2⋅cos⁡(ωt).
The distortion contribution is determined from the linear system [20](41)2Δv=-[(1+a2⋅cos⁡(ωt))2-1-2⋅a2⋅cos⁡(ωt)].
After elementary transformations, we have (42)Δv=-a216-a216⋅cos⁡(2ωt).Thus, we see that the voltage distortion does not contain third harmonics. Here,Δv corresponds to the value Δv(1) from the recursive procedure.For this simple example, it is easy to obtain the solution of (38) in the form(43)v=1+a⋅cos⁡(ωt).This function contains third-order terms in its Taylor expansion (y=a·cos(ωt))(44)1+y=1+y2-y28+y316+O(y4),
which provide third harmonics. Note that if we perform the next step of our computational procedure, we obtain(45)Δv(2)=3a364cos⁡(ωt)+a364⋅cos⁡(3ωt)+O(a4).
This coincides with terms of the explicit Taylor expansion (44).


## 3.1. Distortion Analysis of Periodic Circuits

We now apply the simplified Newton’s method described in Section2 to compute distortion of a periodic circuit simulated with the harmonic balance method. This can be used in place of a multitone harmonic balance simulation in situations where some of the driving signals are small and subject to mild nonlinearity. The quasiperiodic harmonic balance (HB) method [5] transforms (1) to the following nonlinear algebraic system in the frequency domain(20)H(X)=I(X)+Ω⋅Q(X)+Ũ+∑i=1MUi=0,
where vectors X, I(X), Q(X), Ũ, and U are frequency domain representation of vectors x(t), i(x), q(x), ũ(t), and u(t), and Ω is a block-diagonal matrix of combinational frequencies of fundamentals ω̃, ωi.Solving the nonlinear system (20) by the Newton’s method results in the linear system for corrections ΔX(k),(21)JHB(X(k))ΔX(k)=-H(X(k)),
where JHB is the multitone harmonic Jacobian matrix [5].The vector of unknowns contains all harmonics of the linear combinations with integer coefficients of all fundamentalsω=k̃ω̃+ν, where ν denotes the small signal combinational frequency of the form(22)ν=k1⋅ω1+k2⋅ω2+⋯+kM⋅ωM.The number of complex unknowns in (20) or (21) is equal to N·(2K̃·K+K̃+K+1), where N is the circuit dimension, K̃ is the number of large signal harmonics, and K is the number of positive small signal combinational frequencies ν.Now, we consider the approximate solution of (20) by the simplified Newton’s procedure (7) with fixed Jacobian matrix.The initial step (6) is solving the reduced problem that is in our case a single tone HB problem(23)H̃(X̃)=Ĩ(X̃)+Ω̃⋅Q̃(X̃)+Ũ=0.
The nonlinear system (23) is obtained from (20) by applying only the large signal excitation.After determining the initial guessX(0)=X̃, one can apply the iterative process (7) that for system (20) is presented in the form(24)JHB(X̃)ΔX(k)=-H(X(k)).
The Jacobian matrix JHB(X̃) corresponds to (20) linearized about the single tone steady-state solution, and the system (24) presents the problem of the periodic small signal analysis. This problem can be decomposed into the set of low dimensional problems [22] for each small signal combinational frequency(22)(25)JSS(νj)ΔX(k)(νj)=-H(k)(νj).The number of unknowns in each system (25) is equal to N·(2K̃+1), and there are K systems with different νj.The entries of the matrixJSS(νj) can be defined by the Fourier transform of periodically varying conductance G(t) and capacitance C(t) matrices. But there is no need to obtain the small signal matrix in the frequency domain because systems (25) can be solved by efficient Krylov subspace methods with matrix-vector multiplication in the time domain [6].Thus, each step of the simplified Newton’s iterative process includes two numerical procedures:(i)
evaluation of the full rhs vectorH(X(k)) by the same algorithm that is applied in the solving of multitone HB problem,(ii)
solving the linear systems (25) with rhs vectors defined by the frequency decomposition of the full rhs vector for the required small signal frequency set.The evaluation of the full rhs vector for the modified step 3 (16) is performed by the corresponding HB expression(26)H(2)=[JHB(X(1))-JHB(X̃)]ΔX(1),
where JHB(X(1)) is harmonic Jacobian matrix (21) that is computed at the step 2 simultaneously with the computation of H(X(1)).


## 3.2. The Numerical Scheme

Before proceeding with the distortion analysis, it is first necessary to determine the harmonics that must be considered in the computations.The number of large signal harmonics is determined during the solving of the single-tone HB problem (23). To minimize the number of small signal harmonics that provides the accuracy order (10), we take into account the sinusoidal form of small input signals ui=Uiexp(±jωit). In such case, the minimal sets of small signal combinational frequencies at the nth step ν(n)={ν1(n),ν2(n),…} are defined by the following recurrent expressions:(27)ν(1)={…,ωi,…},ν(n+1)={…,|νj(n)±ωi|,…}.In the case of two small excitations, frequency sets (27) of second and third orders result in the forms(28)ν(2)={0,2ω1,2ω2,ω1+ω2,|ω1-ω2|},(29)ν(3)={ω1,ω2,2ω1+ω2,ω1+2ω2,|2ω1-ω2|,|2ω2-ω1|,3ω1,3ω2}.Note that at the last step, it is sufficient to determine only those harmonics of (29) that are required by the user. For example, if only intermodulation distortion is needed, then frequency set (29) is reduced to only one frequency ν(3)={|2ω1-ω2|}.Small signal systems at step 2 are solved only for second-order harmonics (28). It is equivalent to the performing of step 2 by (17), where vector φ represents the rhs vector with zero third-order harmonics (29). So in this case, the third step must be performed by (19), that is, vector (26) must be added by the magnitudes of third-order harmonics obtained from the rhs vector of step 2, H(1)(νj(3)).The computational procedure includes the following steps.Step 1 (Periodic large signal solution).
Solve the system (23) by the single-tone HB method.
Save in memory the solution plus the conductance and capacitance matrices obtained at the last Newton stepx̃(t),G̃(t), C̃(t).Step 2.
This step coincides with the standard periodic small signal analysis [23]. Here, the linear systems corresponding to the fundamentals ω1, ω2 are solved,
(30)JSS(ω1)ΔX(0)(ω1)=-U1,JSS(ω2)ΔX(0)(ω2)=-U2.Step 3.
(a) Determine the first-order solution in time domain by applying inverse Fourier transformΓ-1,
(31)x(1)(t)=x̃(t)+Γ-1ΔX(0).
(b) Compute time domain vectors of circuit charges, currents, and admittance matrices,(32)q(x(1)(t)),i(x(1)(t)),G(x(1)(t)),C(x(1)(t)).
(c) Compute the residual vector in the frequency domain by applying the Fourier transformΓ and using (20),(33)H(X(1))=Γi(x(1)(t))+ΩΓq(x(1)(t))+Ũ+U1+U2.
(d) Decompose the vectorH(X(1)) into H(1)(νj(2)) and H(1)(νj(3)).
(e) Determine corrections of solution by solving the systems(34)JSS(νj(2))ΔX(1)(νj(2))=-H(1)(νj(2)),
for all of 5 small signal combinational frequencies νj(2) from the set (28).Step 4.
(a) Transform the second-order correction to time domain by applying inverse Fourier transform(35)Δx(1)(t)=Γ-1ΔX(1).
(b) Compute vector (26) by matrix-vector multiplication in time domain and Fourier transform to frequency domain(36)H(2)=Γ[(G̃(t)-G(x(1)(t)))Δx(1)(t)]+jΩΓ[(C̃(t)-C(x(1)(t)))Δx(1)(t)].
(c) Decompose vectorH(2) into H(2)(νj(3)).
(d) Determine corrections by solving systems(37)JSS(νj(3))ΔX(2)(νj(3))=-H(2)(νj(3))-H(1)(νj(3)),
for all required third-order harmonics from (29).The conversion of signals from time to frequency domain and vise versa is performed using the multidimensional fast Fourier transform, which results in some aliasing error due to the finite number of sample points. Defining this number in accordance with the Nyquist frequency (≥2n+1), the aliasing error depends on harmonics higher than n. Thus, the required accuracy order (10) is provided for n=3 if the number of sample points is greater or equal to 7 for each small-signal fundamental.


## 3.3. DC Operating Point

The method for distortion analysis of weakly nonlinear circuits described by (2) can be obtained as a special case of the above method for periodically time-varying circuits.In this case, the initial step of the analysis is determining the vector of DC solution. Thus, the vectorX(0) contains N nonzero components corresponding to the DC solution vector, and the N×N Jacobian matrix corresponds to the AC matrix JSS(ν)=G+jνC, where matrices G and C are computed at the DC operating point. The computational procedure then follows the one presented above.


## 3.4. Comparison with Linear-Centric Approach

The linear-centric approach for distortion analysis of time-varying and weakly nonlinear circuits has been described in [20]. The main distinction in the numerical strategy of the approach presented above in comparison with the linear-centric approach is the following: our numerical procedure includes the third step of the recursive process (11)–(13) in the form given by expression (19).The linear-centric technique [20] also avoids the computation of high-order derivatives of circuit nonlinearities. The linear-centric model uses the successive chord iterative method, exploiting the constant Jacobian matrix for solving the nonlinear equations of harmonic balance, and requiring only one linear system solution to estimate distortion components.It can be mentioned that this computational procedure is equivalent in practice to two first steps of our approach. However, as it follows from (10), two steps do not provide the desired accuracy of computations. In the general case, solving the linear problem (13) is required to compute third-order distortion component with required accuracy.The following simple illustrative example demonstrates that the linear-centric model can lead to incorrect results due to limitations of number of solved linear problems.Let a nonlinear resistor with voltage-current relationI=V2 be excited by unit DC current source and a small sinusoidal signal. The corresponding equation has the form(38)v2-1-a⋅cos⁡(ωt)=0.
The DC solution of (38) is v(0)=1.The first-order correction is determined from the linear system obtained by linearization (38) (39)2Δv(0)=a⋅cos⁡(ωt).
Therefore, the first-order solution is expressed as(40)v(1)=v(0)+Δv(0)=1+a2⋅cos⁡(ωt).
The distortion contribution is determined from the linear system [20](41)2Δv=-[(1+a2⋅cos⁡(ωt))2-1-2⋅a2⋅cos⁡(ωt)].
After elementary transformations, we have (42)Δv=-a216-a216⋅cos⁡(2ωt).Thus, we see that the voltage distortion does not contain third harmonics. Here,Δv corresponds to the value Δv(1) from the recursive procedure.For this simple example, it is easy to obtain the solution of (38) in the form(43)v=1+a⋅cos⁡(ωt).This function contains third-order terms in its Taylor expansion (y=a·cos(ωt))(44)1+y=1+y2-y28+y316+O(y4),
which provide third harmonics. Note that if we perform the next step of our computational procedure, we obtain(45)Δv(2)=3a364cos⁡(ωt)+a364⋅cos⁡(3ωt)+O(a4).
This coincides with terms of the explicit Taylor expansion (44).


## 4. Computation of Individual Contributions

### 4.1. Contributions Characterization

An important problem in nonlinear distortion analysis is the computation of the contributions of each nonlinear component to the total circuit distortion. This information allows designers to determine which circuit elements are responsible for the distortion, thus providing guidance in meeting the required specifications.Numerical characterization of contributions is to be achieved if one can evaluate the distortion metric in the vector form(46)D=[D1,D2,…,Dm,…],
with the following properties:(1)
mth component of the vector depends on mth nonlinearity and does not depend on any other nonlinearity,(2)
total distortion is the sum of all components of the vector(47)Dtotal=∑mDm.However, for third-order distortion, this approach is incomplete because various nonlinearities can interact, and this interaction prevents different contributions from being associated with only one nonlinearity [3], as illustrated in the next example.Consider the two-stage amplifier shown in Figure1. Each stage of the amplifier is described by a transfer factor (K) with small nonlinearities of second (a) and third (b) order (48)vout1(vinp1)=K1vinp1+a1vinp12+b1vinp13,vout2(vinp2)=K2vinp2+a2vinp22+b2vinp23.
The circuit contains two nonlinearities defined by parameters a1, b1 for the first nonlinearity and a2, b2 for the second nonlinearity.Figure 1
Two-stage amplifier.The output signal of the amplifier is obtained by(49)vout(v)=vout2(vout1(v)).
After substituting (48) into (49) and neglecting terms higher than third order, this expression can be presented in the form(50)vout(v)=vout(1)+vout(2)+vout(3),
where vout(1)=K1K2v, vout(2)=K2a1v2+K12a2v2, and vout(3)=K2b1v3+K13b2v3+2K1a1a2v3.One can see that the expression for the second-order distortionvout(2) is the sum of two terms, each of them corresponds to only one nonlinearity (a1 or a2) and so can be considered as the contribution of this nonlinearity. Thus, the second-order contributions D(2) can be presented in the form (46) with D1(2)=K2a1v2, D2(2)=K12a2v2.In contrast, the expression for the third-order distortion contains not only terms corresponding to each nonlinearity (K2b1v3,K13b2v3) but also the term 2K1a1a2v3 depending on both nonlinearities.Hence, the influences of circuit nonlinearities on the third-order distortion can be presented as a matrix of contributions(51)D(3)=[D11(3)0D21(3)D22(3)],
where D11(3)=K2b1v3, and D22(3)=K13b2v3, D21(3)=2K1a1a2v3.The meaning of the nondiagonal termD21(3) can be more clear if it is presented as a function of the two output signals of the first stage(52)D21(3)=2a2vout1(1)vout1(2),
where vout1(1)=K1v is linear output signal, and vout1(2)=a1v2=a1vv is the second-order distortion signal produced by the first nonlinearity. Then we can interpret D21(3) (52) as a result of mixing signals vout1(1), vout1(2) by the nonlinearity of the stage 2.Here, the matrixD(3) is triangular because the circuit is unidirectional, but if feedback loop is present, then D12(3)≠0. To see this, consider the next example shown in Figure 2.Figure 2
The amplifier with feedback loop.The transfer functionvout(v) for the circuit is defined by the implicit expression vout(v)=vout1(v-vout2(vout(v))). After substituting (48) and (49) into this expression and equating terms of equal orders, one can obtain the following expressions for output linear signal, second- and third-order distortion contributions:(53)vout1(1)=K1v1+K1K2,(54)D1(2)=a1v2(1+K1K2)3,D2(2)=-a2K13v2(1+K1K2)3,(55)D11(3)=-2a1D1(2)K2v3(1+K1K2)2+b1v3(1+K1K2)4,D22(3)=-2a2D2(2)K12K2v3(1+K1K2)2+-b2K14K23v3(1+K1K2)4,(56)D12(3)=2a1D2(2)v3K1(1+K1K2)2,D21(3)=-2a2D1(2)K12K2v3(1+K1K2)2.The similar generation of third-order distortions by the mixing of linear and second-order signals is present in any nonlinear circuit (as can be proved by Volterra series theory). So the third-order distortion contributions of the circuit can be characterized by the square matrixD(3) with entries Dnm(3).The diagonal entryDnn(3) defines the third-order signal that is produced in nth nonlinearity by the mixing of the linear signal with the second-order signal from the same nonlinearity (first term in (55)) and the mixing of linear signals by the third order coefficient of the nonlinearity (second term in (55)). The nondiagonal entry Dnm(3) defines the third-order signal that is produced in nth nonlinearity by the mixing of linear signals with the second-order signals produced in mth nonlinearity.The following properties of the contribution matrix exist.(1)
The sum of all matrix entries is equal to the total third-order distortion(57)Dtotal(3)=∑n∑mDnm(3).(2)
If parameters of thenth nonlinearity are changed, then only entries of nth column and nth row are modified. All other entries remain unchanged. In particular, if parameters of the nth nonlinearity are set to zero, then all entries of nth column and nth row are zeros. All other entries remain unchanged.(3)
The sum of entries ofnth column and nth row(58)Dn(3)=∑m≠n(Dnm(3)+Dmn(3))+Dnn(3)defines the total influence ofnth nonlinearity on the third-order distortion. It is equal to the difference between the total distortion and the distortion obtained without this nonlinearity.  Note that value Dn(3) cannot be considered as the contribution of nth nonlinearity in the sense (46) because property (47) is not true.(4)
If the circuit contains a group of nonlinearitiesS={i1,i2,…} and it is desired to present the group as one nonlinearity with index k, then the matrix is obtained using the following expressions:(59)Dkm(3)=∑n∈SDnm(3),Dnk(3)=∑m∈SDnm(3),Dkk(3)=∑n∈S∑m∈SDnm(3).Thus, the contribution matrix allows one to obtain various information on the influence of circuit nonlinearities upon the output nonlinear distortion.The utility of the nondiagonal contribution entries can be illustrated by the following way. If the main contribution is defined by diagonal entry, then the overall circuit distortion can be improved by reducing of corresponding nonlinearity. But if such entry is nondiagonal oneDnm(3), then a designer can also change frequency characteristics of the signal path from mth to nth nonlinearity to attenuate the second-order signal through the path.


### 4.2. Evaluation of Contributions by Simplified Newton’s Method

For the determination of individual contributions, the functional description of circuit modelf(z) (5) must be presented as the sum of nonlinear dependencies (60)f(z)=∑mf(m)(z),
corresponding to the defined nonlinearities. The nonlinear dependencies can be, for example, currents of circuit components or individual terminal currents or even separate capacitance and DC currents. The representation (60) is provided in circuit simulators at the step of equation formulation from models of individual nonlinearities, with interconnection equations obtained by using Kirchhoff’s current law.In accordance with (60), the Jacobian matrix of (5) also can be presented as the sum of corresponding matrices of the partial derivatives(61)J(z)=∑mJ(m)(z),J(m)(z)=∂f(m)(z)∂z.For the representations (60), (61), step 2 (12) can be written in the form(62)J(z(0))Δz(1)=∑m(J(m)(z(0))Δz(0)-f(m)(z(1))).Here, we exploit the fact that from (11), we have(63)∑m(J(m)(z(0))Δz(0))=J(z(0))Δz(0)=w.Each term in the residual vector of (62) defines the effect of mth nonlinearity because it presents the difference between the linearized J(m)(z(0))Δz(0) and the nonlinear f(m)(z(1)) dependencies of the nonlinearity. Hence, the second-order distortion can be written as a sum of individual contributions(64)Δz(1)=∑mΔz(1,m),Δz(1,m)=J-1(z(0))(J(m)(z(0))Δz(0)-f(m)(z(1))).To obtain third-order distortion as the sum of individual contributions, we use a modified expression (16) for step 3. After substituting (62), (64) into (16), we obtain(65)J(z(0))Δz(2)=∑n(J(n)(z(0))-J(n)(z(1)))∑mΔz(1,m).Thus, the third-order distortion is(66)Δz(2)=∑n∑mΔz(2,n,m),Δz(2,n,m)=J-1(z(0))(J(n)(z(0))-J(n)(z(1)))Δz(1,m).
Each vector (66) depends only on two nonlinearities m, n.Usually the distortion at the output node is of interest. The components of vectors (64) corresponding to the output node form a vector of second-order contributions Dm(2)=Δzout(1,m), and corresponding components of the vectors (66) form the matrix of third-order contributions Dnm(3)=Δzout(2,n,m), where subscript out denotes the output node.The algorithm for obtaining individual contributions is based on expressions (64) and (66), which may be represented in the frequency domain (in accordance with (34) and (37)) as (67)JSS(νj(2))ΔX(1,m)(νj(2))=-H(1,m)(νj(2)),(68)JSS(νj(3))ΔX(2,n,m)(νj(3))=-H(2,n,m)(νj(3))-δnmH(1,n)(νj(3)),
where δnm is the Kronecker delta. In accordance with (26), H(2,n,m)=[JHB(n)(X(1))-JHB(n)(X̃)]ΔX(1,m).The computations can be performed in the same numerical framework presented in Section3. To obtain vectors H(1,m)(X(1)) and matrices JHB(n)(X(1)), the HB loading procedure must be developed on the per-nonlinearity basis.LetNN be the number of nonlinearities in the circuit, and K2, K3 are the numbers of second- and third-order small signal frequencies, respectively. Then the computation of second-order contributions requires solving K2NN linear systems (67) with only K2 different matrices. For the third-order contribution matrix K3NN2 linear systems (68) should be solved. Fortunately, using the technique of adjoint linear system, this work is reduced to solving K3 linear system with transposed matrix and K3NN2 computations of dot products.


## 4.1. Contributions Characterization

An important problem in nonlinear distortion analysis is the computation of the contributions of each nonlinear component to the total circuit distortion. This information allows designers to determine which circuit elements are responsible for the distortion, thus providing guidance in meeting the required specifications.Numerical characterization of contributions is to be achieved if one can evaluate the distortion metric in the vector form(46)D=[D1,D2,…,Dm,…],
with the following properties:(1)
mth component of the vector depends on mth nonlinearity and does not depend on any other nonlinearity,(2)
total distortion is the sum of all components of the vector(47)Dtotal=∑mDm.However, for third-order distortion, this approach is incomplete because various nonlinearities can interact, and this interaction prevents different contributions from being associated with only one nonlinearity [3], as illustrated in the next example.Consider the two-stage amplifier shown in Figure1. Each stage of the amplifier is described by a transfer factor (K) with small nonlinearities of second (a) and third (b) order (48)vout1(vinp1)=K1vinp1+a1vinp12+b1vinp13,vout2(vinp2)=K2vinp2+a2vinp22+b2vinp23.
The circuit contains two nonlinearities defined by parameters a1, b1 for the first nonlinearity and a2, b2 for the second nonlinearity.Figure 1
Two-stage amplifier.The output signal of the amplifier is obtained by(49)vout(v)=vout2(vout1(v)).
After substituting (48) into (49) and neglecting terms higher than third order, this expression can be presented in the form(50)vout(v)=vout(1)+vout(2)+vout(3),
where vout(1)=K1K2v, vout(2)=K2a1v2+K12a2v2, and vout(3)=K2b1v3+K13b2v3+2K1a1a2v3.One can see that the expression for the second-order distortionvout(2) is the sum of two terms, each of them corresponds to only one nonlinearity (a1 or a2) and so can be considered as the contribution of this nonlinearity. Thus, the second-order contributions D(2) can be presented in the form (46) with D1(2)=K2a1v2, D2(2)=K12a2v2.In contrast, the expression for the third-order distortion contains not only terms corresponding to each nonlinearity (K2b1v3,K13b2v3) but also the term 2K1a1a2v3 depending on both nonlinearities.Hence, the influences of circuit nonlinearities on the third-order distortion can be presented as a matrix of contributions(51)D(3)=[D11(3)0D21(3)D22(3)],
where D11(3)=K2b1v3, and D22(3)=K13b2v3, D21(3)=2K1a1a2v3.The meaning of the nondiagonal termD21(3) can be more clear if it is presented as a function of the two output signals of the first stage(52)D21(3)=2a2vout1(1)vout1(2),
where vout1(1)=K1v is linear output signal, and vout1(2)=a1v2=a1vv is the second-order distortion signal produced by the first nonlinearity. Then we can interpret D21(3) (52) as a result of mixing signals vout1(1), vout1(2) by the nonlinearity of the stage 2.Here, the matrixD(3) is triangular because the circuit is unidirectional, but if feedback loop is present, then D12(3)≠0. To see this, consider the next example shown in Figure 2.Figure 2
The amplifier with feedback loop.The transfer functionvout(v) for the circuit is defined by the implicit expression vout(v)=vout1(v-vout2(vout(v))). After substituting (48) and (49) into this expression and equating terms of equal orders, one can obtain the following expressions for output linear signal, second- and third-order distortion contributions:(53)vout1(1)=K1v1+K1K2,(54)D1(2)=a1v2(1+K1K2)3,D2(2)=-a2K13v2(1+K1K2)3,(55)D11(3)=-2a1D1(2)K2v3(1+K1K2)2+b1v3(1+K1K2)4,D22(3)=-2a2D2(2)K12K2v3(1+K1K2)2+-b2K14K23v3(1+K1K2)4,(56)D12(3)=2a1D2(2)v3K1(1+K1K2)2,D21(3)=-2a2D1(2)K12K2v3(1+K1K2)2.The similar generation of third-order distortions by the mixing of linear and second-order signals is present in any nonlinear circuit (as can be proved by Volterra series theory). So the third-order distortion contributions of the circuit can be characterized by the square matrixD(3) with entries Dnm(3).The diagonal entryDnn(3) defines the third-order signal that is produced in nth nonlinearity by the mixing of the linear signal with the second-order signal from the same nonlinearity (first term in (55)) and the mixing of linear signals by the third order coefficient of the nonlinearity (second term in (55)). The nondiagonal entry Dnm(3) defines the third-order signal that is produced in nth nonlinearity by the mixing of linear signals with the second-order signals produced in mth nonlinearity.The following properties of the contribution matrix exist.(1)
The sum of all matrix entries is equal to the total third-order distortion(57)Dtotal(3)=∑n∑mDnm(3).(2)
If parameters of thenth nonlinearity are changed, then only entries of nth column and nth row are modified. All other entries remain unchanged. In particular, if parameters of the nth nonlinearity are set to zero, then all entries of nth column and nth row are zeros. All other entries remain unchanged.(3)
The sum of entries ofnth column and nth row(58)Dn(3)=∑m≠n(Dnm(3)+Dmn(3))+Dnn(3)defines the total influence ofnth nonlinearity on the third-order distortion. It is equal to the difference between the total distortion and the distortion obtained without this nonlinearity.  Note that value Dn(3) cannot be considered as the contribution of nth nonlinearity in the sense (46) because property (47) is not true.(4)
If the circuit contains a group of nonlinearitiesS={i1,i2,…} and it is desired to present the group as one nonlinearity with index k, then the matrix is obtained using the following expressions:(59)Dkm(3)=∑n∈SDnm(3),Dnk(3)=∑m∈SDnm(3),Dkk(3)=∑n∈S∑m∈SDnm(3).Thus, the contribution matrix allows one to obtain various information on the influence of circuit nonlinearities upon the output nonlinear distortion.The utility of the nondiagonal contribution entries can be illustrated by the following way. If the main contribution is defined by diagonal entry, then the overall circuit distortion can be improved by reducing of corresponding nonlinearity. But if such entry is nondiagonal oneDnm(3), then a designer can also change frequency characteristics of the signal path from mth to nth nonlinearity to attenuate the second-order signal through the path.


## 4.2. Evaluation of Contributions by Simplified Newton’s Method

For the determination of individual contributions, the functional description of circuit modelf(z) (5) must be presented as the sum of nonlinear dependencies (60)f(z)=∑mf(m)(z),
corresponding to the defined nonlinearities. The nonlinear dependencies can be, for example, currents of circuit components or individual terminal currents or even separate capacitance and DC currents. The representation (60) is provided in circuit simulators at the step of equation formulation from models of individual nonlinearities, with interconnection equations obtained by using Kirchhoff’s current law.In accordance with (60), the Jacobian matrix of (5) also can be presented as the sum of corresponding matrices of the partial derivatives(61)J(z)=∑mJ(m)(z),J(m)(z)=∂f(m)(z)∂z.For the representations (60), (61), step 2 (12) can be written in the form(62)J(z(0))Δz(1)=∑m(J(m)(z(0))Δz(0)-f(m)(z(1))).Here, we exploit the fact that from (11), we have(63)∑m(J(m)(z(0))Δz(0))=J(z(0))Δz(0)=w.Each term in the residual vector of (62) defines the effect of mth nonlinearity because it presents the difference between the linearized J(m)(z(0))Δz(0) and the nonlinear f(m)(z(1)) dependencies of the nonlinearity. Hence, the second-order distortion can be written as a sum of individual contributions(64)Δz(1)=∑mΔz(1,m),Δz(1,m)=J-1(z(0))(J(m)(z(0))Δz(0)-f(m)(z(1))).To obtain third-order distortion as the sum of individual contributions, we use a modified expression (16) for step 3. After substituting (62), (64) into (16), we obtain(65)J(z(0))Δz(2)=∑n(J(n)(z(0))-J(n)(z(1)))∑mΔz(1,m).Thus, the third-order distortion is(66)Δz(2)=∑n∑mΔz(2,n,m),Δz(2,n,m)=J-1(z(0))(J(n)(z(0))-J(n)(z(1)))Δz(1,m).
Each vector (66) depends only on two nonlinearities m, n.Usually the distortion at the output node is of interest. The components of vectors (64) corresponding to the output node form a vector of second-order contributions Dm(2)=Δzout(1,m), and corresponding components of the vectors (66) form the matrix of third-order contributions Dnm(3)=Δzout(2,n,m), where subscript out denotes the output node.The algorithm for obtaining individual contributions is based on expressions (64) and (66), which may be represented in the frequency domain (in accordance with (34) and (37)) as (67)JSS(νj(2))ΔX(1,m)(νj(2))=-H(1,m)(νj(2)),(68)JSS(νj(3))ΔX(2,n,m)(νj(3))=-H(2,n,m)(νj(3))-δnmH(1,n)(νj(3)),
where δnm is the Kronecker delta. In accordance with (26), H(2,n,m)=[JHB(n)(X(1))-JHB(n)(X̃)]ΔX(1,m).The computations can be performed in the same numerical framework presented in Section3. To obtain vectors H(1,m)(X(1)) and matrices JHB(n)(X(1)), the HB loading procedure must be developed on the per-nonlinearity basis.LetNN be the number of nonlinearities in the circuit, and K2, K3 are the numbers of second- and third-order small signal frequencies, respectively. Then the computation of second-order contributions requires solving K2NN linear systems (67) with only K2 different matrices. For the third-order contribution matrix K3NN2 linear systems (68) should be solved. Fortunately, using the technique of adjoint linear system, this work is reduced to solving K3 linear system with transposed matrix and K3NN2 computations of dot products.


## 5. Numerical Results

### 5.1. Comparison with Harmonic Balance

We compare our new approach for the periodic distortion analysis based on HB method (PDHB) with complete multitone HB analysis. Numerical experiments are performed for two MOSFET mixers. The first circuit is the single-balanced mixer containing 3 MOS transistors. The first mixer has input parametersALO = 2.5 V, fLO = 375 MHz, f1 = 394.9 MHz, and f2 = 395.1 MHz, with rf input varying from 0 to 0.5 V. The second circuit is the double-balanced mixer [20] containing 6 MOS transistors and biasing circuitry. The second mixer has parameters ALO = 0.75 V, fLO = 985 MHz, f1 = 899 MHz, and f2 = 901 MHz, with rf input varying from 0 to 0.15 V. In both circuits, f1 and f2 are treated as small as signals in the PDHB analysis, whereas the HB analysis must be a full 3-tone simulation.The computed dependencies of the output power of the frequency component (fLO-f1) on the input power are shown in Figure 3 for the first mixer. The obtained curves are in good agreement, especially for small rf input amplitudes as expected. Similar results are obtained for the second mixer. The dependencies of the output power at the third-order intermodulation product (fLO-2f1+f2) are shown in Figure 4 for the first mixer. The obtained curves are in good agreement for practical input amplitudes.Figure 3
Output power of frequency component (fLO-f1) versus input power.Figure 4
Output power of frequency component (fLO-2f1+f2) versus input power.The computational efficiency of the PDHB approach for the two mixer circuits can be seen from Table1. The number of HB harmonics for the three tones fLO, f1, f2 is (15, 3, 3) and (7, 3, 3) for the first and second mixer, respectively. The corresponding orders of linear systems for the two different methods are given in columns 3, 4. The number of linear systems to be solved is 8 and 18 in HB and PDHB analyses, respectively, for the first example. The second mixer requires solving 9 and 19 linear systems. In spite of increased number of linear systems, the CPU time is less in PDHB analysis due to the much smaller order of systems.Table 1CircuitLinear system orderCPU time, secNumberNHBPDHBHBPDHB1309114018603601325079380150039915


### 5.2. Examples of Individual Contributions

The first example is the two-stage bipolar amplifier shown in Figure5. The distortion analysis is applied to compute second- and third-order harmonic distortion. Figure 6 shows the third harmonic together with contributions from different transistors as a function of frequency. It is seen from the figure that the most dominant contribution is due to the transistor of the output stage. The contribution of the first-stage transistor q1 is smaller by a factor of 8–12, and the contribution due to the interaction of transistors q1 and q2 is negligible in this simple example. This relation is changed with introducing feedback loop in the amplifier.Figure 5
Two-stage bipolar amplifier.Figure 6
Third harmonic and its contributions.The second example is an OpAmp UA741 [24]. The transistors q1, q2, q3 belong to the input stage of the amplifier, and transistors q14, q20 are from the output stage. Figure 7 shows the dominant contributions related with these transistors to third harmonic as a function of frequency. Note that contributions of the output stage transistors are dominated over input stage transistors at low frequencies while they become the same order at high frequencies. Also note that contributions due to interaction between transistors q14, q20 are comparable with contributions of q14 and q20 individually and therefore cannot be neglected in this example.Figure 7
Contributions to the third harmonic.


## 5.1. Comparison with Harmonic Balance

We compare our new approach for the periodic distortion analysis based on HB method (PDHB) with complete multitone HB analysis. Numerical experiments are performed for two MOSFET mixers. The first circuit is the single-balanced mixer containing 3 MOS transistors. The first mixer has input parametersALO = 2.5 V, fLO = 375 MHz, f1 = 394.9 MHz, and f2 = 395.1 MHz, with rf input varying from 0 to 0.5 V. The second circuit is the double-balanced mixer [20] containing 6 MOS transistors and biasing circuitry. The second mixer has parameters ALO = 0.75 V, fLO = 985 MHz, f1 = 899 MHz, and f2 = 901 MHz, with rf input varying from 0 to 0.15 V. In both circuits, f1 and f2 are treated as small as signals in the PDHB analysis, whereas the HB analysis must be a full 3-tone simulation.The computed dependencies of the output power of the frequency component (fLO-f1) on the input power are shown in Figure 3 for the first mixer. The obtained curves are in good agreement, especially for small rf input amplitudes as expected. Similar results are obtained for the second mixer. The dependencies of the output power at the third-order intermodulation product (fLO-2f1+f2) are shown in Figure 4 for the first mixer. The obtained curves are in good agreement for practical input amplitudes.Figure 3
Output power of frequency component (fLO-f1) versus input power.Figure 4
Output power of frequency component (fLO-2f1+f2) versus input power.The computational efficiency of the PDHB approach for the two mixer circuits can be seen from Table1. The number of HB harmonics for the three tones fLO, f1, f2 is (15, 3, 3) and (7, 3, 3) for the first and second mixer, respectively. The corresponding orders of linear systems for the two different methods are given in columns 3, 4. The number of linear systems to be solved is 8 and 18 in HB and PDHB analyses, respectively, for the first example. The second mixer requires solving 9 and 19 linear systems. In spite of increased number of linear systems, the CPU time is less in PDHB analysis due to the much smaller order of systems.Table 1CircuitLinear system orderCPU time, secNumberNHBPDHBHBPDHB1309114018603601325079380150039915


## 5.2. Examples of Individual Contributions

The first example is the two-stage bipolar amplifier shown in Figure5. The distortion analysis is applied to compute second- and third-order harmonic distortion. Figure 6 shows the third harmonic together with contributions from different transistors as a function of frequency. It is seen from the figure that the most dominant contribution is due to the transistor of the output stage. The contribution of the first-stage transistor q1 is smaller by a factor of 8–12, and the contribution due to the interaction of transistors q1 and q2 is negligible in this simple example. This relation is changed with introducing feedback loop in the amplifier.Figure 5
Two-stage bipolar amplifier.Figure 6
Third harmonic and its contributions.The second example is an OpAmp UA741 [24]. The transistors q1, q2, q3 belong to the input stage of the amplifier, and transistors q14, q20 are from the output stage. Figure 7 shows the dominant contributions related with these transistors to third harmonic as a function of frequency. Note that contributions of the output stage transistors are dominated over input stage transistors at low frequencies while they become the same order at high frequencies. Also note that contributions due to interaction between transistors q14, q20 are comparable with contributions of q14 and q20 individually and therefore cannot be neglected in this example.Figure 7
Contributions to the third harmonic.


## 6. Conclusion

A new approach for nonlinear distortion analysis has been presented. This approach does not require the computation of second and third derivatives of the semiconductor device models and, thus, has an important advantage over the Volterra series-based distortion analysis. The new technique is applicable to the class of weakly nonlinear and periodic time-varying circuits.A computational algorithm for periodic distortion analysis of communication circuits has been developed. Three steps of the algorithm provide the same order of accuracy that three steps of the distortion procedure based on Volterra series methods.In comparison with complete nonlinear steady-state analysis (such as harmonic balance), the approach provides significant reducing of computational efforts.In addition, the new technique allows computing the contributions of individual circuit components to the distortion. A matrix form for the characterization of third-order individual distortion contributions has been proposed. The numerical procedures for computing individual contributions have been developed, which provide designers with valuable information on the influence of various nonlinearities upon the total distortion.



---

*Source: 540305-2011-08-22.xml*