## Note Essay

### Essay preview

Simulation of RF integrated circuits

2013

2

Contents

1

2

Formulation of Circuit Equations
1.1 Modiﬁed Nodal Analysis . . . . . . . . . . . . . .
1.1.1 Resistor Stamps . . . . . . . . . . . . . . .
1.1.2 Stamps of Independent Current Sources . .
1.1.3 Capacitor Stamp . . . . . . . . . . . . . .
1.1.4 Stamp of Independent Voltage Source (IVS)
1.1.5 Inductor Stamp . . . . . . . . . . . . . . .
1.1.6 Stamp of nonlinear conductance . . . . . .
1.1.7 Stamp of nonlinear capacitors and inductors

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

3
3
6
7
7
11
12
15
17

Harmonic Balance for Circuits with Single-Tone Excitation
2.1 Background on Fourier Transform for Single-Tone Signals . . . . 2.1.1 Discrete Fourier Transformation . . . . . . . . . . . . . . 2.1.2 Fast Fourier Transformation . . . . . . . . . . . . . . . . 2.1.3 Numerical Example . . . . . . . . . . . . . . . . . . . . 2.2 Harmonic Balance for Linear Circuits . . . . . . . . . . . . . . . 2.2.1 Steady-state response for a simple circuit . . . . . . . . . 2.2.2 Steady-state response for general linear circuits using HB 2.3 Harmonic Balance for Nonlinear Circuits . . . . . . . . . . . . . 2.3.1 General nonlinear circuits . . . . . . . . . . . . . . . . . 2.3.2 The N-R method . . . . . . . . . . . . . . . . . . . . . . 2.3.3 Solving the HB equations using N-R method . . . . . . . 2.3.4 Computational complexity of HB for nonlinear circuits . . 2.4 Implementation Issues . . . . . . . . . . . . . . . . . . . . . . . 2.4.1 Convergence of the NR method . . . . . . . . . . . . . . 2.4.2 Source Stepping . . . . . . . . . . . . . . . . . . . . . . 2.4.3 Homotopy Methods . . . . . . . . . . . . . . . . . . . . 2.4.4 Error Dynamics . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

23
24
25
27
28
29
29
31
33
35
41
42
53
56
56
57
57
58

3

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.

CONTENTS

3

4

5

Harmonic Balance for Circuits with Multi-Tone Excitation
3.1 Multi-tone and Almost-Periodic Waveforms . . . . . . . . . . . 3.2 Response of a Nonlinear Circuit to Almost-Periodic Excitations 3.3 Truncating the Spectrum of Almost-Periodic signals . . . . . . . 3.3.1 Box truncation for two-tone waveform . . . . . . . . . . 3.3.2 Diamond truncation for two-tone waveforms . . . . . .

3.3.3 Box and diamond truncation for multi-tone waveforms .
3.4 Harmonic Balance Analysis . . . . . . . . . . . . . . . . . . . 3.5 Generalized Fourier Transform . . . . . . . . . . . . . . . . . . 3.5.1 Multi-Dimensional Fourier Transform (MDFT) . . . . .

3.5.2 Almost-Periodic Fourier Transform (APFT) . . . . . . . 3.5.3 Frequency Mapping Techniques (FMT) . . . . . . . . .
3.5.4 Implementation Procedures for FMT . . . . . . . . . . . Time Domain Methods
4.1 Steady-State Analysis using IVP approaches . . . . . . . . . 4.1.1 Illustrative Examples . . . . . . . . . . . . . . . . . 4.1.2 Transient Analysis of Circuits using BE and TR . . .
4.1.3 Steady-State Analysis via IVP . . . . . . . . . . . .
4.2 Steady-State Analysis using Boundary-Value Problem (BVP) 4.2.1 Shooting Method . . . . . . . . . . . . . . . . . . .
4.2.2 Newton Shooting Method (NSM) . . . . . . . . . .

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.

65
65
67
69
70
72
73
74
77
81
83
84
91

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

95
96
99
106
107
109
109
110

Volterra Series
5.1 Linear Systems and Transfer Function . . . . . . . . . . . . . . . . 5.1.1 Steady-State Response in Linear Systems via (TF) . . . . . 5.2 Volterra-Series and Small-Signal Analysis . . . . . . . . . . . . . . 5.2.1 First-Order Volterra Operator . . . . . . . . . . . . . . . . 5.2.2 Second-Order Volterra Operator . . . . . . . . . . . . . . . 5.2.3 Third-Order Volterra Operator . . . . . . . . . . . . . . . . 5.2.4 General nth -order Volterra Kernel . . . . . . . . . . . . . . 5.3 Constructing a Response to a Sinusoidal input using Volterra Series 5.3.1 Small-Signal Response using Volterra-Series . . . . . . . . 5.4 The Method of Nonlinear Currents . . . . . . . . . . . . . . . . . . 5.4.1 Analytical Expressions for the Frequency-domain kernels . 5.5 Extension to Multi-Tone Analysis . . . . . . . . . . . . . . . . . . 5.5.1 Computing vo1 (t) . . . . . . . . . . . . . . . . . . . . . . . 5.5.2 Discussion on the Method of Nonlinear Currents . . . . . . 5.6 Application of the Method of Nonlinear Currents to General Circuits

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.
.
.
.
.
.
.

