➠➡
Distortion Analysis of Nonlinear Networks Based on
SPICE-Oriented Harmonic Balance Method
Y.Yamagami†, H.Yabe†, Y.Nishio†, A.Ushida‡†Department of Electrical and Electronic Engineering,
Tokushima University, Tokushima, 770-8506 JAPAN‡Department of Mechanical and Electronic Engineering
Tokushima Bunri University, Kagawa, 769-2193 JAPAN
Abstract
Distortion analysis in the frequency domain is very important for designing the analog integrated circuits. We
propose here a new method based on the SPICE-oriented
harmonic balance method, where the Fourier transformations to the bipolar transistors and/or MOSFETs are carried out using new Fourier transfer circuit model. The circuit is composed of the analog behavior models(ABM) ofSPICE, which can be applied to any kind of nonlinear circuit elements described by the exponential and/or piecewise
linear functions. Furthermore, the determining equation of
our harmonic balance method is schematically described bycoupled resistive DC, Cosine and Sine circuits, so that wecan 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 harmonic balance method.

1. Introduction
The distortion analysis in the frequency domain is very
important for designing the nonlinear analog integrated circuits. The Volterra series methods are widely used for theanalysis of weakly nonlinear systems [1-3]. Although the algorithm is theoretically elegant, it is not so easy to derive theVolterra kernels in the higher order distortions analysis [4].
Furthermore, their nonlinear characteristic in the method
must be described in the form of polynomial, so that they
can be obtained by the Taylor expansion in the vicinity of theDC operating point of a nonlinear elements [2]. This task isalso not so easy especially for the complicated circuit models such as the high frequency Gummel-Poon model of bipolar transistors and/or Shichman-Hodges model of MOSFETs[5]. On the other hand, there have been proposed many algorithms for calculating the exact steady-state waveforms ofnonlinear circuits driven by the periodic and/or multi-tone
signals [6-8]. Unfortunately, they are ineﬃcient for calculating the nonlinear frequency response curves.
In this paper, we present a new harmonic balance method
for calculating the frequency response curves in the distortion analysis of nonlinear integrated circuits. The determining equations can be schematically described by the equivalent DC, Cosine and Sine circuits, where Fourier transfer
circuit model can be eﬃciently used for the Fourier expansion of the nonlinear characteristics. The circuit is composedof the resistive analog behavior models (ABM) of SPICE sothat the frequency response curves can be obtained by the0-7803-8251-X/04/$17.00 ©2004 IEEE

application of curve tracing algorithm [12]. Note that theFourier expansion can be eﬃciently applied to any kind ofcircuit elements such as bipolar transistors, MOSFETs andso on. Remark that the method in the reference [9] canbe only applied to the circuit whose elements are describedby the polynomial functions. The Fourier circuit model isshown in section 2, and the equivalent DC, Cosine and Sinecircuits are shown in section 3. Each circuit topology is thesame as the original one. Of course, the technique can beeasily extended to the higher harmonic distortion analyses.
We show an interesting illustrative example in section 4.
Since our method needs not to derive both the circuit equation and the determining equations, it is really user-friendly
algorithm for calculating nonlinear response curves.

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]. Let us discuss how
to calculate the frequency response curves. In this case, weﬁrst need the Fourier expansion to the nonlinear elements
driven by the periodic input.
Let us assume the input and the output waveforms asfollows;
x(t) = A0 +M
k=1y(t) = B0 +M

k=1

⎫⎪
⎪
(A2k−1 cos kωt + A2k sin kωt) ⎪⎪⎬⎪⎪(B2k−1 cos kωt + B2k sin kωt) ⎪⎪⎭(1)

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

B0 = g0 (A0 , A1 , · · · , A2M ) ⎪⎬B1 = g1 (A0 , A1 , · · · , A2M )...............................
⎪⎭
B2M = g2M (A0 , A1 , · · · , A2M )(2)Note that the Fourier coeﬃcients Bk can be given by theexplicit functions of (A0 , A1 , . . . , A2M ), only when the nonlinear function is described by the polynomial function [11].
Therefore, we propose here a new technique using Fouriertransfer circuit model as shown in Fig.1 and 2.

IV - 633

ISCAS 2004

➡

➡
A0
A1
A2AK-1
AKTo investigate the numerical accuracy, we ﬁrst calculate
a modiﬁed Bessel function as follows;

B0
B1
B2FourierTransformationCircuitIN (x) =BK-1
BK

Fig.1 Block diagram realizing the relation (2).

Each Fourier coeﬃcient for the function f (x) can be numerically obtained in the following formula:
B0 =

12π

 2πf (x)dt0

 2π1
B2k−1 =π 0
k = 1, 2, . . . , Mf (x) cos kωtdt, B2k =1π

 2π0⎫⎪⎪⎪⎬f (x) sin kωtdt ⎪⎪⎪⎭(3)

Let us apply the trapezoidal integration formula to (3) asfollows; bh(f0 + fn ) + h(f1 + f2 + · · · + fn−1 )23.0

