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.

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 [511] 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, 1315]. 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. Section 2 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]𝑑𝑑𝑡𝑞(𝑥(𝑡))+𝑖(𝑥(𝑡))+̃𝑢(𝑡)+𝑢(𝑡)=0,(1) where 𝑞 is the vector of charges, 𝑥 is the vectors of nodal voltages, ̃𝑢(𝑡) is a large periodic excitation with frequency 𝜔, and 𝑢(𝑡)=𝑢1(𝑡)+𝑢2(𝑡)++𝑢𝑀(𝑡) is a vector of small periodic excitations corresponding to the set of independent sources with incommensurate frequencies 𝜔𝑖,𝑖=1,,𝑀.

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 ̃𝑢(𝑡)=0, that is,𝑑𝑑𝑡𝑞(𝑥(𝑡))+𝑖(𝑥(𝑡))+𝑢(𝑡)=0.(2)

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𝑥(𝑡)𝑥(𝑛)(𝑡)=̃𝑥(𝑡)+Δ𝑥(0)(𝑡)++Δ𝑥(𝑛1)(𝑡),(3) where ̃𝑥(𝑡) is the periodic solution without small excitations, Δ𝑥(0)(𝑡) is the linear response of the system, that is, the result of periodic AC analysis [1], and Δ𝑥(𝑘)(𝑡) (𝑘>0) is the deviation of (𝑘+1)th order with respect to the norm of small excitations Δ𝑥(𝑘)(𝑡)=𝑂(𝑢(𝑡)𝑘+1). Here, the Landau symbol 𝑂 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𝑥(𝑛)(𝑡)𝑥(𝑡)=𝑂𝑢(𝑡)𝑛+1.(4)

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 solution 𝑥(1)(𝑡) is then used to calculate the equivalent current sources to determine the second-order solution 𝑥(2)(𝑡). Both the first- and second-order solutions are used to compute the third-order solution 𝑥(3)(𝑡). 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) where 𝑤, 𝑧 are 𝑁-dimensional vectors of input and internal signals in frequency or time domain. Let the vector 𝑧(0) be a solution of the problem obtained from (5) by neglecting small excitations, that is,𝑓𝑧(0)=0.(6) The nonlinear algebraic equation (5) are typically solved by the Newton’s method [19]𝐽𝑧(𝑘)Δ𝑧(𝑘)𝑧=𝑤𝑓(𝑘),𝑧(𝑘+1)=𝑧(𝑘)+Δ𝑧(𝑘),𝑘=0,1,2,,(7) where𝐽𝑧(𝑘)=𝜕𝑓1𝑧(𝑘)𝜕𝑧1𝜕𝑓1𝑧(𝑘)𝜕𝑧𝑁𝜕𝑓𝑁𝑧(𝑘)𝜕𝑧1𝜕𝑓𝑁𝑧(𝑘)𝜕𝑧𝑁(8) is the Jacobian matrix of the system (5), and Δ𝑧(𝑘) 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 𝐽(𝑧(0)). Thus, the method does not require any update of the Jacobian matrix.

The numerical scheme contains one Newton step and (𝑘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]: 𝑧(𝑘)𝑧𝑧=𝑂(0)𝑧𝑘+1,𝑘=1,2,,(9) where 𝑧 is the exact solution.

