2.1.1 Representation of Circuits and Systems

In the wireless transmitter example in Figure 2.1 we can recognize a very heterogeneous set of circuits such as digital signal-processing and control, mixed digital-analog, analog IF or baseband, and RF. Considering only this latter subset, we can note circuits that operate in an autonomous way, like the local oscillator, and ones that are driven, and thus present a forced behavior. In addition, we may also divide them into circuits solely composed of lumped elements (like ideal resistors, capacitors, inductors, and controlled sources) or into ones that include distributed elements (such as transmission lines, microstrip discontinuities, waveguides, etc.).

2.1.1.1 Circuits of Lumped Elements

As seen in Chapter 1, the response of dynamic circuits or systems cannot be given as a function of solely the instantaneous input, but depends also on the input past, stored as the circuit’s state. For example, in lumped circuits, i.e., circuits composed of only lumped elements, their state comprises the electric charge stored in the capacitors and the magnetic flux stored in the inductors, the so-called system’s state-variables or the state-vector. Assuming that these capacitor charges and inductor fluxes are memoryless functions of the corresponding voltages and currents (the quasi-static assumption in which accumulated charges, fluxes, or currents are assumed to be instantaneous functions of the applied fields), the state-equation is thus the mathematical rule that describes how the state-vector evolves in time. It is an ordinary differential equation, ODE, in time that expresses the time-rate of change of the state-vector, i.e., its derivative, with respect to time, as a function of the present state and excitation:

dstdt=fstxt(2.1)

where s(t) is state-vector, x(t) is the vector of inputs, or excitation-vector, (independent current and voltage sources) and f[s(t), x(t)]fstxt is a set of functions, each one derived from applying Kirchhoff’s currents law, KCL, to the nodes, and voltages law, KVL, to the branches.

With this state-equation, one can calculate the trajectory of the state-vector over time, s(t). Then, it is easy to compute the output-vector, y(t), i.e., the set of observable currents, voltages, or any algebraic operation on these (such as power, for example) from the following output-equation:

y(t) = g [s(t), x(t)]

yt=gstxt(2.2)

Example 2.1 The State- and Output-Equations

For an illustrative example of this state- and corresponding output-equations, let us use the circuit of Figure 2.2. The application of Kirchhoff currents law to the input and output nodes, leads to

vGt−vItRs+CgdvGtdt+CbddtvGt−vDt=0CbddtvDt−vGt+iDvGvD+iLt+vDtRL=0vDt−VDD=LddiLtdt(2.3)

which can be rewritten as

dvGtdt=1CgvIt−vGtRs−iDvGvD−iLt−vDtRLdvDtdt=1CgvIt−vGtRs−Cb+CgCbCgiDvGvD+iLt+vDtRLdiLtdt=vDt−VDDLd(2.4)

Defining now s₁(t) ≡ v_G(t), s₂(t) ≡ v_D(t), s₃(t) ≡ i_L(t), x₁(t) ≡ v_I(t)s1t≡vGt,s2t≡vDt,s3t≡iLt,x1t≡vIt and x₂(t) ≡ V_DDx2t≡VDD allows us to express these nodal equations in the canonical state-equation form of (2.1) as

dsntdt=fns1t…snt…sNtx1t…xmt…xMt(2.5)

where N = 3, M = 2N=3,M=2 and

f1stxt=1Cgx1t−s1tRs−iDs1ts2t−s3t−s2tRLf2stxt=1Cgx1t−s1tRs−Cb+CgCbCgiDs1ts2t+s3t+s2tRLf3stxt=s2t−x2tLd(2.6)

Assuming that the desired observable is the output voltage, v_O(t)vOt, the canonical output equation can be obtained following a similar procedure, defining y(t) ≡ v_O(t) = v_D(t) = s₂(t)yt≡vOt=vDt=s2t

y(t) = g[s₁(t), s₂(t), s₃(t), x₁(t), x₂(t), ]

yt=gs1ts2ts3tx1tx2t(2.7)

which, in this case, is simply y(t) = s₂(t)yt=s2t.

Figure 2.2 Simplified FET amplifier circuit (a) and its equivalent circuit (b) intended to illustrate the analysis’s algorithms. v_I(t) = V_GG + v_S(t) where v_S(t) is the excitation, V_GG = 3 V, V_DD = 10 V, R_s = 50 Ω, C_g = 50 pF, C_b = 4.4 pF, L_d = 120 nH, and R_L = 6.6 Ω. The MOSFET has a transconductance of g_m = 10 A/V, a threshold voltage of V_T = 3 V, and a knee voltage of V_K = 1 V.

Because this system is a third order one, its behavior is completely described by a three-dimensional state-vector, i.e., involving three state-variables. So, as explained in Section 1.4.4, its phase-portrait should be represented in a volume whose axis are s₁, s₂, and s₃, or v_G, v_D, and i_L, as is shown in Figure 2.3. Please note that the state-variables correspond to the expected circuit’s three independent energy storing elements: s₁(t) ≡ v_G(t)s1t≡vGt is a measure of the accumulated charge (electric field) in the input capacitor C_g, s₂(t) ≡ v_D(t)s2t≡vDt is needed to build v_GD(t) = v_G(t) − v_D(t) = s₁(t) − s₂(t)vGDt=vGt−vDt=s1t−s2t, a measure of the accumulated charge in the feedback capacitor C_b, and s₃(t) ≡ i_L(t)s3t≡iLt stands for the accumulated magnetic flux in the output drain inductor L_d.

Figure 2.3 (a) Illustrative trajectory of the three-dimensional system state, [x y z] ≡ [v_G v_D i_L], of our example circuit. (b) Projection of the 3D state-space trajectory onto the [v_G v_D] plane. (c) Projection of the 3D state-space trajectory onto the [v_G i_L] plane. (d) Projection of the 3D state-space trajectory onto the [v_D i_L] plane.

Although the phase-portrait is a parametric curve in which the evolving parameter is time, it does not allow a direct observation of the trajectory of each of the state-variables over time. In fact, a graph like Figure 2.3 (a) cannot tell us when the system is more active or, on the contrary, when it is almost latent. Such a study would require separate graphs of s₁(t), s₂(t) or s₃(t) versus time as is exemplified in Figure 2.4 for s₁(t) and s₂(t).

Figure 2.4 Time evolution of two of the system’s state-variables over time. (a) s₁(t) = v_G (t) and (b) s₂(t) = v_D(t).

2.1.1.2 Circuits with Distributed Elements

Contrary to circuits that include only lumped elements, circuits or systems that involve distributed components, such as transmission lines or other distributed electromagnetic multi-ports, cannot be described as ODEs in time, requiring a description in both time and space. That is the case of the power divider/combiner, the circulator or the output microwave filter of Figure 2.1, and of the transmission line of Figure 2.5. As an example, the transmission line is described by the following partial differential equations, PDEs, in time, t, and space, x.

∂vxt∂x=−L∂ixt∂t−Rixt∂ixt∂x=−C∂vxt∂t−Gvxt(2.8)

where i(x, t)ixt and v(x, t)vxt are the distributed current and voltage along the line and R, L, C, and G, stand for the usual line’s distributed resistance, inductance, capacitance, and conductance per unit length, respectively.

Figure 2.5 The uniform transmission line as an example of a distributed component. Such distributed and dispersive components are represented by partial differential equations in time and space, impeding the use of the state-equation as a general model.

Please note that, because all time-domain circuit simulators are programmed to solve ODEs, they require that the circuit be represented by the state- and output-equations of (2.1) and (2.2). Consequently, a lumped representation is needed for distributed components. Such lumped representation is usually given in the form of an approximate equivalent circuit, obtained from the frequency-domain behavior observed at the device ports, which mimics the operation of the original multiport [1]. However, alternative behavioral representations, such as the one based on the impulse response function and subsequent convolution with the component’s excitation, are also possible [2].

It should be emphasized that these are behavioral representations, not real equivalent circuits, as distributed elements do not have a true lumped representation. For example, even our uniform transmission line of Figure 2.5 would need an infinite number of elementary lumped RLCG sections. Therefore, if one is willing to simulate such a circuit in the time domain, he or she should be aware that he is not simulating the actual circuit but another one that, hopefully, produces a similar result. For example, if the extracted “equivalent” circuit is not guaranteed to obey certain properties, such as passivity [3], one may end up with nonconvergence results or erroneous solutions such as a nonphysical (and so impossible) oscillatory response from an original circuit composed of only passive components.

2.1.2 Representation of Signals

As far as the signals are concerned, we could group them into baseband – or slowly evolving – information signals, fast-varying RF carriers, and modulated signals, in which the first two are combined. Furthermore, as shown in the spectra included in Figure 2.1, we could also note that, contrary to the RF carriers, which are periodic and thus narrowband, baseband signals are wideband (that span through many octaves or decades) and are necessarily aperiodic, as they carry information (nonrepeatable, by definition). Therefore, while the periodic RF carrier can be represented in either the time or frequency domains, baseband signals must be necessarily treated in the time domain since the discrete Fourier series assumes periodicity in both time and frequency. Actually, because we will be using a digital computer, i.e., a machine with a finite number of states, it will necessarily impose a discrete, or sampled and of confined nature, description of any time- or frequency-domain entities. This means that if you want to represent, in the frequency domain, a certain continuous-time signal, extending from, say, 0 to T seconds, what you will get is the frequency-domain representation of another discrete and periodic signal, of period T, which extends from –∞ to +∞, and whose sampling frequency, f_s, is determined by the computer time resolution. Such a signal has a discrete and periodic spectrum whose frequency resolution is Δf  =  1/T and whose period is f_s. Conversely, any continuous spectrum signal will be treated as one of discrete spectrum, which has a periodic time-domain representation. In addition, any Fourier conversion between these two domains will necessarily be executed via the discrete Fourier transform, DFT, or its fast version, the FFT, in case the number of time and frequency points is a multiple of two.

