IEICE TRANS. FUNDAMENTALS, VOL.E90–A, NO.2 FEBRUARY 2007

406

LETTER

Special Section on Analog Circuit Techniques and Related Topics

Spice-Oriented Frequency-Domain Analysis of Nonlinear
Electronic Circuits
Junji KAWATA†a) , Yousuke TANIGUCHI††b) , Masayoshi ODA††c) , Yoshihiro YAMAGAMI††d) ,Yoshifumi NISHIO††e) , Members, and Akio USHIDA†f) , Fellow

SUMMARY
Distortion analysis of nonlinear circuits is very importantfor designing analog integrated circuits and communication systems. In
this letter, we propose an eﬃcient frequency-domain approach for calculating frequency response curves, which is based on HB (harmonic balance) method combining with ABMs (Analog Behavior Models) of Spice.
Firstly, nonlinear devices such as bipolar transistors and MOSFETs aretransformed into the HB device modules executing the Fourier transformations. Using these modules, the determining equation of the HB method isformed by the equivalent sine-cosine circuit in the schematic form or netlist. It consists of the coupled resistive circuits, so that it can be eﬃciently
solved by the DC analysis of Spice. In our algorithm, we need not to derive
any troublesome circuit equations, and any kinds of the transformations.
key words: frequency-domain analysis, Volterra series method, spiceoriented harmonic balance method, HB device modules

1.

Introduction

Frequency-domain analysis of nonlinear electronic circuits
is very important for designing integrated circuits and communication systems. The Volterra series methods are widelyused for the analysis [1]–[3], especially in Europe, because
the methods give the solutions in analytical forms. Although
their algorithms are theoretically elegant, it is not so easy toderive the higher order Volterra kernels for the large scale
systems containing many nonlinear elements. They can be
only applied to the weakly nonlinear circuits, and the solutions may be erroneous for strong nonlinear circuits, because the methods are based on the bilinear theorem. Furthermore, the nonlinear characteristics should be approximated by the polynomial functions, where Taylor expansions are used to the approximations at the vicinities of DC
operating points [3]. However, these tasks are not easy forthe complicated devices such as Gummel-Poon model ofbipolar transistors, higher frequency models of MOSFETs
and the piecewise linear models [4], [5].
Manuscript received September 8, 2006.
Final manuscript received October 10, 2006.
†
The authors are with the Department of Mechanical and Electronic Engineering, Faculty of Engineering, Tokushima Bunri University, Sanuki-shi, 769-2193 Japan.
††
The authors are with the Department of Electrical and Electronic Engineering, Faculty of Engineering, Tokushima University,
Tokushima-shi, 770-8506 Japan.
a) E-mail: kawata@fe.bunri-u.ac.jpb) E-mail: taniguchi@ee.tokushima-u.ac.jpc) E-mail: masayoshi@ee.tokushima-u.ac.jpd) E-mail: yamagami@ee.tokushima-u.ac.jpe) E-mail: nishio@ee.tokushima-u.ac.jpf) E-mail: ushida@fe.bunri-u.ac.jp
DOI: 10.1093/ietfec/e90–a.2.406

On the other hand, many algorithms in the time-domain
have been proposed for calculating accurate steady-state
waveforms [6], [7]. Unfortunately, they are rather timeconsuming for the cases of calculating frequency response
curves.
HB (harmonic balance) method is well-known, which
gives good results even for choosing a few essential harmonic components [8], [9]. Thus, the HB method, in contrast to Volterra method, can be stably applied to strong nonlinear systems. Two types of HB simulators are proposed
for the steady-state analysis of electronic circuits. First one
is a nodal type in the meaning that HB method is applied
to all the nodal equations [6], [10]–[12]. The second one isbased on circuit partitioning techniques [7], [13]–[15]. They
are using DFT or FFT at the nodal equations and the partitioning ports. The HB determining equations are solved bynumerical method such as Newton-Raphson and/or relaxation methods. Therefore, the solutions are only found in
the numerical forms.
We propose here Spice-oriented HB method for obtaining the frequency response curves of nonlinear circuits.
Firstly, we have developed Fourier transformation circuit
executing the Fourier expansions to nonlinear devices. Using the circuit, all kinds of the nonlinear devices such as
bipolar transistors and MOSFETs are transformed into thecorresponding HB modules. Combining with them, a given
circuit is transformed into the equivalent HB circuit in theschematic form or net-list [16], [17], which corresponds tothe determining equation of HB method and can be solved
by DC analysis of Spice simulator or solution curve tracing method [18], [19]. Thus, our HB method can continuously calculate the characteristic curves like as Volterra series method.
2.

