2005 International Symposium on Nonlinear
Theory and its Applications (NOLTA2005)
Bruges, Belgium, October 18-21, 2005

Intermodulation Analysis of Nonlinear Circuits Using
Two-Dimensional Fourier Transformation
A.Ushida†, J.Kawata†, Y.Yamagami‡, Y.Nishio‡†Department of Mechanical and Electronic Engineering,
Tokushima Bunri University, Kagawa, 769-2193 JAPAN‡Department of Electrical Electronic Engineering,
Tokushima University, Tokushima, 770-8506 JAPAN
Abstract

In this paper, we propose a new Spice-oriented harmonicbalance method for solving the intermodulation phenomenaof the nonlinear circuits. Generally, mixers and modulator
circuits are driven by multiple inputs, so that the steadystate waveforms behave as quasi-periodic functions. In thiscase, the harmonic balance method needs to apply multidimensional Fourier transformation [14] 1 . On the other
hand, ICs consist of many kinds of nonlinear devices suchas diodes, bipolar transistors and MOSFETs, whose device
models are described by special functions and/or piecewise
continuous functions. Then, using ABMs of Spice [12], wehave developed the Fourier transformers which can be applied to any kinds of models, and constructed the Fouriertransfer modules for these devices, where the relations between the terminal voltages and currents are described byfunctions of the Fourier coeﬃcients. We call the modules“packaged device modules” . Thus, in our harmonic balancemethod, all the nonlinear devices in the circuits are replacedby the corresponding packaged device modules, and the linear sub-circuits are transformed into the corresponding Cosine and Sine circuits [11] composed of the linear resistive
elements and controlled sources. Remark that the equivalentcircuit obtained in the above corresponds to the determining
equations of the harmonic balance method. We will simply
call the circuit Fourier transfer circuit, which can be solvedby the DC analysis of Spice, and we obtain the characteristic
curves such as frequency response curves.
Thus, once these packaged device modules for all kinds ofnonlinear devices are installed in our computers as library,
the harmonic balance method can be easily applied to any
electronic circuits. In this way, we have realized user-friendly
simulators for calculating the intermodulation analysis.
We show the Fourier transfer circuit in section 2, and atechnique for getting the packaged device modules in section
3. The interesting illustrative examples are shown in section
4.

Distortion analysis of nonlinear circuits is very importantfor designing analog integrated circuits and communication
systems. In this paper, we propose an eﬃcient frequency domain algorithm for analyzing the intermodulations of mixers and modulators driven by multiple inputs. Firstly, using
ABMs(Analog Behavior Models) of Spice, we have developed
a Fourier transformer executing the two-dimensional Fouriertransformation. Then, the device modules for nonlinear elements such as diodes, bipolar transistors and MOSFETs
are modeled by the Fourier transformer. Applying thesemodules, we can formulate the determining equations of theharmonic balance method in the form of equivalent circuits.
Thus, they can be solved by the DC analysis of Spice, sothat we can easily obtain the characteristics such as the frequency response curves and intermodulation phenomena. In
our algorithm, we need not derive any troublesome circuit
equation and the transformations into the determining equations, so that our simulators are user-friendly for solving
these circuits.

1. Introduction
The frequency response curves of nonlinear electronic circuits are very useful for the investigation of the global circuit
behaviors. The Volterra series methods have been widelyused for these purposes[1-5]. The kernel functions are givenin the analytical forms, so that we can easily understandthe circuit behaviors from the functions. Although the algorithms are theoretically elegant, it is not so easy to derive
the higher order Volterra kernels, especially, for the large
scale systems containing many nonlinear elements [1]. Remark that since the ideas are based on the bilinear theorem,
they can be only applied to the weakly nonlinear circuits andthe convergences are not guaranteed for strong nonlinear circuits. Furthermore, the characteristics of nonlinear devices
in the Volterra series should be approximated by the polynomial forms, which can be done by the Taylor expansions in
the vicinity at each DC operating point of nonlinear devices[2]. These tasks are not easy to the complicated device models such as the Gummel-Poon model of bipolar transistors
and the higher level models of MOSFETs, especially, in thehigh frequency domain[6,7]. On the other hand, many algorithms have been proposed for calculating the exact steadystate waveforms in the time-domain [8-10]. Unfortunately,
they are not eﬃcient for analyzing the characteristics suchas frequency response curves.

