Spice-Oriented Distortion Analysis of Nonlinear Networks
A.Ushida†, Y.Yamagami‡, Y.Nishio‡†Department of Mechanical and Electronic Engineering
Tokushima Bunri University, Kagawa, 769-2193 JAPAN‡Department of Electrical and Electronic Engineering,
Tokushima University, Tokushima, 770-8506 JAPAN
Abstract
Frequency-domain distortion analysis is very importantfor designing the analog integrated circuits. We propose
here a new method based on the Spice-oriented harmonicbalance method, where the Fourier transformations of theresistive circuits driven by a periodic input are carried out
by the Fourier transfer circuit model. The circuit is composed of the analog behavior models(ABM) of Spice, which
can be applied to any kind of nonlinear circuits containing
the exponential and/or piecewise linear characteristics. Furthermore, the determining equation of our harmonic balancemethod is schematically described by coupled resistive DC,
Cosine and Sine circuits, so that we can easily obtain the frequency response curves with the DC analysis of Spice. Thus,
our approach is quite user-friendly, because we need not derive the circuit equation and the determining equations in
our method.

1. Introduction
The frequency-domain distortion analysis is very important for designing the nonlinear analog integrated circuits.
The Volterra series methods are widely used for the analysis of weakly nonlinear systems [1-3]. Although the algorithm is theoretically elegant, it is not so easy to derive
the Volterra kernels for the higher order distortion analysis [4]. Furthermore, their nonlinear characteristics in themethod should be described in the form of power series,
so that we need to apply the Taylor expansion techniquein the vicinity of the DC operating point of the nonlinear
elements [2]. This task is not so easy especially for the complicated circuit models such as the high frequency GummelPoon model of bipolar transistors and/or Shichman-Hodges
model of MOSFETs [5]. On the other hand, there have been
proposed many time-domain algorithms for calculating theexact steady-state waveforms of nonlinear circuits driven by
the periodic and/or multi-tone signals [6-8]. Unfortunately,
they are rather time-consuming for getting the nonlinear frequency response curves.
In this paper, we propose a new harmonic balance method
for calculating the frequency response curves of nonlinear integrated circuits. The determining equations can be
schematically described by the equivalent DC, Cosine andSine circuits, where Fourier transfer circuit model are efficiently used for the Fourier expansion of the nonlinear resistive circuits driven by a periodic input. The circuit isrealized by the resistive analog behavior models (ABM) of

Spice, so that the frequency response curves can be efficiently calculated by the application of curve tracing algorithm [12]. Note that the Fourier expansion can be applied
to any kind of nonlinear circuit elements such as bipolartransistors, MOSFETs and so on. Remark that the methodin the reference [9] can be only applied to the circuit whose
elements are described by the power series. The Fourier circuit model is shown in section 2, and the equivalent DC,
Cosine and Sine circuits are shown in section 3. Note thateach of the equivalent circuit has the same topology as theoriginal one. We show an interesting illustrative example insection 4. Since our method needs not to derive any troublesome circuit equations and the determining equations, it isreally user-friendly algorithm for calculating the frequency
response curves of nonlinear circuits.

2. Fourier transfer circuit model
Analog integrated circuits are usually composed of manykinds of nonlinear elements such as diodes, bipolar transistors and MOSFETs, whose models are described by the
special functions such as the exponential, square-root, piecewise continuous functions and so on [5]. We propose here
an efficient harmonic balance method for calculating the frequency response curves of ICs. In this case, we first need todevelop a simulator carrying out the Fourier expansion of
the nonlinear resistive circuits driven by the periodic input.
Let us assume the input and output waveforms as follows;
v(t) = V0 +

M
X




(V2k−1 cos kωt + V2k sin kωt) 



k=1

i(t) = I0 +

M
X



(I2k−1 cos kωt + I2k sin kωt) 



(1)

k=1

where M denotes the highest harmonic component to takeaccount in the analysis. Thus, the output Fourier coefficientsneed to be described by the functions of input amplitudes
as follows;

I0 = f0 (V0 , V1 , · · · , V2M ) 
I1 = f1 (V0 , V1 , · · · , V2M )...............................

I2M = f2M (V0 , V1 , · · · , V2M )(2)Note that the Fourier coefficients Ik can be given by the
explicit functions of (V0 , V1 , . . . , V2M ) only if the nonlinear function is described by the power series function [11].

Therefore, we need to introduce something a new technique
for the Fourier transformation of general resistive circuits
driven by the periodic input. The circuit model is shown in
Fig.1, where each Fourier coefficient for the function i = f (v)can be numerically calculated in the following formula:
I0 =

1
2π

Z 2π







f (v)dt
0

Z 2π

1
1
f (v) cos kωtdt, I2k =π 0
πk = 1, 2, . . . , M