Fourier Transformation Circuit Model

Analog integrated circuits consist of many kinds of nonlinear devices such as diodes, bipolar transistors and MOSFETs, whose Spice models are usually described by theseveral special functions such as exponential, square-root,
piecewise continuous functions and so on [4], [5]. For thesedevice models, the Fourier coeﬃcients of the output waveforms cannot be described by the analytical forms. Therefore, we will introduce Spice-oriented Fourier expansion
techniques to these device models. Once we have obtained
the HB modules, we can easily carry out the frequency-

c 2007 The Institute of Electronics, Information and Communication Engineers
Copyright 

LETTER

407

then the Fourier transformation circuit model for calculating
the Nth harmonic component is shown in Fig. 1, where
⎫

⎪
1 2π
1
⎪
⎪
⎪
(ĝ0 + ĝ2K )ĝ(v) cos Nθ dθ =I2N−1 =
⎪
⎪
⎪
π 0
2K
⎪
⎪
⎪
⎪
⎪
⎪
1
⎪
⎪
+ (ĝ1 cos Nθ1 + ĝ2 cos Nθ2 +· · ·+ ĝ2K−1 cos Nθ2K−1 ) ⎪
⎬
K
.
⎪
⎪
⎪
 2π
⎪
⎪
⎪
1
1
⎪
⎪
I2N =
ĝ(v) sin Nθ dθ = (ĝ1 sin Nθ1 + ĝ2 sin Nθ2 ⎪
⎪
⎪
⎪
π 0
K
⎪
⎪
⎪
⎪
⎭
+ · · · + ĝ2K−1 sin Nθ2K−1 )(6)
Fig. 1

Fourier transformation circuit model.

domain analysis.
To understand our Fourier transformation circuit
shown in Fig. 1, we consider one-dimensional Fourier expansion of a nonlinear function given by
i = ĝ(v).

Each block in Fig. 1 is realized by the ABMs of Spice [20],and the outputs of this circuit provide the terms I2N−1 and I2Nin (6), where the interval [0, 2π] of the integration is divided
by 2K equal divisions. The value of θk = 2πk/2K is givenby the node voltage at the kth resistor in the resistive circuit.
To investigate the accuracy of our Fourier transformation circuit, we calculate the following Fourier expansion:(1)

Suppose the input and output waveforms as follows:
⎫
M

⎪
⎪
⎪
v(t) = V0 +(V2k−1 cos kωt + V2k sin kωt) ⎪
⎪
⎪
⎪
⎪
⎬
k=1
,
⎪
M
⎪

⎪
⎪
⎪
⎪
i(t) = I0 +{I2k−1 cos kωt + I2k sin kωt) ⎪
⎪
⎭

(2)

k=1

where M denotes the highest harmonic component to betaken into account in our analysis. The coeﬃcients can becalculated by the following formulae:
⎫
 2π
⎪
1
⎪
⎪
⎪
ĝ(v) d(ωt)
I0 =
⎪
⎪
⎪
2π 0
⎪
⎪
⎪
⎪
 2π
⎪
⎪
⎪
1
⎪
⎪
ĝ(v) cos kωt d(ωt), ⎪I2k−1 =
⎬
π 0
.
(3)
⎪
⎪
⎪
⎪
 2π
⎪
⎪
⎪
1
⎪
I2k =
ĝ(v) sin kωt d(ωt) ⎪
⎪
⎪
⎪
π 0
⎪
⎪
⎪
⎪
⎪
⎭
k = 1, 2, . . . , M

e x cos θ = I0 (x) + I1 (x) cos θ + I2 (x) cos 2θ + · · · ,(7)whose Fourier coeﬃcients IN (x), N = 0, 1, 2 . . . are given by
the modified Bessel functions as follows: π1
IN (x) =e x cos θ cos Nθ dθ, N = 0, 1, 2, . . . . (8)2π −πAccording to the simulation result for the case of N = 1,
x = 10 and h = 2π/20, I1 (10) = 2761, which is exactly
equal to the table of Bessel function [21]. Thus, Fouriertransformation circuit model can get the suﬃciently accurate solution with 2K = 20 divisions of the interval [0, 2π].
3.