2. Fourier transfer circuit model
Analog integrated circuits are usually composed of manykinds of nonlinear devices such as diodes, bipolar transistors
and MOSFETs, whose Spice models are described by the
special functions containing the exponential, square-root,
1 If the circuit is driven by N independent frequency components, we need to apply an N-dimensional Fourier transformationthat is really time-consuming for large N.

1282ω1 -component, and (2K + 1)-times expansions for ω2 component2 . Thus, the total number of the Fourier expansions is at least (2N + 1) × (2K + 1).
In our algorithm, we will apply a well-known discrete
Fourier transform for each Fourier expansion as follows;

piecewise continuous functions and so on [7]. For these devices, the Fourier coeﬃcients cannot be described in the analytical forms. Therefore, using ABMs of Spice [12], we
realize the equivalent circuits which execute the Fourier expansions in the harmonic balance methods. Now, consider
a nonlinear function described by
i = ĝ(v)

I0 =

(1)

v(τ1 , τ2 ) = V0 +{V2k−1 cos(m1k ω1 τ1 + m2k ω2 τ2 ) bĝ(v)dt =a

for a given bounded N and K. Let us set{V2k−1 cos(m1k ω1 + m2k ω2 )t

k=1+V2k sin(m1k ω1 + m2k ω2 )t}(4)v(θ) = V0 +Kĝ(v) sin kωtdt ⎪
⎪

⎪
⎭

h
(ĝ0 + ĝn ) + h(ĝ1 + ĝ2 + · · · + ĝn−1 )
2

(10)

M


(V2k−1 cos kθ + V2k sin kθ),

θ = ωt(11)

The Fourier transfer circuit model for calculating the N th higher
harmonic component is shown in Fig.1. Applying integral formula(10) to (9), we have

{α2k−1 (τ1 ) cos kω2 τ2 + α2k (τ1 ) sin kω2 τ2 }

k=1

(5)

where α0 (τ1 ), α1 (τ1 ), . . . , α2K (τ1 ) are periodic functions
with the period T1 = 2π/ω1 , and K is the highest harmonic component of τ2 .
Here, if we execute the Fourier transformation of (5) for ﬁxed τ1 atω1 τ1 ={0, ∆θ1 , 2∆θ1 , . . . , 2π} for ∆θ1 = π/N , then we canobtain αk (0),αk (∆θ1 ),αk (2∆θ1 ),. . . . , αk (2π). Thus, αk (τ1 )can be also expanded into a Fourier series;
αk (τ1 ) = αk,0 +0

k=1Hence, the Fourier expansion of (1) for the input (4) can beobtained by the use of the response for the input (2) insteadof (4). The output response i(τ1 , τ2 ) is given byi(τ1 , τ2 ) = α0 (τ1 ) + 2π

where the step size of the integration is h = (a − b)/n. Then,the truncation error is given by ĝ (2) h2 /12n, where (2) shows thesecond derivative. Replacing the integrations in (9) by (10), we
can realize the equivalent circuit model satisfying the relations(9). To understand the circuit model, we assume the input(3)

in (2). Then, it is reduced to
M


0 2π

⎫
⎪
⎪
⎪
⎬

(9)
Now, let us apply a trapezoidal integration formula to (9) as
follows;+V2k sin(m1k ω1 τ1 + m2k ω2 τ2 )}, for |m1k | < N, |m2k | < K(2)v(t) = V0 +ĝ(v)dt

I2k1 =

k=1

