The 15th IEEE International Workshop on Nonlinear Dynamics of Electronic Systems
NDES'07, Tokushima, Japan, July 23-26, 2007

Comparisons Between Spice-Oriented Harmonic Balance
Method and Volterra Series Method
Junji Kawata†, Takaaki Kinouchi‡, Yoshihiro Yamagami‡, Yoshifumi Nishio‡, Akio Ushida††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
It is very important to analyze the mixer and modulator
circuits driven by multiple inputs frequencies. Both Volterra
series and HB (harmonic balance) methods are widely used
for the frequency-domain analysis of nonlinear circuits. We
compared the errors of the above two methods using illustrative examples. In our HB method, the Fourier expansion
is executed by MATLAB in symbolic form, and it is transformed into the net-list with Spice. We found that it is sometime happened to exceed the capacities of Spice, when wechoose higher order polynomial approximations of the nonlinear devises, because the Fourier components are rapidly
increased. Hence, we propose here a new perturbation like
Fourier expansion method to the approximations.

1. Introduction
Frequency-domain analysis of nonlinear electronic circuits
driven by multiple inputs are very important for designing RF circuits and communication systems. Volterra series
methods are widely used for the analysis [1-3], because theygive the solutions in analytical forms. There have been published many papers about computer-aided HB(harmonic balance) methods [4-6] during 1990-2000 which are using DFT,FFT and/or interpolation methods for the Fourier expansion. They apply the Fourier expansions to all the nodal
equations of the circuit, so that the resultant systems become usually very large. Hence, they propose to solve themwith the applications of a method like the Krylov-subspace
GMRES algorithm instead of Newton method which will reduce the computation time. It seems that their methods arestill time-consuming for large scale circuits.
On the other hand, the classical HB methods [7-8] can beapplied to relatively small scale circuits, which find many
interesting nonlinear phenomena such as bifurcations andchaos phenomena. Now, we are going to apply it to electronic circuits containing transistors and diodes. These devices have strong nonlinear characteristics and they are modeled by special functions such as exponential in Spice. In
order to apply the HB method, their characteristics must
be approximated by the polynomial functions. In this way,
the Fourier expansions of the nonlinear devices are stored as
modules in the net-lists of Spice with MATLAB. In this case,
we have sometime happened a serious problem such thatwhen we apply high order polynomial approximations to thedevices, the Fourier expansion contains too many terms andit exceeds the capacity of Spice. So, we propose here a new
perturbation like Fourier expansion method using MATLAB
which reduces the resultant terms and the size of Spice netlist. The ideas are shown in section 2.
Using the above Fourier expansion method, nonlinear devices are replaced by HB modules as follows: Firstly, bipolartransistors and MOSFET are approximated by the Taylor
expansions in the polynomial forms. The Fourier expansions can be executed by MATLAB in the symbolic forms
and they are stored as the modules of bipolar transistor in
our Spice library. Thus, it is possible to formulate the determining equation of HB in the forms of net-list with Spice.
It can be easily solved by the DC analysis of Spice and thefrequency-domain solutions are easily obtained.
Using an example, we will compare the results of Volttera
series method and our HB method in section 2. In section
3, we show an interesting example of a mixer circuit.

2. Fourier expansion of nonlinear
devices
2.1 Taylor approximation of nonlinear devices In thissub-section, we study the errors of HB method when we applied the Fourier expansion to the Taylor approximation ofan exponential function. Transistors are usually modeled by
Ebers-Moll model consisted of two diodes, whose nonlinear
characteristic (id , vd ) is described by exponential function as
follows;
id = Is exp vd /VT ,VT ' 0.025.
(1)
Firstly, let us calculate the Fourier coefficients DC, cos ωt,
cos 2ωt, cos 3ωt of (1) driven by
vd = Vd0 + Vd cos ωt,for Vd0 = 0.5822,(2)

