Spice-Oriented Intermodulation Analysis Combining with
MATLAB
Junji Kawata∗ , Yoshihiro Yamagami† , Yoshifumi Nishio† and Akio Ushida∗∗ Department of Mechanical and Electronic Engineering

Tokushima Bunri University
Email: {kawata,ushida}@fe.bunri-u.ac.jp† Department of Electrical and Electronic Engineering
Tokushima University
Email: {yamagami,nishio}@ee.tokushima-u.ac.jp

Abstract— It is very important to analyze the mixer and modulator
circuits which are driven by multiple input frequencies. Volterra series
methods are widely used for the frequency-domain analysis of nonlinear
circuits, which can be applied to relatively strong nonlinear circuits.
We compare properties of the above two methods using illustrative
examples. In this paper, we propose a Spice-oriented HB method
combining with MATLAB. Firstly, for the circuits containing nonlinear
devices characterized by the special function, they are approximatedby the Taylor series. Then, the Fourier coefficients can be calculated
by MATLAB in the symbolic forms. Thus, the determining equations of
HB method can be formulated by net-list in Spice simulator. It can be
efficiently solved by the DC analysis of Spice, and the frequency-domain
solutions are obtained. We found from examples that, although Volterra
series method can be efficiently applied to weakly nonlinear circuits, it
becomes erroneous to strong nonlinear circuits.

I. I NTRODUCTION
Frequency-domain analysis of nonlinear electronic circuits driven
by multiple inputs is very important for designing integrated circuits
and communication systems. Volterra series methods are widely used
for the analysis [1-3] because they give the solutions in analyticalforms. The algorithms are based on bilinear theorem [1], and they
can be effectively applied to the intermodulation analysis, where
nonlinear elements must be described by the polynomial functions.
HB (harmonic balance) method is also well-known in thefrequency-domain analysis, which will give good results even for
relatively strong nonlinear circuits. The classical HB methods [46] can be applied to relatively small scale circuits and find many
interesting nonlinear phenomena such as bifurcations and chaos. Note
that, in order to obtain the Fourier coefficients, nonlinear elements
must be also described by polynomial functions in HB method. Itcomplicate the applications of HB method for large scale circuits.
For example, if the circuit equation is composed of “n” functionsand “K” frequency components are considered in the analysis, thenthe determining equation consist of a number of “n(2K+1)” functions.
Therefore, we propose here another type of Spice-oriented HB
method combining with MATLAB. Firstly, bipolar transistors and
MOSFET modeled by the special functions are approximated by the
Taylor expansions in the polynomial forms. Thus, in this case, theFourier coefficients can be described in the symbolic forms if we use
MATLAB. It is further possible to formulate the determining equationof HB method in the form of net-list of Spice. It can be easily
solved by DC analysis of Spice, and we can get the frequency-domain
solutions, efficiently. Thus, we need not any troublesome tasks suchas formulation of the circuit equations and many transformations
required to get the determining equations of HB method.
Using an example, we will compare the Volterra series methodwith our HB method in section 2. In section 3, we show interesting
examples of mixer circuit.

1-4244-1474-1/07/$25.00 ©2007 IEEE.

II. C OMPARISON BETWEEN VOLTERRA SERIES AND HBMETHODS

A. Taylor Approximation of Nonlinear Devices
Transistors are usually modeled by Ebers-Moll model consisted of
two diodes, where the nonlinear diode (id , vd ) is described by
id = Is exp vd /VT .

(1)

where Is = 10−12 , VT = 0.026. Firstly, we estimate the Fourier
coefficients DC, cos ωt, cos 2ωt, cos 3ωt driven by
vd = Vd0 + Vd cos ωt,for Vd0 = 0.6,

(2)

and compare them with those obtained from 3rd order approximation
with Taylor expansion. Then, its approximation is described asfollows:
id = k0 + k1 vd + k2 vd2 + k3 vd3 ,(3)Then, the Fourier coefficients of DC, cos ωt, cos 2ωt, cos 3ωt are
shown in Fig.1, where the curves are the coefficients obtained fromthe 3rd order approximation given by (3), and the bar graphs arethose of the exponential function (1), where horizontal axis is Vd in
(2).
0.6
Id

