2007 International Symposium on Nonlinear Theory and itsApplications
NOLTA'07, Vancouver, Canada, September 16-19, 2007

Spice-Oriented Harmonic Balance
and Volterra Series Methods
Takaaki Kinouchi†, Yoshihiro Yamagami†, Yoshifumi Nishio†, Akio 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

response curves. We propose here another type of Spiceoriented HB method which is based on the classical HB. It
is largely improved such that that it can be easily applied tolarge scale nonlinear circuits combining with MATLAB andSpice, where the determining equation called Sine-Cosine
circuit [11,12] can be obtained in the form of schematic
and/or net-list of Spice. The circuit can be solved by DC
analysis of Spice simulator and easily get the frequency characteristic.
We briefly show the Volttera series method in Section 2,
and our HB method in Section 3. Comparisons between our
HB and the Volterra methods are given in Section 4.

There are typical two types of frequency-domain analysis
of nonlinear circuits. One is Volterra series method which
is widely used in European countries, because the algorithm
is based on the bilinear theorem, and that the higher order solutions can be obtained in the analytical forms. Unfortunately, it can be only applied to the weakly nonlinear
systems. The other is the harmonic balance (HB) method.
In this paper, we have greatly improved the HB by combining with MATLAB and Spice, such that it can be efficiently applied to any kind of nonlinear electronic circuits.
Namely, the Fourier coefficients can be obtained in the symbolic forms with MATLAB. The determining equation ofharmonic balance method is given by the net-lists, which
can be solved by the DC analysis of Spice. We will compare
the Volterra series with HB methods in the point of theiraccuracies.

2. Volterra series method
To understand the Volterra series method [3-5], we consider a simple circuit shown in Fig. 1, where nonlinear resistor is described by the following polynomial function;

1. Introduction23
iG = k1 vG + k2 vG+ k3 v G.

For designing integrated circuits and communication systems, it is very important to analyze the frequency-domaincharacteristics of nonlinear circuits. There are two kinds ofthe time-domain and frequency-domain methods. Although,the first one can get the accurate steady-state waveform, it
is inefficient to calculate the frequency response curve[1,2].
The latter is useful to obtain the frequency response curve ofnonlinear circuits. There are two methods of Volterra series
and harmonic balance. The Volterra series method is widelyused for the purposes [3-5] which gives the solution in analytical form, and the algorithm can be easily used combining with MATLAB [6]. The algorithms are based on bilinear
theorem, and can be also effectively applied to the intermodulation’s analysis driven by multiple inputs. However, they
can be only applied to the weakly nonlinear circuits, and thesolutions for strong nonlinear circuits become erroneous.
The harmonic balance (HB) is also a well-known method
for the frequency domain analysis, which gives good results
even for relatively strong nonlinear circuits. There are two
types of classical HB [7-9] and computer-aided HB methods [10,11]. The classical HB can be applied to relatively
small scale circuits, but they have found many interesting
nonlinear phenomena such as bifurcations and chaos. The
latter HBs are rather numerical methods, and find the solutions with the Newton-Raphson and/or relaxation methods. Therefore, they are inefficient to calculate the frequency

(1)

iG

+
vG
-

LinearEm

NetworkFig. 1 Nonlinear resistive circuit.
Now, let us caluculate the kernels in complex frequencydomain “s”. The schematic diagram is shown in Fig. 2. The
first order kernel corresponds to the output voltage for thenormalized unit inputs U (s).
Linear+V2 (s1 ,s2 ) (b)
2nd order
I2 (s1 ,s2 )rLinearU(s)rNetwork+
V1 (s)-

Network

Linear(a)r

1st order

Network+
V3 (s1 ,s2 ,s3 )(c)
3rd order
I3 (s1 ,s2 ,s3 )Fig. 2 Volterra kernels.
and the 1st order kernel is given by
H1 (s) = V1 (s) = K(s)U (s)where K(s) is the transfer function.