117
117
118
122
123
123
123
124
124
126
130
130
140
141
143
144

1

.
.
.
.
.
.
.

.
.
.
.
.
.
.

.
.
.
.
.
.
.

CONTENTS

2

Chapter

1

Formulation of Circuit Equations
The main objective in this chapter is to introduce the basic idea used to automatically represent general nonlinear circuits in the mathematical domain. This representation is the very ﬁrst step needed to apply almost all of the simulation techniques currently implemented in commercial circuit simulators (e.g., frequency-domain, time-domain, steady-state, distortion analysis, etc.). The techniques of representing a general circuit in the mathematical domain have matured over two or three decades ago, and all of modern circuit simulators such as HSPICE, ADS have adopted them. The basic approach used to represent circuits mathematically is known as the Modiﬁed Nodal Admittance (MNA) approach is the main subject of this chapter.

1.1

Modiﬁed Nodal Analysis

Consider the circuit shown in Figure 1.1. We will ﬁrst derive the nodal equations for this circuit. Using these derivations we will be able to deduce the roles of each element in the general formulation. The knowledge of each element’s role will then enable us to derive the systematic approach employed in the MNA technique to represent general circuits. By writing Kirchhoff Current Law (KCL) at each node in the circuit, we obtain 4 differ3

1.1. MODIFIED NODAL ANALYSIS

2

3
4

1

0

Figure 1.1: A simple resistive circuit.
ent equations,
Currents @ node 1 = 0
Currents @ node 2 = 0
Currents @ node 3 = 0
Currents @ node 4 = 0
Note that we did not include a nodal equation at the reference since the voltage at this node is assumed to be zero by convention. Using the currents designated at each node and substituting in the above equations yield,

2

i1,k = 0
k=1
3

i2,k = 0
k=1
2

i3,k = 0
k=1
2

i4,k = 0
k=1

The currents in the above equations can be expressed in terms of the nodes voltages using 4

CHAPTER 1. FORMULATION OF CIRCUIT EQUATIONS

simple Ohm’s law. For example,
i1,2 =

1
(v1 − v2 )
R1

while i1,1 is related to the current from the independent current source, i1,1 = −iS1
Substituting for the currents in terms of the nodes voltages gives us 1
(v1 − v2 ) − iS1
R2
1
v2
v2 − v3
(v2 − v1 ) +
+
R1
R2
R3
1
(v3 − v2 ) − iS2
R3
1
v4 + iS2
R4

= 0
= 0
= 0
= 0

The above equations can be arranged in matrix/vector form as follows, 
 1
1

 
− R1
0
0
R1
iS1
 v1

1
1
1
1
1
+ R2 + R3
− R3 0   v2   0 
 − R1
R1


= i
1
1
S2
 0
− R3
0  v3
R3
v4
−iS2
1
0
0
0
R4

(1.1)

or more succinctly,
Gx = b
where

G=

1
R1
1
− R1

0
0

(1.2)

1
− R1
1
R1

+

1
R2
1
R3

0

5

+

0
1
R3

1
− R3
1
− R3
0

0

0 

0 

1
R4

(1.3)

1.1. MODIFIED NODAL ANALYSIS

and

−iS1
 0 
b= i

S2
−iS2

v1
 v 
x =  v2 
3
v4

(1.4)

Obviously, solving for the unknown variables x will give us all nodes voltages, which in turn will enable computing all branch currents.
We would like now to examine the role of each element in the ﬁnal formation of the circuit equations in (1.1). This will enable us to infer the “stamp” that each element leaves on the equations formulations.

1.1.1

Resistor Stamps

Consider the resistor R1 in the circuit of Figure 1.1. Such a resistor is connected between nodes 1 and 2. The effect that this resistor has on the formulation is having 1/R1 added to the entries at (1, 1) and (2, 2) of the G matrix. Also the existence of this resistor caused a value of 1/R1 to be subtracted from the entries at (1, 2) and (2, 1). Similarly resistor R3 (connected between nodes 2 and 3) inﬂuenced the G matrix by adding 1/R3 to the entries at (2, 2) and (3, 3) and subtracting 1/R3 from the entries at (2, 3) and (3, 2).

This observation can be further generalized to any resistor having a resistance R Ω and connected between nodes i and j. The stamp that results from this general resistor is the addition of 1/R to the entries at (i, i) and (j, j) and the subtraction of 1/R from the entries of (i, j) and (j, i), respectively. Figure 1.2 illustrates this notion graphically.

i

j

Figure 1.2: Resistor stamp.
However, if one of the nodes happened to be the reference node, such as in the case of the resistor R4 , then the stamp reduces to only adding 1/R to the (k, k) entry of the G matrix, as shown in Figure 1.3

6

CHAPTER 1. FORMULATION OF CIRCUIT EQUATIONS

2

0

Figure 1.3: Resistor stamp, where one node is the reference node. 1.1.2

Stamps of Independent Current Sources

By examining Eq. (1.1), we can also deduce the “stamp” that an independent current source leaves on the MNA formulation. This stamp appears only on the right-side vector b, and therefore this vector is usually referred to as the source vector. As can be observed from Eq. (1.1), the current source iS2 , which is connected between nodes 3 and 4 has left its Ampere value of iS2 at the 3rd and 4th entries of the source vector b. Figure 1.4 demonstrates the general rule for including a general independent current source connected between two general nodes. Figure 1.5 also shows the stamp of an independent current source when one node is the reference node.