Harmonic components

0.5
0.4

:DC component:Fundamental:2nd harmonic:3rd harmonic

0.3
0.2
0.1
0

0.02

0

0.05

0.07

Vd 0.10

Fig.1 Fourier coefficients obtained from (1) and (3).

Observe that when we apply 3rd order approximation of Taylor
expansion to (2), the errors increase rapidly beyond Vd ≈ 0.05.
B. Volterra Series Method
In order to understand the Volterra series method driven by 2 inputs[1-3], we consider a simple circuit shown in Fig.2, where it has one
nonlinear resistive element given by (3).

1107

iG
vin1(ω t)Linearvin2 (ω 2 t)

Network

1

+
vG
-

v(t)

Fig.2 Nonlinear circuit driven by two inputs.

Now, let us derive the Volterra kernels in the complex frequency
domain “s”[3]. The 1st order kernel corresponds to the output
voltages at the resistor “r” for the unit inputs U1 (s1 ) and U2 (s1 )as shown in the linear circuits of Fig.3(a). Thus, 1st order kernel isgiven byH1k (s1 ) = H1 (s1 ),(4)
for two inputs. The nonlinear current source for 2nd order is given
by [3]
I2k (s1 , s2 ) = k2 H1k (s1 )H1k (s2 ).
(5)

C. Spice-Oriented HB Method
We propose here a Spice-oriented HB (harmonic balance) method
combining with MATLAB, such that the determining equation is
given by the equivalent circuit and/or net-list of Spice. To use
MATLAB for the Fourier expansion of nonlinear devices, we need todescribe the nonlinear characteristics in the polynomial forms. The
schematic diagram of diode is shown in Fig.4(a), where it is driven
by the two inputs
v1 (t) = V1 cos ω1 t,

v2 (t) = V2 cos ω2 t.

(12)

Let the admittance from the output terminal be Y (s). Then, the 2ndorder kernel H2k (s1 , s2 ) is derived byFirstly, it is approximated by 3rd order Taylor expansion as follows:

Y (s1 + s2 )H2k (s1 , s2 ) = −I2k (s1 , s2 ).

with an unknown operating point Vd0 . Now, we consider 13 intermodulation frequency components:(6)

id = k0 (Vd0 ) + k1 (Vd0 )vd + k2 (Vd0 )vd2 + k3 (Vd0 )vd3

In the same way, the nonlinear current source for the 3rd order isgiven byI3k (s1 , s2 , s3 ) = k3 H1k (s1 )H1k (s2 )H3k (s3 )+2
k (H1 (s1 )H2 (s2 , s3 ) + H1 (s2 )H2 (s1 , s3 ) + H1 (s3 )H2 (s1 , s2 ))
3 2(7)

U (s1)

Linearr
Network+H (s1 )
-11

Linearr
Network

DC, ω1 , ω2 , 2ω1 , 2ω2 , 3ω1 , 3ω2
ω1 + ω2 , |ω1 − ω2 |, 2ω1 + ω2 , |2ω1 − ω2 |
ω1 + 2ω2 , |ω1 − 2ω2 |r
U (s1)

Network+
H12(s1 )-Linear+H3k(s ,s ,s )- 1 2 3I (s ,s s3 ) (c) 3rd order

NetworkKX

vd (t) = Vd0 ++H (s ,s )
1 2
-2k
(b)
2nd order
I (s1 ,s2 )

r

.

(14)




(Vd,2k−1 cos νk t + Vd,2k sin νk t) 



k=1

i(t) = Id0 +

K
X



{Id,2k−1 cos νk t + Id,2k sin νk t) 



,

(15)

k=1

Substituting the first equation of (15) into (13), and neglecting theterms except for (14), we have second relation of (15).

3k 1 2(a) 1st order

Vd0
Vd1 cosν1 t
Vd2 sinν 1t

Fig.3 Computations of Volterra kernels.

idand the 3rd order kernel H3k (s1 , s2 , s3 ) is derived by
Y (s1 + s2 + s3 )H3k (s1 , s2 , s3 ) = −I3k (s1 , s2 , s3 ).

(8)

The output voltage for the following 2 inputs
vin2 (t) = A2 sin ω2 t