Using the estimate 𝑧(0)𝑧=𝑂(𝑤), we have from (9) for 𝑘=𝑛 that𝑧(𝑛)𝑧=𝑂𝑤𝑛+1.(10)

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 Δ𝑧(𝑘) up to the third order, the computational procedure after the determination of the initial guess by (6) includes the following three steps:(1)𝐽𝑧(0)Δ𝑧(0)=𝑤,(11)(2)𝐽𝑧(0)Δ𝑧(1)𝑧=𝑤𝑓(1),(12)(3)𝐽𝑧(0)Δ𝑧(2)𝑧=𝑤𝑓(2).(13) 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 Δ𝑧(0), Δ𝑧(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 of 𝑓(𝑧(2)) yields𝑓𝑧(2)𝑧=𝑓(1)+Δ𝑧(1)𝑧=𝑓(1)𝑧+𝐽(1)Δ𝑧(1)+𝑂Δ𝑧(1)2.(14)

Taking into account that Δ𝑧(1)=𝑂(𝑤2), we can evaluate 𝑂(Δ𝑧(1)2)=𝑂(𝑤4). So neglecting this term does not change the order of error at step 3 (see (10) for 𝑛=3). The substitution (14) into (13) leads to the modified expression of step 3𝐽𝑧(0)Δ𝑧(2)𝑧=𝑤𝑓(1)𝑧𝐽(1)Δ𝑧(1).(15)

Subtracting (12) from (15), a more convenient form of this expression is obtained 𝐽𝑧(0)Δ𝑧(2)=𝐽𝑧(0)𝑧𝐽(1)Δ𝑧(1).(16) In many cases, it is efficient to perform step 2 by the equation 𝐽𝑧(0)Δ𝑧(1)=𝑤𝜑,(17) where vector 𝜑 approximates 𝑓(𝑧(1)) with the accuracy𝑧𝜑𝑓(1)=𝑂𝑤3.(18)

If step 2 is performed by (17) instead of (12), then previous considerations lead to the following expression instead of (16):𝐽𝑧(0)Δ𝑧(2)=𝐽𝑧(0)𝑧𝐽(1)Δ𝑧(1)𝑧+𝜑𝑓(1).(19)

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 Section 2 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𝐻(𝑋)=𝐼(𝑋)+Ω𝑄(𝑋)+𝑈+𝑀𝑖=1𝑈𝑖=0,(20) where vectors 𝑋, 𝐼(𝑋), 𝑄(𝑋), 𝑈, and 𝑈 are frequency domain representation of vectors 𝑥(𝑡), 𝑖(𝑥), 𝑞(𝑥), ̃𝑢(𝑡), and 𝑢(𝑡), and Ω is a block-diagonal matrix of combinational frequencies of fundamentals 𝜔, 𝜔𝑖.

Solving the nonlinear system (20) by the Newton’s method results in the linear system for corrections Δ𝑋(𝑘),𝐽HB𝑋(𝑘)Δ𝑋(𝑘)𝑋=𝐻(𝑘),(21) where 𝐽HB is the multitone harmonic Jacobian matrix [5].

The vector of unknowns contains all harmonics of the linear combinations with integer coefficients of all fundamentals ̃𝜔=𝑘𝜔+𝜈, where 𝜈 denotes the small signal combinational frequency of the form𝜈=𝑘1𝜔1+𝑘2𝜔2++𝑘𝑀𝜔𝑀.(22)

The number of complex unknowns in (20) or (21) is equal to 𝑁(2𝐾𝐾+𝐾+𝐾+1), where 𝑁 is the circuit dimension, 𝐾 is the number of large signal harmonics, and 𝐾 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𝐻𝑋=𝐼𝑋+𝑄𝑋+Ω𝑈=0.(23) The nonlinear system (23) is obtained from (20) by applying only the large signal excitation.

After determining the initial guess 𝑋(0)=𝑋, one can apply the iterative process (7) that for system (20) is presented in the form𝐽HB𝑋Δ𝑋(𝑘)𝑋=𝐻(𝑘).(24) The Jacobian matrix 𝐽HB(𝑋) 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)𝐽SS𝜈𝑗Δ𝑋(𝑘)𝜈𝑗=𝐻(𝑘)𝜈𝑗.(25)

The number of unknowns in each system (25) is equal to 𝑁(2𝐾+1), and there are 𝐾 systems with different 𝜈𝑗.