2.2.1 Time-Step Integration Engines

Time-step integration is used in all SPICE-like transient simulators [4], [5] but also in many system simulators [6–8]. Since it aims at the evolution of the circuit in time, it can be seen as the most “natural” circuit analysis technique because it closely follows the conceptual interpretation we make of the circuit operation. Actually, assuming one starts from a known excitation, x(t), discretized in a number of appropriate time-points, x(t_n), and an also known, initial state, s(t₀), one can easily compute the response at all future time samples, s(t_n), from the discretization of the state-equation. This discretization reflects the knowledge that the derivative of the state at the time, t_n, is the limit when the time-step length, h_n, tends to zero, of the ratio [s(t_n + 1) − s(t_n)]/h_nstn+1−stn/hn. Therefore, the state-equation (2.1) can be approximated by the following difference equation

stn+1−stnhn=fstnxtn(2.9)

which allows the calculation of all s(t_n) using the following recursive scheme known as the forward Euler method [9]

in a time-step by time-step basis. It is this time-step integration of the state-equation that constitutes the root of the name of this family of methods.

Alternatively, we could have also approximated the state derivative at t_n by

stn+1−stnhn=fstn+1xtn+1(2.11)

which would then lead to the alternative recursive scheme of

which is known as the backward Euler method.

Despite the obvious similarities between (2.10) and (2.12) (originated in their common genesis), they correspond to very distinct numerical methods. First of all, (2.10) gives the next state, s(t_n+1), in an explicit way, since everything in its right-hand side is a priori known. On the contrary, (2.12) has an implicit formulation since the next state s(t_n+1) depends on the known s(t_n) and x(t_n+1), but also on itself. Hence, s(t_n+1) can be obtained through simple and fast algebraic operations from (2.10), while the solution of (2.12) demands for an iterative nonlinear solver, for which it is usually formulated as

Actually, (2.13) can be solved using, e.g., a gradient method like the Newton-Raphson iteration [10].

A Brief Note on the Newton-Raphson Iteration

The Newton-Raphson iteration plays such an important role in nonlinear circuit simulation that we will outline it here. Based on a first order Taylor series approximation, it finds the zero of a nonlinear function, or vector of functions, through the solution of a successively approximated linear system of equations. In the present case, the vector of functions whose zero is to be calculated is φ[s(t_n + 1)] ≡ s(t_n + 1) − s(t_n) − h_nf[s(t_n + 1), x(t_n + 1)]φstn+1≡stn+1−stn−hnfstn+1xtn+1 and so its first-order (linear) Taylor approximation in the vicinity of a certain predetermined estimate, s(t_n + 1)_istn+1i, becomes

φ[s(t_n + 1)]≈φ[s(t_n + 1)_i] + J_φ[s(t_n + 1)_i][s(t_n + 1)_i + 1 − s(t_n + 1)_i] = 0

φstn+1≈φstn+1i+Jφstn+1istn+1i+1−stn+1i=0(2.14)

where J_φ[s(t_n + 1)_i]Jφstn+1i is the Jacobian matrix of φ[s(t_n + 1)]φstn+1 whose element of row m and column l is calculated by

Jφstn+1iml≡∂φmstn+1∂sltn+1stn+1i(2.15)

i.e., it is the derivative of the m’th component of the vector of functions φ[s(t_n + 1)]φstn+1 with respect to the l’th component of the state-vector s(t_n + 1)stn+1, evaluated at the previous estimate s(t_n + 1)_istn+1i.

Since (2.14) is now a linear system in the unknown s(t_n + 1)_i + 1stn+1i+1, we can get a (hopefully) closer estimate of the desired zero by the following iteration:

s(t_n + 1)_i + 1 = s(t_n + 1)_i − [J_φ[s(t_n + 1)_i]]⁻¹φ[s(t_n + 1)_i]

stn+1i+1=stn+1i−Jφstn+1i−1φstn+1i(2.16)

This is graphically illustrated in Figure 2.6 (a) for the simplest one-dimensional case. Please note that, if φ[s(t_n + 1)]φstn+1 were linear, its Jacobian matrix would be constant, i.e., independent of s(t_n + 1)_istn+1i, and the Newton-Raphson algorithm would converge in only one iteration. In general, the Newton-Raphson method is faster in presence of mildly nonlinear problems and slower – i.e., it needs more iterations – when dealing with strongly nonlinear problems (higher curvature in the function of Figure 2.6). Furthermore, it requires an initial estimate of the solution, and converges if the function φ[s(t_n + 1)]φstn+1 is well behaved, i.e., if its Jacobian matrix has no points of zero determinant. In fact, if s(t_n + 1)_istn+1i is a point in the state-space in which the determinant of the Jacobian is zero, or close to zero, the next estimate will suffer a big jump, possibly to a point very far from the desired solution φ[s(t_n + 1)] = 0φstn+1=0, even if it were already close to it at s(t_n + 1)_istn+1i. As illustrated now in Figure 2.6 (b), the Newton-Raphson iteration of (2.16) may diverge from the target solution even when it is already quite close to it. Alternatively, since in the vicinity of a zero determinant J_φ[s(t_n + 1)_i]Jφstn+1i changes sign, the iteration may be stalled, endlessly jumping from two distinct points as depicted in Figure 2.6 (c). In this case, (2.16) does not diverge, but it cannot also converge to the solution.

Figure 2.6 Graphical illustration of the Newton-Raphson iteration in the simplest one-dimensional case. (a) Normal convergence of the Newton-Raphson iteration. (b) When the function of which we want to calculate the zero, φ(s)φs, has a zero in its derivative (the Jacobian), the iteration may diverge from the desired solution. (c) Situation in which the iteration does not either diverge or converge but is kept wandering around the vicinity of a zero derivative point.

Unfortunately, the simplicity of (2.10) has a high price in numerical stability. Since it is a predictive method, its associated error tends to build up very rapidly, a problem that must be obviated by reducing the time-step, and thus increasing the overall simulation time. In practice, this increase in simulation time and the potential numerical instability are so significant, when compared to (2.12), that (2.10) is hardly ever used in practical circuit simulators.

Actually, (2.10) and (2.12) do not constitute the only way that we could have used to integrate the state-equation. They are just examples of a zero-order method in which the derivative is assumed constant during the time-step interval. For example, if a first-order integration (known as the trapezoidal integration rule) would have been used, we would obtain [8]

stn+1=stn+hnfstn+1xtn+1+fstnxtn2(2.17)

In general, the time-step integration scheme is formulated as

stn+1=stn+∫tntn+1fstxtdt(2.18)

from which the 4th-order Runge-Kutta method [8] is one of its most popular implementations (see Exercise 2.1).

Finally, it should be added that there is no need to keep the time-step, h_n, constant. In fact, all modern time-domain simulators use dynamic time-steps, measuring the variation in time of s(t). Larger time-steps are used when s(t) is slowly varying, while smaller ones are adopted in zones of fast change of state. This way, a good compromise between accuracy and computation time is obtained.

Example 2.2 The Time-Step Integration Formulation

As an example of time-step integration, the trapezoidal rule discretization of our circuit example described by (2.5) and (2.6) would lead to

s1tn+1=s1tn+hnf1tn+1+f1tn2s2tn+1=s2tn+hnf2tn+1+f2tn2s3tn+1=s3tn+hnf3tn+1+f3tn2(2.19)

where the three f(t_n)ftn were defined in (2.6).

For each time step, t_n, (2.19) is a nonlinear system of three equations in three unknowns, s₁(t_n+1), s₂(t_n+1) and s₃(t_n+1), which can be solved by the Newton-Raphson iteration of (2.16) making

φstn+1=s1tn+1−s1tn−hn2f1tn+1+f1tns2tn+1−s2tn−hn2f2tn+1+f2tns3tn+1−s3tn−hn2f3tn+1+f3tn=0(2.20a)

or in a compacted form

φstn+1=stn+1−stn−hn2ftn+1+ftn=0(2.20b)

and

Jφstn+1=1−hn2∂f1tn+1∂s1tn+1−hn2∂f1tn+1∂s2tn+1−hn2∂f1tn+1∂s3tn+1−hn2∂f2tn+1∂s1tn+11−hn2∂f2tn+1∂s2tn+1−hn2∂f2tn+1∂s3tn+1−hn2∂f3tn+1∂s1tn+1−hn2∂f3tn+1∂s2tn+11−hn2∂f3tn+1∂s3tn+1

=1+hn2Cg1Rs+gmtn+1hn2Cg1RL+gdstn+1hn2Cghn2Cg1Rs+Cb+CgCbgmtn+11+hnCb+Cg2CgCb1RL+gdstn+1hnCb+Cg2CgCb0−hn2Ld1(2.21a)