HB Modules of Nonlinear Devices

We will derive HB modules of the nonlinear devices such as
bipolar transistors and MOSFETs using Fourier transformation circuit shown in Fig. 1.
3.1 HB Module of Bipolar Transistor

The coeﬃcients cannot be given in the analytical forms forthe special functions ĝ(v). Therefore, we apply the trapezoidal formula to the integration (3) as follows: bhĝ(v) d(ωt) = (ĝ0 + ĝ2K ) + h(ĝ1 + ĝ2 + · · · + ĝ2K−1 ),2aĝi = ĝ(v(ti )), for ωti = 0, π/K, · · · , (2K − 1)π/K, 2π(4)

Let us obtain HB module to Gummel-Poon model of bipolartransistor. We assume the input waveforms of collector, baseand emitter (resp. vC (t), vB (t) and vE (t)) as follows:

where the step size of the integration is h = (b − a)/2K(=
π/K) for 2K divisions. Then, the truncation error is givenby ĝ(2) h3 K/6, where (2) denotes the second derivative. We
assume the input as follows:

We transform them into the voltage diﬀerences of the basecollector and base-emitterv(θ) = V0 +

M

k=1

(V2k−1 cos kθ + V2k sin kθ), θ = ωt (5)vi (t) = Vi,0 +M(Vi,2k−1 cos kωt + Vi,2k sin kωt),

(9.1)

k=1

i = C, B, E

vBC (t) = vB (t) − vC (t), vBE (t) = vB (t) − vE (t).

(9.2)

Firstly, using the Fourier transformation circuit ofFig. 1, the Fourier coeﬃcients of nonlinear resistive currents {iB (t), iC (t)} and charges of parasitic capacitors

IEICE TRANS. FUNDAMENTALS, VOL.E90–A, NO.2 FEBRUARY 2007

408

Fig. 2

HB module of bipolar transistor.
Fig. 4

Diode circuit. (R = 1 [kΩ], C = 1 [µF], E = 1 [V])(a) DC circuit
Fig. 3

HB module of MOSFET.

{qBC (t), qBE (t)} are calculated. Combining these currents
and charges, the output currents are estimated from the following relations:
dqBE (t) dqBC (t) ⎫⎪⎪⎪+iB (t) = iB (t) +
⎪
⎬
dt
dt ⎪
,
(9.3)
⎪
⎪
⎪
(t)
dq
⎪
BC
⎪
⎭
iC (t) = iC (t) −dtThey are also described by the Fourier expansions(b) cos ωt circuit(c) sin ωt circuit(d) cos 2ωt circuit(e) sin 2ωt circuit(f) cos 3ωt circuit(g) sin 3ωt circuitM

ii (t) = Ii,0 +(Ii,2k−1 cos kωt+Ii,2k sin kωt), i =C, B,

k=1

(9.4)
where the coeﬃcients are given as follows:IB,0 = IB,0
, IB,2k−1 = IB,2k−1+ kω(QBE,2k + QBC,2k ),
IB,2k = IB,2k − kω(QBE,2k−1 + QBC,2k−1 )Fig. 5Equivalent HB circuit of Fig. 4.

(9.5)


, IC,2k−1 = IC,2k−1− kωQBC,2k ,IC,0 = IC,0

IC,2k = IC,2k + kωQBC,2k−1

4.

Illustrative Examples(9.6)

k = 1, 2, . . . , M
Thus, we have the HB module for bipolar transistor as
shown in Fig. 2. The blocks {“cos kωt, sin kωt”, k =1. . . . , M} correspond to the Fourier coeﬃcients (9.5) and(9.6), which are constructed by ABMs of Spice and theFourier transformation circuit of Fig. 1. Note that they arecoupled resistively and do not contain any dynamical element.
3.2 HB Module of MOSFET
There are many types of MOS models depending on the levels of accuracy [4], [5]. Their HB modules are also built by
the same ways as the bipolar transistor. The three terminal
HB module is shown in Fig. 3.

4.1 Example of HB Method and Comparison with Volterra
SeriesTo understand our HB method, consider a simple diode circuit as shown in Fig. 4. We assume the input voltage and thediode current as follows:e(t) = Em cos ωt,