(11)

are calculated by setting { si = {±ω1 , ±ω2 }, i = 1, 2 }in (9). Thus,
the output amplitudes are given in the Table 1 [3]. Note that thesolutions from the Volterra series method can be easily solved by
MATLAB.
Table 1 The output voltage in complex domain [3]Frequency
Amplitude of responseDC| 12 A21 H2k (jω1 , −jω1 ) + 21 A22 H2k (jω2 , −jω2 )|ω1|A1 H1k (jω1 ) + 23 A1 A22 H3k (ω1 , jω2 , −jω2 )+ 43 A31 H3k (ω1 , jω1 , −jω1 )|ω2|A2 H1k (jω2 ) + 23 A21 A2 H3k (ω1 , −jω1 , jω2 )+ 43 A32 H3k (ω2 , jω2 , −jω2 )|ω1 + ω2|A1 A2 H2k (jω1 , jω2 )||ω1 − ω2 ||A1 A2 H2k (jω1 , −jω2 )|2ω1| 12 A21 H2k (jω1 , jω1 )|2ω2| 12 A22 H2k (jω2 , jω2 )|3 2
2ω1 + ω2
A A H (ω1 , jω1 , jω2 )|4 1 2 3k
3 2
|2ω1 − ω2 |A A H (ω1 , jω1 , −jω2 )|4 1 2 3k3
ω1 + 2ω2
A A2 H (ω1 , jω2 , jω2 )|4 1 2 3k3
|ω1 − 2ω2 |A A2 H (ω1 , −jω2 , −jω2 )|4 1 2 3k3ω1
| 14 A31 H3k (ω1 , jω1 , jω1 )|3ω2| 14 A32 H3k (ω2 , jω2 , jω2 )|vd(a)Hout (s1 , s2 , s3 ) = H1k (s1 ) + H2k (s1 , s2 ) + H3k (s1 , s2 , s3 ). (9)Diode HB

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

model

+

Thus, the output in the complex frequency domain is given byvin1 (t) = A1 sin ω1 t,)

Thus, the input voltage and output current waveforms are assumed
as follows:

2k

Linear(13)-

I d2K sinνKt

Vd2K sinνKt(b)

Fig.4 (a) Diode model, (b) HB module with MATLAB.

Note that the number of Fourier coefficients arisen by the above
arithmetic calculations is
1 + 27 + 272 + 273 = 20440.

(16)

We choose only the frequency components corresponding to frequencies (14) using MATLAB, and formulate the net-list in the Spice asshown in Fig.4(b). Thus, we can formulate “Sine-Cosine circuit”[7,8],
which corresponds to the determining equation of HB. The circuit
is resistive, and it can be easily solved with DC analysis of Spice.
Thus, we can get the frequency-domain solutions such as frequency
response curves and so on.
It is convenient to develop the HB modules for bipolar transistors.
In this case, if we model the transistor with Ebers-Moll model, theemitter, base and corrector currents are given as follows:

iE = Is exp( vVBE
) − αR Is exp( vVBC
)