or

Jφstn+1=I−hn2∂ftn+1∂stn+1(2.21b)

where I stands for the identity matrix, and g_m(t) and g_ds(t) are, respectively, the FET’s transcoductance and output conductance, defined by gmt=∂iDvGvD∂vGvGt,vDt and gdst=∂iDvGvD∂vDvGt,vDt.

Figures 2.3 and 2.4 depict the results of the simulation of (2.16) for an excitation of v_I(t) = V_GG + v_S(t) = 3 + 20 cos (2πf_ct)vIt=VGG+vSt=3+20cos2πfct, f_c = 1 GHz, and with s₁(0) = 3.0 V, s₂(0) = 10.0 V and s₃(0) = −10/6.6 A initial conditions. This time-step integration used the trapezoidal rule of (2.19) and a fixed time step of h = 1/(100f_c).

If the frequency-domain response representation is desired, as is quite common in RF and microwave circuit analysis and design, then a DFT must be performed on the time-domain response, which requires that the signal is uniformly sampled and periodic.

Unfortunately, and as we have seen, simulation speed determines the use of dynamic time-steps, which produces a non-uniformly sampled response. So, either one forces the simulator to use a fixed time-step, e.g., using the simulator’s time-step ceiling option – turning the time-domain simulation highly inefficient – or relies on the hidden interpolation and resampling functions operated before the DFT is applied, which are not error free – turning the frequency-domain representation potentially inaccurate.

If the signal is not periodic, then the DFT will make it periodic, with a period equal to the simulated time, producing an error that is visible as spectral leakage [10], i.e., in which the energy of individual spectral lines is distributed onto (leaks to) their vicinity. This error is caused by the time-domain signal discontinuities that may exist in the beginning and end points of the constructed periodic signal and is illustrated in Figure 2.7. It is particularly harmful in distortion analysis where the distortion components may be totally immersed in the spectral leakage error. Hence, it is quite common to adopt measures to mitigate spectral leakage, increasing the simulation time (in order to reduce the relative weight of these discontinuities in the overall signal) or to simply eliminate these discontinuities, forcing the output signal to start and end in the same zero value – y(0) = y(T) = 0 – multiplying y(t) by an appropriate window function [10].

Figure 2.7 Illustrative example of the spectral leakage generated when an aperiodic time-domain function is (discrete) Fourier transformed to the frequency domain. Please note the spectral leakage manifested as erroneous spectral power around the desired signal and its mitigation by windowing. In this example, VdTω is the frequency-domain representation of five periods of the result of our transient simulation, v_D(t), collected between 50 ns and 55 ns, while VdHω is the corresponding signal passed through a Hanning window [10]. The dotted spectrum, VdSω, is the same five periods but now collected after the steady-state was reached, i.e., collected between 380 ns and 385 ns.

This lack of periodicity of the output signal can be due to the nature of the problem at hand. For example, a two-tone signal of f₁ and f₂ is only periodic if its central frequency f_c = (f₁ + f₂)/2 is a multiple of its frequency separation Δf = f₂ − f₁. But it may also be due to the presence of the transient response. In such a case, it is necessary to truncate the signal to its last few periods. This process is usually quite frustrating when simulating RF circuits since the length of this annoying transient response can be many orders of magnitude higher than the length of the desired period, or of a few periods. Indeed, the length of these long transients is usually determined by the time-constants associated to the large dc decoupling capacitors and RF chokes of bias networks.

Therefore, this transient time-domain simulation of microwave and RF circuits tends to be either inefficient, inaccurate, or both, which motivated the development of analysis techniques specifically devoted to calculate the periodic steady-state.

2.2.2 Periodic Steady-State Analysis in the Time-Domain

Although one might think that the integration of the transient response is unavoidable, as it corresponds to the way the real circuit behaves, that is not necessarily the case. To understand why, let us recall the previously discussed problem that consists in determining the response of a forced system to a periodic excitation of period T. To guarantee that this problem is solvable, let us also consider that the steady-state response is, indeed, periodic, of the same period T of the excitation.

Obviously, the complete response starts with the aperiodic transient, which is then followed by the periodic steady-state. This transient is naturally governed by the system, the excitation and the assumed initial conditions, or initial state, while the periodic steady-state is only dependent on the system and the excitation (in fact, its period is solely determined by the excitation). Therefore, different initial states lead to distinct transient responses, but always to the same steady-state. Actually, a close observation of Figures 2.3 and 2.4 indicates that different initial states may lead to transients of very different length, depending on whether we start from an initial state that is close, or very far, from the final cyclic orbit that describes the periodic steady-state. In fact, the notion of state tells us that all states of that cyclic orbit keep the system in the same orbit, i.e., lead to a transient response of zero length or, in other words, a response without transient.

Now that we are convinced that our problem of bypassing the transient response may have a solution, it becomes restricted to find a way to determine these appropriate initial states (see Exercise 2.2). There are, at least, three known techniques to do that.

The first method we can imagine to determine the proper initial states (i.e., the ones that are already contained in the cyclic orbit) was already discussed. Indeed, although the time-step integration starts from an initial state that may be very far from the desired one, the state it reaches after one period, is closer. Then, this new updated state is inherently used as the initial state for the time-step integration of the second period, from which a new update is found, in a process that we know is convergent to the desired proper initial state. Our problem is that this process is usually too slow due to the large dc decoupling bias network capacitors and RF chokes. In an attempt to speed up the transient calculation, a simple dc analysis (in which all capacitors and inductors are treated as open or short circuits, respectively) is performed prior to any transient simulation. Then, the dc open circuit voltages across the capacitors and the short-circuit currents across the inductors are used as the initial conditions of the corresponding voltages and currents in the actual transient simulation, unless other initial conditions are imposed by the user. These dc values constitute, therefore, the default initial state that the simulator uses. When these dc quiescent values are already very close to the actual dc bias values, as is the case of all linear or quasi-linear circuits, such as small-signal or class A power amplifiers, the transient response is short and this method of finding the default initial conditions can be very effective. Otherwise, i.e., in case strongly nonlinear circuits such as class C amplifiers or rectifiers are to be simulated, these default values are still very far from the actual bias point and a large transient response subsists.

The second and third methods find the proper initial condition using a completely different strategy: they change the nature of the numerical problem.

2.2.2.2 Finite-Differences in Time-Domain

Having this boundary-value problem in mind, the second method is known as the finite-differences in time-domain and is directly formulated from the discretization of the state-equation and the periodic boundary condition s(t₀ + T) = s(t₀). Using again a zero-order backward Euler discretization in N + 1 time-samples, we can define the following (incomplete) system (for each state-variable or each circuit node) of N equations in N + 1 unknowns

st1=st0+h0fst1xt1⋮stn+1=stn+hnfstn+1xtn+1⋮stN=stN−1+hNfstNxtN(2.22)

which, augmented by the boundary equation s(t_N) = s(t₀), constitutes a (now complete) system of N + 1 equations in N + 1 unknowns that can be readily solved through a N + 1-dimensional Newton-Raphson iteration scheme.

2.2.2.3 Shooting Method

The third alternative to find the proper initial condition that leads to the desired steady-state is known as the shooting or periodic steady-state, PSS, method. Conceptually, it can be described as a trial-and-error process in which we start by simulating the circuit for one period with some initial state guess; for example, the one given by an a priori dc analysis, s(t₀)₁st01. In our circuit example, this attempt would lead to the first set of trajectories depicted in the phase-space plots of Figure 2.8, in which s(t_N)₁≠s(t₀)₁stN1≠st01, indicating, as expected, that we failed the periodic steady-state. By slightly changing the tested initial condition by some Δs(t₀)Δst0, say s(t₀)₂ = s(t₀)₁ + Δs(t₀)st02=st01+Δst0, we can try the alternative attempt identified as the second set of trajectories of Figure 2.8. Although this one led to another failure, since s(t_N)₂stN2 is still different from s(t₀)₂st02, it may provide us an idea of the sensitivity of the final condition to the initial condition, Δs(t_N)/Δs(t₀)ΔstN/Δst0, which we then may use to direct our search for better trial s(t₀)₃st03, s(t₀)₃st03, …, s(t₀)_Kst0K, until the desired steady-state is achieved, i.e., in which s(t_N)_K = s(t₀)_KstNK=st0K.

Figure 2.8 Illustration of the shooting-Newton method of our circuit example excited by v_I(t) = 3 + 20 cos (2π10⁹t)vIt=3+20cos2π109t, where all the required five Newton iterations are shown. As done for the previous transient analysis, we started from the initial conditions derived from dc analysis, v_G(0) = 3.0 V, v_D(0) = 10.0 V and i_L(0) = −10/6.6 A, and ended up on v_G(0) = 3.118026121336410 V, v_D(0) = 11.183189430567937 V and i_L(0) = −2.625492967125298 A.

Actually, using a first order Taylor series approximation of the desired null error, s(t_N) − s(t₀) = 0stN−st0=0, we may write [11], [12]:

stNi+1−st0i+1≈stNi−st0i+ΔstNΔst0−Ist0i+1−st0i=0(2.23)

which leads to the following Newton-Raphson iteration:

st0i+1=st0i−ΔstNΔst0−I−1stNi−st0i(2.24)

known as the shooting-Newton PSS method. It should be noted that this sensitivity matrix, ΔstNΔst0, can be evaluated along the time-step integration of the period, using the chain differentiation rule as

ΔstNΔst0≡∂stN∂st0=∂stN∂stN−1∂stN−1∂stN−2⋯∂stn+1∂st0(2.25a)

where

∂stn+1∂st0=∂stn+1∂stn∂stn∂stn−1⋯∂st1∂st0(2.25b)

directly from the time-step integration formula of (2.18).

Example 2.4 The Shooting-Newton Formulation

In our circuit example, where we used the trapezoidal integration rule, this recurrent formula can be derived from (2.20) and (2.21) by noting that, if φ[s(t_n + 1)] = 0φstn+1=0, then

∂φstN∂st0=0(2.26)

and so

∂stn+1∂st0−∂stn∂st0−hn2∂ftn+1∂st0+∂ftn∂st0

=∂stn+1∂st0−∂stn∂st0−hn2∂ftn+1∂stn+1∂stn+1∂st0+∂ftn∂stn∂stn∂st0

=I−hn2∂ftn+1∂stn+1∂stn+1∂st0−I+hn2∂ftn∂stn∂stn∂st0

=Jφstn+1∂stn+1∂st0−2I−Jφstn∂stn∂st0=0

and finally

∂stn+1∂st0=Jφstn+1−12I−Jφstn∂stn∂st0(2.27)

Although the shooting-Newton PSS is based on an iterative process to find the appropriate initial state, it is, by far, much more efficient than the other two preceding methods for at least three reasons [12]. The first one is that, while (2.22) attempts to solve a N + 1-dimensional problem at once, (2.24) solves it resolving N + 1 one-dimensional problems, one at a time. The second is that the sensitivity, or the Jacobian, of (2.24) can be computed along with the transient simulation. Finally, the third is that although the circuit may be strongly nonlinear, the dependence of the final condition on the initial condition seems to be only mildly nonlinear in most practical problems, so that the Newton-Raphson solver converges in just a few iteration steps. For example, as is illustrated in the phase-state plots of Figure 2.8, the application of the shooting-Newton algorithm to our circuit example needed only five iterations to reach the desired periodic steady-state, even if we started from the same initial condition used for the transient analysis.

These properties make the PSS the technique of choice for efficiently computing the periodic steady-state regimes in time-domain of RF and microwave circuits. However, the inexistence of time-domain models for many RF and microwave components (mostly passive components such as filters and distributed elements, which are usually described by scattering parameter matrices) and the fact that, many times, RF/microwave designers operate on stimuli and responses represented as spectra, make frequency-domain methods, such as the harmonic-balance, the most widely used simulation techniques in almost all RF/microwave circuit designs.

2.3.1 Harmonic-Balance Engines

The application of the Fourier series representation to the analysis of nonlinear systems assumes that, for a broad set of systems (namely the ones that are stable and of fading memory [13]), their steady-state response to a sinusoid, or to a periodic excitation, of period T, is still a periodic function of the same period. However, as we already know from the harmonic generation property of nonlinear systems, this does not mean that the frequency content of the response is the same of the input. It simply means that it is located on the same frequency grid, ω_k = kω₀ωk=kω0. In general, the harmonic content of the output is richer than the one of the input, although there is, at least, one counter-intuitive example of this: linearization by pre- or post-distortion, in which a nonlinear function and its inverse are cascaded. Please note that this assumption of periodicity of the steady-state response of our system to a periodic excitation (in the same frequency grid) is a very strong one. If it fails (for example, if the system is chaotic or is unstable, oscillating in a frequency that falls outside the original frequency grid) frequency-domain techniques may not converge to any solution or, worse, may converge to wrong solutions without any warning.

The extension of the frequency-domain analysis of linear systems to nonlinear ones is known as the harmonic-balance, HB, since it tries to find an equilibrium between the Fourier components that describe the currents and voltages in all circuit nodes or branches. It is quite easy to formulate, although, as we will see, not so easy to solve. To formulate it, we simply state again the Fourier expansions of (2.29) and (2.30) and substitute them in the state-equation (2.1), to obtain

∑k=−KKjkω0Skejkω0t=f∑k=−KKSkejkω0t∑k=−KKXkejkω0t(2.34)

Within the stable and fading memory system restriction, we do know that the right-hand side of (2.34) will be periodic of fundamental frequency ω₀, and thus admits the following Fourier expansion:

f∑k=−KKSkejkω0t∑k=−KKXkejkω0t=∑k=−KKFkSωXωejkω0t(2.35)

in which S(ω) and X(ω) stand for the Fourier components vectors of s(t) and x(t), respectively. Substituting (2.35) in (2.34) results in

∑k=−KKjkω0Skejkω0t=∑k=−KKFkSωXωejkω0t(2.36)

which would not be of much use if it weren’t for the orthogonality property of the Fourier series. In fact, since the complex exponentials of the summations of (2.36) are orthogonal – which means that no component of e^jk₁ω₀tejk1ω0t can contribute to e^jk₂ω₀tejk2ω0t unless k₁ = k₂ − (2.36) can actually be unfolded as a (nonlinear) system of 2K + 1 equations (for each state-variable or circuit node), of the form

in 2K + 1 unknowns (the S_k). Written in vector-matrix form, (2.37) can be expressed as

in which jΩ is a diagonal matrix whose diagonal values are the jkω₀.

Please note that since there is no known method to compute the nonlinear mapping [S(ω),X(ω)]→F(ω) for any general function f[s(t),x(t)] or F[S(ω),X(ω)] directly in the frequency domain, one has to rely on the time-domain calculation, f[s(t),x(t)], converting first S(ω) and X(ω) to the time domain, and then converting the resulting f(t) = f[s(t),x(t)] back to the frequency domain. These frequency-time-frequency domain transforms constitute a heavy burden to the HB algorithm, which has motivated the quest for a pure frequency-domain alternative, nowadays known as spectral-balance [14].

The idea is to express f[s(t),x(t)] with only arithmetic operations, for which the frequency-domain mapping is known. For example, if f[s(t),x(t)] is expressed as a polynomial – which is nothing else but a summation of products – the time-domain addition corresponds to a frequency-domain addition, while the time-domain product corresponds to a frequency-domain spectral convolution (which can be efficiently implemented as the product of a Toplitz matrix by a vector [15]). Actually, this same idea can be extended to more powerful rational-function approximators (i.e., ratios of polynomials) since the time-domain division corresponds to a frequency-domain deconvolution, which is implemented as the product of the inverse of the Toplitz matrix just referred by a vector [15]. Although the spectral-balance has not seen the same widespread implementation of the traditional harmonic-balance, it is still interesting as a pure frequency-domain nonlinear simulation technique.

To solve the nonlinear system of equations of (2.38), we can again make use of the multidimensional Newton-Raphson algorithm, for which we need to rewrite (2.38) in its canonical form as

Then, the Newton-Raphson iteration would be

in which J_F[S(ω), X(ω)] stands for the Jacobian matrix of the frequency-domain function F[S(ω), X(ω)].

Unfortunately, the previous phrase may be misunderstood because of its apparent simplicity. In fact, if we said that we did not know any direct way to compute F[S(ω), X(ω)] in the frequency domain, we certainly do not know of any way to compute its Jacobian either. However, this Jacobian plays such an important role in nonlinear RF circuit analysis that we cannot avoid digging into its calculation. Therefore, to retain the essential, without being distracted by accessory numerical details, let us imagine we have a simple one-dimensional function of a one-dimensional state-variable: f(t) = f[s(t)] or F(ω) = F[S(ω)]:

FkSω=1T∫−T/2T/2fste−jkω0tdt(2.41)

The row m and column l of the wanted Jacobian matrix, J_F[S(ω)]_ml, is thus the partial derivative of the m’th Fourier component of F(ω), F_m, with respect to the l’th Fourier component of S(ω), S_l:

JFSωml≡∂Fm∂Sl=∂∂Sl1T∫−T/2T/2f∑k=−KKSkejkω0te−jmω0tdt(2.42)

Now, as the order of the derivative and the integral can be interchanged, using the derivative chain rule that allows us to write

∂∂Slfst=dfsds∂s∂Sl=gs∂∂Sl∑k=−KKSkejkω0t=gsejlω0t(2.43)

and thus

JFSωml≡∂Fm∂Sl=1T∫−T/2T/2gste−jm−lω0tdt(2.44)

which is the (m − l)’th Fourier component of the derivative of f(s) with respect to s, g(s).

A physical interpretation of this result can be given if we consider that, for example, f(s) represents a FET’s drain-source current dependence on v_GS, i_DS[v_GS(t)]. In this case, g[s(t)] is the FET’s transconductance and so (2.44) stands for the (m − l)’th Fourier component of the time-variant transconductance, g_m(t) = g_m[v_GS(t)], in the large-signal operating point, LSOP, [v_GS(t), i_DS(t)]. It constitutes, thus, the frequency-domain representation of the device operation linearized in the vicinity of its LSOP, an extension of the linear gain, or transfer function, when the device is no longer biased with a dc quiescent point, but with one that is time-variant and periodic.