i

j

Figure 1.4: Stamp of independent current sources.

1.1.3

Capacitor Stamp

We now derive the stamp of a general capacitor in a similar manner by ﬁrst formulating the circuit equations for the circuit shown in Figure 1.6. Writing KCL at each node gives, 7

1.1. MODIFIED NODAL ANALYSIS

k

0

Figure 1.5: Stamps of independent current sources, where one of the nodes is a reference node.

2
1

3

0

Figure 1.6: A circuit with a capacitive element.

8

CHAPTER 1. FORMULATION OF CIRCUIT EQUATIONS

i1,1 + i1,2 = 0
i2,1 + i2,2 + i2,3 = 0
i3,1 + i3,2 = 0

(1.5)

Substituting for each current in terms of the nodes voltages, we get i1,1 = −Is1
1
(v1 − v2 )
i1,2 =
R1
1
i2,1 =
(v2 − v1 )
R1
dv2
i2,2 = C1
dt
dv2 dv3
i2,3 = C2 (

)
dt
dt
dv3 dv2
i3,1 = C2 (

)
dt
dt
v3
i3,2 =
R2
Substituting from the above equations into Eqs. (1.5) and arranging in matrix form, we get  1

 dv 
1
1
− R1 0
is1
v1
0
0
0
R1
dt
1
1
 −R
0  v2 + 0 C1 + C2 −C2  dv2  =
0
(1.6)
R1
dt
1
1
dv3
v3
0
−C2
C2
0
0
0
R2
dt
or in a more compact form
Gx(t) + C

dx(t)
= b(t)
dt

(1.7)

where
x(t) =

C =

v1
v2
v3

G=

1
R1
1
− R1

1
− R1

0

0

0
0
0
0 C1 + C2 −C2
0
−C2
C2

1
R1

0
0 

(1.8)

1
R2

b(t) =

i s1
0
0

(1.9)

One obvious conclusion from Eq. (1.6), is that capacitors do not have any effect on the 9

1.1. MODIFIED NODAL ANALYSIS

G matrix. Instead they appear in another matrix, denoted by C matrix, which is multiplied by the derivatives of the nodes voltages. It is also easy to see from the structure of the C matrix that capacitor C2 , which is connected between nodes 2 and 3, was added to the entries (2, 2) and (3, 3) and subtracted from the entries at (2, 3) and (3, 2) of the C matrix. We can generalize this observation to any general capacitor as shown in Figures 1.7 and 1.8.

j

i

Figure 1.7: Illustrating the stamp of a capacitor on the MNA formulation.

k

Figure 1.8: Illustrating the stamp of a capacitor on the MNA formulation (One node is the reference node).

10

CHAPTER 1. FORMULATION OF CIRCUIT EQUATIONS

1.1.4

Stamp of Independent Voltage Source (IVS)

Independent Voltage Sources (IVS) require special attention in the MNA formulation. To illustrate how IVS’s are handled in the MNA formulation, consider the circuit shown in 1.9. Writing KCL at each node gives us

1

2

Figure 1.9: A circuit example with a independent voltage source. i1,1 + i1,2 = 0
i2,1 + i2,2 = 0

(1.10)

Substitution using Ohm’s law,
i1,1 = Is
1
R
1
= (v2 − v1 ) ×
R
dv2
= C1
dt

(1.11)

i1,2 = (v1 − v2 ) ×

(1.12)

i2,1

(1.13)

i2,2

(1.14)

Hence,
v1 v2

= 0
R
R
v1 v2
dv2
− +
+ C1
= 0
R
R
dt
i1,1 +

11

(1.15)
(1.16)

1.1. MODIFIED NODAL ANALYSIS

Arranging in matrix form gives us
1
R
1
−R

1
−R 1
1
0
R

v1
v2
Is

0 0 0
0 C1 0

+

dv1 /dt
dv2 /dt
dIs /dt

=0

(1.17)

The above system of equations is unbalanced, in the sense that there are more unknowns (v1 , v2 , i1,1 ) than equations. To overcome this problem, we need one more equation, which can be obtained by noting that

v1 = VS
(1.18)
Adding Eq. (1.18) to the system of equations in (1.17) yields 

1
1
−R 1
dv1 /dt
v1
0 0 0
R
1
 −1
 v2 + 0 C1 0
dv2 /dt
0
R
R
Is
0 0 0
dIs /dt
1
0 0

0
0
VS

=

(1.19)

which is a square system, i.e., a system in which the number of unknowns is equal to the number of equations.
Using the above example, we would like to derive the effect that an IVS has on the general MNA formulation. The ﬁrst observation is that having an IVS increments the number of equations by 1. In additions, the value of the voltage source appears in the right hand side vector. Figures 1.10 and 1.11 summarizes the inﬂuences of IVS on the MNA formulation. Note the +1 and −1 entries in the G matrix.

Extra Column to account for
the current in the source
i

j

Extra Row added for the extra equation

Figure 1.10: Incorporating independent voltage source in the MNA formulation. The voltage source is connected between two general nodes.

1.1.5

Inductor Stamp

To derive the inductor’s stamp, consider the circuit shown in Figure 1.12. Applying KCL 12

CHAPTER 1. FORMULATION OF CIRCUIT EQUATIONS