T
T
vBC
iC = αF Is exp( vVBE
)
−
I
)
(17)
s exp( VT
T
vBC
iB = (1 − αF )Is exp( vVBE
)
+
(1
−
α
)
R )Is exp( VT
T
for IS = 10−12 , VT = 0.026, αF = 0.99, αR = 0.3.
They have some exponential functions, so that HB transistor module
can be developed by the use of the diode HB module shown by
Fig.4(b).

1108

16
[V]

D. Comparisons between Volterra Series and HB Methods

Output waveform

14

We consider a simple base-modulator circuit as shown in Fig.5(a)[11]. Firstly, let us solve it with the Volterra series method mentioned
in subsection 2.1.

12
10
8
6

C

R

L

Q

E

e (t)
1

R2

e (t)

αI d

C0

R1
C0

V1 (jω 1 )Id

V2 (jω 2)

g

d

0

R

L Vout
-

2

(a)

10

15

20

25

[µsec]

30

Fig.5(e) Steady-state waveform of base modulator circuit.

+
C

5

(b)

Fig.5(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]

Note that although the response curves obtained by Volterra series
and HB methods are similar, the amplitudes are largely different
around the resonant point. We have also calculated the transient
responses, where transistors are modeled by (17), instead of 3rd order
Taylor approximation in Volterra and HB methods. The ω1 and ω2
frequency-components are calculated with DFT from the steady-state
waveforms and shown in Table 2.
Table 2 Comparisons of amplitudes.
(ω1 = 107 [rad/s], ω2 = 8.2 × 106 [rad/s])
Volterra
HB
DFT
V1 = V2 ω1
ω2
ω1
ω2
ω1
ω2
0.01
0.950 0.693 0.929 0.691 0.913 0.698
0.02
2.20
1.41
1.72
1.31
1.68
1.29
0.03
3.91
2.32
2.50
1.92
2.17
1.72
0.04
6.25
3.35
3.35
2.59
2.50
2.11
0.05
9.80
4.55
4.31
3.34
2.75
2.38

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

(18)

Thus, we have
k1 = 0.52, k2 = 10.40, k3 = 138.60(19)in (3). Note that we have vBC < 0 and Is exp(vBC /VT ) = 0 forthese parameters, so that 1st order Volterra kernel is given by Fig.5(b).
The 2nd and 3rd Volterra kernels are calculated by methods shown insubsection 2.2. The frequency-domain intermodulation phenomena as
shown in Fig.5(c) is easily obtained by the application of MATLAB
to the Volterra kernels.

Note that the results obtained from transient analysis will be mostaccurate as the results of the base modulation circuit, because thetransistor is only modeled by Ebers-Moll model. On the other
hand, the Volterra method approximates the Ebers-Moll model atfixed operating point. Our HB method approximates it at unknown
operating point decided by HB analysis. Other differences are comefrom the analytical methods. For smaller input, the intermodulation
phenomena obtained by Volterra and HB methods have almost samefeatures. However, the amplitudes from the Volterra method are muchlarger than the others, especially around the resonant point. Therefore,
it seems from the results in Table 2 that the Volterra series method
becomes erroneous for strong nonlinear circuits with large inputs.

[V]

2.5
[V]

Intermodulation phenomena(Volterra)

2.0

ω1 component

III. I LLUSTRATIVE EXAMPLES

ω 2 component

1.2

ω 2 component

1.5

Intermodulation phenomena

1.6

ω1 component

0.8

1.0
|ω1 −ω2 |component

2ω1 component

2ω1 component ω1 +ω2 component

0.5

0.4

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

0

0.0
0

5

10

15

|2ω1 −ω2 |component

20
(c)

25
6
[x10 rad/s]

30

0

5

10

| 2ω1 −ω2 | component

15

20
(d)

25
30
6
[x10 rad/s]

We consider a mixer circuit containing 3 transistors as shown in
Fig.6(a) [11]. Note that it becomes 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 modeledby the Ebers-Moll model, where
id = Is eλvd , for Is = 10−12 [A], λ = 40, and α = 0.99

Fig.5(c),(d) The intermodulation phenomena
obtained by Volterra series method and HB method.

(20)

and they are approximated by Taylor series as follows:
Next, we analyze the circuit with HB method. Firstly, the transistor
is modeled by Ebers-Moll model and two diodes are replaced by HB
diode modules shown in Fig.4(b), where the diodes are approximatedby the 3rd order Taylor expansion. And then, we formulated the“Sine-Cosine circuits”[7,8] corresponding to the frequency components in (14). These sub-circuits are coupled each other with
controlled sources which correspond to the determining equationof Spice. Therefore, the frequency-domain intermodulation can becalculated by the DC analysis of Spice as shown in Fig.5(d).

≡ k0 (vd0 ) + k1 (vd0 )v + k2 (vd0 )v 2 + k3 (vd0 )v 3

for vd = vd0 + v(21)
In this case, we consider the 13 intermodulation frequencies givenby (14). Firstly, all the transistors are replaced by the corresponding
HB transistor modules, and RLC circuit is transformed into “SineCosine circuits” [7-8]. Thus, the circuit is transformed into the 13
sub-circuits coupled with controlled sources, it can be solved by the
DC analysis of Spice.

1109

Note that we can get good result from our HB method.

+
v

V
out C

R

CC

L

IV. C ONCLUSIONS AND REMARKS

Q

Q

1

2

v1 (t)Q3v 2 (t)

R2
V

V

1

2

Fig.6(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 phenomena as
shown in Fig.6(b) and the steady-state waveforms of Fig.6(c)
1
[V]

13
[V]

ω 2 component

Output waveform

12

0.8

ω component
1

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

0.6

10
ω1 +ω2 component[x20]0.49|ω1 −ω2 | component[x20]

0.2

8

0
3

4

5

6
(b)

7

8

R EFERENCES

7910
ω1 [x107 rad/s]

20

21

22

(c)

23

24

[µsec]

25

Fig.6(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) Steady-state waveform for ω1 = 50 × 106 .

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]
Vout,2 = 0.9603[V](HBmethod) V2 = 0.9585[V](DFT)for component ω2 = 50 × 106 [rad/s]







The intermodulation components of ω1 + ω2 and |ω1 − ω2 | are smallenough, and the |ω1 − ω2 |-component becomes zero at ω1 = ω2 . The
computation time is 102.05[s].
Next, we analyze the circuit as modulator. In this case, the resonant
frequency is equal to 107 [rad/s]. Figure 6 (d) and (e) shows the intermodulation phenomena and the steady-state waveform, respectively.
4[V]

Intermodulation phenomena6ω2=0.2[x10 rad/s]

3

ω 1 component

15[V]

Output waveform

13
11

2

9
ω1 +ω2 component[x20]1

In this paper, we compared Volterra series method with our Spiceoriented HB method for the analysis of intermodulation phenomena,
and we found the following results:
1) Both methods belong to the symbolical analysis, where thenonlinear devices must be described by the power series. In
this case, Volterra method must found the operating point usedin Taylor expansion. Our HB method need only to describe it
as Taylor expansion with unknown operating point which is
decided by HB analysis.
2) Volterra series method finds the solution by solving the kernels
with MATLAB, efficiently.
3) Our HB method transforms the circuit into “Sine-Cosine circuit” corresponding to the determining equation of HB method,where we use HB modules. The circuit is solved by DC analysis
of Spice.
4) Volterra series method can be only applied to weakly nonlinear
circuits. Our HB method can be applied to relatively strong
nonlinear circuits with a few harmonics.
As future works, we need to apply our HB method to strong nonlinear
elements, and develop of HB modules of MOSFETs.