Therefore, this Jacobian is the frequency-domain time-variant transfer function of a mixer whose LSOP is determined by the local-oscillator (maybe with also a dc bias), known in the mixer design context as the conversion matrix [16], [17]. Beyond this, the position of the poles of the Jacobian’s determinant provide valuable information on the large-signal stability of the circuit, in much the same way as the poles of the frequency-domain gain of a time-invariant linear amplifier – biased at a dc quiescent point – determine its small-signal stability [18].

But, this Jacobian can represent much more than this if a convenient linear transformation of variables is performed in the circuits’ currents and voltages to build incident and scattered waves. If F[S(ω)] is defined so that it maps frequency-domain incident waves, A(ω) = S(ω), into their corresponding scattered waves, B(ω) = F[A(ω)], the Jacobian would correspond to a form of large-signal S-parameters (available in some commercial simulators), or, more rigorously, to the S and T parameters of the poly-harmonic distortion, PHD, or X-Parameters model, [19], [20] (as will be discussed in Chapter 4).

As always, an important issue that arises when one is trying to find the solution of a nonlinear equation like (2.39) with the Newton-Raphson iteration of (2.40), is the selection of an appropriate initial condition. From the various possible alternatives we might think of, there is one that usually gives good results. It assumes that most practical systems are mildly nonlinear, i.e., that their actual response is not substantially distinct from the one obtained if the circuit is considered linear in the vicinity of its quiescent point. Therefore, the simulator starts by determining this dc quiescent point, setting to zero all ac excitations, and then estimates an initial condition from the response of this linearized circuit.

If this does not work, i.e., if it is found that the Newton-Raphson iteration cannot converge to a solution, then the simulator reduces the excitation level down to a point where this linear or, quasi-linear, approximation holds. Having reached convergence for this reduced stimulus, it uses this response as the initial condition for another simulation of the same circuit, but now with a slightly higher excitation level. Then, this source-stepping process is repeated up to the excitation level originally desired.

Example 2.5 The Harmonic-Balance Formulation

Continuing to use our circuit example of Figure 2.2, with the same periodic stimulus v_I(t) = 3 + 20 cos (2π10⁹t)vIt=3+20cos2π109t, we will now solve it for its periodic steady-state using the harmonic-balance formulation. For that, we will use the time-domain state equations of (2.5) and (2.6) to identify F[S(ω), X(ω)]FSωXω in (2.39) and (2.40) as

FSωXω=−1CgS1ω−X1ωRs+IdS1S2+S3ω+S2ωRLX1ω−S1ωCgRs−Cb+CgCbCgIdS1S2+S3ω+S2ωRLS2ω−X2ωLb(2.45)

The application of the Newton-Raphson algorithm to (2.39) requires the construction of the Jacobian matrix J_F[S(ω), X(ω)]. In this example of three state-variables, the Jacobian matrix will be a 3 × 3 block matrix whose elements are given by:

J_F[S(ω), X(ω)]_ml

JFSωXωml

=−−1Cg1Rs+Gmωm−l−1CgGdsωm−l+1RL−1Cg1CgRs−Cb+CgCgCbGmωm−l−Cb+CgCgCbGdsωm−l+1RL−Cb+CgCgCb01Lb0(2.46)

in which G_m(ω_m − l)Gmωm−l and G_ds(ω_m − l)Gdsωm−l are the (m − l)m−l frequency component of the device’s time-varying transconductance, g_m(t)gmt, and output conductance, g_ds(t)gdst, respectively.

Figure 2.9 illustrates the evolution of the harmonic content of the output voltage of our example circuit, |V_d(kω₀)| = |S₃(kω₀)|, during the Newton-Raphson solution of (2.40), while Figure 2.10 provides this same result but now seen in the time-domain, v_D(t) = s₃(t). Although the 3-D plot of Figure 2.9 presents only the first 15 harmonics plus dc, 50 harmonics were actually used in the simulation.

Figure 2.9 Illustration of the evolution of the first 15 harmonics, plus dc, of the amplitude spectrum of the output voltage, |V_d (kω₀)| = |S₃(kω₀)|, during the first 50 Newton-Raphson iterations of the HB simulation of our example circuit. Note how the simulation starts by assuming that the output voltage is a sinusoid (at the fundamental frequency, k = 1) plus a dc value (k = 0), and then evolves to build its full harmonic content.

Figure 2.10 Illustration of the evolution of the time-domain steady-state harmonic-balance solution of v_D(t) = s₃(t), during the whole Newton-Raphson iterations. Note how the simulation starts by assuming that the output voltage is a sinusoid plus a dc value, and then evolves to build the correct waveform shape.

2.3.2 Piece-Wise Harmonic-Balance

In this section we introduce what became known as the piece-wise harmonic-balance method, to distinguish it from the nodal-based harmonic-balance analysis technique just described. In fact, in piece-wise HB we do not start by a nodal analysis – from which we derived the state-equation – but begin by separating all of the circuit’s current, charge, or flux nonlinear sources from the remainder of the linear subcircuit network, as is illustrated in Figure 2.11. Then, we only analyze the exposed circuit nodes, i.e., the ones to which nonlinear voltage-dependent charge or current sources are connected. Let us, for the sake of language simplicity, designate these exposed nodes as nonlinear nodes and linear nodes the ones contained in the linear subcircuit network. (Please note that, despite not mentioning them here, nonlinear current-dependent flux sources could have also been considered, although requiring a slightly different formulation).

Figure 2.11 Decomposition of a nonlinear microwave circuit into a subcircuit linear block to which the nonlinear voltage dependent current or charge sources are connected, as needed for the application of the piece-wise harmonic-balance technique.

In piece-wise HB, the variables we will determine are no longer the circuit’s state-variables, but only the voltages of these exposed nodes, V(ω). Because most microwave circuits have typically many more linear nodes than nonlinear ones, the dimension of V(ω) is usually much smaller than the one of S(ω), which constitutes an enormous computational advantage of the piece-wise HB with respect to its nodal-based counterpart. Now, we apply frequency-domain nodal analysis to these nonlinear nodes obtaining the following desired HB equation:

In (2.47), I_nl[V(ω)]InlVω is the frequency-domain current vector representing the currents entering the nonlinear current sources, which, as before, must be evaluated in the time-domain as

InlkVω=1T∫−T2T2inlvte−jkω0tdt=1T∫−T/2T/2inl∑m=−KKVmejmω0te−jkω0tdt(2.48)

in the case of voltage-dependent nonlinear current sources, or as

InlkVω=jkω0Qnl,kVω=jkω01T∫−T2T2qnlvte−jkω0tdt=jkω01T∫−T/2T/2qnl∑m=−KKVmejmω0te−jkω0tdt(2.49)

in the case of voltage-dependent nonlinear charge sources (nonlinear capacitances).

I_l[V(ω)]IlVω stands for the current vector entering the linear subcircuit network, which is easily computed in the frequency domain by

in which Y(ω)Yω is the linear subcircuit network admittance matrix.

Finally, I_x[X(ω)]IxXω are the excitation currents that can be computed as

where the transfer matrix A_x(ω) can have either current gain or transadmittance gain entries, whether they refer to actual independent current source or independent voltage source excitations, respectively.

Please note that, while the nonlinear current sources, or the nonlinear subcircuit network, which is composed of only quasi-static elements, is treated in the time domain, the linear subcircuit, which contains all dynamic components (even the time derivatives of nonlinear charges), is directly treated in the frequency-domain.

This has an important advantage in terms of computational efficiency. Because static nonlinearities are treated in the time domain, they tend to require fast product relations, while they would need much more expensive convolution relations if dealt with in the frequency domain (as was seen when we discussed spectral-balance). On the other hand, linear dynamic elements necessitate convolution in the time domain, but only require product relations when treated in the frequency domain.

Finally, we should say that the piece-wise HB equation of (2.47) can be solved for the desired V(ω) vectors through the following Newton-Raphson iteration:

Vi+1ω=Viω−JInlViω+Yω−1.InlViω+YωViω+AxωXω(2.52)

For an application example of the piece-wise harmonic-balance method, see Exercise 2.4.

To summarize this description of the basics of the harmonic-balance – formulated in both its nodal and piece-wise form – Figure 2.12 depicts its algorithm in the form of a flowchart.

Figure 2.12 An illustrative flowchart of the harmonic-balance algorithm. Although this flowchart was formulated for HB applied to each node or state-variable, in its nodal form, it could be easily extended to be applied to only the nonlinear nodes, as is required for the piece-wise HB formulation.

2.3.3 Harmonic-Balance in Oscillator Analysis

Because oscillators are autonomous systems, i.e., systems in which there is no excitation, we would not think of them in the context of harmonic-balance. Indeed, the whole idea of HB is based on the a priori knowledge of the state-vector frequency set, something we are looking for in the oscillator analysis. However, a clever modification of HB allows it to be used for the determination of the steady-state regime of this kind of autonomous systems.

To understand how it works, let us consider the piece-wise HB formulation of (2.47), in which all but the dc components of the excitation vector, I_x[X(ω)]IxXω, are zero.

Please note that, dividing both sides of (2.53) evaluated at the fundamental by the corresponding (non-null) node voltage vector, would lead to

InlVω0Vω0+IlVω0Vω0=Ynlω0+Ylω0=0(2.54)