Z 2π

f (v) sin kωtdt I2k−1 =




0

(3)

Let us apply the trapezoidal integration formula to (3) as
follows;

Z b
f (v)dt =ah(f0 + fn ) + h(f1 + f2 + · · · + fn−1 )
2

(4)

M
X

1
IN (x) =2π

The simulation results with h = 2π/20 is shown in Fig.2.
The value I1 (10) = 2761 at N = 1, x = 10 is exactly equal to
the result from the Table of Bessel function [10]. Note that
this kind of Fourier expansion for the exponential function
is very important to the analysis of the circuit containing
diodes and bipolar transistors [5].
(x10 3 )
3.0
N=0
N=1
N=2
N=3

2.5

2.0

1.5

1.0

0.0
0

2

0.5f(V,θ2K)cos Nθ 2K

6Σ/2Κ
Coefficient of CosNθ2π

Next, we apply it to the Fourier expansion of MOSFET,
whose characteristic in the Spice model is described by apiecewise continuous function [5] as follows:
1. Linear region: (VGS > VT , VGS − VT ≥ VDS > 0)

h

ID =

i

KW
VDS
(VGS − VT ) −VDS (1 + λVDS )
L
2
(7.1)

2. Saturation region: (VGS > VT , VDS > VGS − VT )KW(VGS − VT )2 (1 + λVDS )L

R

0.5f(V,θ0)cosNθ 0
0.5f(V,θ0) sin Nθ 0

(7.2)

30
ID[µ]
20

f(V,θ1 )cosNθ 1f(V,θ1 ) sin Nθ 1

10

Fig.2 Fourier transformation for modified Bessel function.

f(V,θ2K-1 )sin Nθ 2K-1

R

8

The result of Fourier expansions for the input VGS cos ωt isshown in Fig.3.

0.5 f(V,θ2K) sin Nθ 2K
f(V,θ2K-1 )cosNθ 2K-1

4

x

ID =

V0

(6)

−π

(5)

We set θ = ωt. The Fourier transfer circuit model for
calculating the N th higher harmonic components is shownby Fig.1. A number of 2K + 1 blocks calculate the N thcomponents using (3) which is composed of the nonlinear resistive circuits and ABMs of Spice. On the other
hand, the integration interval [0, 2π] is divided by 2K sections using 2K resistors, so that each input node voltageof the ABM blocks is given by θk = 2πk/2K at the kth
node. Thus, the resultant current sources are given byf (V, θk ) cos N θk , f (V, θk ) sin N θk , k = 1, 2, . . . , 2K. Summing them, their outputs correspond to the coefficients ofcos N θ and sin N θ.

R

ex cos θ cos N θdθ

0.5(V2k−1 cos kωt + V2k sin kωt)

k=1

R

Z π

Modified Bessel FunctionIN

where the step-size of the integration is h = (a−b)/n. Then,the truncation error is given by f (2) h2 /12n. Using this formula, we will realize the equivalent circuit model satisfying(3) with ABMs of Spice. To understand the circuit model,
we assume that the input isv(t) = V0 +To investigate the numerical accuracy, we first calculate
a modified Bessel function as follows;

Fourier Expansion of MOSFET

N=1
N=2
N=3
N=4
N=5

Σ/2Κ
Coefficient of SinNθ

V2M

Fig.1 Fourier transfer circuit model.

The fact that the circuit model is only composed of the resistive circuits is important for the calculation of the frequency
response curves with the curve tracing algorithm [12].

10

0
0.0

2.0

4.0

6.0

8.0

VGS[V]

Fig.3 Fourier transformation for MOSFET.
VDS = 3[V ], K = 3.87[µA], λ = 0.01605,
W = 2[µm], L = 2[µ], VT = 0.827[V ]10.0

Thus, the Fourier transfer circuit model shown in Fig.1 canbe efficiently applied to any kind of circuit elements contained in analog integrated circuits.

3. Frequency response curves ofnonlinear circuits
Nowadays, the frequency response curves of the fundamental frequency (H1 ), the second harmonic distortion(HD2 ), the third harmonic distortion(HD3 ) and so on areusually calculated by the use of the Volterra series method[3], where the nonlinear characteristics must be described by
the power series functions. On the other hand, our distortion
analysis using the Fourier transfer circuit model can be efficiently applied to the circuits containing any type of the elements described by exponential, piecewise linear functions
and others. To understand our harmonic balance method,
we consider the following circuit equation;f (v̇, v, w, ωt) = 0, f : R2n+m 7→ Rn+m

(8)

Although the steady-state waveforms may contain manyhigher harmonic components, for simplicity, we consider
only the DC and fundamental frequency components as follows;
v(t) = V0 + V1 cos ωt + V2 sin ωtw(t) = W0 + W1 cos ωt + W2 sin ωt