Extra Column to account for
the current in the source
k

the extra equation

Figure 1.11: Incorporating independent voltage source in the MNA formulation. The voltage source is connected between a general node and the reference node.

2
1

Figure 1.12: A circuit with an inductor.

13

1.1. MODIFIED NODAL ANALYSIS

at nodes 1, and 2 results in the following two equations
i1,1 + i1,2 = 0
i2,1 + i2,2 + i2,3 = 0

(1.20)
(1.21)

Given that using Laplace elements
I1,1 = −IS
1
(V1 − V2 )
I1,2 =
sL
1
I2,1 =
(V2 − V1 )
sL
1
I2,2 =
V2
R
I2,3 = sCV2
Substituting in (1.20), rearranging in matrix form results in the following system of equations 0 0
0 1/R

V1
V2

0 0
0 C

+s

V1
V2

+

1
s

1/L −1/L
−1/L 1/L

V1
V2

=

Is
0
(1.22)

The above formulation presents a particular problem when transformed to the timedomain, t
dx(t)
+L
x(τ )dτ = 0,
(1.23)
Gx(t) + C
dt
0
since the resulting system is of the integro-differential equations, which is more difﬁcult to deal with than a system of differential equations. Therefore, it would be much better if we can ﬁnd an alternative method to incorporate the inductor while keeping the differential equations form.

To this end, we add the current in the inductor as an additional unknown variable. With this modiﬁcation we have
I2,1 = −IL

I1,2 = IL ,

(1.24)

Substituting from (1.24) into (1.20), rearranging in matrix form gives us 0
0
1
0 1/R1 −1

V1
V2
IL

+s
14

0 0 0
0 C1 0

V1
V2
IL

=

Is
0
0

(1.25)

CHAPTER 1. FORMULATION OF CIRCUIT EQUATIONS

Extra Column to account for
the current in the inductor

i

j

the extra equation

Figure 1.13: Inductor’s stamp. The inductor is connected between two general nodes in the circuit.
Since we have one less equation than the number of unknowns, we will need one more equation. That equation can be obtained through applying Ohm’s law on the inductor as follows,
V1 − V2 = sIL L
(1.26)
Adding Eq. (1.26) to the system of Eqs. in (1.25) results in the following system of equations
0
0
1
0 1/R1 −1
1 −1
0

V1
V2
IL

0 0 0
0 C 0
0 0 −L

+s

V1
V2
IL

=

Is
0
0

(1.27)

in which the number of equations is equal to the number of unknowns. It can also be observed that transforming the system in (1.27) to the time-domain results in the following system of differential equations.

We are now in a position to deduce the general technique that enables automatic incorporation of the inductor in the general MNA formulation. First, it is clear that having the inductor led to adding one more unknown (the inductor’s current) with an additional equation. Furthermore, the inductance value appeared in the extra entry added to the C matrix. Figures illustrate the general stamp of an inductor.

1.1.6

Stamp of nonlinear conductance

Figure 1.15 shows a circuit with a nonlinear conductance whose current is controlled by a quadratic nonlinearity function of the applied voltage. Writing KCL at nodes 1 and 2 gives i1,1 + i1,2 = 0

i2,1 + i2,2 + i2,3 = 0
15

1.1. MODIFIED NODAL ANALYSIS

Extra Column to account for
the current in the inductor
k

Extra Row added for the extra equation

Figure 1.14: Inductor’s stamp. The inductor is connected between one node and the reference node.

Substituting for the currents using
i1,1 = −IS
i1,2 = −i = −f (v2 − v1 ) = − (v2 − v1 )2
i2,1 = i = f (v2 − v1 ) = (v2 − v1 )2
v2
i2,2 =
R
dv2
i2,3 = C
dt
and rearranging in matrix form produces the following system of equations 0 0
0 1/R

v1
v2

+

0 0
0 C

dv1
dt
dv2
dt

+

−(v2 − v1 )2
(v2 − v1 )2

=

Is
0

(1.28)

It is straightforward to see from Eq. (1.28) to see that the presence of a nonlinear conductance can be accounted for by a adding a new vector whose entries are dedicated to nonlinear functions. This vector will be denoted henceforth by f (x(t)), where x(t) is the vector of the MNA variables. Figure 1.16 shows the stamp that a nonlinear conductance leaves on the MNA formulation. In general having a nonlinear conductance characterized by an analytic function g(· · · ) connected between nodes i and j results in having g(·) added to/subtracted from the i-th entry and subtracted from/added to the j-th entry of the nonlinear function vector f (x(t)). Notice that the stamp developed here can be used to represent a large class of two-terminal devices. One notable example is the diode whose 16

CHAPTER 1. FORMULATION OF CIRCUIT EQUATIONS

2
1

Figure 1.15: A circuit with nonlinear conductance.
current is controlled by the applied voltage through the exponential function, i = Io ev/VT − 1

(1.29)

Figure 1.16: The MNA stamp of nonlinear conductance.
The general formulation developed so far for s general nonlinear circuit takes the following form dx(t)
Gx(t) + C
+ f (x(t)) = b(t)
(1.30)
dt
1.1.7

Stamp of nonlinear capacitors and inductors

Figure 1.17 shows the circuit symbol that is often used to represent a nonlinear capacitor. A nonlinear capacitor is a capacitor in which the charge is a nonlinear function of the applied 17