which can be recognized as the Kurokawa oscillator condition [21]. So, using HB for oscillator analysis sums up to solve (2.53) for V(ω)Vω. The problem is that not only (2.53) is insufficient since it involves, at least, one more unknown than the available number of equations, as we already know one possible solution of it, that we do not want. In fact, if (2.47) was a determined system, (2.53) can no longer be because the oscillation base frequency ω₀ is now also unknown. In addition, because I_x[X(ω)]IxXω is now zero, except for ω = kω₀ = 0, one possible solution of (2.53) is the so-called trivial, or degenerate, nonoscillatory condition where all but the dc components of V(ω)Vω are zero. In the following, we will deal with these two problems separately.

The first of these problems can be solved noting that if it is true that because I_x[X(ω)]IxXω is now zero (except for ω = 0), and so there is no way to a priori know the basis of the frequency set, it also true that the phases of V(ω)Vω are undetermined, as we lost the phase reference. Well, rigorously speaking, only the phase of the fundamental is undetermined, as the phases of the harmonics are referenced to it. Actually, what is undetermined is the (arbitrary) starting time, t₀, at which we define the periodic steady-state. In fact, since the system is time-invariant, it should be insensitive to any time shift, which means that if v(t) is a periodic solution of (2.53), than v(t − t₀) should also be an identical solution of that same equation. In the frequency domain, this means that the phase of the fundamental frequency ϕ₁ = −ω₀t₀ can be made arbitrary, while the phase of all the other higher-order harmonics are related to it as ϕ_k = −kω₀t₀ = −kϕ₁. So, (2.53) is still a system of (2K + 1) equations in 2K + 1 unknown vectors: the vector of dc values V(0), K vectors of amplitudes of the fundamental and its harmonics, (K − 1) vectors of their phases (ϕ₂…ϕ_K), plus the unknown frequency basis ω₀. Therefore, it can be solved in much the same way as before, using, for example, the Newton-Raphson iteration, except that now there is a column vector in the Jacobian that is no longer calculated as ∂F_m/∂S_l but as ∂F_m/∂ω₀.

The second problem consists in preventing the Newton-Raphson iteration on (2.53) to converge to the trivial, or degenerate, solution. As explained, this solution is the dc solution in which all but the 0’th harmonic components of V(ω)Vω are zero. One way to prevent this is to maximize the error for this particular trivial solution. This can be done normalizing the error of (2.53) by the root-mean-square voltage magnitude of non-dc V(ω)Vω:

InlVω+IlVω∑k=1KVk2=0(2.55)

In commercial HB simulators, the adopted approach uses the so-called oscillator-probe. Its underlying idea is to analyze the oscillator, not as an autonomous system, but as a forced one. Hence, the user should select a particular node in which it is known that a small perturbation can significantly disturb the oscillator and apply a special source (the oscillator-probe) to it, as is represented in Figure 2.13.

Figure 2.13 Oscillator analysis via harmonic-balance of a nonautonomous system using the oscillator-probe.

Such an appropriate node is one that is in the oscillator feedback loop or close to the resonator. The oscillator-probe is a voltage source whose frequency ω₀ is left as an optimization variable, and whose source impedance is infinite at all harmonics but the fundamental at which it is zero. So, this Z_S(ω) behaves as a filter which is a short-circuit to the fundamental frequency, ω₀, and a nonperturbing open-circuit to all the other harmonics kω₀, where k≠1. Then, the oscillator analysis is based on finding the so-called nonperturbation condition, i.e., to determine ω₀ and |V_S| that make I_S(ω₀) = 0. Under this condition, the circuit behavior is the same whether it has the oscillator-probe attached to it or not. So, although this solution was found as the solution of a nonautonomous problem, it coincides with the solution of the original autonomous circuit.

2.3.4 Multitone Harmonic-Balance

The frequency-domain analysis technique presented above assumes a periodic regime for both the excitation and the state-vector, so that these quantities can be represented using the discrete Fourier transform. Unfortunately, there are a significant set of steady-state simulation problems that do not obey this restriction. Some of them, as chaotic systems, are so obviously aperiodic that we would not even dare to submit them to any frequency-domain analysis. Their spectra are continuous and so there is no point in trying to find the amplitude and phase of a set of finite discrete spectral lines. There are, however, other situations that are still aperiodic and thus cannot be treated using the DFT, but that can still be represented with a finite set of discrete spectral lines. They are the so-called multitone or almost-periodic, regimes. For example, consider a signal composed by the addition of two sinusoidal functions:

This excitation can be either periodic or aperiodic, depending on whether or not ω₁ and ω₂ are two harmonics of the same fundamental frequency. That is, if ω₁ = k₁ω₀ and ω₂ = k₂ω₀, then x(t) is periodic, and it is aperiodic if there is no ω₀ that verifies both of these conditions. An equivalent way of stating x(t) periodicity is to stipulate that ω₂/ω₁ = k₂/k₁ is a rational number. Conversely, stating that x(t) is aperiodic is equivalent to stipulating that there are no two nonzero integer numbers, m₁ and m₂, that satisfy m₁ω₁ + m₂ω₂ = 0. Therefore, periodicity requires that ω₁ and ω₂ are commensurate (ω₂/ω₁ is a rational number), while we say that the aperiodic two-tone excitation is composed of two incommensurate frequencies (ω₂/ω₁ is an irrational number). In general, an excitation composed by a set of frequencies, ω₁, …, ω_Kω1,…,ωK is aperiodic if, and only if, there is no set of nonzero integers that verify

In this context, we will restrict our multitone excitation to the case of incommensurate tones, as the other one can always be treated with the above described harmonic-balance algorithm for periodic regimes.

Let us then consider the two-tone excitation of (2.56), in which ω₁ and ω₂ are incommensurate. If it is true that neither x(t) nor s(t) can be expressed as a Fourier series, it is as well true that we may write x(t) as

xt=∑k=−KKXkejωkt=∑k1=−KK∑k2=−KKXk1,k2ejk1ω1+k2ω2t(2.58)

in which, X_{k1, k2} = 0Xk1,k2=0 when |k₁| + |k₂|>1, for a pure two-tone excitation, and

st=∑k=−KKSkejωkt=∑k1=−KK∑k2=−KKSk1,k2ejk1ω1+k2ω2t(2.59)

for s(t). In fact, experience tells us that the vast majority of nonlinear systems responds to a set of discrete spectral lines at ω₁ and ω₂, with a set of discrete lines located at all harmonics of ω₁, all harmonics of ω₂ and all possible mixing products of ω₁ and ω₂. So, although we cannot affirm that ω_k = kω₀, because the response is aperiodic, and so neither (2.58) nor (2.59) are Fourier series, we can state that these signals can still be expressed as a sum of discrete spectral lines whose location is a priori known.

Naturally, as it happened with periodic harmonic-balance, a certain spectrum truncation must be established. However, while in the periodic regime this was easily set in a single axis of frequencies as |k|ω₀ ≤ Kω₀, in the almost-periodic case, since we have two harmonic limits, K₁ and K₂, there is no unique obvious choice. On the contrary, there are at least two possibilities often used.

The first one establishes that |k₁| ≤ K₁ and that |k₂| ≤ K₂, in which K₁ and K₂ need not to be equal. For example, in mixer analysis, the highest harmonic order for the large-signal local oscillator must be significantly higher than the one used for the small RF signal. If the mixing products, k₁ω₁ + k₂ω₂, are represented as points in a plane defined by the orthogonal axes k₁ω₁ and k₂ω₂, then these conditions define a rectangular figure, or box, limited by K₁ω₁ and K₂ω₂. This is illustrated in Figure 2.14 for K₁ = 2 and K₂ = 3. That is the reason why this spectrum truncation scheme is known as box truncation.

Figure 2.14 Box truncation frequency set for |k₁| ≤ 2 and |k₂| ≤ 3 (35 mixing products). “X” indicates a used position of the bidimensional spectrum truncation.

The second truncation scheme derives from the practical knowledge that typical spectra tend to have smaller and smaller amplitudes when the order of nonlinearity is increased. This means that, if it can be accepted that the highest harmonic order for ω₁ is, say, K₁, and the maximum order for ω₂ is K₂, then a mixing product at K₁ω₁ + K₂ω₂ will certainly have negligible amplitude and box truncation may waste precious computational resources. So, the alternative is to consider a maximum harmonic order K and then state that |k₁| + |k₂| ≤ K. As illustrated in Figure 2.15, this condition leads to a diamond-shaped coverage of the frequency plane above described, which is the reason why it is known as the diamond truncation.

Figure 2.15 Diamond truncation frequency set for |k₁| + |k₂| ≤ 3 (25 mixing products). “X” indicates a used position of the bidimensional spectrum truncation.

After stating that the response to be found is given by (2.59), under an appropriate spectrum truncation scheme, multitone HB determines the coefficients S_kSk, or S_{k1, k2}Sk1,k2, knowing that they can be related by

∑k1=−KK∑k2=−KKjk1ω1+kω2Sk1,k2ejk1ω1+k2ω2t=f∑k1=−KK∑k2=−KKSk1,k2ejk1ω1+k2ω2t∑k1=−KK∑k2=−KKXk1,k2ejk1ω1+k2ω2t(2.60)

