6
IEEE
TRANSACTIONS ON CIRCUITS AND SYSTEMS-I: FUNDAMENTAL THEORY AND APPLICATIONS, VOL.
42,
NO.
I,
JANUARY
1995
Pole-Zero Computation
in
Microwave Circuits
Using Multipoint Pad6 Approximation
Mustafa
Celik,
Ogan
Ocali, Mehmet
A.
Tan,
Senior Member,
IEEE,
and
Abdullah Atalar,
Senior
Member,
IEEE
Abstruct-
A
new method is proposed for dominant pole-
zero
(or
pole-residue) analysis
of
large linear microwave circuits
containing both lumped and distributed elements. The method is
based on a multipoint Pad6 approximation. It finds a reduced-
order rational s-domain transfer function using a data set ob-
tained by solving the circuit at only a few frequency points. We
propose two techniques in order to obtain the coefficients of
the transfer function from the data set. The proposed method
provides a more efficient computation of both transient and
frequency domain responses than conventional simulators and
more accurate results than the techniques based on single-point
PadC approximation such as asymptotic waveform evaluation.
I.
INTRODUCTION
OLE-ZERO computation in a linear circuit is alge-
P
braically equivalent to the computation of the eigenvalues
of
a
circuit matrix. There exist many numerical eigenvalue
algorithms including QR, QZ, deflation-QZ, Muller, MD-QR
algorithms. A detailed comparison of them is given in [l].
Another approach in pole-zero computation is to obtain
a
rational network function in the s-domain and find the poles
and zeros from the polynomials of this rational function by
means of
a
standard root finding algorithm. A review of
symbolic frequency domain network analysis methods can be
found in
[6].
Among these methods the most popular one is
the numerical interpolation method. It is based on polynomial
interpolation with arbitrary selection of frequency points. It
is shown in
[7]
that the best result is obtained when the
interpolation points are uniformly distributed on the unit circle,
which is known
as
interpolation using FFT algorithm. Most
of the symbolic analysis methods, including the interpolation
algorithm, try
to
compute an exact form of the network
functions in the frequency domain.
The methods mentioned above are not practical for mi-
crowave circuits for two reasons:
1)
Usually, practical microwave circuits are of large size,
therefore difficult to analyze. Even
a
simple circuit
may have
a
very large equivalent circuit due to highly
complex device models and parasitic elements that may
be obtained by means of layout extractors.
2)
Circuits which contain distributed elements are infinite-
dimensional systems and have an infinite number of
poles. Therefore, the methods which attempt to find an
exact solution would not be successful.
Manuscript received June
7,
1993; revised March
1,
1994 and September
15,
1994. This work was supported by NATO-SFS project TU-MIMIC. This
paper was recommended by Associate Editor Amedeo Premoli.
The authors are with the Department
of
Electrical and Electronics Engi-
neering, Bilkent University, Bilkent 06533, Ankara, Turkey.
1057-7
122/95$04
One solution to these problems is the dominant pole-zero (or
pole-residue) approximation using the Asymptotic Waveform
Evaluation (AWE)
[2]
method. The AWE technique employs
a
form of PadC approximation to approximate the behavior
of the higher order linear circuit with
a
reduced order model.
The moments of the circuit, which result from
a
Taylor series
expansion of the circuit response about
s
=
0,
are matched
to
a
reduced order rational function. Since the moments
convey information about the low-frequency characteristics
of
the circuit, the AWE technique can only extract the low-
frequency poles. However, for some applications, e.g., the
interconnect circuits, the mid and high frequency ranges
are
more important.
In order
to
improve the accuracy and generality of the
AWE method, many techniques have been proposed. It has
been extended to handle lossy coupled transmission lines
131,
1141, [lo]. In addition to the moments, the Markov parameters,
which are the coefficients of the Taylor series expansion at
s
=
m,
are used to improve the accuracy of the transient response
near
t
=
0
[5]. Moment matching techniques have been refined
in order to obtain accurate and stable low-frequency poles
[11]-[13].
Recently, Chiprout and Nakhla have introduced the complex
frequency hopping (CFH) technique [9] in order to find all
of the dominant poles within the frequency range of interest.
In
this technique,
a
number of single point expansions is
performed at different frequency points. The expansion points
are chosen on the
jw
axis using
a
binary search technique
and then the poles which are considered to be accurate under
some criteria are collected.
This paper proposes
a
new pole-zero (or pole-residue) ap-
proximation technique for the analysis of large linear circuits.
The novelty of this method over AWE based methods is
that the approximation holds for the entire frequency range
under consideration rather than for the low frequencies only.
This method also provides a better approximation than the
previously proposed work which uses multipoint moment
matching methods such
as
the CFH technique,
as
will be
shown to reader in the following sections.
The proposed approach requires the solution of the circuit
matrix at
a
few frequency points. The derivatives of the
network function with respect to
s
are obtained efficiently
from these solutions. By using these derivatives
at
different
complex frequency points,
a
multipoint Pad6 approximation
is used in order to obtain
a
reduced order s-domain network
function. Poles and zeros (or poles and residues) can be found
from this rational network function using standard techniques.
.OO
0
1995 IEEE
CELIK
er
al.:
POLE-ZERO
COMPUTATION
IN
MICROWAVE
CIRCUITS
7
In the next section, the order reduction technique is ex-
plained briefly. Frequency-shifted moments are defined in
Section
111.
We propose two methods for the calculation
of
the coefficients of the network function in Section IV. Practical
and numerical considerations of the proposed method are given
in Section V. In Section VI we present some examples to show
the performance of the proposed method.
11.
ORDER
REDUCTION
Consider a linear system modeled by a coupled set of linear
algebraic equations in Laplace domain,
T(s)x
=
w
(1)
where
T
is the system matrix, the vector
x
is the system
response and the vector
w
is
the
system excitation. In gen-
eral, the system matrix
T
is an arbitrary function of
s.
Let
the system output be any linear combination of the system
response,
H
=
dTx.
(2)
Using Cramer's rule one can obtain
(3)
If the elements of the system matrix
are
polynomials in
s
(e.g., in lumped networks
T
=
T1
+
ST^),
the expansion of
the determinants in
(3)
leads to polynomials in
s,
(4)
In this paper, we consider the circuits containing distributed
components as well as lumped elements. Those circuits can
be regarded as infinite-dimensional systems and the network
functions for an infinite-dimensional system cannot be ex-
pressed as a ratio of two polynomials of finite-degree. Our
aim is to approximate the network function
H
(s)
-
regard-
less of whether it is a rational
or
irrational function of
s,
with a rational function
B(s)
which has approximately the
same frequency characteristics as the original circuit. Let the
approximate function be of the form
bo
+
bls
+
. . .
+
bq-1sq-'
H(s)
=
(5)
1
+
a1s
+
. .
.
+
a&
Since there are
29
parameters to compute in the reduced model,
we need
2q
constraints from the actual circuit. In the AWE
technique
2q
unknowns are calculated by matching the first
T
moments and the first
(2q
-
T)
Markov parameters of the
original circuit to the approximate rational function
[5].
In this work, we propose a method which uses a data
set obtained from the circui! to construct the approximate
s-domain rational function,
H(s).
This data set contains the
frequency-shifed moments
obtained at different complex fre-
quency points. In the following section, we present the eval-
uation of the frequency-shifted moments.
111. FREQUENCY-SHIFTED MOMENTS
The system response
x(s),
can be expanded into a Taylor
series at
s
=
Sk
as:
m
i=O
provided that
x(s)
is analytic at
s
=
Sk.
The coefficient
xki
in
(6)
is called the vector of ith
frequency-shifted moments'
at
s
=
Sk
and
(7)
The first moment vector is the solution of the circuit at
s
=
sk,
X~O
=
T-'(s~)w.
(8)
It can be
shown that the higher order moments can be
evaluated recursively as,
..
r=l
where superscript
(T-)
indicates the rth derivative with respect
to
s.
If the circuit contains only lumped components, then
T(')
=
0
for
T-
>
1.
If it contains distributed elements,
then the derivatives can be found efficiently using either the
eigenvalue moment method [3],[ 141 or the matrix exponential
method
[lo].
are obtained from the moment vectors
xki's
as,
The moments for a particular output of the circuit at
s
=
Sk
.
where
nk
is the number of the moments at
s
=
sk.
We denote
the point
s
=
0
by
SO
and the moments at
s
=
0
are represented
with
moi.
Note that, the frequepcy-shifted moments at
s
=
sz
are complex conjugates of the frequency-shifted moments at
s
=
sk.
We represent the frequency point
SE
by
s-k
and the
moments at
s
=
s-k
by
m-ki.
Let
N
be the total number of
moments, then we have
n
n
k=-n
k=l
where
n
is the number of expansion points in the upper half
s-plane.
IV. MULTIPOINT MOMENT MATCHING
We match the moment set, which contains
N
frequency-
shifted moments obtained from
2n
+
1
points, to a qth order
proper rational function
(q
=
N/2,
which implies
N
should
be even), which is denoted by
[(q
-
l)/q]:
k=-n
,..,
0
,..,
n
We propose two methods for the calculation of the coefficients
of
the rational function. In the first method a set of linear
'
Hereafter,
the
frequency-shifted
moments
will
be
referred
to
as
moments.
8
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS-1: FUNDAMENTAL THEORY AND APPLICATIONS,
VOL.
42,
NO.
1,
JANUARY
1995
equations is obtained for the coefficients and solved directly.
The second method finds the coefficients in
a
recursive manner
without requiring
a
matrix inversion.
matrix equation is formed:
A.
Method
1
In this part we will show that finding the coefficients of
the rational function is equivalent to solving
a
set of linear
equations.
Let us first consider the one-point moment matching case
at
s
=
sk.
If we write the polynomials of the rational function
in terms of
i
=
s
-
Sk
rather than
s,
we obtain:
60
+
b1i
+
...
+
bq-$-l
- -
&I
+
ii1i
+
...
+
iiqiq
mkO
+
mkli
+
...
+
mk(nk-l)ink-l
+
0(i7'n),
(12)
where,
60
=
1
+
E;-
a&
6
t-
-E"-'
l=z
bl(;)s;-z,
(13)
iiz
=
E;='=,
al(Z)S;-Z,
-11
i
=
1,2,
...q
i=0,1,
...
9-1.
Multiplying the denominator polynomial with right hand side
in (12) and equating the coefficients gives
bo
=
mkoao
-
bnb-1
-
mk06~~~-1
+
...
+
mk(n~-l)~o
or in the matrix form,
Mkc
=
mk.
(15)
In
(151,
c
=
[bo...b,-lal
...a,]T
is the coefficient vector
and,
mk
=
[mko
.
.
.
mk(r~~-l)]*
and the
nk
by
N
matrix
Mk
is equivalent to
[C1
i
CZ BC2 BCa].
where
...
I
mkO
B=
[mkl
mkO
:
mk(nk-1) mk(nk-2) mkO
and
C1, Cz
and
C3
are defined
as:
[CI
I
cz
I
Ca]=
Note that,
M-k
=
M;
and
m-k
=
m;,
Finally, collecting
the equations obtained from all points, the following
N
by
N
The solution of this matrix equation yields the coefficients of
the rational transfer function we
are
seeking.
B.
Method
2
Alternatively, the coefficients of the rational function can be
found by means of
a
recursive computation scheme starting
from
a
polynomial which interpolates the given data set.
This method corresponds to the computation of
a
cross-
diagonal sequence in Pad6 table, i.e., it gives all
[m/Z]
Pad6
approximations such that
m
+
1
=
N
-
1.
We denote the rational function that corresponds to the
[(N-l-Z)/Z]
Pad6 approximation by
pl/ql.
In other words, for
a
denominator polynomial the degree is its subscript and for
a
numerator polynomial the degree is the difference between
its subscript and
N
-
1.
Now, let us suppose that the rational
function
pl/ql
interpolates the given data set. That is,
lc
=
-n,
..,
0,
..n
Let us also define the polynomial
g(
s)
of degree
N
as
71
k=-n
Now, assuming
pl(s)
and
q(s)
are co-prime with
g(s),
we
claim that
if
Pl(s)qm(s)
-
Pm(Skl(S)
=
ds)r(s),
(22)
where
~(s)
is
a
polynomial such that
deg[r(s)]
=
max(l
-
m
-
1,
m
-
1
-
l),
then the rational function
p,/q,
also
interpolates the given set. We can prove this claim
as
follows.
Dividing both sides
of
(22) by
ql(s)qm(s)
we obtain,
(23)
Pl(s)
-
P,o
=
g(s)r(s)
ql(s)
4m(s)
41(.5)4m(S).
From the definition of
g(s)
we have
IC
=
-n,
..,0,
..n.
(24)
Therefore,
i
=
0,1,
..,
nk
-
1,
IC
=
-72,
..,0,
..TI,.
(25)
The converse of the claim is also true. That is, if
pl/ql
and
p,
/qm
are two rational polynomials interpolating the moment
set, then, (22) holds.
CELIK
et
al.:
POLE-ZERO
COMPUTATION
IN MICROWAVE
CIRCUITS
9
Now, we can construct our method as follows. Let
pl/ql
and
plPl
/ql-
be two consecutive solutions, then (22) becomes:
Pl(s)ql-l(s)
-
Pl-l(s)ql(s)
=
s(s)ro,
(26)
where
TO
is any real number. Now if we divide the polynomial
pl-l(s)
by
pl(s)
we obtain
Pl-l(S)
=
Pl(S)C(S)
-
4s),
(27)
where,
c(s)
is the quotient polynomial with degree
1
and
d(s)
is the remainder polynomial with degree
N
-
1
-
2. If we add
and subtract the polynomial
pl(
s)ql
(s)c(
s)
to the left hand
side of (26) we obtain
(28)
where
e(s)
=
ql(s)c(s)-ql-I(s)
with degree
1+1.
Now since
d(s)
and
e(s)
have degrees of
N
-
1
-
2
and
1
+
1,
respectively
and since they satisfy (22), we conclude that
d/e
is nothing
but
[(N
-
1
-
2)/(1
+
l)] Pad6 approximation:
d(s)ql(s)
-
Pl(s)e(s)
=
s(s)ro,
Pl+l(S)
=
4s)
(29)
ql+l(s)
=
4s).
(30)
As a summary, given
pl-l/ql-I
and
pl/ql,
we can find
pl+l/ql+1
simply by evaluating one polynomial division, one
polynomial multiplication and one polynomial addition. The
first two approximations
are
found as follows: Let
f(s)
be
a polynomial of degree
N
-
1
interpolating the computed
moment set,
IC
=
-n
,..,
O,..n.
(31)
The polynomial
f(s)
is, in fact,
[(N
-
1)/0] Pad6 approxi-
mation and
to
obtain
[(N
-
2)/1] approximation we rewrite
(26),
Pl(S)
-
f(s)q1(s)
=
g(s)ro,
(32)
which means that the quotient of the division of
g(s)
by
f(s)
gives
q1(s)
and the remainder is
pl(s).
Therefore, having
found the data interpolating polynomial
f
(s),
we can compute
recursively the cross diagonal sequence of the Pad6 table.
V.
PRACTICAL CONSIDERATIONS
In the previous sections, we treated the subject theoret-
ically. Now, we will discuss some topics on the practical
implementation of the proposed method.
A.
Calculation
of
Frequency-Shijled Moments
A recursive scheme for computing the frequency-shifted
moments at a point is given in
(9).
Since the LU factorization
of the circuit matrix
T(sk)
is known from the solution of the
first moment vector, each higher order moment vector can be
obtained by one forward and one back substitution
(FBS)
only.
Totally, we need
n
+
1
LU decompositions of the circuit matrix
which
is
equivalent to obtaining the ac response of the circuit
at dc and
n
points. In addition to LU decompositions we also
need
ni
FBS’s
in order to calculate the
N
moments.
B.
Selection
of
Frequency Points
A crucial step in our method is the selection of frequency
points. Only the poles close to the
jw
axis
are
important in
both time and frequency analyses. Therefore, we choose the
expansion points on the
jw
axis. Once the frequency region
of interest is specified, which is generally between dc and a
maximum frequency, the location and the number of expansion
points can be chosen using the complex frequency hopping
(CFH) algorithm [9]. This algorithm first performs one-point
expansions at the points
s
=
0
and
s
=
jwmaz.
The poles
are
calculated separately from these expansions and if there exists
any common pole, then the search is completed. Otherwise,
more frequencies are selected using a binary search and, an
expansion is performed at every new point until every two suc-
cessive expansions result in at least one common pole. In the
CFH method, the poles about an expansion point
are
calculated
independently from the other expansions, and this conse-
quently decreases the accuracy. In contrast, in our algorithm,
all of the poles and the corresponding residues are obtained
considering all the expansion points simultaneously. This
approach yields a more compact and accurate approximation.
The CFH technique gives the number of accurate poles
about an expansion point. We can choose the number of
moments at each frequency point using this information but not
less than
8
moments at one point. Having chosen the expansion
points and the number of moments at each point, we can find
the coefficients
of
the rational transfer function using one of
the two methods presented in Section
IV.
C.
Finding the Coeflcients Using Method
1
The complex conjugate of every row also exists in (19).
Therefore, this
N
by
N
complex equation set is equivalent to
an
N
by
N
real equation set and can be solved using ordinary
elimination algorithms.
Since our method is proposed for relatively complex cir-
cuits, the orders of approximations
(10
N
50)
are generally
large compared to the orders of typical approximations
(<
12)
seen in the AWE technique. Hence, very large numbers can
appear in the entries of the matrix given in (19) because of
the powers of the expansion points, and consequently, the
matrix can become ill-conditioned. Therefore, we use double
precision arithmetic, and, we also perform a frequency scaling
such that the absolute values of the expansion points,
I
Sk
l’s,
are reduced to around unity.
Another important topic is the stability of the approxi-
mated poles. Similar to the AWE technique, multipoint Pad6
approximations may also result in spurious right hand side
poles. When some unstable poles appear in an approximation,
we discard them and solve a new matrix equation for the
coefficients of the numerator polynomial. The new matrix
equation is obtained as follows. First, let
q’
be the number of
stable poles of the
q
approximated poles. Now, let us rewrite
(19)
in the form,
(33)
or
10
IEEE TRANSACTIONS ON CIRCUITS AND SYSTEMS-I: FUNDAMENTAL THEORY AND APPLICATIONS,
VOL.
42,
NO.
I,
JANUARY
1995
where the vector
b’
contains the coefficients of the numerator
polynomial to be solved, the vector
a’
corresponds to the
coefficients of the stable denominator polynomial and the other
primed quantities are obtained in a similar manner defined
in Section IV, but this time the number of moments are
reduced such that
2q’
=
7th
+
2
E;=,
nk.
Equation
(34)
is
overdetermined and a solution
in
the least square sense can be
obtained. However, this solution does not exactly match the
first
q’
moments. Therefore, we form the primed quantities
in
(34)
by taking the first
nk/2
moments from each point which
results in a
q’
by
q’
matrix equation.
D.
Finding the CoefJicients Using Method
2
In contrast to the first one, this method does not require
any matrix inversion. It is more efficient than the first method.
Moreover, this method allows a search on the Pad6 table in
order to choose stable and accurate approximations. The major
difficulty with this method is the construction of the data
interpolation polynomial,
f(
s).
This polynomial is obtained
using the method of divided differences
[8]
and, generally,
it
is numerically hard to obtain
the
coefficients for higher order
interpolations. Our experiments show that approximations
having orders up to about
20-25
can be found. However, the
first method yields accurate results up to
40-50
poles.
E.
Finding the Poles and the Residues
After having obtained the coefficients of the network func-
tion, to find the poles and residues, a partial fraction de-
composition routine is employed which requires a polynomial
factorization with an associated extra CPU time. This extra
cost becomes important only when the order of approximation
must be increased to about
50.
For this size, this task is about
1.6
million floating point operations. Even in this case, this is
less than the CPU time required for the moment computation.
As the circuit size grows, moment computation gets more
costly, but the cost of the polynomial factorization remains
the same.
The obtained set of poles and residues may be inaccu-
rate due to round-off errors both in the computation of the
coefficients of the network function and in the polynomial
factorization. This set of poles and residues can be verified for
accuracy by using an error criterion defined in the following.
Consider the approximated transfer function:
where the
pj
are the
q
approximate poles, and the
kj
their
corresponding residues. Then, the approximate moments of
this transfer function can be computed as,
IC
=
-n
,..,
O,..n.
(36)
In the absence of round-off errors, we should have
i
=
O,l,
..,nk-l,
riZki
=
mki,
IC
=
-n,
..,O,
..n
(37)
0,
1
exact
Methcd
1
AWE
___
--I
C
-60
-70
-80
Fig.
1.
Frequency response
of
the lowpass filter.
where
mki
are the exact moments obtained from the circuit.
A
normalized error can be defined as
The error increases as the order of approximation increases and
an error beyond a tolerance limit indicates that the approxima-
tion is inaccurate. In this case, a new lower order multipoint
Pad6 approximation should be performed. However, our ex-
periments have revealed that, even for an approximation of
order
50,
this error is less than
1
x
which corresponds
to
a good accuracy.
VI.
EXAMPLES
The following examples demonstrate the performance of the
proposed method on both transient and frequency analyses.
Example
I:
The first example is a lowpass filter imple-
mented with transmission lines. It has a cutoff frequency of
4
GHz. Because of the repetitive properties of the transmission
lines, the response is significant also at higher frequencies.
The actual frequency response of the filter is shown in
Fig.
1.
Applying the enhanced moment matching techniques
[ll],
we found that the best AWE result is the
[6/9]
Pad6
approximation which
is
plotted in Fig.
1.
We have chosen the
maximum frequency as
50
GHz, and, found the expansion
points to be
0,
12.5,
25,
37.5, and
50
GHz using CFH
technique
(so
=
0,
s1
=
j27r12.5
x
lo9,
s-1
=
-j27r12.5
x
lo9,
...,
s4
=
j27r5Ox1O9,
sV4
=
-j27r5Ox1O9).
Solving
(19),
we have found the 47th order approximation
(no
=
14,ni
=
10,
for
i
=
1, -1,
...,
4,
-4)
which resulted in
41
stable
poles. Then solving
(34),
we have obtained the stable
[40/41]
approximation whose frequency response is
also
shown in
Fig.
1.
The step response of the filter is also computed from
the approximations mentioned above. The results are shown
in Fig.
2
together with the HSPICE simulation result for
comparison.
Using this circuit again, we have compared the accuracies of
the CFH technique and our first method. Let
0
and
11
GHz be