Rf(A,θ 2 )cosNθ 2f(A,θ 2 )sinNθ 2(x10 3 )
N=0
N=1
N=2
N=3

Modified Bessel FunctionIN

1.00.5f(A,θ 1)cosNθ 10.5f(A,θ1 )sinNθ 1(5)

−π

2.0

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 operation, we simply assume that the input is x = A sin θ. TheFourier transfer circuit model shown by Fig.2 is composedby ABM blocks, whose operations such as multiplications,exponential, trigonometric functions and so on can be doneby the functions written by the Fortran languages. On theother hand, the integration interval [0, 2π] is divided by nsections using a number of n resistors, so that each input
node voltage of the ABM block is given by θk = 2πk/n at
kth node. Thus, the resultant current sources are given byf (A, θk ) cos N θk , f (A, θk ) sin N θk , k = 1, 2, . . . , n. Summing
them, the outputs correspond to the coeﬃcient of cos N θand sin N θ (13). In this case, if the input contains the higher
harmonic components, we need to set the number of m terminals (A1 , A2 , . . . , Am ) instead of the single terminal A inFig.2.

Rex cos θ cos N θdθ

2.5

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

(4)f (x)dx =12π0.50.00246810xFig.3 Fourier transformation for modiﬁed Bessel function.

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

ID =

KW
VDS(VGS − VT ) −VDS (1 + λVDS )L2(6.1)

2. Saturation region: (VGS > VT , VDS > VGS − VT )ID =KW(VGS − VT )2 (1 + λVDS )L(6.2)The result of Fourier expansions for the input VGS cos ωt isshown in Fig.4.
30

Σ/nCoefficient of CosNθ

ID[µ]
202π

Fourier Expansion of MOSFET

N=1N=2
N=3
N=4
N=5f(A,θn-1)cosNθ n-1Rf(A,θn-1)sinNθ n-10.5f(A,θn )cosNθ nR

Σ/nCoefficient of SinNθ

100.5f(A,θn )sinNθ nA0N

0.0

Fig.2 Fourier transfer circuit model.

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

IV - 634

2.0

4.0

6.0

8.0

VGS[V]

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

➡
For the voltage vC = VC0 + VC1 cos ωt + VC2 sin ωt, wehaveThus, the Fourier transfer circuit model shown in Fig.2 canbe eﬃciently applied to any kind of circuit elements contained in analog integrated circuits.

iC (t)  −ωQC1 sin ωt + ωQC2 cos ωt

3. Frequency response curve ofnonlinear circuits

where QC1 and QC2 are function of {VL0 , VL1 , VL2 },and they are also estimated by the Fourier transfer
circuit shown in Fig.2.
For the linear capacitor C, the Fourier coeﬃcients aregiven by the linear voltage-controlled current sources
ωCVV 2 in the Cosine-circuit and −ωCVC1 in the Sinecircuit, respectively [9].

Nowadays, the frequency response curves of the fundamental
frequency (H1 ), the second harmonic distortion(HD2 ) andthe third harmonic distortion(HD3 ) are usually calculatedby the use of the Volterra series method [3], where the nonlinear characteristics must be described by the polynomialfunctions. On the other hand, our distortion analysis using
the Fourier transfer circuit model can be eﬃciently applied
to the circuits containing any kind of the elements describedby 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 → Rn+mc. Resistive elements: Many of integrated circuit
elements are described by the nonlinear voltagecontrolled current sources as follows;
iG = îG (vG )

(7)

(8)

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

In generally, the derivation of (9) from the circuit equation,and the analysis of the determining equations are troublesome. We propose a SPICE-oriented algorithm for getting
the DC-circuit, Cosine-circuit and Sine-circuit corresponding to the above the relations (9).
a. Inductive elements: Assume that the nonlinear inductor is described by current-controlled characteristic as follows;(10.1)
φL = φ̂L (iL )iG (t)  IG0 + IG1 cos ωt + IG2 sin ωt(10.2)

where ΦL1 and ΦL2 are function of {IL0 , IL1 , IL2 }.
They are calculated by the Fourier transfer circuit
shown by Fig.2.
For the linear inductor L, the coeﬃcient are simply
given by the linear current-controlled voltage sources
ωLIL2 in the Cosine-circuit and −ωLIL1 in the Sinecircuit, respectively [9].
b. Capacitive elements: High frequency bipolar transistor and MOFET models [5] usually contain nonlinear depletion and diﬀusion capacitors which aredescribed by the voltage-controlled characteristics asfollow;(11.1)
qC = q̂C (vC )(12.2)

where IG0 , IG1 and IG2 are described by the functionof {VG0 , VG1 , VG2 }. They are estimated by the use of
the Fourier transfer circuit model shown by Fig.2. For
linear resistors, we need not any transformation [9].
Note that the inductive elements in the DC-circuit are removed by the shorted-circuits, and the capacitive elements
by the opened-circuits.
Thus, the circuit topologies corresponding to the DC, Cosine and Sine circuits and their higher harmonic components
are equal to the original circuit, and they are coupled in eachother with the controlled sources.

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