t = τ 1 = τ2

 2π

1
1
ĝ(v) cos kωtdt, I2k =π 0
πk = 1, 2, . . . , N, or K

Assume the variable v is consisted of the combination tonesof ω1 , ω2 , and is described by two variables τ1 and τ2 ;namely,
M


1
2π

N


{αn,2k−1 cos nω1 τ1 + αk,2n sin nω1 τ1 }

 2π

n=1

k = 0, 1, 2, . . . , 2K

(6)

0.5^g(V,θ2K)cos Nθ 2K^ ,θ ) sin Nθ
0.5g(VR

2K

Substituting (6) into (5), we obtaini(τ1 , τ2 ) = α0,0 +

+

+

N






k=1

n=1
× cos kω2 τ2N{α2k−1,0 +K{α2k,0 +

k=1



Σ/2Κ
Coefficient of CosNθ

2π[α2k−1,2n−1 cos nω1 τ1 + α2k−1,2n sin nω1 τ1 ]}[α2k,2n−1 cos nω1 τ1 + α2k,2n sin nω1 τ1 ]}n=1
× sin kω2 τ2R

0.5^g(V,θ0)cosNθ 0
0.5^g(V,θ0) sin Nθ 0Σ/2Κ
Coefficient of SinNθV2M

To investigate the accuracy of our Fourier transfer circuit model,
we calculate the following Fourier expansion:M{I2k−1 cos(m1k ω1 + m2k ω2 )t^ ,θ1 )cosNθ 1g(V^ ,θ ) sin Nθ 1
g(V1

Fig.1 Fourier transfer circuit model.

Now, setting t = τ1 = τ2 in (7) and applying some trigonometric arithmetics, we obtain it in the form of

k=1+I2k sin(m1k ω1 + m2k ω2 )t}

R

V0

(7)

i(t) = I0 +2K^ ,θ2K-1 )cosNθ
g(V2K-1^ ,θ )sin Nθ
g(V
2K-1
2K-1

R

{α0,2k−1 cos nω1 τ1 + α0,2n sin nω1 τ1 }

k=n
N

K

(8)

ex cos θ = I0 (x) + I1 (x) cos θ + I2 (x) cos 2θ + · · ·2 The least number comes from the sampling theorem.

(13)

For
getting exact solution, we need the more number of Fourier expansions.

Observe that the above two-dimensional Fourier expansion takes at least (2N + 1)-times Fourier expansions for

⎫

1
1
2
⎪
(ĝ0 + ĝ2K ) +(ĝ1 cos N θ1 ⎪ĝ(v) cos N θdθ =I2N −1 =
⎪
⎪
π 0
K
K
⎪
⎬
+ĝ2 cos N θ2 + · · · + ĝ2K−1 cos N θ2K−1 )
 2π
1
2
⎪
⎪
(ĝ1 sin N θ1 + ĝ2 sin N θ2
I2N =ĝ(v) sin N θdθ =
⎪
⎪
π 0
K
⎪
⎭
+ · · · + ĝ2K−1 sin N θ2K−1 )(12)
The blocks in Fig.1 are constructed by the ABMs of Spice whichcalculates the each term of the relation (12), where the interval[0, 2π] of the integration is divided by 2K equal divisions. In Fig.
1, the value of θk = 2π/2K is obtained by the node voltage at
the kth resistor in the resistive circuit.

2283

whose Fourier coeﬃcients are given by modiﬁed Bessel functions[13] as follows:1
IN (x) =2π

 π

ex cos θ cos N θdθ,

N = 0, 1, 2, . . .

On the other hand, the charge qBE and qBC are given by
Cj0 Anτ IS An
vi
qi =
+
exp
vi , for (vi < Fc φB )[1 − (vi /φB )]mnVt
nVt(17.1)(14)−π

The simulation result for h = 2π/20 in Fig.1 is I1 (10) = 2761
at N = 1, x = 10 which is exactly equal to the value from thetable of Bessel function [13]. Thus, we found that the Fouriertransfer circuit model can get the suﬃciently exact solution evenwith 2K = 10 to 20 divisions of the interval 2π.
Next, we consider the two-dimensional Fourier transformation
for the input waveform (4). The algorithm is as follows;

qi =



Cj0 Anmvi
φB[1 − Fc ](1+m)
τ IS Anvi+
exp) vi , for (vi ≥ Fc φB )(17.2)
nVt
nVt
i = BC, BE
Now, we assume the input voltage waveforms at the collector,
base and emitter vC , vB , vE as follows;