1.1. MODIFIED NODAL ANALYSIS

voltage,
q = g(v)

(1.31)

where g(v) can be any arbitrary nonlinear function. This is to be contrasted to a linear capacitor whose charge is linearly related to the applied voltage via the capacitance. q = Cv

(1.32)

Figure 1.17: Circuit symbol for nonlinear capacitor
Since the current in the capacitor is equal to the rate of ﬂow charge, the i-v relation for the nonlinear capacitor should be given by
i=

d
dq
= g(v)
dt
dt

(1.33)

To show how a nonlinear capacitor is incorporated in the general MNA formulation, we consider the circuit shown in Figure 1.18.
Applying KCL at each node results in the following two equations, Is + i = 0
v2
dv2
+C
+ f (v2 ) − i = 0
R
dt

(1.34)
(1.35)

where i is the current ﬂowing in the nonlinear capacitor. Substituting i with in the above two equations and arranging them in matrix form, we obtain 0 0
0 1/R

v1
v2

+

0 0
0 C

dv1
dt
dv2
dt

+

0
f (v2 )

d
+
dt
18

g(v1 − v2 )
−g(v1 − v2 )

=

Is
0

d
dt

(g(v1 − v2 ))

(1.36)

CHAPTER 1. FORMULATION OF CIRCUIT EQUATIONS

or in a more concise form,
Gx(t) + C

dx(t)
d
+ f (x(t)) + g(x(t)) = b
dt
dt

(1.37)

Figure 1.18: A circuit with nonlinear capacitor
It can be seen from the above equations that the presence of a nonlinear capacitor has forced one more vector of nonlinear functions, whose components are to be differentiated with respect to time t. That may not be a desirable feature, especially for some of the techniques used in Chapter 3, where the nonlinearities will be required to appear in pure algebraic rather than the differential form of the fourth term (1.37). Overcoming this difﬁculty is typically done by treating the charge in the nonlinear capacitor as an additional unknown variable in the circuit, and adding one more equation to the system. The addition equation is basically the constitutive equation for the nonlinear capacitor given in (1.31). Thus, for the example circuit of Figure 1.18, we should have three equations given by

dq
= 0
dt
v2
dv2
dq
+C
+ f (v2 ) −
= 0
R
dt
dt
q − g(v2 ) = 0
Is +

Arranging in matrix form results in the following system of equations  dv1 
0
0 0 0
v1
0 0 +1
dt
f (v2 )
v2 + 0 C −1  dv2  +
0 1/R 0
=
dt
dq
q
0 0 0
0 0 1
−g(v2 )
dt

19

(1.38)
(1.39)
(1.40)

Is
0
0

(1.41)

1.1. MODIFIED NODAL ANALYSIS

which is equivalent to
dx(t)
+ f (x(t)) = b
(1.42)
dt
where x(t) is the vector of the MNA variables which contains the charge waveform in the nonlinear capacitor in addition to other unknown variables. Based on the above analysis, it is possible to deduce the stamp that a nonlinear capacitor will leave on the general MNA formulation as shown in Figures 1.21 and 1.22.

Gx(t) + C

Extra Column to account for
the charge in the capacitor

i

j

the extra equation

Figure 1.19: Stamp of a nonlinear capacitor connected between nodes i and j.

Extra Column to account for
the charge in the capacitor

k

the extra equation

Figure 1.20: Stamp of a nonlinear capacitor connected between node k and the reference node.
The above analysis can be used in an analogous manner to deduce the MNA stamp of a nonlinear inductor. In a nonlinear inductor, the ﬂux is a nonlinear function of the applied voltage,
φ = u(v)
(1.43)
20

CHAPTER 1. FORMULATION OF CIRCUIT EQUATIONS

where u can be any arbitrary nonlinear function. Using the fact that the current in the inductor is the rate of ﬂux change, the i-v relation in the inductor is given by i=

d
u(v)
dt

(1.44)

The MNA stamp for the inductor is shown in Figures 1.21 and 1.22. Extra Column to account for
the flux in the inductor

i

j

the extra equation

Figure 1.21: Stamp of a nonlinear inductor connected between nodes i and j.

Extra Column to account for
the charge in the capacitor

k

the extra equation

Figure 1.22: Stamp of a nonlinear inductor connected between node k and the reference node.

21

1.1. MODIFIED NODAL ANALYSIS

22

Chapter

2

Harmonic Balance for Circuits with
Single-Tone Excitation
This chapter introduces the Harmonic Balance (HB) approach for ﬁnding the steady-state response of circuits with single-tone excitations. The term “single-tone” excitation here refers to the fact that the stimulus excitation waveforms considered in this chapter contain a single angular frequency and integer multiples of it. It is well known that such a waveform can be represented by an inﬁnite Fourier series expansion of the form, ∞

Uk ekω0 t

u(t) =

(2.1)

k=−∞

where ω0 is known as the fundamental (angular) frequency, and Uk is the coefﬁcient of the k-th harmonic. The fundamental period T is the period of time for which u(t) = u(t + T ), 2π
and is given in terms of the fundamental angular frequency by T = ω0 . If u(t) is a real signal, it is easy to show that Um = U−m , where represents the complex conjugate operation.
Another equivalent representation for single-tone periodic waveforms can be provided using sinusoidal waveforms as follows,