id (t) = 10−8 exp(40vd )(10)Then we set the capacitor and diode voltages as:⎫⎪⎪⎪(VC,2k−1 cos kωt + VC,2k sin kωt) ⎪vC (t) = VC,0 +
⎪
⎪
⎪
⎪
⎬
k=1
,
⎪
3
⎪

⎪
⎪
⎪
⎪
(Vd,2k−1 cos kωt + Vd,2k sin kωt) ⎪vd (t) = Vd,0 +
⎪
⎭
3


k=1

(11)
and the Fourier expansion of diode as:
id (t) = Id,0 +3(Id,2k−1 cos kωt + Id,2k sin kωt).
k=1

(12)

LETTER

409

(a)
(a)

(b)
Fig. 6(a) 2nd and (b) 3rd order harmonic distortion. (Em = 1.0 [V])

Table 1
Em
0.1
0.5
0.7
1.0
2.0

Comparison of accuracy. (ω = 20 × 103 )Volterra series
HD2
HD3
0.01608 0.00031
0.08044 0.00773
0.11262 0.01515
0.16088 0.03092
0.32176 0.12367

Harmonic balance
HD2
HD3
0.01602 0.00031
0.07526 0.00679
0.09898 0.01179
0.12436 0.01886
0.14042 0.03438

Transient analysis
HD2
HD3
0.01276 0.00184
0.07113 0.00491
0.09803 0.01002
0.12978 0.01873
0.16978 0.03546

We have the equivalent HB circuit as shown in Fig. 5. The
nonlinear controlled currents { Id,k , k = 0, 1, . . . , 6 } are thefunctions of { Vd,k , k = 0, 1, . . . , 6 }, which are obtained by
the Fourier transformation circuit of Fig. 1. Note that theequivalent HB circuit of Fig. 5 corresponds to the determining equation of the HB method, and can be solved by
DC analysis of Spice. We also solved the circuit of Fig. 4
by Volterra series method, and compared both results. The
2nd and 3rd order harmonic distortion curves are shown in
Fig. 6.
We also compared them with the results obtained bytransient analysis and FFT† . Table 1 shows the results fordiﬀerent input voltages. Note that there are slight diﬀerences between HB and Volterra methods for small input,
but the results are completely diﬀerent for large input. Although we neglected harmonics larger than 3rd order in our
HB method, HB method is more accurate than the Volterra
series method.

(b)
Fig. 7 (a) RF power amplifier, (b) Frequency response curve. L1 = L2 =
0.1 [mH], Lg = 10 [nH], C1 = C2 = C3 = C4 = 1 [nF], VCC = 3 [V],Vbb = 3 [V], Rb1 = Rb2 = 1 [Ω], Rb = 1 [Ω], ICC = 1 [mA], RL = 36 [Ω],vin = 0.1 cos ωt.

4.2 RF Power Amplifier
We consider a high frequency RF power amplifier shown in
Fig. 7(a), which is used in a mobile-phone [23], where Lg ’s
are parasitic inductances. Firstly, the transistors are modeledby the Gummel-Poon model, whose typical parameters are
given in [5]. We considered up to 3rd harmonic components
and implemented the equivalent HB circuit with transistor
modules of Fig. 3.
We have obtained the frequency response curves forω = 108 –1011 [rad/sec], as shown in Fig. 7(b), where only
fundamental frequency component is shown, because thewaveform is nearly sinusoidal. The waveform is suﬃciently
accurate in comparison to the result of the transient analysis by Spice. It takes 41 [sec] to get the frequency response
curve with PSpice†† .

†The results may have truncation errors in transient analysis,
but they can be assumed to be accurate.
††
The simulation is carried out by a PC with the following spec;
CPU: PentiumM 733 (1.10 GHz), Main memory: 256 MB, OS:
Windows XP Professional SP2.

IEICE TRANS. FUNDAMENTALS, VOL.E90–A, NO.2 FEBRUARY 2007

410

5.

Conclusions and Remarks