Two-Dimensional Fourier Transformationvi (t) = Vi,0 +Step 1 Set τ1 and τ2 as follows;M{Vi,2k−1 cos(m1k ω1 + m2k ω2 )t

k=1+Vi,2k sin(m1k ω1 + m2k ω2 )t},⎫ω1 τ1 = 0, ∆θ1 , 2∆θ1 , . . . , (2N − 1)∆θ1 , 2π, ⎬for ∆θ1 = π/Nω2 τ2 = 0, ∆θ2 , 2∆θ2 , . . . , (2K − 1)∆θ2 , 2π, ⎭for ∆θ2 = π/K

1 − Fc (1 + m) +i = C, B, E,
Firstly, using two-dimensional Fourier expansion given in section
2, we execute the Fourier expansion to the base and corrector
currents {iB (t), iC (t)} given by (16), and the base-collector andbase-emitter charges {qBC (t), qBE (t)} given by (17). After then,we use the following relations;⎫dqBC (t) ⎬dqBE (t)+iB (t) = iB (t) +
dt
dt
(19)
dqBC (t)⎭iC (t) = iC (t) −dtThus, we have the two-dimensional Fourier expansions of thebase and collector currents as follows;(15)

Step 2 Calculate i(ω1 τ1 , ω2 τ2 ) at the points {i(n∆θ1 , k∆θ2 ), n =0, 1, . . . , 2N, k = 0, 1, . . . , 2K}.
Step 3 For ﬁxed ω1 τ1 = n∆θ1 , n = 0, 1, . . . , 2N , execute theFourier transformation for θ2 . Thus, we get {αk (n∆θ1 ), k =0, 1, 2 . . . , 2K} from the relations (5).
Step 4 Next, for ﬁxed k, αk (τ1 ) is expanded into the Fourier series
using the values {αk (n∆θ1 ), n = 0, 1, 2, . . . , 2N } . Thus, wehave the relations (6) and (7).

ii (t) = Ii,0 +Step 5 Finally, from some trigonometric arithmetics, we obtain thetwo dimensional Fourier expansion (9).

M{Ii,2k−1 cos(m1k ω1 + m2k ω2 )t

k=1+Ii,2k sin(m1k ω1 + m2k ω2 )t},

The ﬂowchart for a three terminal element such as bipolar transistors and MOSFETs is shown in Fig. 2, where the inputs are
given as the Fourier coeﬃcients of 3 terminal voltages, and thecorresponding output current coeﬃcients are obtained from two
terminals.

iD
+
vD
-

C

vBC
B
ib

i1 (τ 1 , τ 2)vα (τ1 , τ 2)
Input
v2

Response ofi 2 (τ1 , τ 2)

Nonlinear Element

Waveformvi (t) , i=1,2,3vβ (τ1 , τ 2)

at (τ 1, τ 2)

q (τ 1 , τ 2 )Transformanionq (τ1 , τ 2)i i(t) , i=1,2

i2

Analog ICs contain many kind of nonlinear devices such as diodes,
bipolar transistors and MOSFETs. In this section, we showa technique to construct the packaged device modules with theFourier transfer circuit given in section 2.These customizations
are very useful to develop user-friendly simulators in our harmonic balance method. Now, we consider npn bipolar transistor
as shown in Fig. 3. The DC currents with Gummel-Poon model
is given as follows[6-7];