C
S
Uk cos(kω0 t) + Uk sin(kω0 t)

u(t) = U0 +

(2.2)

k=1
C
S
C
where U0 , Uk and Uk are all real coefﬁcients. It is well-known that Uk = 2 (Uk ), S
Uk = −2 (Uk ), and U0 = U0 . Typically in ﬁnding the steady-state response of a gen-

23

2.1. BACKGROUND ON FOURIER TRANSFORM FOR SINGLE-TONE SIGNALS

eral circuit, one always assumes that all independent sources in the circuit exhibit a certain periodic waveform with a certain period, say T , and the main objective is then focused on computing the time-domain circuit response at steady-state, i.e., as t → ∞, where by “response” here we mean node voltages, branch currents, capacitor charges, inductor ﬂuxes or a combination thereof. Given that the excitation sources are periodic at steady-state, it is intuitive to expect that the response will also be periodic, with the same period, T . This basic knowledge about the general nature of the response at steady-state represents the main key idea in the HB technique and in other steady-state analysis approaches. This chapter starts ﬁrst by offering a brief background material in Section 2.1 on the Fourier transformation techniques used to convert a periodic signal, such as u(t), from the time-domain to its Fourier domain representation, and vice-versa. The notation and mathematical operations laid out in this section are central to the derivation of the HB approach in the following sections.

2.1

Background on Fourier Transform for Single-Tone Signals

Let x(t) ba an arbitrary periodic signal with period T having a Fourier series representation similar to (2.2). To make the process of computing its Fourier series representation tractable, we will have to truncate the inﬁnite series to a limited number of terms, say K terms, and rewrite (2.2) as follows

K
C
S
Xk cos(kω0 t) + Xk sin(kω0 t)

x(t) = X0 +

(2.3)

k=1

The above truncation implicitly assumes that the effect of coefﬁcients corresponding to terms of order higher than K are negligibly small, and would therefore add no further accuracy should x(t) be reconstructed using the Fourier series expansion. The value of integer K deemed sufﬁcient for a good approximation usually depends on the shape of the waveform x(t). For example, in the extreme case where x(t) is a smooth cosinusoid, we C

S
C
S
would have K = 1 with X 0 = Xm = Xm = 0 for m > 1, while, X1 = 1 and X1 = 0. On the other hand, if x(t) is a square-wave with sharp bends and turning points, K will range from few hundreds to few thousands of terms.

Assume that x(t) is sampled at equally spaced H time points, where H = 2K + 1, at the time instants t0 , t1 , · · · , t2K in the interval [0, T ] as shown in Figure 2.1, with tn = T
n H , n = 0, · · · , 2K. Generally, there are two main techniques that are typically used to compute the Fourier Coefﬁcients for a periodic waveform given its time-domain samples (or vice versa, obtain the time-domain samples from a given Fourier series representation). We start with most basic one of them, the Discrete Fourier Transform , DFT (or its inverse known as Inverse Fourier Transform, the IDFT) in the following subsection. 24

CHAPTER 2. HARMONIC BALANCE FOR CIRCUITS WITH SINGLE-TONE EXCITATION

Figure 2.1: A periodic waveform sampled at 2K + 1 time points. 2.1.1

Discrete Fourier Transformation

Writing the Fourier series at an arbitrary time instant ti gives us, K
C
S
Xk cos(kω0 ti ) + Xk sin(kωo ti )

x(ti ) = X0 +

(2.4)

k=1

Repeating the above equation for all time points enables us to write the following matrix equations,

x0

 X 
0
1 cos(ω0 t0 ) sin(ω0 t0 ) · · · cos(Kω0 t0 ) sin(Kω0 t0 )  x1 
C
 .   1 cos(ω0 t1 ) sin(ω0 t1 ) · · · cos(Kω0 t1 ) sin(Kω0 t1 )   X1   .  
  XS 
.
.
.
.
..
 .   .
.
.
.
.
.
 1 
.
.
.
.
.

= .
 . 

 
 . 
.

 C 
 .  
 X 
K
 . 
.
S
1 cos(ω0 t2K ) sin(ω0 t2K ) · · · cos(Kω0 t2K ) sin(Kω0 t2K ) XK
x2K
Γ−1

(2.5)
where xi = x(ti ). Eq. (2.5) can be represented more compactly as follows, x = Γ−1 X

(2.6)
T

C
S
C
S
where x = [x0 , · · · , x2K ]T and X = X0 , X1 , X1 , · · · , XK , XK , and Γ−1 is the matrix shown in Eq. (2.5).
A closer inspection of the elements in the matrix Γ−1 reveals that they are independent of the fundamental angular frequency ω0 as can be seen from the argument of any co- or

25

2.1. BACKGROUND ON FOURIER TRANSFORM FOR SINGLE-TONE SIGNALS

sinusoidal function,
kω0 ti = k

T

i

T
H

= ki

,
H

(1 ≤ k ≤ K, 0 ≤ i ≤ H − 1)

(2.7)

Thus the only condition to make Γ−1 independent of the fundamental angular frequency is to have the time samples equally spaced in the period [0, T ]. In that sense, we can see that the matrix Γ−1 acts as an operator that converts a vector of 2K + 1 Fourier-domain coefﬁcients to its time-domain representation. The essence of the operation that needs to be performed is a matrix-vector multiplication and results in H = 2K + 1 time-domain samples of x(t) at equally spaced time points. The matrix Γ−1 is therefore termed the Inverse Discrete Fourier Transformation (IDFT) operator.