The entries of the matrix 𝐽SS(𝜈𝑗) can be defined by the Fourier transform of periodically varying conductance 𝐺(𝑡) and capacitance 𝐶(𝑡) 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 vector 𝐻(𝑋(𝑘)) 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𝐻(2)=𝐽HB𝑋(1)𝐽HB𝑋Δ𝑋(1),(26) where 𝐽HB(𝑋(1)) is harmonic Jacobian matrix (21) that is computed at the step 2 simultaneously with the computation of 𝐻(𝑋(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 𝑢𝑖=𝑈𝑖exp(±𝑗𝜔𝑖𝑡). In such case, the minimal sets of small signal combinational frequencies at the 𝑛th step 𝜈(𝑛)={𝜈1(𝑛),𝜈2(𝑛),} are defined by the following recurrent expressions:𝜈(1)=,𝜔𝑖,𝜈,(𝑛+1)=||𝜈,𝑗(𝑛)±𝜔𝑖||.,(27)

In the case of two small excitations, frequency sets (27) of second and third orders result in the forms𝜈(2)=0,2𝜔1,2𝜔2,𝜔1+𝜔2,||𝜔1𝜔2||𝜈,(28)(3)=𝜔1,𝜔2,2𝜔1+𝜔2,𝜔1+2𝜔2,||2𝜔1𝜔2||,||2𝜔2𝜔1||,3𝜔1,3𝜔2.(29)

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, 𝐻(1)(𝜈𝑗(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 step ̃𝑥(𝑡),𝐺(𝑡), 𝐶(𝑡).

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, 𝐽SS𝜔1Δ𝑋(0)𝜔1=𝑈1,𝐽SS𝜔2Δ𝑋(0)𝜔2=𝑈2.(30)

Step 3. (a) Determine the first-order solution in time domain by applying inverse Fourier transform Γ1, 𝑥(1)(𝑡)=̃𝑥(𝑡)+Γ1Δ𝑋(0).(31)
(b) Compute time domain vectors of circuit charges, currents, and admittance matrices,𝑞𝑥(1)𝑥(𝑡),𝑖(1)𝑥(𝑡),𝐺(1)𝑥(𝑡),𝐶(1)(𝑡).(32)
(c) Compute the residual vector in the frequency domain by applying the Fourier transform Γ and using (20),𝐻𝑋(1)𝑥=Γ𝑖(1)𝑥(𝑡)+ΩΓ𝑞(1)+(𝑡)𝑈+𝑈1+𝑈2.(33)
(d) Decompose the vector 𝐻(𝑋(1)) into 𝐻(1)(𝜈𝑗(2)) and 𝐻(1)(𝜈𝑗(3)).
(e) Determine corrections of solution by solving the systems𝐽SS𝜈𝑗(2)Δ𝑋(1)𝜈𝑗(2)=𝐻(1)𝜈𝑗(2),(34) for all of 5 small signal combinational frequencies 𝜈𝑗(2) from the set (28).

Step 4. (a) Transform the second-order correction to time domain by applying inverse Fourier transform Δ𝑥(1)(𝑡)=Γ1Δ𝑋(1).(35)
(b) Compute vector (26) by matrix-vector multiplication in time domain and Fourier transform to frequency domain𝐻(2)𝑥=Γ𝐺(𝑡)𝐺(1)(𝑡)Δ𝑥(1)𝑥(𝑡)+𝑗ΩΓ𝐶(𝑡)𝐶(1)(𝑡)Δ𝑥(1).(𝑡)(36)
(c) Decompose vector 𝐻(2) into 𝐻(2)(𝜈𝑗(3)).
(d) Determine corrections by solving systems𝐽SS𝜈𝑗(3)Δ𝑋(2)𝜈𝑗(3)=𝐻(2)𝜈𝑗(3)𝐻(1)𝜈𝑗(3),(37) 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 (2𝑛+1), the aliasing error depends on harmonics higher than 𝑛. Thus, the required accuracy order (10) is provided for 𝑛=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 vector 𝑋(0) contains 𝑁 nonzero components corresponding to the DC solution vector, and the 𝑁×𝑁 Jacobian matrix corresponds to the AC matrix 𝐽SS(𝜈)=𝐺+𝑗𝜈𝐶, where matrices 𝐺 and 𝐶 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 relation 𝐼=𝑉2 be excited by unit DC current source and a small sinusoidal signal. The corresponding equation has the form𝑣21𝑎cos(𝜔𝑡)=0.(38) The DC solution of (38) is 𝑣(0)=1.

The first-order correction is determined from the linear system obtained by linearization (38) 2Δ𝑣(0)=𝑎cos(𝜔𝑡).(39) Therefore, the first-order solution is expressed as𝑣(1)=𝑣(0)+Δ𝑣(0)𝑎=1+2cos(𝜔𝑡).(40) The distortion contribution is determined from the linear system [20]𝑎2Δ𝑣=1+2cos(𝜔𝑡)2𝑎122cos(𝜔𝑡).(41) After elementary transformations, we have 𝑎Δ𝑣=2𝑎16216cos(2𝜔𝑡).(42)

Thus, we see that the voltage distortion does not contain third harmonics. Here, Δ𝑣 corresponds to the value Δ𝑣(1) from the recursive procedure.

For this simple example, it is easy to obtain the solution of (38) in the form𝑣=1+𝑎cos(𝜔𝑡).(43)

This function contains third-order terms in its Taylor expansion (𝑦=𝑎cos(𝜔𝑡))𝑦1+𝑦=1+2𝑦28+𝑦3𝑦16+𝑂4,(44) which provide third harmonics. Note that if we perform the next step of our computational procedure, we obtainΔ𝑣(2)=3𝑎3𝑎64cos(𝜔𝑡)+3𝑎64cos(3𝜔𝑡)+𝑂4.(45) 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𝐷𝐷=1,𝐷2,,𝐷𝑚,,(46) with the following properties:(1)𝑚th component of the vector depends on 𝑚th nonlinearity and does not depend on any other nonlinearity,(2)total distortion is the sum of all components of the vector𝐷total=𝑚𝐷𝑚.(47)

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 Figure 1. Each stage of the amplifier is described by a transfer factor (𝐾) with small nonlinearities of second (𝑎) and third (𝑏) order 𝑣out1𝑣inp1=𝐾1𝑣inp1+𝑎1𝑣2inp1+𝑏1𝑣3inp1,𝑣out2𝑣inp2=𝐾2𝑣inp2+𝑎2𝑣2inp2+𝑏2𝑣3inp2.(48) The circuit contains two nonlinearities defined by parameters 𝑎1, 𝑏1 for the first nonlinearity and 𝑎2, 𝑏2 for the second nonlinearity.

The output signal of the amplifier is obtained by 𝑣out(𝑣)=𝑣out2𝑣out1(𝑣).(49) After substituting (48) into (49) and neglecting terms higher than third order, this expression can be presented in the form𝑣out(𝑣)=𝑣(1)out+𝑣(2)out+𝑣(3)out,(50) where 𝑣(1)out=𝐾1𝐾2𝑣, 𝑣(2)out=𝐾2𝑎1𝑣2+𝐾21𝑎2𝑣2, and 𝑣(3)out=𝐾2𝑏1𝑣3+𝐾31𝑏2𝑣3+2𝐾1𝑎1𝑎2𝑣3.

One can see that the expression for the second-order distortion 𝑣(2)out is the sum of two terms, each of them corresponds to only one nonlinearity (𝑎1 or 𝑎2) and so can be considered as the contribution of this nonlinearity. Thus, the second-order contributions 𝐷(2) can be presented in the form (46) with 𝐷1(2)=𝐾2𝑎1𝑣2, 𝐷2(2)=𝐾21𝑎2𝑣2.

In contrast, the expression for the third-order distortion contains not only terms corresponding to each nonlinearity (𝐾2𝑏1𝑣3,𝐾31𝑏2𝑣3) but also the term 2𝐾1𝑎1𝑎2𝑣3 depending on both nonlinearities.

Hence, the influences of circuit nonlinearities on the third-order distortion can be presented as a matrix of contributions 𝐷(3)=𝐷(3)110𝐷(3)21𝐷(3)22,(51) where 𝐷(3)11=𝐾2𝑏1𝑣3, and 𝐷(3)22=𝐾31𝑏2𝑣3, 𝐷(3)21=2𝐾1𝑎1𝑎2𝑣3.

The meaning of the nondiagonal term 𝐷(3)21 can be more clear if it is presented as a function of the two output signals of the first stage𝐷(3)21=2𝑎2𝑣(1)out1𝑣(2)out1,(52) where 𝑣(1)out1=𝐾1𝑣 is linear output signal, and 𝑣(2)out1=𝑎1𝑣2=𝑎1𝑣𝑣 is the second-order distortion signal produced by the first nonlinearity. Then we can interpret 𝐷(3)21 (52) as a result of mixing signals 𝑣(1)out1, 𝑣(2)out1 by the nonlinearity of the stage 2.

Here, the matrix 𝐷(3) is triangular because the circuit is unidirectional, but if feedback loop is present, then 𝐷(3)120. To see this, consider the next example shown in Figure 2.

The transfer function 𝑣out(𝑣) for the circuit is defined by the implicit expression 𝑣out(𝑣)=𝑣out1(𝑣𝑣out2(𝑣out(𝑣))). 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:𝑣(1)out1=𝐾1𝑣1+𝐾1𝐾2,𝐷(53)1(2)=𝑎1𝑣21+𝐾1𝐾23,𝐷2(2)𝑎=2𝐾31𝑣21+𝐾1𝐾23𝐷,(54)(3)11=2𝑎1𝐷1(2)𝐾2𝑣31+𝐾1𝐾22+𝑏1𝑣31+𝐾1𝐾24,𝐷(3)22=2𝑎2𝐷2(2)𝐾21𝐾2𝑣31+𝐾1𝐾22+𝑏2𝐾41𝐾32𝑣31+𝐾1𝐾24,𝐷(55)(3)12=2𝑎1𝐷2(2)𝑣3𝐾11+𝐾1𝐾22,𝐷(3)21=2𝑎2𝐷1(2)𝐾21𝐾2𝑣31+𝐾1𝐾22.(56)

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 matrix 𝐷(3) with entries 𝐷(3)𝑛𝑚.

The diagonal entry 𝐷(3)𝑛𝑛 defines the third-order signal that is produced in 𝑛th 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 𝐷(3)𝑛𝑚 defines the third-order signal that is produced in 𝑛th nonlinearity by the mixing of linear signals with the second-order signals produced in 𝑚th nonlinearity.

The following properties of the contribution matrix exist.(1)The sum of all matrix entries is equal to the total third-order distortion 𝐷(3)total=𝑛𝑚𝐷(3)𝑛𝑚.(57)(2)If parameters of the 𝑛th nonlinearity are changed, then only entries of 𝑛th column and 𝑛th row are modified. All other entries remain unchanged. In particular, if parameters of the 𝑛th nonlinearity are set to zero, then all entries of 𝑛th column and 𝑛th row are zeros. All other entries remain unchanged.(3)The sum of entries of 𝑛th column and 𝑛th row 𝐷𝑛(3)=𝑚𝑛𝐷(3)𝑛𝑚+𝐷(3)𝑚𝑛+𝐷(3)𝑛𝑛(58) defines the total influence of 𝑛th 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 𝐷𝑛(3) cannot be considered as the contribution of 𝑛th nonlinearity in the sense (46) because property (47) is not true.(4)If the circuit contains a group of nonlinearities 𝑆={𝑖1,𝑖2,} and it is desired to present the group as one nonlinearity with index 𝑘, then the matrix is obtained using the following expressions:𝐷(3)𝑘𝑚=𝑛𝑆𝐷(3)𝑛𝑚,𝐷(3)𝑛𝑘=𝑚𝑆𝐷(3)𝑛𝑚,𝐷(3)𝑘𝑘=𝑛𝑆𝑚𝑆𝐷(3)𝑛𝑚.(59)

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 one 𝐷(3)𝑛𝑚, then a designer can also change frequency characteristics of the signal path from 𝑚th to 𝑛th 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 model 𝑓(𝑧) (5) must be presented as the sum of nonlinear dependencies 𝑓(𝑧)=𝑚𝑓(𝑚)(𝑧),(60) 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)

For the representations (60), (61), step 2 (12) can be written in the form𝐽𝑧(0)Δ𝑧(1)=𝑚𝐽(𝑚)𝑧(0)Δ𝑧(0)𝑓(𝑚)𝑧(1).(62)

Here, we exploit the fact that from (11), we have𝑚𝐽(𝑚)𝑧(0)Δ𝑧(0)𝑧=𝐽(0)Δ𝑧(0)=𝑤.(63)

Each term in the residual vector of (62) defines the effect of 𝑚th nonlinearity because it presents the difference between the linearized 𝐽(𝑚)(𝑧(0))Δ𝑧(0) and the nonlinear 𝑓(𝑚)(𝑧(1)) dependencies of the nonlinearity. Hence, the second-order distortion can be written as a sum of individual contributionsΔ𝑧(1)=𝑚Δ𝑧(1,𝑚),Δ𝑧(1,𝑚)=𝐽1𝑧(0)𝐽(𝑚)𝑧(0)Δ𝑧(0)𝑓(𝑚)𝑧(1).(64)

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𝐽𝑧(0)Δ𝑧(2)=𝑛𝐽(𝑛)𝑧(0)𝐽(𝑛)𝑧(1)𝑚Δ𝑧(1,𝑚).(65)

Thus, the third-order distortion isΔ𝑧(2)=𝑛𝑚Δ𝑧(2,𝑛,𝑚),Δ𝑧(2,𝑛,𝑚)=𝐽1𝑧(0)𝐽(𝑛)𝑧(0)𝐽(𝑛)𝑧(1)Δ𝑧(1,𝑚).(66) Each vector (66) depends only on two nonlinearities 𝑚, 𝑛.

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 𝐷𝑚(2)=Δ𝑧(1,𝑚)out, and corresponding components of the vectors (66) form the matrix of third-order contributions 𝐷(3)𝑛𝑚=Δ𝑧(2,𝑛,𝑚)out, 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 𝐽SS𝜈𝑗(2)Δ𝑋(1,𝑚)𝜈𝑗(2)=𝐻(1,𝑚)𝜈𝑗(2)𝐽,(67)SS𝜈𝑗(3)Δ𝑋(2,𝑛,𝑚)𝜈𝑗(3)=𝐻(2,𝑛,𝑚)𝜈𝑗(3)𝛿𝑚𝑛𝐻(1,𝑛)𝜈𝑗(3),(68) where 𝛿𝑚𝑛 is the Kronecker delta. In accordance with (26), 𝐻(2,𝑛,𝑚)=[𝐽HB(𝑛)(𝑋(1))𝐽HB(𝑛)(𝑋)]Δ𝑋(1,𝑚).

The computations can be performed in the same numerical framework presented in Section 3. To obtain vectors 𝐻(1,𝑚)(𝑋(1)) and matrices 𝐽HB(𝑛)(𝑋(1)), the HB loading procedure must be developed on the per-nonlinearity basis.

Let 𝑁𝑁 be the number of nonlinearities in the circuit, and 𝐾2, 𝐾3 are the numbers of second- and third-order small signal frequencies, respectively. Then the computation of second-order contributions requires solving 𝐾2𝑁𝑁 linear systems (67) with only 𝐾2 different matrices. For the third-order contribution matrix 𝐾3𝑁2𝑁 linear systems (68) should be solved. Fortunately, using the technique of adjoint linear system, this work is reduced to solving 𝐾3 linear system with transposed matrix and 𝐾3𝑁2𝑁 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 parameters 𝐴LO = 2.5 V, 𝑓LO = 375 MHz, 𝑓1 = 394.9 MHz, and 𝑓2 = 395.1 MHz, with 𝑟𝑓 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 𝐴LO = 0.75 V, 𝑓LO = 985 MHz, 𝑓1 = 899 MHz, and 𝑓2 = 901 MHz, with 𝑟𝑓 input varying from 0 to 0.15 V. In both circuits, 𝑓1 and 𝑓2 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 (𝑓LO𝑓1) on the input power are shown in Figure 3 for the first mixer. The obtained curves are in good agreement, especially for small 𝑟𝑓 input amplitudes as expected. Similar results are obtained for the second mixer. The dependencies of the output power at the third-order intermodulation product (𝑓LO2𝑓1+𝑓2) are shown in Figure 4 for the first mixer. The obtained curves are in good agreement for practical input amplitudes.

The computational efficiency of the PDHB approach for the two mixer circuits can be seen from Table 1. The number of HB harmonics for the three tones 𝑓LO, 𝑓1, 𝑓2 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.

5.2. Examples of Individual Contributions

The first example is the two-stage bipolar amplifier shown in Figure 5. 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 𝑞1 is smaller by a factor of 8–12, and the contribution due to the interaction of transistors 𝑞1 and 𝑞2 is negligible in this simple example. This relation is changed with introducing feedback loop in the amplifier.

The second example is an OpAmp UA741 [24]. The transistors 𝑞1, 𝑞2, 𝑞3 belong to the input stage of the amplifier, and transistors 𝑞14, 𝑞20 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 𝑞14, 𝑞20 are comparable with contributions of 𝑞14 and 𝑞20 individually and therefore cannot be neglected in this example.

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.