We also compare them with the results when we applied it
to the 3rd order Taylor approximation as follows:
id = k0 + k1 vd + k2 vd2 + k3 vd3 ,(3)for Is = 10−12 , VT = 0.026 in (1), which are shown in Table
1.
Table 1.
Exact Fourier coefficients and those of 3rd order Taylor expansion.
Exact solutions
3rd order Taylor series
Vd
DC Fund. 2nd 3rd
DC Fund. 2nd 3rd
0.02 0.0151 0.0113 0.0022 0.0003 0.0151 0.0112 0.0021 0.0003
0.05 0.0295 0.0413 0.0179 0.0055 0.0259 0.0390 0.0130 0.0043
0.07 0.0538 0.0859 0.0469 0.0191 0.0383 0.0720 0.0255 0.0119
0.10 0.1460 0.2542 0.1624 0.0870 0.0647 0.1560 0.0520 0.0347- 14 -

We can see that the errors for the small input Vd are smallenough, but they are rapidly increased for the large input.
Next, when the higher order Taylor approximations are used,
the results for the input such as Vd = 0.1 are shown together
with the exact ones in Table 2.

Firstly, it is approximated by Mth order Taylor expansion
as follow:
id = k0 (Vd0 ) + k1 (Vd0 )vd + · · · + kM (Vd0 )vdM

with an unknown operating point Vd0 . Now, we rewrite (9)into the following Horner’s scheme

Table 2.Solutions of order of the Taylor expansions and theexact solutions, for Vd = 0.1
DC
Fund.
2nd
3rd
Exact
0.14590.2538 0.1670
0.0870
3rd order 0.0646
0.1559 0.0519
0.03475th order 0.1162
0.2250 0.1216
0.0696
7th order 0.1397
0.2483 0.1559
0.0832

id1 (vd ) = KM −1 (Vd0 ) + KM (Vd0 )vd 
id2 (vd ) = kM −2 (Vd0 ) + id1 (vd )vd,..........................................

idM (vd ) = k0 (Vd0 ) + idM −1 (vd )vd

Observe that if we apply it to 7th order Taylor approximation, the errors is again reduced to small enough. Thus,
we found that the higher order Taylor approximation of the
exponential function is needed for larger input.
2.2 Perturbation like Fourier expansion algorithm In
practical transistor circuits, the input signals are sufficient
small compared to the DC bias voltages. Let us describe it
as follows:vin (t) = ε(V1 cos ω1 t + V2 cos ω2 t).

XK

vd (t) = Vd0 +
k=1

εk

" 2k
X

ε order ν1m
ε2 order ν2m
ε3 order ν3m

···
···
···
···

ω1 , omega2
2ω1 , 2ω2 , ω1 + ω2 , |ω1 − ω2 |
3ω1 , 3ω2 , 2ω1 + ω2 , ω1 + 2ω2 ,|2ω1 − ω2 |, |ω1 − 2ω2 |(6)

In this case, two cross term v1 (t) × v2 (t) having the formsof (5.1) is given by' V30 +k=1

εk

v3 (t) = v1 (t) × v2 (t)

" 2k
X

#

(V3,k,2m−1 cos νkm t + V3,k,2m sin νkm t) .
m=1

(7)

If we consider all the terms lower than Kth order about ε,they can be estimated by the combinations of the terms of(5.1) where the higher order terms than ε( K + 1) are neglected.
Thus, we call the algorithm perturbation like Fourier expansion in the meaning that our method contain all the terms
less than order Kth of ε. Thus, we can reduce the numberof cross terms. Observe that if choose a larger K, we can
expect much more exact Fourier expansion.
2.3 Spice-oriented HB method We propose here a
Spice-oriented HB (harmonic balance) method combining
with MATLAB, where the determining equation is given by
the equivalent circuit and/or net-list of Spice [9]. Firstly,
we describe the nonlinear characteristics in the polynomial
forms. The schematic diagram of diode example is shown in
Fig.1, where it is driven by the two inputs
v1 (t) = εV1 cos ω1 t,

v2 (t) = εV2 cos ω2 t.

-

vd

I d2K sinνKt

Vd2K sinνKt

(a)

(b)

Fig.1 (a) a diode model, (b) HB module with MATLAB.

The circuit is resistive, and it can be easily solved with
the DC analysis of Spice. Thus, we can get the frequencydomain solutions such as frequency response curves and soon.
It is convenient to develop the HB modules for bipolartransistors. In this case, if we model the transistor with
Ebers-Moll model, the emitter, base and corrector currents
are given as follows:



etc.

K
X

Id0
I d1 cosν1 t
I d2 sinν 1t

Diode HB
model

+

#

(5.1)