Let Γ denote the matrix inverse of Γ−1 , assuming it is nonsingular. Then Eq. (2.6) can be written as
X = Γx
(2.8)
Hence, the matrix Γ functions as an operator that, when presented with a vector of timedomain samples of a periodic waveform, produces the corresponding coefﬁcients of its harmonic content. Based on this notion, Γ is termed the Discrete Fourier Transformation (DFT) operator.

The IDFT operator matrix Γ−1 has another nice property; its columns are orthogonal and have a norm equal to H/2, except for the ﬁrst column which has a norm equal to H. Thus if γ m and γ n are two columns in Γ−1 , then their inner-product is given by 

m=n
 0

H
γT γn =
(2.9)
m=n=1
m
 2

H m=n=1
This property makes the matrix product Γ−1

T

Γ−1 a diagonal matrix,

Γ−1

T

Γ−1

H
0 ···
 0 H/2 0
= .
..
 .
. ···
.
0 ···

···

0
0
0
H/2

(2.10)

Post-multiplying both sides of (2.10) by the DFT operator matrix Γ shows that it is possible to write the rows of Γ in terms of the columns of Γ−1 . If the row vectors of Γ are denoted 26

CHAPTER 2. HARMONIC BALANCE FOR CIRCUITS WITH SINGLE-TONE EXCITATION

by ρ1 , ρ2 , · · · , ρH , then
1 T
γ
H m

m=1

ρm =

2 T
γ
H m

m = 2, 3, · · · , H

(2.11)

That shows that there is no need to perform the cumbersome matrix inversion to obtain Γ.

2.1.2

Fast Fourier Transformation

As can be seen from Eq.(2.5) the matrices Γ−1 and Γ are full matrices with dimensions H × H, where H is the number of time-domain sampling points. It is well-known that multiplying an H × H matrix by an H × 1 vector requires performing H 2 multiplications and H 2 additions. This fact is expressed by stating that the computational complexity associated with a DFT (or IDFT) is of the order of H 2 , or O (H 2 ). The second order growth of the computational complexity with the number of time-domain samples can represent a signiﬁcant burden on the computational resources for a large number of time samples, especially if the DFT (IDFT) operation has to be performed repeatedly, such as in the course of HB analysis as will be shown in this chapter.

Nevertheless, the very structures of the matrices Γ−1 and Γ offer a special advantage that can be used to reduce the computational complexity. In fact, the characterization of this structure led to an efﬁcient implementation of the DFT (IDFT) known as Fast Fourier Transform or FFT (or IFFT for the inverse implementation). A complete description of the FFT (IFFT) is beyond the scope of this text. However, it would be sufﬁcient to state some basic facts about the FFT (or IFFT). The most important fact is related to the computational complexity which is equivalent to O (H log2 H) if the number of time samples, H, is a power of 2. Other implementations can also be performed in O (H log H) for all H. The reader is referred to the link in  for a C library that implements the state-of-the-art in FFT. In fact, this is the version upon which Matlab functions fft and ifft are built. For the purpose of this course, we use the DFT (IDFT) for its convenience in the mathematical description of the HB analysis, where it will be implicitly understood that the FFT (IFFT) can be used in the actual implementation.

It is worth noting that the above discussion implicitly assumed that the number of timedomain samples is always odd (H = 2K + 1). Yet, this should by no means be always the case, because this number can be made even (or order of 2) by including the cosinusoidal coefﬁcient of the next high order while neglecting the associated sinusoidal coefﬁcient. In this case the Fourier expansion takes the following unbalanced form K

C
S
C
Xk cos(kω0 t) + Xk sin(kω0 t) + XK+1 cos((K + 1)ω0 t)

x(t) = X0 +
k=1

27

(2.12)

2.1. BACKGROUND ON FOURIER TRANSFORM FOR SINGLE-TONE SIGNALS

which yields an even number of time-domain samples (2K + 2). 2.1.3

Numerical Example

We consider here a simple numerical example to illustrate the usage of FFT and DFT techniques. For that purpose, we consider a simple sinusoidal signal given by, x(t) = sin(10t)

(2.13)

The Fourier coefﬁcients of this signal are trivial and can be written as follows, S
C
C
S
X1 = 1, X1 = 0, Xm = Xm = X0 = 0

for m > 1

(2.14)

Another illustrative example can be obtained by considering a square wave with 50% duty cycle. Program 1 is a Matlab script that generates a square wave and computes its Fourier series coefﬁcients using fft function. Also note that because Matlab fft function uses the complex form of the Fourier transformation (i.e., (2.1)), some additional steps had to be done to obtain the corresponding real form.