In fact, there is no substantial difference between the multitone HB algorithm of (2.60) and its original periodic HB except that, now, we are impeded to evaluate the right-hand side of (2.60) in the time domain as we lack a mathematical tool to jump between time and frequency domains. Hence, what we will outline in the following are three possible alternatives to overcome this limitation.

2.3.4.1 Almost-Periodic Fourier Transform

The first multitone HB technique we are going to address takes advantage of the fact that, although neither x(t) nor s(t) are periodic, they can still be represented by a finite set of discrete spectral lines. They are, thus, named almost-periodic and its frequency-domain coefficients can be determined in a way that is very similar to the DFT and is named the almost-periodic Fourier transform, APFT [22]. To understand how, let us assume that, for simplicity, s(t) is a one-dimensional vector, s(t), defined in a discrete time grid of t = nT_s. If it can be expressed as (2.59) we may write

s1⋮sn⋮sN=e−jωKTs…1…ejωkTs…ejωKTs⋮⋮⋮⋮e−jωKnTs…1…ejωknTs…ejωKnTs⋮⋮⋮⋮e−jωKNTs…1…ejωkNTs…ejωKNTs.S−K⋮S0⋮Sk⋮SK(2.61a)

or

s_n = S_k

sn=ΓSk(2.61b)

In the periodic regime, NT_s can be made one entire period so that (2.61) corresponds to a square system of full rank where N = 2K + 1. In that case, the coefficients S_k can be determined by

S_k=⁻¹s_n

Sk=Γ−1sn(2.62)

which corresponds to the orthogonal DFT.

In the almost-periodic regime, the best we can do is to get an estimation of the S_k. As always, the estimation is improved whenever we provide it more information, which, in this case, corresponds to increasing the sampled time, NT_s, and thus the size of s_n. Actually, if we would make a test with only periodic regimes, but progressively longer periods (i.e., gradually tending to the almost-periodic regime), S_k would smoothly increase in size (because the frequency resolution is Δf = 1/T) and the result of (2.62) would always be exact. Unfortunately, we cannot afford to do this and have to limit the size of S_k. This leads to a nonsquare system in which N>>2K + 1, and the S_k can be estimated in the least-squares sense as

S_k = (Γ^†Γ)⁻¹Γ^†s_n

Sk=Γ†Γ−1Γ†sn(2.63)

in which Γ^† stands for the conjugate transpose or Hermitian transpose of Γ. Actually, because the DFT matrix is a unitary matrix, i.e., Γ^† = Γ⁻¹, the APFT of (2.63) reduces to the DFT of (2.62), under periodic regimes.

Because this almost-periodic Fourier transform constitutes an approximate way to handle the aperiodic regime, as if it were periodic, it is prone to errors. Therefore, many times, the following two alternative methods are preferred.

2.3.4.2 Multidimensional Fourier Transform

To treat the two-tone problem exactly, we need to start by accepting that it cannot be described by any conventional Fourier series. Instead, it must be represented in a frequency set supported by a bidimensional basis defined by the incommensurate frequencies ω₁ and ω₂. In fact, stating the incommensurate nature of ω₁ and ω₂ is similar to stating their independence, in the sense that they do not share any common fundamental. But, if ω₁ and ω₂ are independent, then θ₁ = ω₁t and θ₂ = ω₂t are also two independent variables, which means that (2.59) could be rewritten as a bidimensional Fourier series expansion [23].

Actually,

sn1,n2=∑k1=−K1K1∑k2=−K2K2Sk1,k2ejk1ω1n1Ts1ejk2ω2n2Ts2(2.64)

and

Sk1,k2=12K1+12K2+1∑n1=02K1+1∑n2=02K2+1sn1,n2ejk1ω1n1Ts1ejk2ω2n2Ts2(2.65)

define a bidimensional DFT pair, in which the time-domain data is sampled according to

Ts1=2π2K1+1ω1andTs2=2π2K2+1ω2(2.66)

Since s_{n1, n2}sn1,n2 is periodic in both 2π/ω₁ and 2π/ω₂, the bidimensional DFT coefficients S_{k1, k2}Sk1,k2can be computed as (2K₂ + 1) one-dimensional DFTs followed by another set of (2K₁ + 1) DFTs (actually, much more computationally efficient FFTs if 2K₁ + 1 and 2K₂ + 1 are made integer multiples of 2).

Since this bidimensional 2D – DFT (or, in general, the multidimensional DFT, MDFT) is benefiting from the periodicity of the waveforms in the two bases ω₁ and ω₂, it does not suffer from the approximation error of the almost-periodic Fourier transform. However, it is much more computationally demanding which limits its use to a small number of incommensurate tones. Moreover, it requires that the ω₁ and ω₂ are, indeed, two incommensurate frequencies. If they are not, then there will be certain mixing products, at ω_k1 = m₁ω₁ + m₂ω₂ and, say, ω_k2 = m₃ω₁ + m₄ω₂, that will coincide in frequency, and so whose addition should be made in phase, while the MDFT based harmonic-balance will add them in power, as if they were uncorrelated. This may naturally introduce simulation errors, whose significance depends on the amplitudes, and thus on the orders of the referred mixing products (higher order products, i.e., with higher |m₁| + |m₂|, usually have smaller amplitudes).

2.3.4.3 Artificial Frequency Mapping

To understand the underlying idea of artificial frequency mapping techniques, we first need to recognize two important harmonic-balance features.

The first one is that the role played by the Fourier transform in HB is restricted to solve the problem of calculating the spectral mapping imposed by a certain static nonlinearity, Y(ω) = F[X(ω)]. So, as long as this spectral mapping is correctly obtained, we do not care about what numerical technique, or trick, was used to get it. Actually, the difference between the traditional DFT-based HB and the spectral-balance is the most obvious proof of this.

The second feature worth discussing is less obvious. It refers to the fact that, since f[x(t)] is assumed to be static, the actual positions of the spectral lines of X(ω) are irrelevant for the calculation of Y(ω); only their relative positions are important. In fact, note that, if you have no doubt that the amplitudes of the output spectral lines of, e.g., a squarer, a₀ + a₁x(t) + a₂x²(t), subject to a tone at ω₁ of amplitude A₁, and another one at ω₂ of amplitude A₂, are a₀ + (1/2)a₂A₁² + (1/2)a₂A₂² at dc, a₁A₁ and a₁A₂ at ω₁ and ω₂, (1/2)a₂A₁² and (1/2)a₂A₂² at 2ω₁ and 2ω₂, a₂A₁A₂ at ω₁ − ω₂ and ω₁ + ω₂, it is because you know that the response shape (the amplitude of each frequency component) is independent of the frequencies ω₁ and ω₂. Otherwise, you would not have even dared to think of a result without first asking for the actual values of ω₁ and ω₂. Indeed, you know that you can calculate these amplitudes independently of whether the frequencies are, say, 13 Hz and 14 Hz, 13 GHz and 14 GHz, or even 13 GHz and 102 GHz. However, there is a fundamental difference in HB when you consider 13 GHz and 14 GHz or 13 GHz and 102 GHz as the former pair defines a periodic regime, allowing the use of the DFT, while the latter does not. Artificial frequency mapping techniques are algorithms that map an almost-periodic real time and frequency grid, (n₁T_s1,n₂T_s2) and (k₁ω₁,k₂ω₂), onto one that is periodic, in the artificial time, rT_s, and artificial frequency, sλ₀, before applying the IDFT. Then, after the f(rT_s,) = f[x(rT_s,)] and Y(λ) = DFT[f(rT_s,)] are calculated, they apply the inverse mapping to Y(λ) to obtain the desired Y(ω).

For example, if the original real frequency grid were defined over f₁ = 13f1=13 GHz and f2=102 GHz, and we were assuming box truncation up to K₁ = 2 and K₂ = 3, then, the artificial frequency mapped spectrum could be set on a normalized one-dimensional grid of λ₀ = 1, where k₁f₁ would be transformed onto k₁λ₀ and k₂f₂ onto (1 + 2K₁)k₂λ₀.

This illustrative artificial frequency mapping for box truncation is depicted in Figure 2.16 but many other similar maps have been proposed [24], [25].

Figure 2.16 Artificial frequency mapping technique for box truncated spectra defined by k₁| ≤ 2 and |k₂| ≤ 3.

2.4.1 Representation of Modulated Signals

The signals we are willing to address can be represented as an amplitude, a_x(t)axt, and phase, ϕ_x(t)ϕxt, modulated carrier of frequency ω₀, of the form

In communication systems, the entity x˜t≡axtejϕxt is the so-called low-pass equivalent complex envelope because it is a generalization of the real amplitude envelope a(t)at, but now including also the phase information, ϕ_x(t)ϕxt. Please note that, being a complex signal (i.e., whose imaginary part is not zero), it is not constrained to obey the complex conjugate symmetry of real signals. So, not only X˜ω≠X˜∗−ω, as there may be no relation between X˜ω and X˜−ω. However, as x(t) is still a real signal X(ω) = X^∗(−ω)Xω=X∗−ω.

Because the system is nonlinear, it will distort both the carrier and its modulation. Therefore, not only it will generate harmonics of ω₀, but also the envelope of each of these harmonics will carry distorted versions of the original amplitude, a_x(t)axt, and phase, ϕ_x(t)ϕxt. Hence, in general, we will have to assume that, even if the excitation could be restricted to (2.67), the state-variables (and thus any generalized modulated periodic excitation) will share the form of