νkm = p1,km ω1 + p2,km ω2 , for integers |p1,km |, |p2,km | ≤ k(5.2)
We take an account the following terms for the input (4);Vd0
Vd1 cosν1 t
Vd2 sinν 1t

id(Vd,k,2m−1 cos νkm t + Vd,k,2m sin νkm t)
m=1

(10)

so that we have id = idM (vd ). Each cross product idk (vd ) ×
vd in (10) is calculated by the method of (7) using MATLAB. Thus, we can execute the Fourier expansion and get
the module of the diode as shown in Fig.1(b). Now, we formulate “Sine-Cosine circuit”[9,10], which corresponds to thedetermining equation of HB.

(4)Then, the diode voltages together with the DC component
Vd0 are described by

(9)

(8)

iE = Is exp( vVBE
) − αR Is exp( vVBC
)

T
T
vBE
)
iC = αF Is exp( VT ) − Is exp( vVBCTiB = (1 − αF )Is exp( vVBE
) + (1 − αR )Is exp( vVBC
)
T
T
(11)
for IS − 10−12 , VT = 0.026, αF = 0.99, αR = 0.3
They consist of the exponential functions, so that HB transistor module can be developed with the diode HB module
shown by Fig.1(b).
2.4 Comparisons between Volterra series and HB
methods We consider a simple base-modulator circuit asshown in Fig.2(a) [13]. Firstly, let us solve it with theVolterra series method.

C
R1
C0
e1 (t)
e (t)

R

Q
R2

L

αI d

C0
E

V1 (jω 1 )Id

V2 (jω 2)

d

g

+
C

R

L Vout
-

2

(a)

(b)

Fig.2(a) Base modulation circuit, (b) 1st order Volterra
equivalent circuit.
R1 = 50[kΩ], R2 = 10[kΩ], R = 0.2[kΩ]C = 1[nF ], C0 = 1[nF ], C0 = 1[nF ], L = 10[µF ], E = 10[V ]v1 = 0.02 sin ω1 t, v2 = 0.02 sin ω1 t, for ω2 = 8.2 × 106 [rad/sec]- 15 -

The transistor is replaced by the Ebers-Moll model, andthe operating point is firstly calculated and we have vd0 =0.5822[V ]. Next, the nonlinear characteristic is approximated by the 3rd order Taylor expansion at the point vd0 as
follow:
id = 0.0130 + 0.52vd + 10.40vd2 + 138.60vd3 .

(12)Note that, for these parameters, we have vBC < 0 andIs exp(vBC /VT ) = 0, so that 1st order Volterra kernel isgiven by Fig.2(b). Next, the 2nd and 3rd Volterra kernels
are calculated according to the method of ref. [3]. The
frequency-domain inter-modulation phenomena as shown in
Fig.2(c) is easily obtained by the application of MATLAB.
[V]

2.5
[V]
2.0

Intermodulation phenomena(Volterra)

1.6

ω1 component

1.2

1.5

ω 2 component

Intermodulation phenomena
ω 2 component

Table 4 Modulation results for R = 1kΩ.
ω1omega2 |ω1 − ω2 | ω1 + ω2 |3rd order 20.33
4.92
7.15
3. 37
7th order 10.38
2.81
9.73
4.58
Exact
8.45
2.21
7.58
3.86

ω1 component

0.8

1.0
|ω1 −ω2 |component

2ω1 component ω1 +ω2 component

0.5
0.0

Note that the results from transient analysis may be mostaccurate, because the transistor is modeled by Ebers-Mollmodel. On the other hand, the Volterra method approximate the Ebers-Moll model in 3rd order polynomial at fixed
operating point. Our HB method approximate it as an unknown operating point decided by HB analysis. Other differences are come from the analytical methods as follows;
For smaller input, the intermodulation phenomena both theVolterra and HB methods are almost same features. However, the amplitudes of the Volterra method become muchlarger than those of HB, especially around the resonant point
as shown in Table 3. Therefore, it seems from the resultsthat the Volterra series method becomes erroneous for thestrong nonlinear circuits having large inputs. Next, we will
investigate the modulations and get the following results:

2ω1 component

0.4|ω1 −ω2 |componentω1 +ω component
2

0

5

10

15

|2ω1 −ω2 |component

20
(a)

25
[x10 rad/s]

30

0
5

0

10

| 2ω1 −ω2 | component

15

20
(b)

25
30
6
[x10 rad/s]

Fig.2(c),(d) The intermodulation phenomena
with Volterra series method and HB method.

Next, we analyze the circuit with HB method. Firstly, thetransistor is modeled by Ebers-Moll model and two diodes
are replaced by diode HB modules as shown in Fig.1(b),where the diodes are approximated by the 3rd order Taylor expansion and it is expanded in the Fourier series with
MATLAB. After then, we formulated the “Sine-Cosine circuits”[9,10] corresponding to the frequency components in(6). Their sub-circuits are coupled with controlled sourcesin each others that corresponds to the determining equationof Spice. Therefore, the frequency-domain intermodulation
can be calculated by the DC analysis of Spice as shown in
Fig.6(d).
Note that although the above two response curves obtained by Volterra series and HB methods are similar forms
in each other, the amplitudes are largely different aroundthe resonant point. We have also calculated the transient
responses, where transistors are modeled by the exponentialfunctions. The ω1 and ω2 frequency-components are calculated with DFT from the steady-state waveforms as shownin Table 3.
16[V]

Although both the frequency response curves have almost
the same phenomena, the peak values to the 3rd order Taylor
approximations are larger than those of 7th order.

3. Illustrative examples
We consider a mixer circuit containing 3 transistors asshown in Fig.3(a) [13]. Note that it works as a modulator
if we set v1 (t) carrier signal, v2 (t) input signal and choose
the LC resonant frequency equal to the carrier frequency.
Now, we analyze the intermodulation phenomena with our
HB method. Firstly, the transistors are modeled by theEbers-Moll model, where
id = Is eλvd , for Is = 10−12 [A], λ = 40, and α = 0.99 (13)and they are approximated by 3rd order Taylor series asfollows:≡ k0 (vd0 )+k1 (vd0 )v+k2 (vd0 )v 2 +k3 (vd0 )v 3 for vd = vd0 +v(14)
In this case, we applied 3rd order Taylor approximation
method. Firstly, all the transistors are replaced by the corresponding HB transistor modules, and RLC circuit is transformed into “Sine-Cosine circuits” [9,10]. Thus, the circuit
is transformed into the 13 HB sub-circuits coupled with controlled sources, which can solved by the DC analysis of Spice.

Output waveform

14
+

12

vout C

10
Q

8
6

V

CC

L

Q2

1

v1 (t)

0

5

10

15

20

25

[µsec]

30
Q3

Fig.2(e) Steady-state wave form of base modulator circuit.

V1 = V2
0.01
0.02
0.03
0.04
0.05

R

-

Table 3 Comparisons of amplitudes(ω1 = 107 [dad/s], ω2 = 8.2 × 106 [rad/s])
Volterra
Harmonic balance
DFT
ω1
ω2
ω1
ω2
ω1
ω2
0.950 0.693 0.929 0.691
0.913
0.698
2.20
1.41
1.72
1.31
1.68
1.29
3.91
2.32
2.50
1.92
2.17
1.72
6.25
3.35
3.35
2.59
2.50
2.11
9.80
4.55
4.31
3.34
2.75
2.38

R2
V1

v 2 (t)V2

Fig.3(a) Mixer circuit.
R2 = 10[kΩ], R = 100[Ω], L = 10[µH],
VCC = 10[V ], V1 = −10[V ], V2 = −5[V ]For the mixer circuit, we had the intermodulation result
as shown in Fig.3(b), and the steady-state waveforms in
Fig.3(c)- 16 -

1
[V]

3.0
[V]

ω 2 component

0.8

Intermodulation phenomena[ω2 =57x106 rad/s]

0.6

1.0
0.0

ω1 +ω2 component[x20]0.4

3. Volterra series method can be only applied to weakly
nonlinear circuits. Our HB method can be applied torelatively strong nonlinear circuits with a few harmonics.

Output waveform

2.0

ω component
1

-1.0

0.2
0

|ω1 −ω2 | component[x20]

3

4

5

6
(b)

7

8

-2.0
-3.0

9
10
ω1 [x107 rad/s]

20

21

22

(c)

24

23

25

[µsec]

Fig.3(b) Frequency response of mixer circuit,
C = 1[pF ], v1 = 0.002 sin ω1 t, v2 = 0.2 sin ω2 t, for ω2 =57 × 106 [rad/s].
(c) Inputs v1 (t) = 0.002 sin 50 × 106 t, v2 (t) = 0.2 sin 57 × 106 t

From DFT analysis of the waveform, we have the following
results:
Vout,1 = 0.7919[V ](HBmethod) V1 = 0.8311[V ](DFT)for component ω1 = 57 × 106 [rad/s]
V0ut,2 = 0.9603[V ](HBmethod) V2 = 0.9585[V ](DFT)for component ω2 = 50 × 106 [rad/s]3[3] P.Wambacq and W.Sansen, Distortion Analysis of Analog
Integrated Circuits, Kluwer Academic Pub., 1998.

2.0
0.0

1

0

| ω1−ω2 | component[x20]

6

8

10

12
(d)

-2.0
-4.0

ω2 component[x20]

4

14

16

18
20
ω1 [x10 6 rad/s][5] K.Kundart “Accurate Fourier analysis for circuit simulators,” CICC94, pp.25-32, 1994.
[6] R.Telichevesky, K.S.Kundart,I El-Fadel and J.K.White,“Fast simulation algorithms for RF circuits,” CICC96, pp.18, 1996.
[8] 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. Jour. of Microwave
and Millimeter-Wave Computer-Aided Eng. vol.1, pp.22-37,
1991.

ω 1 componentω1 +ω2 component[x20][4] R.Telichevesky, K.S.Kundart and J.K.White, “Efficient
steady-state analysis based on matrix-free Krylov-subspacemethods,” ACM, pp.480-485, 1995.

[7] Y.Ueda, The Road to Chaos-II, Aerial Press. Inc., 2001.

4.0

2[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.

[V]

Intermodulation phenomena
ω2=0.2[x106 rad/s]

References

The mixer components are ω1 + ω2 and |ω1 − ω2 | , and the|ω1 − ω2 |-component become zero at ω1 = ω2 . The computation time is 102.05[s].
Next, we analyze the circuit as modulator. In this case, theresonant frequency is equal to 1107 [rad/s]. The inptermodulation frequency phenomena is shown by Fig.3(d) and thesteady-state waveform by Fig.3(e).
4[V]As the future work, we need to apply our HB method tomuch more strong nonlinear elements, and develop of HB
modules of MOSFETs.

10

12

14

(e)

16

18

[µsec]

20

Fig.3(d) Frequency response of modulator.
C = 1[nF ], v1 = 0.001 sin ω1 t, v2 = 0.1 sin ω2 t, for ω2 =0.2 × 106 [rad/s].
(e) Input v1 (t) = 0.001 sin 107 t, v2 (t) = 0.1 sin 0.2 × 106 t

We have the following results:
Vout,1 = 3.987[V ](HBmethod) V1 = 4.218[V ](DFT)for component ω1 = 107 [rad/s]
Vout,2 = 0.307[V ](HBmethod) V2 = 0.337[V ](DFT)for component ω1 = 0.2 × 106 [rad/s]
Vout = 85.5[mV ] at ω1 − ω2 , ω1 − ω2









[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.
[10] A.Ushida, J.Kawata, Y.Yamagami and Y.Nishio, “ Intermodulation analysis of nonlinear circuits using twodimensional Fourier transformation,” NOLTA’05, Budges,
pp.282-285, 2005.
[11] R.L.Geiger, P.E.Allen and N.R.Strader, VLSI: Design Techniques for Analog and Digital Circuits, McGraw-Hill, 1990.
[12] A.S.Sedra and K.C.Smith, Microelectronic Circuits, Oxford
Univ. Press, 2004.
[13] K.K.Clarke and D.T.Hess, Communication Circuits: Analysis and Design, Addison-Wesley Pub. Co., 1971.

The computation time is 209.30[s].

4. Conclusions and remarks
In this paper, we compared Volterra series method with
our Spice-oriented HB method for the analysis of intermodulation phenomena, and we found the following results:
1. Both methods belong to the symbolical analysis, in themeaning that the nonlinear devices must be describedby the power series. Volterra method must be firstly
found the operating point, and use by Taylor expansion. Our HB method also use Taylor expansion insymbolical form. The determining equation is given by
the equivalent circuit so that our HB method is a kindof symbolic method.
2. Both method give good results for small input.

- 17 -