In this letter, we have proposed Spice-oriented HB (harmonic balance) method for calculating frequency response
curves. At first, the nonlinear devices such as bipolar transistors and MOSFETs are transformed into the HB modules.
Using these device modules, we transformed the nonlinear
circuit into the equivalent HB circuit which corresponds tothe determining equations of the HB method. It consists of(2M + 1)-coupled sub-circuits when we consider M higher
harmonic components, and it can be easily solved with the
DC analysis of Spice.
References[1] M. Schetzen, The Volterra and Wiener Theorems of Nonlinear Systems, John Wiley and Sons, 1978.
[2] J. Wood and D.E. Root, Fundamentals of Nonlinear Behavioral
Modeling for RF and Microwave Design, Artech House, 2005.
[3] P. Wambacq and W. Sansen, Distortion Analysis of Analog Integrated Circuits, Kluwer Academic Pub., 1998.
[4] M. Miura, T. Shiino, and K. Mori, Circuit-Simulation Technique and
MOSFET Modeling, Riaraizu Rikou Senta, 2003.
[5] A.S. Sedra and K.C. Smith, Microelectronic Circuit, Oxford Univ.
Press, 2004.
[6] K.S. Kundert, J.K. White, and A. Sangiovanni-Vincentelli, SteadyState Methods for Simulating Analog and Microwave Circuits,
Kluwer Academic Pub., 1990.
[7] A. Ushida, T. Adachi, and L.O. Chua, “Steady-state analysis of nonlinear circuits based on hybrid methods,” IEEE Trans. Circuits Syst.
I, Fundam. Theory Appl., vol.39, no.9, pp.649–661, 1992.
[8] C. Hayashi, Nonlinear Oscillations in Physical Systems, McGrawHill, 1964.
[9] Y. Ueda, The Road to Chaos-II, Aerial Press. Inc., 2001.
[10] R. Telichevesky, K.S. Kundart, and J.K. White, “Eﬃcient steadystate analysis based on matrix-free Krylov-subspace methods,” Proc.
ACM/IEEE Design Automation Conference, pp.480–485, 1995.

[11] R.J. Gilmore and M.B. Steer, “Nonlinear circuit analysis using themethod of harmonic balance: A review of the Art. Part I. Introductory concepts,” Int. J. Microw. Millim. Wave Comput.-Aided Eng.,
vol.1, pp.22–37, 1991.
[12] R.J. Gilmore and M.B. Steer, “Nonlinear circuit analysis using themethod of harmonic balance: A review of the Art. Part II. Advanced
concepts,” Int. J. Microw. Millim. Wave Comput.-Aided Eng., vol.1,
pp.159–180, 1991.
[13] H.G. Brachtendorf, G. Welsh, R. Laur, and A. Bunse-Gerstner,“Numerical steady state analysis of electronic circuits driven bymulti-tone signals,” Electrical Engineering, vol.79, no.2, pp.103–
112, Springer-Verlag, 1996.
[14] Y. Yamagami, Y. Nishio, A. Ushida, M. Takahashi, and K. Ogawa,“Analysis of communication circuits based on multidimensional
Fourier transformation,” IEEE Trans. Comput.-Aided Des. Integr.
Circuits Syst., vol.18, no.8, pp.1165–1177, 1999.
[15] S.A. Mass, Nonlinear Microwave Circuits, Artech House, 1988.
[16] A. Ushida, Y. Yamagami, and Y. Nishio, “Frequency responses ofnonlinear networks using curve tracing algorithm,” ISCAS 2002,
vol.I, pp.641–644, 2002.
[17] Y. Yamagami, H. Yabe, Y. Nishio, and A. Ushida, “Distortion analysis of nonlinear networks based on Spice-oriented harmonic balancemethod,” ISCAS04, pp.633–636, 2004.
[18] A. Ushida, Y. Yamagami, Y. Nishio, I. Kinouchi, and Y. Inoue, “Aneﬃcient algorithm for finding multiple DC solutions on the Spiceoriented Newton homotopy method,” IEEE Trans. Comput.-Aided
Des. Integr. Circuits Syst., vol.21, pp.337–348, 2002.
[19] A. Ushida and M. Tanaka, Electronic Circuit Simulation, pp.166–
186, Corona Co., 2002.
[20] MicroSim: PSpice A/D Circuit Analysis User’s Guides, MicroSimCo., 1995.
[21] K.K. Clarke and D.T. Hess, Communication Circuits: Analysis andDesign, Addison-Wesley Pub. Co., 1971.
[22] J.E. Meyer, “MOS models and circuit simulations,” RCA Review,
vol.32, no.1, pp.42–63, 1971.
[23] D. Leenaerts, J.V. Tang, and C. Vaucher, Circuit Design for RF
Transverse, Kluwer Academic Pub., 2001.
[24] R. Gregorian, Introduction to CMOS OP-Amps and Comparators,
John Wiley and Sons, 1999.