exp



vBE
Vt



−1

vCN



QN

vEN

Fig.4 Nonlinear electronic
circuit.

D+

vC1,2kkth vB1,2k-1
Cos-circuit

vE1,2k-1

i D,2k-1

vD,2k-1-

i D,2k-1

vD,2k-1
vC1,2k-1



iCN
iBN

vBN

D





Q1

Network

+

1
vBC
1
+
exp
−1
(16.1)
QB
BR
V
 
 
 t
 
IS
IS
vBE
vBC
exp
−1 +
exp
− 1 (16.2)
iB =
BF
Vt
BR
Vt
where
vBC
vBE
1
1−
−
QB
VAF
VAF
−IS

iB1

where

,
IB,0 = IB,0
+ (m1k ω1 + m2k ω2 )(QBE,2k−1 + QBC,2k−1 )IB,2k−1 = IB,2k−1
(21.1)
IB,2k = IB,2k − (m1k ω1 + m2k ω2 )(QBE,2k + QBC,2k ) ,IC,0 = IC,0− (m1k ω1 + m2k ω2 )QBC,2k−1
IC,2k−1 = IC,2k−1
IC,2k = IC,2k+ (m1k ω1 + m2k ω2 )QBC,2k(21.2)

3. Fourier transformation of nonlinear
devicesISQB

ie

Fig.3 Nonlinear electronic
circuit.

Fig.2 Flowchart of the two-dimensional Fourier transformation
for a three terminal element.

D
iC1

vE1

E

2

iC =

vC1
vB1

LinearvBE

i1

Fourier1

v3

Two Dimensional(20)

i = C, B

ic

v1

(18)

i B1,2k-1

iC1,2k-1
Q1

iBN,2k-1

i CN,2k-1
QN

vCN,2k-1
vBN,2k-1
vEN,2k-1(a)kthvB1,2kSin-Circuit

v E1,2ki B1,2k

i C1,2k
Q1

iBN,2kiCN,2kQN

vCN,2kvBN,2kvEN,2k(b)

Fig.5 Intermodulation circuit of Fig. 4 at m1k ω1 + m2k ω2 .

3284

Remark that “I  s” in (21) show the currents evaluated from(16) and ”Qs” from (17) using two-dimensional Fourier transformation.
In the same way, we can calculate the above input-output twodimensional Fourier transformations to all the nonlinear devices
such diodes and MOSFETs. We also transform all the reactance elements into the linear equivalent Cosine and Sine elements [11]. Thus, the circuit Fig.4 is transformed the equivalent
Cosine-Sine circuit as shown in Fig.5 at each frequency component m1k ω1 +m2k ω2 . It corresponds to the determining equationsof the harmonic balance method . Thus, we can solve the circuitwith DC analysis of Spice, eﬃciently.

Thus, we have ω1 = 0.18[M Hz] and ω2 = 0.2[M Gz] at x = 1.
Note that we take 137 [sec] to get the netlist and 463[sec] for thecomputation with Spice. The mixing ratio is deﬁned by
M = A1 /A2for the two output amplitudes A1 and A2 .

5. Conclusions and remarks
The distortion analysis is very important for designing the highfrequency communication systems.
In this paper, we have proposed an eﬃcient technique for calculating the frequency response curves and mixing ratio of nonlinear electronic circuits. At ﬁrst, the nonlinear devices such as
bipolar transistors and MOSFETs are transformed into the corresponding modules which carry out the two-dimensional Fouriertransformations. Using these device modules, the nonlinear circuit is transformed into the Fourier circuit corresponding to thedetermining equations of the harmonic balance method, and they
are solved with the DC analysis of Spice.
We have found that it will take a lot of computational time forlarge scale circuits containing many transistors when we apply
two-dimensional Fourier transformation to them. Therefore, weneed to improve the algorithm in future.