st=a0t+∑k=−Kk≠0Kaktejkω0t+ϕkt(2.68)

and the objective of the simulation is to obtain the corresponding complex envelope amplitudes, a_k(t)akt, and phases, ϕ_k(t)ϕkt.

2.4.2 Envelope-Transient Harmonic-Balance

The underlying idea behind envelop-following techniques is to admit that, because of their distinct nature, the envelope and the carrier behave as if they were evolving in two independent, or orthogonal, time scales, τ₁ and τ₂. This way, any modulated excitation-vector, x(t), and state-vector, s(t), could be rewritten as

xτ1τ2=ax0τ1+∑k=−Kk≠0Kaxkτ1ejkω0τ2+ϕxkτ1(2.69)

and

sτ1τ2=a0τ1+∑k=−Kk≠0Kakτ1ejkω0τ2+ϕkτ1(2.70)

and the ordinary differential equation, ODE, in time that we used to model our system, the state-equation of (2.1), becomes the following multirate partial differential equation, MPDE, in τ₁ and τ₂ [27]

∂sτ1τ2∂τ1+∂sτ1τ2∂τ2=fsτ1τ2xτ1τ2(2.71)

The solution of this MPDE is a trajectory of s(τ₁,τ₂) whose physical significance is restricted to τ₁ = τ₂ = t.

At this point, the reader may be wondering where we are heading, since, to solve an already difficult problem, it seems we have complicated it even more. And we have to admit he is right. However, there is an important advantage hidden in (2.71): the dynamics with respect to the envelope and the carrier, i.e., with respect to τ₁ and τ₂, are decoupled so that we can now take advantage of the periodicity of s(τ₁,τ₂) in τ₂.

This means that we could try to solve (2.71) in the frequency domain for the periodic regime determined by the RF carrier and keep the time-domain representation for the aperiodic complex envelope. This way, we could rewrite (2.70) as

sτ1τ2=∑k=−KKSkτ1ejkω0τ2(2.72)

which, substituted in the MPDE of (2.71), would lead to

∑k=−KKdSkτ1dτ1ejkω0τ2+∑k=−KKjkω0Skτ1ejkω0τ2=f∑k=−KKSkτ1ejkω0τ2∑k=−KKXkτ1ejkω0τ2(2.73)

or

dSτ1ωdτ1+jΩSτ1ω=fSτ1ωXτ1ω(2.74)

(2.73) and (2.74) can be seen as a τ₁-varying HB equation, and also a 2K + 1 set of ODEs, which describe the evolution in the time τ₁ of the state-vector Fourier coefficients S(ω) or S_k. It constitutes, therefore, a time and frequency mixed-mode simulation technique.

Using, for example, a backward Euler discretization,

the ODE in τ₁ can be solved as a succession of harmonic-balance problems, in which the Fourier coefficients to be determined are the S(τ_1n + 1, ω)Sτ1n+1ω, while, because in each time-step the S(τ_1n, ω)Sτ1nω are known, they are treated as excitations.

This formulation – for obvious reasons known as the envelope-transient harmonic-balance, ETHB – can be interpreted as the solution of a set of periodic steady-state problems in which the excitation is evolving at the much slower pace of the envelope [28–30]. Actually, this process is what is illustrated in Figures 2.17 and 2.18. (As an example of application of ETHB to a circuit, see Exercise 2.6).

Figure 2.17 A simplified flowchart of the envelope-transient harmonic-balance algorithm.

Figure 2.18 An illustrative view of the envelope-transient harmonic-balance process. (a) Representation of the state-variable harmonic components in the carrier frequency and envelope time domains. (b) and (c) depict the dc and fundamental components of the excitation and the response, respectively.

Because it is assumed that, at each envelope sample, the system is kept in a periodic steady-state, the envelope must be slowly varying in comparison to the carrier, so that it can be considered frozen during, at least, one RF period. Hence, ETHB requires that the modulated excitations and the state-variables can be expressed as narrowband signals whose bands, centered at dc and each of the carrier harmonics, do not overlap. Furthermore, contrary to HB, which can handle devices with only a frequency-domain representation (such as a scattering matrix), this ETHB formulation requires a time-domain description for all system elements. If those time-domain models do not exist, they have to be created, which is prone to errors or even to numerical instabilities (namely, when the extracted time-domain models of passive blocks are not guaranteed to be passive). On the other hand, it correctly handles the nonlinearity and the system’s dynamics imposed on both the carrier and the envelope. For example, it can cope with the transients observed in radar pulsed amplifiers or in any other RF pulsed systems, the initial transients of oscillators or even the capture or hold-in processes of phase-locked-loops, PLLs. In fact, ETHB should be the simulation option for all moderately strong nonlinear systems subject to multirate stimuli that only exhibit periodicity in the fast time-scale and are aperiodic in the slow one. This is the case of envelope-modulated RF carriers but could be also applied to phase-noise problems, electro-thermal transient behavior, etc.

2.4.3 Envelope-Transient over Shooting

In the same way we could apply either time-domain shooting or frequency-domain harmonic-balance to determine the system’s RF periodic steady-state, we can now conceive a time-frequency envelope-transient harmonic-balance algorithm or a similar time-time envelope-transient over shooting, EToS [31]. This EToS takes the MPDE of (2.71) discretizes it in τ₁ and τ₂ (using, for example, a backward Euler scheme for both τ₁ and τ₂)

sτ1n1+1τ2n2+1−sτ1n1τ2n2+1hn1+sτ1n1+1τ2n2+1−sτ1n1+1τ2n2hn2=fsτ1n1+1τ2n2+1xτ1n1+1τ2n2+1(2.76)

and carries a transient analysis on τ₁ over successive periodic steady-state analyses using shooting-Newton on τ₂. This envelope following PSS should be preferred to ETHB whenever the periodic regime involves strong nonlinearities and so requires a very large number of harmonics. This is the case, for example, of many new digital assisted, or digital-like, RF circuits based on switching mode operation, such as envelope tracking amplifiers or polar transmitters [32], RF pulse-width modulators [33], and many other software defined radio circuits.

2.4.4 Another View of MDFT-Based Harmonic-Balance

For completeness, we could ask ourselves what would result from using the MPDE formulation in presence of a periodic envelope, i.e., what would result from the MPDE if both time scales were treated in the frequency domain. In that case, the solution would have to be expressed by

sτ1τ2=∑k1=−K1K1∑k2=−K2K2Sk1,k2ejk1ω1τ1ejk2ω2τ2(2.77)

which is essentially equivalent to the MDFT of (2.64). Now, substituting this solution, and the corresponding bidimensional excitation, into the MPDE of (2.71) results in

∑k1=−K1K1∑k2=−K2K2jk1ω1+k2ω2Sk1,k2ejk1ω1τ1+k2ω2τ2=f∑k1=−K1K1∑k2=−K2K2Sk1,k2ejk1ω1τ1+k2ω2τ2∑k1=−K1K1∑k2=−K2K2Xk1,k2ejk1ω1τ1+k2ω2τ2(2.78)

the MDFT-based harmonic-balance equation.

2.4.5 Circuit-Level and System-Level Co-Simulation

In wireless communications, we are also interested in simulating how the modulated data is treated by all blocks of the wireless system, conducting a so-called system-level simulation. Unfortunately, system-level simulators operate on modulated data and behavioral models of the RF or baseband circuits, and so do not support simulation at the circuit level. Envelope-following techniques are especially useful to bridge this gap, incorporating circuit-level simulation in these system-level simulators. While current and voltage waveforms, or spectra, are the desired results in circuit-level simulation, in system-level simulators we aim to predict how the information is processed, regardless of how it is represented by any particular physical magnitude. For that, we extract the information from the original real signal and then use it as the excitation of system-level models especially conceived to disregard the details of the physical magnitudes carrying that information. The process by which we extract the information from the real signal was already addressed in this text when, in Section 2.4.1, we introduced the concept of the complex envelope. The system-level models we are thinking of are low-pass equivalent behavioral models [34], whose excitation and response are their input and output complex envelopes, respectively.

System-level simulation is, therefore, much more computationally efficient than circuit-level simulation as it does not have to simulate each time point of the RF carrier but simply the time samples of the complex envelope, a much slower signal. Unfortunately, as in any other simulation environment, system-level simulators are limited in accuracy by their low-pass equivalent behavioral models. And if these are widely accepted as very accurate when used to represent linear blocks, as filters, they still constitute hot topics of research when they are supposed to represent nonlinear dynamic blocks such as RF power amplifiers [35] or I/Q modulators [36]. In these cases, we would like to simulate all the other blocks at the system-level, and keep, e.g., the RF power amplifier at the circuit-level. This form of system-level/circuit-level co-simulation platform is depicted in Figure 2.19 and is enabled by our studied envelope-following techniques.

Figure 2.19 System-level/circuit-level co-simulation platform of a wireless transmitter in which all blocks except the RF power amplifier are represented as low-pass behavioral models for the complex envelope, whereas the power amplifier is represented at the circuit-level by its equivalent-circuit model for the real RF modulated signal.