For the current iL = IL0 + IL1 cos ωt + IL2 sin ωt, wehavevL (t)  −ωΦL1 sin ωt + ΦL2 cos ωt(12.1)

For the voltage vG = VG0 + VG1 cos ωt + VG2 sin ωt, wehave

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(11.2)

RL = 10[kΩ],αR = 0.3Vs = 800[mV ]33C1 = 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 +3(Vbe,2k−1 cos kωt + Vbe,2k sin kωt)

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

k=1

Applying the harmonic balance method, we obtained theDC, Cosine and Sine circuit as shown in Fig.5 (c),(d) and (e),respectively. The circuit models for calculating the higher
order distortions are obtained in the same manner, whereIV - 635

➡

➠IC0 , IE0 are the DC voltage-controlled current sources, andthe voltage-controlled current sources (IC,k , IE,k , k = 1, 2, 3)
are given by the function of the input frequency ω. These
controlled sources can be calculated by the use of Fouriertransfer circuit model shown in Fig.2.

VD
CRLv outv cbib
RSv inC1ic

VSie

α F iD1BC2vbei D2iD1

α RiD2E(a)(b)

where the nonlinear characteristics must be described by
the polynomial function in the vicinity of the DC operating
point. Thus, the method can be only applied to weakly nonlinear circuits. Furthermore, the higher order kernels of theVolterra series are complicate, so that it is usually restricted
to the analysis of relatively low order distortion analysis.
In this paper, we have proposed an eﬃcient SPICEoriented distortion analysis method based on the harmonicbalance method, where we uses Fourier transfer circuit
model. The model can be eﬃciently applied to any kindof nonlinear elements described by such as exponential andpiecewise continuous functions. 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 and userfriendly. We are going to apply the algorithm to the complicated circuits such as the modulators.

References[1] M.Schetzen, The Volterra and Wiener Theorems ofRL
Nonlinear Systems, John Wiley and Sons, 1978.
RLRL[2]D.D.Weinerand J.E.Spina, Sinusoidal Analysis andI C0I C1IC2Modeling of Weakly Nonlinear Circuits with Application to Nonlinear Interference Eﬀects, Van Nostrand
Reinhold, 1980.
I E2RS
IE0 R SI E1 RS[3] P.Wambacq and W.Sansen, Distortion Analysis of AnaVS
log Integrated Circuits, Kluwer Academic Pub., 1998.
v in[4] L.O.Chua and C-Y.Ng, “Frequency-domain analysis ofnonlinear systems: formulation of transfer functions,”(e)(c)
IEE Trans. Electronic Circuits and Systems, Vol.3,(d)no.6, pp.257-269, 1979.
Fig.5 (a) Ampliﬁer circuit, (b) High frequency Ebers-Moll[5]R.L.Geiger,
P.E.Allen and N.R.Strader, VLSI: Designmodel, (c) DC-circuit, (d) Cosine-circuit, (e) Sine-circuit
Techniques for Analog and Digital Circuits, McGrawContinuously changing the input frequency with the DC
Hill, 1990.
analysis of SPICE, we can obtain the frequency response[6] T.J.Aprille and T.N.Trick, “Steady-state analysis ofcurves for the distortion analysis as shown in Fig.6, which
nonlinear circuits with periodic input,” Proc. IEEE,is almost same as those of reference [3].
vol.60, pp.108-114, 1972.
[7] K.S.Kundert,
J.K.WhiteandA.Sangiovanni1.0e+1
Vincentelli, Steady-State Methods for Simulating Analog and Microwave Circuits, Kluwer Academic Pub.,v m, in =10mV
1.0e+01990.
Fundamental
HD2[8]
A.Ushida,
T.Adachi and L.O.Chua, “Steady-state anal1.0e-1
HD3
ysis of Nonlinear Circuits based on hybrid methods,”
IEEE Trans. on Circuits and Systems-I: Fundamen1.0e-2
v outtal Theory and Applications, vol.39, no.9, pp.649-661,
1992.
1.0e-3[9] A.Ushida, Y.Yamagami and Y.Nishio, “Frequency responses of nonlinear networks using curve tracing algo1.0e-4
rithm,”ISCAS 2002, vol.I, pp.641-644, 2002.
1.0e-5[10] K.K.Clarke and D.T.Hess, Communication Circuits:
Analysis and Design, Addison-Wesley Pub. Co., 1971.
1.0e-6
1.0e+31.0e+41.0e+51.0e+6
1.0e+71.0e+81.0e+9 [11] S.Wolfram, Mathematica, Addison-Wesley Pub. Co.,
(rad/sec)
Omega
1988
Fig.6 Distortion analysis.
[12] A.Ushida, Y.Yamagami, Y.Nishio, I.Kinouchi and
Y.Inoue, ”An eﬃcient algorithm for ﬁnding multiple
5. Conclusions and remarks
DC solutions based on the SPICE-oriented Newton homotopy method,” IEEE Trans. on CAD. vol.21, no.3,
The distortion analysis in frequency domain is very imporpp.337-348, 2002.
tant for designing the high frequency analog integrated circuits. It is usually carried out by the Volterra series method,IV - 636