4. An illustrative example
Let us calculate the frequency response curves of a mixer circuit [13] shown in Fig. 6 which has two independent frequencies
ω1 = 0.18[M Hz] and ω2 = 0.2[M Hz] . Although the circuit may
happen to arise many frequency components, we consider DC, ω1
ω2 , 2ω1 , 2ω2 , ω1 ± ω2 components. We also set N = K = 3 in
(2). The transistor is modeled by a Gummel-Poon model, whoseparameters are given by [6]IS = 10−14 [A]BF = 470
BR = 1

τ = 10−10 [sec]Cj0 = 0.1[pF ]Vt = 0.026[V ]VAF = 150[V ]φB = 0.75
Fc = 0.74

References[1] M.Schetzen, The Volterra and Wiener Theorems of Nonlinear Systems, John Wiley and Sons, 1978.

11[2] D.D.Weiner and J.E.Spina, Sinusoidal Analysis and Modeling of Weakly Nonlinear Circuits with Application to Nonlinear Interference Eﬀects, Van Nostrand Reinhold, 1980.

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

12
V1

V2

L

T1

R1

C

T2

[4] B.J.Leon and D.J.Schaefer, “Volterra series and Picard iteration for nonlinear circuits and systems,” IEEE Trans. on
Circuits and Systems, vol.CAS-25, no.9, pp.789-793, 1078.

9v 1 (t)

T3

8
R2
V3

0

5

10

15

30
25
[x10 -5 sec]

20

v2 (t)

Fig.7 Transient responseof the diﬀerential mixer.

V4

Fig.6 A diﬀerential mixer.

[6] R.L.Geiger, P.E.Allen and N.R.Strader, VLSI: Design Techniques for Analog and Digital Circuits, McGraw-Hill, 1990.
[7] M.Miura, T.Shiino and K.Mori, Circuit-Simulation Technique and MOSFET modeling, Riaraizu Rikou Senta(Japanese), 2003

V1 = V2 = 10[V ], V3 = −10[V ], V4 = 5[V ], L = 10[µH],C = 1[pF ] R1 = 100[Ω], R2 = 10[Ω]
v1 (t) = 0.01 cos 0.18 × 106 t, v2 (t) = 0.2 cos 0.2 × 106 t[8] T.J.Aprille and T.N.Trick, “Steady-state analysis of nonlinear circuits with periodic input,” Proc. IEEE, vol.60, pp.108114, 1972.

3[9] K.S.Kundert, J.K.White and A.Sangiovanni-Vincentelli,
Steady-State Methods for Simulating Analog and Microwave
Circuits, Kluwer Academic Pub., 1990.

xω1 =0.018x10 [MHz]xω2 =0.02x10 [MHz]Mixing ratio[10] A.Ushida, T.Adachi and L.O.Chua, “Steady-state analysis of
Nonlinear Circuits based on hybrid methods,” IEEE Trans.
on Circuits and Systems-I: Fundamental Theory and Applications, vol.39, no.9, pp.649-661, 1992.

2

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

1

ω2[12] MicroSim: PSpice A/D Circuit Analysis User’s Guides, MicroSim Co. 1995.

0
2

2.5

3

3.5

4

x

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

Fig.8 Frequency response curves and mixing ratio of thediﬀerential mixer.

[14] A.Ushida, L.O.chua and T.Sugawara, “A substitution algorithm for solving nonlinear circuits with multi-frequency
components,” Int. Jour. of Circuit Theory and Applications,, vol.15, pp.327-355, 1987.

In the numerical computation in Fig.7, we set two frequencies asfollows;
ω1 = 0.018 × 106 × 10x ,[5] L.O.Chua and C-Y.Ng, “Frequency-domain analysis of nonlinear systems: formulation of transfer functions,” IEE
Trans. Electronic Circuits and Systems, Vol.3, no.6, pp.257269, 1979.

ω2 = 0.02 × 106 × 10x

4285