[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 Modelingfor RF and Microwave Design, Artech House, 2005.
[3] P.Wambacq and W.Sansen, Distortion Analysis of Analog Integrated
Circuits, Kluwer Academic Pub., 1998.
[4] Y.Ueda, The Road to Chaos-II, Aerial Press. Inc., 2001.
[5] R.Telichevesky, K.S.Kundart and J.K.White, “Efficient steady-state analysis based on matrix-free Krylov-subspace methods,” ACM, pp.480-485,
1995.
[6] R.J.Gilmore and M.B.Steer, “Nonlinear circuit analysis using the methodof 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.
[7] A.Ushida, Y.Yamagami and Y.Nishio, “Frequency responses of nonlinear
networks using curve tracing algorithm,”ISCAS 2002, vol.I, pp.641-644,
2002.
[8] A.Ushida, J.Kawata, Y.Yamagami and Y.Nishio, “ Intermodulation analysis of nonlinear circuits using two-dimensional Fourier transformation,”
NOLTA’05, Budges, pp.282-285, 2005.
[9] R.L.Geiger, P.E.Allen and N.R.Strader, VLSI: Design Techniques for
Analog and Digital Circuits, McGraw-Hill, 1990.
[10] A.S.Sedra and K.C.Smith, Microelectronic Circuits, Oxford Univ. Press,
2004.
[11] K.K.Clarke and D.T.Hess, Communication Circuits: Analysis and Design, Addison-Wesley Pub. Co., 1971.

| ω1−ω2 | component[x20]7ω2 component[x20]

0
4

6

8

10
(d)

12

14

16

18
20
ω1 [x10 6 rad/s]

5
10

12

14
(e)

16

18

20
[µsec]

Fig.6(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) Steady-state waveform for ω1 = 107 .

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











The computation time is 209.30[s].

1110