1- 513 -

(2.1)

Next, we consider the 2nd order kernel. The current isobtained by the 1st order kernel.
I2 (s1 , s2 ) = k2 H1 (s1 )H1 (s2 )(2.2)

for integers m1k , m2k , . . . , mpk . The Fourier expansion i(t)for the input v(t) can be carried out symbolically with MATLAB.
Now, we apply Fourier expansion to the capacitive element described by
qC = q̂C (vC ),
(6.1)

(2.3)

and get the current as follow:

and the 2nd order kernel is obtained as follows.
H2 (s1 , s2 ) = Y (s1 + s2 )I2 (s1 , s2 )where Y (s1 + s2 ) is the admittance at the output port.
Finally, the 3rd order current is obtained as follows;iC =I3 (s1 , s2 , s3 ) = k3 (H1 (s1 )jH1 (s2 )H1 (s2 )+ 32 k2 (H1 (s1 )H2 (s2 , s3 ) + H1 (s2 )H2 (s1 , s3 ) + H1 (s3 )H2 (s1 , s2 ))(2.4)

IC,2k−1 = νk QC,2k , IC,2k = −νk QC,2k−1 , k = 1, 2, . . . , M.
(6.3)
Thus, the HB models of capacitors are replaced by thevoltage-controlled current sources which are function of thecapacitor voltages.
Next, we apply Fourier expansion to the inductive element
described by
φL = Φ̂L (iL ),
(7.1)

(2.5)

Thus, the output is given by
Hout (s1 , s2 , s3 ) = H1 (s1 )+H2 (s1 , s2 )+H3 (s1 , s2 , s3 ). (2.6)
For the input(3)

and get the voltage as follow:

by setting si = ±jω, (i=1,2,3), the output amplitudes canbe estimated by the following Table 1.

vL =Table 1 Frequency response by Volterra series method.
Frequency Amplitude of response2
DC| 21 Em
H2 (jω, −jω)|3
ω|Em H1 (jω) + 43 Em
H3 (ω, jω, −jω)|
1 22ω
| 2 Em H2 (jω, jω)|33ω| 41 Em
H3 (ω, jω, jω)|

Cosine component
Inductor
Capacitor= ^(i L )i(t) = I0 +(I2k−1 cos νk t + I2k sin νk t) ResistorvGiG =i^G (vG )



L,2k


VL,2k= − kL,2k-1VC,2k-1

I C,2k-1 = k Q C,2k


VC,2k

I C,2kIC,2k = − k Q C,2k-1
^^

I G,2k-1

I G,2k
VG,2k

VG,2k-1^

I G,2k-1 = I G,2k-1

I G,2k = ^I G,2k

Fig. 3 Sine-Cosine equivalent circuit models.
Thus, the HB models of inductors are replaced by thecurrent-controlled voltage sources, which are function of theinductor currents. All the nonlinear devices such as bipolartransistors and MOSFETs can be replaced by the equivalent
HB modules with MATLAB.
Replacing each element by the corresponding HB model,
we can formulate the Sine-Cosine circuit that is equal to the
HB determining equation. The circuit is efficiently solved by
the DC analysis of Spice.

(5.1)

k=1

where
νk ≡ m1k ω1 + m2k ω2 + · · · + mpk ωpVL,2k

I C,2k-1

vC

iG

(4)

,

VL,2k-1 = k

I L,2kiCq (v C )
q =^

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

k=1

M
X

VL,2k-1vL

Sine component

I L,2k-1

iL

Suppose the circuit is driven by the independent “p” frequency components ω1 , ω2 , . . . , ωp . The input and output
waveforms of (4) is assumed as follows;
v(t) = V0 +(7.2)

VC,2k−1 = νk ΦL,2k , VC,2k = −νk ΦL,2k−1 , k = 1, 2, . . . , M.
(7.3)

In this Section, we propose our Spice-oriented HB
method. Analog integrated circuits are usually composedof many kinds of nonlinear devices such as diodes, bipolar transistors and MOSFETs, whose Spice models are described by the several special functions such as exponential,
square-root, piecewise continuous functions and so on [13].
For the special functions, the Fourier expansions can not be
described by analytical forms. Therefore, we need to approximate them with the Taylor expansions as follows;MX

νk (−ΦL,2k−1 sin νk t + ΦL,2k cos νk t)

where Φ1 , Φ2 , . . . , Φ2M −1 , Φ2M are Fourier coefficient from
(7.1). Thus, we have

3. Spice-oriented HB methodf (v)k0 + k1 v + k2 v 2 + · · · + kK v K .

M
X
k=1

Note that the solution of Volterra series method can be easily
solved by the use of MATLAB.

i=
=

(6.2)

where Q1 , Q2 , . . . , Q2M −1 , Q2M are Fourier coefficient from
(6.1). Thus, we have

V3 (s1 , s2 , s3 ) = H3 (s1 , s2 , s3 ),e(t) = Em sin(ωt),νk (−QC,2k−1 sin νk t + QQ,2k cos νk t)

k=1

and the 3rd order voltage is obtained as follows;= Y (s1 + s2 + s3 )I3 (s1 , s2 , s3 ).

M
X

(5.2)

2

- 514 -

4. Comparisons between Volterra and
HB

Note that the equivalent HB circuit of Fig. 5 corresponds
to the HB determining equation, and can be solved by DC
analysis of Spice. We have calculated the 2nd-3rd orderharmonic distortions (HD2, HD3) to the fundamental harmonic. Thus, we obtained the results shown by Table 2.
We have compare the results with those of Volterra method
and the results from the transient analysis1 . When the input is small, the values of both methods are similar. For
the large inputs, the results from Volterra method may beerroneous. They are coming from the philosophy of Volterra
series method.

Consider a simple diode circuit driven by a single input asshown in Fig. 4(a). In this Section, we compare the Volterra
series method with our HB method in the point of accuracy.

R
e(t)

C

E

+ id
vc
-

+
vd
-

Table 2 Comparisons (ω = 20, 000).

Fig. 4 Diode circuit. R=1[kΩ], C=1[µF], E=1[V],e(t)=Em cos ωt.

Em
0.1
0.5
0.7
1.0
2.0
5.0

Assume the diode characteristic as follow;
id = 10−8 exp(40vd ).

(8.1)

It is approximated by the 3rd order Taylor expansion at vd0
as follow;
id = k0 (Vd0 ) + k1 (Vd0 )vd + k2 (Vd0 )vd2 + k3 (Vd0 )vd3

3
X

(Vd,2k−1 cos kωt + Vd,2k sin kωt).

(9)

[V]

[V]

0.3

Substituting it into (8.2), the Fourier series is described asfollow;

X

2nd harmonic distortion

0.12

E in = 1.0[V]
0.2

3rd harmonic distortion
Ein=1.0[V]

0.08

Volterra series

3(Id,2k−1 cos kωt + Id,2k sin kωt).

Transient
HD2
HD3
0.0128 0.0018
0.0711 0.0049
0.0980 0.0100
0.1298 0.0187
0.1698 0.0355
0.1280 0.0371

We have also solved the circuit in Fig. 4 with Volterra series
and our HB methods in the frequency-domain. The 2nd and3rd order harmonic distortion curves are shown in Figs. 6(a)and (b)2 .

k=1

id = Id,0 +

HB
HD2
HD3
0.0160 0.0003
0.0753 0.0068
0.0990 0.0118
0.1244 0.0189
0.1404 0.0344
0.0872 0.0709(8.2)

for an unknown value Vd0 . We assume the waveform up to
3rd higher harmonic as follows;
vd = Vd,0 +Volterra
HD2
HD3
0.0161 0.0003
0.0804 0.0077
0.1126 0.0151
0.1609 0.0309
0.3218 0.1237
0.8043 0.7729

0.1

(10)

Harmonic balance

k=1

0
0

In the arithmetic computaion, we have

Volterra series

0.04

20

40

0

100
80
3
[x10 rad/sec ]60

Harmonic balance
20

0

3

1 + 7 + 7 + 7 = 400

5. An illustrative example
We consider a base module circuit driven by two inputs asshown in Fig. 7(a) [14]. Firstly, let us solve it with Volterra
series method.

R(a) DC circuit
I d,0

E

C

R

C V2

R
+

+

V1

Vd,1

-

-

C V1

C0

+

+

V2

Vd,2

-

I d,2

Q
e1 (t)
e (t)

-

R2

I

L

E


V (j

V (j
1

1

2

2

C0

d

Id
g

d

+
C

R

L Vout-2(c) sin t circuit

R

R
(a)

+

+

V3

Vd,3
I d,3 --+
- 2 C V3

V4
-

V
I d,4

(d) cos 2 t circuit
R
+

+

V5

Vd,5

(f) cos3 t circuit

I d,5 -d,4-

R1 =50[kΩ], R2 =10[kΩ], R=1[kΩ], C=1[nF], C0 =1[nF],C0 =1[nF], L=10[µF], E=10[V], v1 =0.01sin ω1 t, v2 =0.01sin ω2 t,
for ω1 =8.2×106 [rad/sec].
1 We used Ebers-Moll model for the bipolar transistor in thetransient analysis. Note that the approximation error of HB
method may be arisen from the Taylor expansion and neglecting the higher harmonic component than 3rd order.
2 At the lower frequency, the results from Volterra are erroneous because the diode voltage is large at the frequencies.

(e) sinRt circuit++- 3 C V5

V6

Vd,6

-

I d,6

(b)

Fig. 7 (a) Base modulation circuit. (b) 1st order Volterra
equivalent circuit.

+

3 C V6

RR1

I d,1 -(b) cos t circuit

2 C V480
100[x10 rad/sec](b)Fig. 6 (a) 2nd order harmonic distortion. (b) 3rd orderharmonic distortion.

terms. It is not easy to estimate all the Fourier coefficientsby hand calculations. Hence, we apply MATLAB and getthem in the formula symbolic forms. Thus, we have theequivalent HB circuit as shown in Fig. 5.

Em

60

40

(a)

2

-

(g) sin3 t circuit

Fig. 5 Equivalent HB circuit of Fig. 4.

3- 515 -

Since the input voltages in this example are small, the bothresults of Volterra series method and HB method are almost
the same.

We consider here 8 intermodulation frequencies given in Table 3.
Table 3 The frequency components to take into theconsideration.
ν1 = ω 1
ν5 = ω1 − ω2ν2 = ω 2
ν6 = ω1 + ω2
ν3 = 2ω1 ν7 = 2ω1 − ω2
ν4 = 2ω2 ν8 = 2ω1 + ω2

6. Conclusions and remarks
In this paper, we proposed Spice-oriented HB method
which can efficiently trace the characteristic curves such as
frequency responses. For the nonlinear elements such as
bipolar transistor, the Fourier coefficients are calculated by
MATLAB in the symbolic forms. Thus, we can formulate
the HB determining equation in the form of net-list and solve
it with DC analysis of Spice. We have compared the methodwith Volterra series method, and found that our HB method
can get much accurate solutions compared to Volterra series
method.
As the future work, in order to apply our method to large
scale networks, we need to develop device HB modules forbipolar transistors and MOSFETs.

We have modeled the transistor by the Ebers-Moll model
where
id = Is eλvd , for Is = 10−12 [A], λ = 40, and α = 0.99. (11)Firstly, we calculate the operating point of diode voltage,and have vd0 = 0.5822[V ]. It is approximated by the Taylor
expansion at the point vd0 as follow;
id = 0.0130 + 0.52vd + 10.40vd2 + 138.60vd3 .

(12)Thus, we have

References

gd = 0.52, k2 = 10.40, k3 = 138.60.

(13)[1] T.J.Aprille and T.N.Trick, “Steady-state analysis of nonlinear circuits with periodic input, Proc. IEEE, vol.60, pp.108114, 1972.

7

We have changed the frequency ω2 in the region [0, 3 × 10 ],and get the frequency characteristics as shown Figs. 8(a)and (b), where we consider up to the 3rd order Volterra
kernels, and solved with MATLAB.
5
[V]

0.08
2

Frequency characteristic
Volterra 1st and 2nd order

4

[V]
0.06

1

2

0

10

1

1

2

0

[4] J.Wood and D.E.Root, Fundamentals of Nonlinear Behavioral Modeling for RF and Microwave Design, Artech
House, 2005.

2

0.04

   1[3] M.Schetzen, The Volterra and Wiener Theorems of Nonlinear Systems, John Wiley and Sons, 1978.

characteristic      Frequency
Volterra 3rd order1212
    32[2] R.Telichevesky and K.Kundert, SpectreRF Primer, Cadence
Design Systems, San Jose, 1996.

2

0.02


1



[5] P.Wambacq and W.Sansen, Distortion Analysis of Analog
Integrated Circuits, Kluwer Academic Pub., 1998.



2

1

2

2
20

30

[6] S.Oishi, Numerical Ananlysis using MATLAB (Japanese),
2001.

0
10

0

[MHz]

20

[MHz] 30

Fig. 8 (a) 1st and 2nd order solutions. (b) 3rd ordersolutions.

[7] C.Hayashi, Nonlinear Oscillations in Physical Systems,
McGraw-Hill, 1964.
[8] Y.Ueda, The Road to Chaos-II, Aerial Press. Inc., 2001.

Next, we approximate the nonlinear diode by Taylor expansion in our HB method as follow;
id = k0 (vd0 ) + k1 (vd0 )v + k2 (vd0 )v 2 + k3 (vd0 )v 3 .

[9] R.J.Gilmore and M.B.Steer, “Nonlinear circuit analysis using the method of harmonic balance-A review of the Art.
Part I. Introductory concepts,” Int. J. of Microwave andMillimeter-Wave Computer-Aided Eng., vol.1, pp.22-37,
1991.

(14)

Remark that vd0 is unknown variable in HB. The waveform
is approximated by the frequency components in Table 3;v(t) = V0 +8
X

(V2k−1 cos νk t + V2k sin νk t).

[10] R.Telichevesky, K.S.Kundart and J.K.White, “Efficient
steady-state analysis based on matrix-free Krylov-subspacemethods,” ACM, pp.480-485, 1995.

(15)[11] A.Ushida, Y.Yamagami and Y.Nishio, “Frequency responses
of nonlinear networks using curve tracing algorithm,”Proc.
ISCAS 2002, vol.1, pp.641-644, 2002.

k=1

The Sine-Cosine circuit can be formulated by MATLAB, andis solved by the DC analysis of Spice.

[12] J.Kawata, Y.Taniguchi, M.Oda,Y.Yamagami, Y.Nishio and
A.Ushida, “Spice-oriented frequency-domain analysis ofnonlinear electronic circuits,” IEICE Trans. Fundamentals,
vol.E90-A, no. 2, pp. 406-410, Feb. 2007.

0.08
6
[V]

Frequency characteristic
Harmonic balance[V]

Frequency characteristic
Harmonic balance
2

5
4

0.06



21

12[13] A.S.Sedra and K.C.Smith, Microelectronic Circuits, Oxford
Univ. Press, 2004.



1

2

0.04

3
2

 

1

1

0
0

2

  

0.02

1

1



2

1

2

2

0

2
10

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



2

20

[MHz]

30

0

10

20

[MHz]

30

Fig. 9 (a) Frequency characteristic by HB. (b) Other
frequency components.

4- 516 -