¾
(9)

)
(10)

In generally, although the equation (10) can be obtainedby the application of the harmonic balance method to thecircuit equation, it is really troublesome task. We propose a
Spice-oriented algorithm for getting the relations (10) fromthe DC-circuit, Cosine-circuit and Sine-circuit, directly.
a. Inductive elements: Assume that the nonlinear inductor is described by current-controlled characteristic as follows;
φL = φ̂L (iL )(11.1)For the inductor current,
M
X

iL = I0,L +(I2k−1,L cos kωt + I2k,L sin kωt), (11.2)
k=1

we have






!
Ã


M

X



(Φ2k−1,L cos kωt + Φ2k,L sin kωt)
= Φ0,L +



k=1
ÃM
!
X


×
kω(−I2k−1,L sin kωt + I2k,L sin kωt)




k=1


M

X



(V2k−1,L cos kωt + V2k,L sin kωt)
' V0,L +vL (t) =∂ φ̂L diL∂iL dtk=1

b. Capacitive elements: High frequency bipolar transistor
and MOFET models [5] usually contain nonlinear depletion and diffusion capacitors which are described by thevoltage-controlled characteristics as follow;
qC = q̂C (vC )(12.1)For the capacitor voltage,
M
X(V2k−1,L cos kωt + V2k,C sin kωt),vC = V0,C +k=1(12.2)

we have






Ã
!

M

X



= Q0,C +(Q2k−1,C cos kωt + Q2k,C sin kωt)




k=1
ÃM
!
X

×
kω(−V2k−1,C sin kωt + V2k,C sin kωt)




k=1


M

X



' I0,C +(I2k−1,C cos kωt + I2k,C sin kωt)iC (t) =∂ q̂C dvC∂vC dtk=1

Substituting (9) into (8) and applying the harmonic balancemethod, we have the following determining equation;
F0 (V0 , V1 , V2 , W0 , W1 , W2 , ω) = 0 · · · DC
Fc (V0 , V1 , V2 , W0 , W1 , W2 , ω) = 0 · · · cosωt
Fs (V0 , V1 , V2 , W0 , W1 , W2 , ω) = 0 · · · sinωt

For the linear inductor L, the coefficient are simply givenby the linear current-controlled voltage sources kωLI2k,L
in the Cosine-circuit and −kωLI2k−1,L in the Sine-circuit,respectively [9].

(11.3)
where Φ0,L , Φ2k−1,L and Φ2k,L are function of{I0,L , I2k−1,L , I2k,L }. They are calculated by the Fouriertransfer circuit shown by Fig.1.

(12.3)
where Q0,C , Q2k−1,C and Q2k,C are function of{V0,C , V2k−1,C , V2k,C }, and they are also estimated by the
Fourier transfer circuit shown in Fig.1.
For the linear capacitor C, the Fourier coefficients are
given by the linear voltage-controlled current sourceskωCV2k,C in the Cosine-circuit and −kωCV2k−1,C in the
Sine-circuit, respectively [9].

Note that the inductive elements in the DC-circuit areremoved by the shorted-circuits, and the capacitive elements
by the opened-circuits.
Thus, all the circuit topologies corresponding to the DC,
Cosine-circuits and Sine-circuits for the higher harmoniccomponents are equal to the original circuit, but they arecoupled in each other with the controlled sources.

4. An illustrative example
Now, consider distortion analysis of a simple high frequency amplifier [3] shown by Fig.4(a). The Ebers-Mollmodel of transistor [5] is shown by Fig.4(b), where the diodemodels are given by
iD1 = 10−14 {exp(40vbe − 1)}[A],iD2 = 10−14 {exp(40vbc − 1)}[A]andRs = 10[kΩ],αF = 0.99,

RL = 10[kΩ],αR = 0.3

Vs = 800[mV ]33
C1 = 10+10vbc[pC], C2 = 10+10vbe[pC], vin = Vm,in cos ωt

For calculating the frequency response of the fundamental, the second order distortion (HD2 ) and the third orderdistortion(HD3 ), we assumed the waveforms as follows;
vbe = Vbe,0 +3X(Vbe,2k−1 cos kωt + Vbe,2k sin kωt)

k=1
3vbc = Vbc,0 +X(Vbc,2k−1 cos kωt + Vbc,2k sin kωt)

k=1

Applying the harmonic balance method, we obtained theDC, Cosine-circuit and Sine-circuit as shown in Fig.4 (c),(d)and (e), respectively. The circuit models for calculating thehigher order distortions can be obtained in the same manner. Note that IC0 , IE0 in the figure are the DC i voltagecontrolled current sources, and (IC,k , IE,k , k = 1, 2, 3) are
given by the functions of the input frequency ω. These
controlled sources can be calculated by the use of Fouriertransfer circuit model shown in Fig.1.
VD
C