Program 1 A Matlab script that generates a square-wave and compute its Fourier series. clear all, close all
nTimePoints = 512;
T = 10;
Time = linspace(0,T-T/nTimePoints,nTimePoints);
x = 2*((Time then
comment: Error is higher than an acceptable threshold.
;
Convergence = false
else
comment: Error is lower than an acceptable threshold.
;
Convergence = true;
x(·) ← x(0)
j ← 0;
comment: Repeat until converging to a solution
while Convergence = false do
comment: Compute the Jacobian matrix and the nonlinear function vector at x = x(j)
;
J ← ∂Ψ @ x = x(j) ;
∂x
∆(j) ← Ψ(x) @ x = x(j) ;
comment: Compute the correction step for x(j)
;
δx ← J −1 ∆(j) ;
comment: Obtain the next trial vector.
x(j+1) ← xj − δ x ;
comment: Compute the error at the next trial vector.
∆(j) ← Ψ(x) @ x = x(j) ;
if ∆(j) < then
comment: Error is higher than an acceptable threshold.
;
Convergence = false
else
comment: Error is lower than an acceptable threshold.
;
Convergence = true;
x(·) ← x(j)
j ← j + 1;
43
return x(·)
Figure 2.8: A pseudocode description of the N-R method.

2.3. HARMONIC BALANCE FOR NONLINEAR CIRCUITS

leave us with no other option but to use numerical computation schemes. We consider this computation next, where we summarize the basic ideas and then give a detailed description ¯ ¯ (j) .
for a step-by-step process that results in computing F X
¯ ¯ (j) is assigning some values for the components
The starting point in computing F X
¯ (j)
of the j-th trial vector, X . Those values are typically taken from the values of the previous (j − 1) trial vector, or they could be based on some educated guess, if we happened to be at the very ﬁrst trial vector, (j = 0). More on choosing the initial trial vector will be ¯ ¯ (j) comprises the Fourier coefﬁdiscussed later. Next, recall from Eq. (2.48) that F X cients of the nonlinear functions in f (x(t)). Therefore, obtaining those Fourier coefﬁcients requires the knowledge of the periodic time-domain waveform of each one of those components. However, obtaining the time-domain waveform for a given component in f (x(t)) requires computing the waveform of its arguments, which are typically some selected components of x(t). Thus we need to compute the time-domain waveform for each variable in the MNA formulation. Given that we have the Fourier coefﬁcients at the j-th trial vector, we can make use of those Fourier coefﬁcients to obtain the waveforms for each component ¯ ¯ (j) implementing the

in x(t). We can now describe a ﬁve-step process to compute F X above procedural description but in the reverse order.
• STEP 1. Preprocessing of the j-th trial vector. This step is done in preparation for the following steps. In the next step, the steady-state periodic waveform for each variable in the MNA formulation is computed through IFFT. To perform the IFFT, ¯ (j)

we ﬁrst need to collect the harmonics of each variable from X . As shown in Eq. (2.38), and more clearly in Figure 2.9, the harmonics of each variable are scattered in ¯ (j)
the vector X .
The ordering conﬁguration shown in Figure 2.9 is known as Harmonic-major/Nodeminor ordering in an indication to the fact the components are ordered according to their Harmonic index ﬁrst and then according to their MNA index. Figure 2.10 illustrates these two indices.

To collect the harmonics of each variable, we will need to reorder the components in ¯ (j)
X in Node-major/Harmonic-minor, as shown in Figure 2.11. We denote this new ¯ (j)
vector by X node to emphasize the ordering mode. Vector reordering is usually implemented through swapping memory pointers. A simpler ordering procedure, which is more convenient for the sake of mathematical derivation, is based on the multiplication with a suitable permutation matrix. A permutation matrix is an orthonormal matrix with at most one entry equal to (1) in each row, and (0)’s in the rest of row’s ¯ (j)

entries. Assume that the required permutation matrix needed to obtain X node from ¯ (j)
X is denoted by P, then,
¯ (j)
¯ (j)
X node = P X
(2.60)
44

CHAPTER 2. HARMONIC BALANCE FOR CIRCUITS WITH SINGLE-TONE EXCITATION

The vector of the DC components of all MNA variables

The vector of the Co-sinusoidal coefficients of the first harmonic for all MNA variables

The vector of the sinusoidal coefficients of the first harmonic for all MNA variables

The vector of the sinusoidal coefficients of the Mth harmonic for all MNA variables

¯ (j)
Figure 2.9: The initial ordering of the j-th trial vector X .

Harmonic Index

MNA Index

¯ (j)
Figure 2.10: Illustrating the two main indices used to index the components of X . and by the orthonormality of P we have
P T P = identity matrix

(2.61)

• STEP 2. Performing IFFT. In this step the harmonics of every MNA variable is passed through IFFT to obtain the corresponding periodic time-domain waveform. From the discussion of Section 2, this process is equivalent to multiplying each one of ¯ (j)

the boxed subvectors (shown in Figure 2.11) of X node by the IDFT operator Γ−1 . The resulting time-domain waveforms are represented by (2M + 1) sample time points at the time instants t0 , t1 , · · · , t2M . We can collect those time samples for each MNA 45

2.3. HARMONIC BALANCE FOR NONLINEAR CIRCUITS

¯
Figure 2.11: Reordering the components of X

(j)

in Node-major/Harmonic-minor.
(j)

variable in a single vector and denote that vector by xnode I . Figure 2.12 displays the ¯
(j)
structure of xnode . Note that the operation performed in this step can be described ¯
mathematically as follows
(j)
¯ −1 ¯ (j)
xnode = Γ X node
¯
(2.62)
−1

¯
where Γ is an N (2M + 1) × N (2M + 1) block diagonal matrix having the matrix −1
Γ as its diagonal blocks.
• STEP 3. Computing the no...