RL
v out
v cb
ib

RL
i D2

C1

ic

B
C2

RS

vbe

v in

I C0

α F iD1

iD1

ie

α Ri

D2

RS

IE0

VS

VS

References

E

DC

(b)

(a)
RL

RL

IC1C

IC1S

RS
I E1C

Vm,in

cos ω tto the weakly nonlinear circuits. Furthermore, the higher
order kernels of the Volterra series are complicate, so that
it is usually restricted to the analysis of relatively low orderdistortion analysis [3].
In this paper, we proposed an efficient Spice-oriented distortion analysis based on the harmonic balance method,where we proposed the Fourier transfer circuit model. The
model can be efficiently applied to any kind of nonlinear elements described by such as exponential and piecewise continuous functions and/or resistive circuits. We also proposed
the DC, Cosine and Sine circuits corresponding to the determining equation of the harmonic balance method, whose
circuit topologies are equal to the original circuit. Thus,
the algorithm of our distortion analysis is quite simple anduser-friendly. We are going to apply the algorithm to thecomplicated circuits such as the modulators.

RS

(d)

I E1S

RS

I C2C

I C2S

I E2C R S

I E2S

cos 2ω tsin ω t

(c)

(e)

RL

RL

RL

[1] M.Schetzen, The Volterra and Wiener Theorems of
Nonlinear Systems, John Wiley and Sons, 1978.
RL

I C3C

I E3C

RS

sin2ω tcos 3ω t

(f)

I C3S

I E3S

RS

sin3 ω t[2] D.D.Weiner and J.E.Spina, Sinusoidal Analysis andModeling of Weakly Nonlinear Circuits with Application to Nonlinear Interference Effects, Van Nostrand
Reinhold, 1980.
[3] P.Wambacq and W.Sansen, Distortion Analysis of Analog Integrated Circuits, Kluwer Academic Pub., 1998.

Fig.4 (a) Amplifier circuit, (b) High frequency Ebers-Mollmodel, (c) DC-circuit, (d) Cosine-circuit, (e) Sine-circuit[4] L.O.Chua and C-Y.Ng, “Frequency-domain analysis ofnonlinear systems: formulation of transfer functions,”
IEE Trans. Electronic Circuits and Systems, Vol.3,
no.6, pp.257-269, 1979.

Continuously changing the input frequency ω, we can obtainthe frequency response curves for the distortion analyses forthe fundamental, HD2 and HD3 as shown in Fig.5, which isalmost the same as those of reference [3].

[5] R.L.Geiger, P.E.Allen and N.R.Strader, VLSI: Design
Techniques for Analog and Digital Circuits, McGrawHill, 1990.

1.0e+1v in =10mV

1.0e+0Fundamental
HD2
HD3

1.0e-1[6] T.J.Aprille and T.N.Trick, “Steady-state analysis ofnonlinear circuits with periodic input,” Proc. IEEE,
vol.60, pp.108-114, 1972.
[7] K.S.Kundert,
J.K.White
andA.SangiovanniVincentelli, Steady-State Methods for Simulating Analog and Microwave Circuits, Kluwer Academic Pub.,
1990.
[8] A.Ushida, T.Adachi and L.O.Chua, “Steady-state analysis of Nonlinear Circuits based on hybrid methods,”
IEEE Trans. on Circuits and Systems-I: Fundamental Theory and Applications, vol.39, no.9, pp.649-661,
1992.

1.0e-2v out

1.0e-31.0e-4[9] A.Ushida, Y.Yamagami and Y.Nishio, “Frequency responses of nonlinear networks using curve tracing algorithm,”ISCAS 2002, vol.I, pp.641-644, 2002.

1.0e-51.0e-6
1.0e+3(c)

1.0e+41.0e+51.0e+6Omega

1.0e+71.0e+81.0e+9[10](rad/sec)

Fig.5 Distortion analysis.

5. Conclusions and remarks
The frequency domain distortion analysis is very important for designing the high frequency analog integrated circuits. It has been usually carried out by the Volterra series
method, where the nonlinear characteristics should be described by the power series function in the vicinity of the
DC operating point. Thus, the method can be only applied

K.K.Clarke and D.T.Hess, Communication Circuits:
Analysis and Design, Addison-Wesley Pub. Co., 1971.

[11] S.Wolfram, Mathematica, Addison-Wesley Pub. Co.,
1988.
[12] A.Ushida, Y.Yamagami, Y.Nishio, I.Kinouchi and
Y.Inoue, “An efficient algorithm for finding multiple
DC solutions based on the SPICE-oriented Newton homotopy method,”IEEE Trans. on CAD. vol.21, no.3,
pp.337-348, 2002.