- Email: [email protected]
Singular root distribution problem for delta-operator based real polynomials
Automatica 35 (1999) 791}807
Singular root distribution problem for delta-operator based real polynomials H. Fan* Department of Electrical & Computer Engineering and Computer Science, University of Cincinnati, Cincinnati, OH 45221-0030, USA Received 14 November 1996; revised 17 October 1997; received in "nal form 5 November 1998
Delta operator based Schur}Cohn tests for root distribution of a polynomial in the singular cases I and II are developed. These singular case tests connect and unify the traditional discrete time Schur}Cohn tests and their continuous-time counterparts in the singular cases.
Abstract E$cient zero location tests for the delta-operator-based polynomials have been a recent topic of interest. These tests not only provide the needed e$cient tests for the delta-operator formulated systems, control, and signal processing problems, which are numerically better behaved than the traditional shift-operator based tests for fast sampling, they also provide missing links between classical tests in the z-domain and the s-domain such as the Schur}Cohn test (SCT) and the Routh test (RT). These will also be accomplished in this paper, but for the singular cases. Speci"cally, the recently developed singular case I Schur}Cohn test by Pal and Kailath (1994) is modi"ed and transformed into the delta-domain. As the sampling interval vanishes this new algorithm converges to an s-domain singular case I test associated with the recently developed s-domain SCT. Singular case II is also considered and the di!erence between the traditional approaches in the z-domain and the s-domain reconciled. Numerical examples are given to illustrate better numerical properties of the new test over Pal-Kailath recursion for fast sampling. 1999 Elsevier Science Ltd. All rights reserved. Keywords: Stability tests; Inertia tests; Delta operators; Singularities
1. Introduction The root distribution problem, or in other words the zero location tests, has been studied for over 100 years. Classical e$cient tests include the Schur}Cohn test (SCT) in the discrete-time shift operator case, denoted the z-domain, and the Routh test (RT) in the continuoustime case, denoted the s-domain (Jury, 1982). Recently, the d operator, i.e., d"(q!1)/* where qx(t)"x(t#1) 22222 * Corresponding author. Tel.: (513) 556-4765; fax: (513) 556-7326; e-mail: [email protected]. This paper was not presented at any IFAC meeting. This paper was recommended for publication in revised form by Associate Editor C. V. Hollot under the direction of Editor R. Tempo. This work is supported by the O$ce of Naval Research under Grant N00014-96-1-0241, and by the National Science Foundation under Grant INT-9321821.
and * is the sampling interval, has been used to develop such tests for the d-operator based discrete-time systems (Fan, 1997a, b). These new tests shall be referred to as in the d-domain or in terms of its frequency variable the c-domain. The new c-domain tests are computationally e$cient with O(N) computations where N is the system order. They are numerically better behaved than the traditional shift-operator based tests such as the SCT for fast sampling. In addition, they provide smooth transitions between the z-domain tests and the s-domain ones. For example, the limit (as *P0) of the d-SCT gives an s-domain SCT, called p-SCT (Fan, 1997a). However, these tests are all developed for the strongly regular case. They all fail in the singular cases. First let us review the strongly regular case of SCT, d-SCT, and p-SCT (Fan, 1997a) and show why they would fail in the singular cases. The formulation of this section will also be useful for the subsequent
0005-1098/99/$*see front matter 1999 Elsevier Science Ltd. All rights reserved PII: S 0 0 0 5 - 1 0 9 8 ( 9 8 ) 0 0 2 1 8 - 0
792
H. Fan/Automatica 35 (1999) 791}807
development. Consider an Nth-order real polynomial in the z-domain (assuming a O0): f (z)Oa z,#a z,\#2#a z#a . (1.1) ,\ , Using the vectors z O[z,, z,\, 2, z, 1]2 ,
Using the relationship c"(z!1)/* (Goodwin et al., 1992), it has been shown (Vijayan et al., 1991) that c "¹2 z where ¹ is a lower-triangular matrix whose , , , , l, kth entry is given by
(!1)J\I N!k (¹ ) " , 04l, k4N, l5k. , J I *,\I l!k
(1.10)
(1.2) aO[a , a , 2, a , a ]2 ,\ , Eq. (1.1) can be written as f (z)"z2 a. Here and through, out this paper, boldface letters will denote vectors whereas unboldfaced letters denote scalars. Furthermore, de"ne
To keep the identity f (z)"z2 a"*,c2 bO*,F(c), one , , needs to de"ne a"*,¹ b as in Vijayan et al. (1991). , C De"ne b "[b b 2b ]2 and b "[bC bC 2bC ]2, L L L LL L L L LL and
a OJ a, , J Oantidiag.+121,, (1.3) , and the lower-order parameter vectors a "[a a 2 L L L a ]2 and a "[a 2a a ]2 for all n"0, 1, 2, N!1. LL L LL L L Then the SCT is given by
Left multiply *\,¹\ and *\L\¹\ to the SCT, and , L> after some algebra, one obtains
a
b "*\L¹\ a , bC"*\L¹\ a . L L L L L L
0
b
(1.4)
a L "a #! a , n"N!2, N!3, 2, 0, L> L> L> 0 (1.5)
#*
,\ 0 0
,\ "a#! a , , 0
(1.11)
b ,\
b ,\ "b#! bC, , 0 0
(1.12)
#* b "! ,\ 0
b ,\ 0
(1#! ! ) (! #! ) C ,\ , b# ,\ , b # * *
(1.13)
and for n4N!3, a a ! "! LL , ! "! , , L , a a L
(1.6)
n>"var 1, 1!! , (1!! ) (1!! ), 2, , , ,\ , (1!!)G G G
(1.7)
or n\"C of entries in +a , ,a , ,\ 2 having the same sign as a , (1.8) where n> is the number of outside unit circle roots and n\ that of inside ones, and var+x, denotes the variation of signs in the sequence +x,. The singular cases occur when a "0 for some n. From Eq. (1.6) it is clear that L when this happens ! cannot be computed, therefore the L SCT cannot be continued. The same situations occur with the d-SCT of Fan (1997a). We "rst brie#y review the d-SCT. It turns out that the d-SCT can be easily obtained by a linear transformation of a scaled version of the SCT in Eqs. (1.4)}(1.8). First de"ne the d-operator quantities. c O[c,, c,\, 2, c, 1]2, , bO[b , b , 2, b , b ]2. ,\ ,
(1.9)
0
0
1 ! 0 #* b " 1# L> L * ! L> b 0 L
0
b L>
b ! ! b # L> L> ! L> L> b , L> ! 0 ! b L> L> L> where
(1.14)
1 !" [(!1)L>, (!1)L*, (!1)L\*, 2, !*L]b , L b L L (1.15) bC"*\,¹\ a "¹\J ¹ b. , , , , See Fan (1997a) for details. The inertia test is
(1.16)
b b b ,\ , ,\ ,\ , 2, n> b "var 1, b b b ,\ , b L L\ , (1.17) b L L n\"C of entries in +b , ,b , ,\ 2 having the same sign as b . (1.18) The singular cases of the d-SCT occur when b "0 for L some n. Similar to the SCT, from Eq. (1.15), it is seen that
H. Fan/Automatica 35 (1999) 791}807
in the singular cases ! cannot be computed, thus the L d-SCT cannot be continued. The following example illustrates one singular case for both the SCT and the d-SCT. Example 1. f (z)"z!4.999z#10.001z!10.002z#5.001z!1, (1.19) F(c)"c#0.004c#0.08c#0.448c#1.024c #1.024, *"0.25.
(1.20)
These two characteristic functions correspond to the same system/process but formulated using the shiftoperator and the d-operator, respectively. Now perform the SCT on f (z) and the d-SCT on F(c). SC¹: a"[1 !4.999 10.001 !10.002 5.001 !1]2, ! "1, a "[0 0.002 !0.001 !0.001 0.002]2, ! "?
a "0, d-SC¹:
b"[1 0.004 0.08 0.448 1.024 1.024]2, bC"[!1 0.004 0.032 0.064 0.256 1.024]2, ! "1, b "[0 0.032 0.32 0.768]2, b "0,
! "?.
The same situation may occur to the p-SCT in the s-domain as well. The p-SCT in the s-domain can be obtained by simply taking the limit *P0 of the d-SCT. By doing so ! P(!1)L> and Eq. (1.12) becomes L b Pn "[2m , 0, 2m , 0, 2]2, ,\ ,\
(1.21)
where n"[m , m , 2, m ]2 is the continuous-time sys , tem parameter vector which is the limit of b as *P0. Furthermore, if one de"nes n "[2m , 0, 2m , 0, 2]2, ,
(1.22)
then Eqs. (1.13) and (1.14) give the same limit 0
0
m n 0 P 0 " L> n ! L> L> m 0 L> b n L L
(1.23)
793
for n4N!2. Therefore, Eqs. (1.21)}(1.23) constitute the p-SCT. The inertia test is Eqs. (1.17) and (1.18) with b's replaced by m's. Note that m "0 for j odd and for all L H n. It is also clear from Eq. (1.23) that singularity occurs when m "0 for some n, for which case Eq. (1.23) cannot L be continued at stage n!2. Note that the above kind of singularities is not to be confused with the numerical singularity of the SCT due to fast sampling. To be speci"c, assume that the discretetime characteristic polynomial is obtained by sampling a continuous-time system or process whose characteristic polynomial is strongly regular in the sense of the p-SCT. When sampling is fast (*P0) the SCT gives "! "P1 and L a P0 for all n (Fan, 1997a). This is the numerical L singularity of the SCT due to fast sampling, which is now overcome by the d-SCT of Fan (1997a). There are two essential steps in the development of the d-SCT: (1) Scaling by 1/* and 1/* to eliminate such singularity; and (2) transformation to provide a meaningful limit as *P0. In other words, numerical singularity of the SCT due to fast sampling can be eliminated by scaling, i.e., if scaled properly as in Eqs. (21)}(25) of Fan (1997a), a as L obtained by Eqs. (1.4) and (1.5) do not vanish as *P0 when the underlying s-domain characteristic polynomial is strongly regular in the sense of p-SCT. Therefore, the kind of singularities that will be dealt within this paper is such that a "0 regardless of the L speed of sampling (without taking the limit *P0), and can be called &&algorithmic'', as opposed to &&numerical''. Speci"cally, two types of algorithmic singularities will be considered. In type I (or case I) a "a "2"a "0, L L LH b "b "2"b "0, and m "m "2"m "0 L L LH L L LH for some n and some j such that 0(j(n. Clearly Example 1 yields type I singularity. Type II (or case II) is type I with j"n, i.e., a "0, b "0, and n "0 for some L L L n. The strongly regular SCT, d-SCT, and p-SCT of Fan (1997a) as summarized above all fail in case of either Type I or Type II algorithmic singularities. In this paper the algorithmic singular cases will be investigated in the framework of the SCT, d-SCT, and p-SCT. The next three sections will be devoted to algorithmic singularity type I. In Section 2 an order recursive solution for the SCT due to Pal and Kailath (1994) will be scaled and modi"ed. It will then be transformed into the d-domain in Section 3, resulting in a new singular case I d-SCT. Its limit as *P0 naturally yields a new singular case I p-SCT in Section 4. For the singular case II, the traditional approach of taking derivative in the z-domain will be used to yield a similar method in the c-domain in Section 5. The traditional approach in the s-domain, however, is to use the original polynomial plus its derivative. This seeming discrepancy with the z-domain approach is easily reconciled by taking the limit of the new c-domain method. Numerical examples will be given to demonstrate better numerical behavior of the new tests over the traditional singular case SCT that
794
H. Fan/Automatica 35 (1999) 791}807
is based on the shift-operator, for fast sampling. Real polynomials (those with real coe$cients) will be assumed throughout this paper. Generalization to the complex case should not be too di$cult.
singularity. At the last stage of Eqs. (2.4) and (2.5) when there are no zeros at either end of g ,a and (\R( (\R( a are recovered by (\R(
2. Modi5ed singularity type I SCT
1 1 a " +g #h ,, a " +g !h ,. (\R( 2 (\R( (\R( (\R( 2! (\R( (\R( ( (2.8)
The singularity type I approach of Pal and Kailath (1994) for the (unscaled) SCT is "rst presented here. Assume that in either Eqs. (1.4) or (1.5) singularity type I is encountered at n#1"J, i.e.,
Then the regular SCT can resume. The zero distribution for a is counted by (
g Oa #! a , h Oa !! a , ( ( ( ( ( ( ( (
where n is the number of roots that are on the unit circle. If J(N, i.e., a is not at the beginning stage, then ( Theorem 2 of Stoica and Moses (1992) should be used together with Eq. (2.9) in determining zero distribution of a. To eliminate numerical singularity due to fast sampling and to facilitate later development, the above Pal}Kailath recursion needs to be scaled and modi"ed. The new version is as follows:
(2.1)
where the (J#1)-dimensional vector g would have ( t zeros at the top, i.e., g "0 for j"0, 1, 2, t !1. ( ( H ( When J"N, a is used in place of a per Eq. (1.4). By the ( de"nition of ! in Eq. (1.6), zero entries on top of g are ( ( possible only when ! "$1. It then follows that ( g "$g . In other words, g must also have t zeros at ( ( ( ( the bottom, i.e., g "0 for j"J!t #1, J!t # ( H ( ( 2, 2, J. Since only singularity type I is considered here, one then must have J J even: t 4 , ( 2
J!1 J odd: t 4 . ( 2
(2.2)
The Pal}Kailath recursion (Pal and Kailath, 1994) can then be written in the following vector form. For j"J, J!2, 2, J!2t #2: ( h h q " HH "! H , H g g H H\RH H RH
(2.3)
0 "g , g H H\ 0 0
(2.4)
+n\, n, n>,"+t , 0, t ,#+those of a ,, ( ( (\R(
1 g " g if J(N, g "g if J"N, ( * ( ( ( h "h . ( (
(2.10)
For j"J, J!2, 2, J!2t #2, ( hM q "! H , H gN H RH
t "t !1, H\ H
(2.5)
(2.6)
where 0 denotes a 2t -dimensional zero vector, and RH H RH g is a ( j!2t #1)-dimensional vector de"ned as RH H H RH g "[g , g , ,g ]2. RH H H RH H RH> 2 H H\RH
(2.7)
Obviously RH g "R( g for all j, due to Eqs. (2.4) and RH H R( ( (2.6). The second equality of Eq. (2.3) holds since ! "$1 and hence g "$g and h "Gh . The reH H H H H cursions of Eqs. (2.4) and (2.5) eliminate all the top zero entries of g , one at a time, thus overcoming the (
(2.11)
0
1 g " g , H\ * H 0 0
RH g 0 RH , h "h #q RH H !q H\ H H 0 H RH g H H H R R 0
(2.9)
(2.12)
1 q h " h # H H\ * H * 0
RH g t RH H ! H 0 1 RH
#2#(!1)RH
q ! H *
t H g H t H
0 t RH ! H RH g 1 RH H
#2#(!1)RH
0 0
RH g RH H RH\
0
RH\ RH g RH H 0
t H g . H t H
(2.13)
795
H. Fan/Automatica 35 (1999) 791}807
At the last stage a
(\R(
and a are recovered by (\R(
1 a " +*g #h ,, (\R( (\R( (\R( 2 1 a " +*g !h ,, if J(N, (\R( 2! (\R( (\R( (
1 h 2" [8 !1 1 8], *
q "8,
1 h 2" [7 7] *
E out of singularity,
1 " +g !h ,, a (\R( (\R( (\R( 2! (
if J"N.
(2.14)
Otherwise the new version is the same as the Pal}Kailath recursion. Similar to the scaled SCT of Fan (1997a), it is easily seen that the scaling factors 1/* and 1/* in Eqs. (2.10)}(2.14) only scale g by 1/*J> or 1/*J, (\J h by 1/*J, and q by * for J(N from their unscaled (\J H versions. Thus there is no essential change from the Pal}Kailath recursion in this kind of scaling except that the numerical singularity due to fast sampling is eliminated. Note that Eqs. (2.10)}(2.14) have no meaningful limit form as *P0. Thus, a linear transformation will be performed in Section 3. The additional terms of the g vector in Eq. (2.13) would facilitate such a transformaH tion. However, it is due to these additional terms in Eq. (2.13) that the equivalence of the new version (2.10)}(2.14) for *"1 with the Pal}Kailath recursion (2.1)}(2.8) is not obvious. The following lemma establishes such equivalence. Lemma 1. ¸et *"1 in Eqs. (2.10)}(2.14). ¹hen in the modi,ed Pal}Kailath recursions g "g for all j"J, H H J!2, 2, J!2t and h "h where g and h are ( (\R( (\R( H H given in Eq. (2.1). See Appendix A for a proof. This lemma establishes that, other than the scaling, t , q , and a are all ( ( (\R( unchanged from the Pal}Kailath recursion. Therefore the newly modi"ed version is equivalent to Pal}Kailath recursion. The following example is taken from Pal and Kailath (1994). Example 2. f (z)"z#z#z#2z#z#1, ! "!1,
g 2"[0 0 !1 1 0 0], t "2, E singularity type I, h 2"[2 2 3 3 2 2],
t "1,
1 g 2" [!1 1], *
1 a " +g #h ,, (\R( 2 (\R( (\R(
a2"[1 1 1 2 1 1],
1 g 2" [0 !1 1 0], *
q "2,
1 a2" [3 4], *
1 a 2" [4 3]. *
Since t "2 and "! "'1, n\"2 and n>"3. This is identical (other than the scaling) to Pal and Kailath (1994). Before leaving this section, we introduce the following interesting properties of the Pal}Kailath recursion which justi"es the use of ! in Eqs. (2.8) and (2.14) and will also ( be useful in Section 4. Lemma 2. In the Pal}Kailath recursions (2.1)}(2.7), de,ne a g "a #! a , h "a !! a , ! "! HH H H H H H H H H H a H
(2.15)
for all j. ¹hen 1 a " +g #h ,, H 2 H H
1 a " +g !h ,, H 2! H H H
! "! , j"J!2, 2, J!2t #2, H ( ( *P0
! && ! . (\R( (
(2.16) (2.17) (2.18)
See Appendix B for a proof. Note that Eqs. (2.15) and (2.16) are valid for any j including both those in the singularity I block and in the regular block.
3. Singularity type I d-SCT The linear transformation described in Section 1 can be applied to the modi"ed Pal}Kailath recursion of Section 2 to obtain a singularity type I recursion for d-SCT. First de"ne u "*\H¹\ g , w "*\H¹\ h H H H H H H
(3.1)
which are consistent with Eq. (1.11). Note that there are t zeros at the top of u for each j, but not at the bottom of H H u since ¹\ is lower-triangular. Left multiplying H H *\H¹\ to Eq. (2.12) and (2.13) one obtains the following H result.
796
H. Fan/Automatica 35 (1999) 791}807
Theorem 1 (Singular I d-SC¹). ¸et J be the order index at which the d-SC¹ becomes singular (¹ype I), i.e., b "0. If J"N, de,ne (\ u "b#! bC, w "b!! bC. ( ( ( (
(3.2)
If J(N, de,ne
1 ! 0 " 1# ( u * ! ( (>
*! b 0 ( ! ( b ! 0 (> (
*! b ( b # ( , (> ! b (> (>
(3.4)
0 #* u
0 u H\ 0
H\ 0
"u , H
#* w
w H\
0
0
b
#*
(\R(\
(3.10)
b 1#! ! (\R(\ " ( (\R( u (\R( 2 0
1!! ! ( (\R( w # , (\R( 2
RH Hu H 0 H\RH
0 RH RHu H 0 H R
,
0
b
(3.11)
(\R(\
#*
(3.12)
b 1#! ! (\R(\ " ( (\R( u (\R( 0 2
1!! ! ( (\R( w # . (\R( 2*
(3.13)
¹hen the regular d-SC¹ resumes. ¹he zero distribution of b is given by ( +n\, n, n>,"+t , 0, t ,#+those of b ,. ( ( (\R(
(3.6)
where lu is u cutting out the ,rst l entries, i.e., H H (3.7)
and t q "! H , H u H RH
b " +u #w ,, (\R( (\R( (\R(
0 0 RH RH #*j RH>u #* j RH>u H H H H 0 0 H H R > R >
u2"[u l, u l , 2, u ] H H H > H H
(3.9)
=hen u and w are obtained at the end of the (\R( (\R( recursions, b and b are recovered by, for (\R( (\R(\ J"N,
b " +*u #w ,, (\R( (\R( (\R(
RHu H !(!1)RH q "w #q H H 0 H RH
l
and for J(N,
H\ 0
#2#*H\RH j H H\RH
t J j "(!1)J H # (!1) G KG m G HJ 1 G GKGIGJ O KGIG J t K t K t H ; H 2 ! H . k k l
(3.5)
0
0
t t H ! H , 2 3
where u has t zeros at its top, and Eq. (2.2) is satis,ed ( ( for t . ¹he singular I d-SC¹ recursions are given by, for ( j"J, J!2, 2, J!2t #2, ( 0
t t j " H ! H , H 1 2
$ (3.3)
! 0 " 1! ( w ! (> (
t j "! H , H 1
t t j "! H #2 H H 1 1
! 0 b ( # ( b ! 0 ( (>
! b ( b ! ( , (> ! b (> (>
and the j l's are pre-calculated constants and are given by H
(3.8)
(3.14)
See Appendix C for a proof. When J(N, Theorem 2 of Stoica and Moses (1992) with a replaced by b and a replaced by b should L L be used together with Eq. (3.14) to determine the zero distribution of b. It is also important to recognize that Eq. (3.3) has exactly the same form as Eq. (1.14). In other words, in executing Eq. (1.14) for the regular case if the RHS gives more than one zero at the top, then singular type I is naturally entered as in Eq. (3.3). The total amount of computation for the singular block is roughly O((J!2t ) t #2Jt ) #ops. The following example illus( ( ( trates the numerical advantage gained using the new
797
H. Fan/Automatica 35 (1999) 791}807
d-version against the original Pal}Kailath recursion for fast sampling. Example 1 (Continued).
"ve decimal places would the correct conclusion be reached. This kind of improvement in numerical behavior is similar to that of Fan (1997a) and is expected to become more dramatic for a smaller *.
f (z)"z!4.999z#10.001z!10.002z#5.001z!1, (3.15) F(c)"c#0.004c#0.08c#0.448c#1.024c #1.024, *"0.25.
(3.16)
Firstly the new test in the c-domain is illustrated. Only four decimal places after the decimal point are saved in the calculation. b2"[1 0.004 0.08 0.448 1.024 1.024], bC2"[!1 0.004 0.032 0.064 0.256 1.024], ! "1E singularity, J"5"N, u2"b2#! bC2 "[0 0.008 0.112 0.512 1.28 2.048], t "1 w2"[2 0 0.048 0.384 0.768 0], q "!250 u2"[0.032 0.32 0.768 2.048], [0 0 w2 ]#*[0 w2 0] "[0 !36 !208 !512 !1024 0], w2"[!144 !256 !1024 0], E out of singularity, b2"[!71.984 !127.84 !511.616 1.024], ! "1.0004, [0 b2]#*[b2 0] "[0.0608 0.3713 0.9730 2.0484], b2"[0.2561 0.4609 2.0484], ! "!1.05, [0 0 b2]#*[0 b2 0] "[!0.0025 !0.0615 !0.3484 !0.41], b2"[!0.2459 !0.41], ! "0.5832, [0 0 b2]#*[0 b2 0] "[0 !0.4554 !1.8215], b2"!1.8215. Applying Eq. (1.18) to b , it is seen that in the sequence +b , b , b , there are two entries having the same sign as b . Thus b has n\"2 and n>"1. Finally, Eq. (3.14) yields that F(c) or b have n\"3 and n>"2. Now if one uses the original (unscaled) Pal}Kailath recursion on f (z) of Eq. (3.15) and still saves four decimal places after the decimal point, an erroneous conclusion will be reached (n\"1, n>"4). Only by saving at least
4. Singularity type I p-SCT Similar to the regular case, an s-domain version of singularity type I can be obtained by taking the limit of the singularity type I d-SCT as *P0. This limit shall be called the singularity type I p-SCT. It needs to be pointed out that in the limit as *P0 the rule of counting top zeros of Eq. (2.6) is to be changed. This can be seen by taking the limit of (C.1). The new rule is t "t !2 and t "t !2l. (4.1) H\ H H\J H Taking into consideration that t will always be odd (Pal ( and Kailath, 1992), it follows that t will be odd for all j. H Thus assuming b Pn , u Pl , and w Pm as *P0, H H H H H H taking the limit of Eqs. (3.5) and (3.6) yields the following result. Theorem 2 (Singular I p-SCT). ¸et J be the order index at which the p-SC¹ becomes singular (¹ype I), i.e., n "0. If J"N, i.e., if m "0, de,ne (\ l "[0, 2m , 0, 2m , 0, 2]2, (4.2) ( m "[2m , 0, 2m , 0, 2]2. (4.3) ( If J(N, de"ne
0 n m " ( n ! ( , (4.4) (> l m 0 ( (> 0 0 "2 , (4.5) m n ( ( where l has t zeros at its top, and t satis,es ( ( ( J even: t !1(J, J odd: t (J. (4.6) ( ( ¹he singular I p-SC¹ recursions are given by, for j"J, J!2, J!4, 2, J!t #3, ( 0 0 l H\ 0
"l , H
(4.7)
0 RH RHl H #q RHl , 0 "m #q H H 0 H H RH m 0 H\ RH l q "! H , H k H RH t "t !2, H\ H
(4.8)
(4.9) (4.10)
798
H. Fan/Automatica 35 (1999) 791}807
where RHl is l cutting out the ,rst 2t entries. At the end H H H of the above recursions, execute Eqs. (4.8) and (4.9) one more time for j"J!t #1 to obtain m . ¹hen ( (\R(\ n and n are recovered by (\R( (\R(\ n " m , (\R(\ (\R(\
(4.11)
(4.12)
0
n
(\R(
"l . (\R(>
¹hen the regular p-SC¹ can resume. ¹he zero distribution of n is given by: ( for q (0 and (t !1)/2 even, or q '0 and (t !1)/2 ( ( ( ( odd,
t !1 t !1 #1, 0, ( #+those of n ,; +n\, n, n>," ( (\R( 2 2 (4.13) for q '0 and (t !1)/2 even, or q (0 and (t !1)/2 ( ( ( ( odd,
t !1 t !1 +n\, n, n>," ( , 0, ( #1 #+those of n ,, (\R( 2 2 (4.14)
where n\, n, and n> are the numbers of the left-side plane zeros, the jw-axis zeros, and the right-side plane zeros, respectively. Note that there are t zeros at the top of l with a H H reduction of two zeros for each iteration. Thus (t !1)/2 ( replaces many functions of t in the previous sections. ( Also note that in case J(N, Theorem 2 of Stoica and Moses (1992) with a replaced by m and a replaced by L L m should again be used to count the zero distribution of n. Furthermore, t "1 is not a singular case, i.e., singu( larity is declared only for t "3, 5, 7, etc. The following ( example illustrates the utility of this new algorithm. Example 4. F(s)"s#s#s#s#s#s#s#2s#1, (4.15) n2"[1 0 1 0 1 0 1 0 1], n2"[1 0 1 0 1 0 2 0], m n2![n2, 0]"[0 0 0 0 0 0 !1 0 1 ], m GFFFHFFFI l2 singularity type I, t "5: m2"[2 0 2 0 2 0 4 0], l2"[0 0 0!1 0 1], m2"[4 0 2 0 4 0],
q "2,
t "3, q "4,
l2"[0 !1 0 1], t "1, m2"[6 0 4 0], q "6, m2"[16 0], E out of singularity, n2"[!1 0 1], n2"[8 0], regular: n "!8. Since there is no sign change from m to m , n\"1 and n>"0 for that part. For the singular block, since q '0 and (t !1)/2 even, Eq. (4.14) gives n\"2 and n>"3. Finally, the sequence +m , m , m , yields n\"1 and n>"1. Thus the root distribution of n is n\"4 and n>"4.
5. Singularity type II When a "b "n "0, singularity type II is encoun( ( ( tered. All regular as well as singular type I methods fail in this case. The traditional solution in the z-domain is to replace f (z) by its derivative f (z) since n> is the (> (> same for both f (z) and f (z) (Pal and Kailath, 1994; (> (> Stoica and Moses, 1992). On the other hand, the traditional solution in the s-domain is to replace F (s) by (> F (s)#F (s) since they have the same n> (Pal and (> (> Kailath, 1992). This seeming discrepancy can be reconciled by the d-domain version as a bridge. First de"ne the parameters of f (z) by a and note the relationship (> ( a "D a , (5.1) ( (> (> where
0
J#1
$
J
D " (>
\
0
0 $
1
.
(5.2)
0
Next transform Eq. (5.1) into the d-domain by left multiplying *\(¹\ to both sides and de"ne b "*\(¹\a : ( ( ( ( b "+*¹\D ¹ , +*\(\¹\ a ,"D b , ( ( (> (> (> (> (> (> (5.3) where the identity *¹\D ¹ "D is easily veri( (> (> (> "ed by straight forward multiplication using the de"nition (1.10). Identity (5.3) implies that b is precisely the ( coe$cients of F (c), the derivative of F (c). There(> (> fore for both the z-domain and the c-domain, if singularity type II is encountered one would then work with a or b , respectively, and continue the regular case ( ( recursion and count only n> a , the number of the outside (Y unit circle roots of a , or n> b , respectively. The root (Y ( distribution of a is then (Stoica and Moses, 1992) n\ a "n> a #s\, (Y
n> a "n> a #s>, (Y
na "J#1!2n> a , (Y (5.4)
799
H. Fan/Automatica 35 (1999) 791}807
where s\(s>) is the number of entries in +a , ,\ a , 2, a , having the same sign as a (!a ). ,\ (> The same holds true for b by replacing all quantities of a by those of b. However, since the regular case d-SCT recursion depends on two previous stages, this aspect needs further examination. Let the parameter vectors of the subsequent stages after b be b , b , etc. To calculate b one cannot ( (\ (\ (\ use (1.14) since the relationship between b and b (> ( is not due to the regular d-SCT. Instead the d-SCT has to be restarted based on b using a scaled version of ( Eq. (1.12) since J(N (see Fan, 1997a for details), i.e.,
0
#*
1 b (\ " +b #! b C,. ( ( 0 * (
b (\ To calculate b C, de"ne ( 0 1
$ DC " (> $ 0
0
2
.
\
0
(5.5)
(5.6)
J#1
Then b C can be expressed by ( b C"*\(¹\a "*\(¹\DC a ( ( ( ( (> (> "+*¹\DC ¹ , +*\(\¹\ a , ( (> (> (> (> (5.7) "+!D #*DC ,bC , (> (> (> where the equality of the "rst brackets in the last step is again easily veri"ed by straight forward multiplication. Substitute Eq. (5.7) into Eq. (5.5) and observe that, since b #! bC "0, bC "(!1/! )b , (> (> (> (> (> (> 0 b #* (\ b 0 (\ 1 ! ( +D !*DC , b " b # ( (> (> (> * *! (> 1 ! ! , (5.8) " 1# ( b ! ( DC b ( ! (> (> * ! (> (> where Eq. (5.3) is used. Furthermore, using Eq. (5.3) and writing the vectors componentwise, it can be veri"ed that
b b 0 ( b ! ( " C . (5.9) (> b 0 D b (> (> (> Substitute Eq. (5.9) into Eq. (5.8) to "nally obtain 0 0 b (\
0
1 ! #* b " 1# ( (\ * ! (> 0
! b ! b ( ! ( ( b # ( . (> ! 0 ! b (> (> (>
0 b (
(5.10)
It is interesting to see that Eq. (5.10) has exactly the same form as Eq. (1.14). Thus if one uses b to replace b "0, ( ( the regular case d-SCT can continue without any interruption. Eqs. (5.3) and (5.4) with a replaced by b, Eq. (5.10), and subsequently Eq. (1.14) comprise the singular type II d-SCT. The partial inertia count n> b is given ( by Eq. (1.17) with b replacing b and b replacing b. G G ( Therefore we have proven the following. Theorem 3 (Singular II d-SCT). In the regular d-SC¹ if b "0 at some stage J, then take derivative of the ( c-polynomial associated with b ,F (c)"b c(> (> (> (> #b c(#2#b c#b , and call the co(> (> ( (> (> e.cient vector of dF (c)/dc as b . ¹he regular d-SC¹ (> ( can then continue with b replacing b . ¹he zero distribu( ( tion of b is given by Eq. (5.4) with all a's replaced by b's, in which n> b is given by Eq. (1.17) with b replaced by b and ( G G b replaced by b . ( To make the singular type II d-SCT a bridge between the SCT and the p-SCT, observe that since b and (> ,a and a "(J#1)a b "(J#1) b ( (> (> ( (> have the same signs, the value of the RHS of Eq. (1.17) with b replaced by b is unchanged by the following G G modi"cation:
n> b "var 1, (
b b b ( , ( (\ , 2, b b b (> (> (
b (> ( b L ( L\ b b (> L L
(5.11)
"n> b
(>
such that Eq. (5.4) becomes n\ b "n> b #s\, (>
n> b "n> b #s>, (>
nb"J#1!2n> . b (>
(5.12) A corresponding version for a is by simply substituting b with a. Note that although n> is as expected, the b "n> b ( (> particular form of n> in Eq. (5.11) is not trivial. This b (> form states that in counting the zero distributions one can treat b , b , 2, etc., as subsequent stages of b , ( (\ (> similar to Eq. (5.10). This point of view is essential in connecting to the p-SCT. The singular type II p-SCT is easily obtained by taking the limit of the d-version as *P0. Since b Pn , it (> (> follows that b Pn , the derivative of n . Eq. (5.3) ( ( (> becomes n "D n , and Eq. (5.10) becomes ( (> (> 0 0 m n 0 P 0 " ( n ! ( , (5.13) (> m 0 (> b n (\ (\
800
H. Fan/Automatica 35 (1999) 791}807
which is the same as Eq. (1.23). In other words, if one replaces n "0 by n , then the regular case p-SCT recur( ( sion continues without any interruption. The zero distribution count would be the same as Eqs. (5.11) and (5.12) with n replacing b. Thus again we have proven the following result. Theorem 4 (Singular II p-SCT). In the regular p-SC¹ if n "0 at some stage J, then take derivative of the ( s-polynomial associated with n , F (s)"m s(> (> (> (> #m s(#2#m s#m , and call the coef(> (> ( (> (> ,cient vector of dF (s)/ds as n . ¹he regular p-SC¹ can (> ( then continue with n replacing n . ¹he zero distribution of ( ( n is given by Eq. (5.4) with all a's replaced by n's, in which n> n is given by Eq. (5.11) with all b's replaced by m's. (
It is important to point out that, since m "0 and L J m "0 for l odd, at any stage j the p-SCT can be viewed K J as a recursive procedure to test zero distribution of an s-polynomial with the coe$cient vector
1 0 n# 2 H n H\
1 " [m , m ,m ,m , ]2. 2 H H\ H H\ 2 (5.14)
This is certainly true for j"N in light of Eqs. (1.21) and (1.22). For j(N, if one applies the rules of Eqs. (1.21) and (1.22) to test the zero distribution of the vector on the RHS of Eq. (5.14), this vector will be separated into n and H n . Thus, the p-SCT is unchanged whether n and H\ H n are obtained as intermediate stages from the p-SCT H\ or as the beginning stages with Eq. (5.14). In other words, the statement leading to Eq. (5.14) is true for all j. In particular, for j"J#1 the singular-type II p-SCT developed above, due to the continuity in Eqs. (5.10) and (5.11), is then to test the zero distribution of an s-polynomial with the parameter vector +n #[n ], and the ( (> zero distribution count (5.12) with s\"s>"0. Clearly, this polynomial is +F (s)#F (s),, which is identi (> (> cal (except a factor of 2) to the traditional solution to singularity-type II (Pal and Kailath, 1992). Thus, the seemingly di!erent approaches in the z-domain and the s-domain for singularity-type II are reconciled. The following example demonstrates the singularitytype II algorithms in all three domains, thus showing the numerical advantages of the d-SCT and p-SCT.
a"[1 !4 5.9999 !3.9998 0.9999]2, ! "!0.9999, a " [0.00019999 !0.00059998 0.00059999 !0.0002]2. Even without going into the Pal/Kailath recursion, it is clear from a that unless one keeps 8 digits after the decimal point, a would be rounded to a "[0.0002 !0.0006 0.0006 !0.0002]2 which corresponds to a triple root at 1. This erroneous result is due to the sensitivity of the SCT to fast sampling. Now perform the d-SCT and keep only four digits after the decimal point. d-SC¹: b"[1 0 !1 0 0]2, bC"[0.9999 !0.02 !1 0 0]2, ! "!0.9999, [0 b2]#*[b2 0] "[0.0002 0.02 !0.0001 0 0]2, b "[0.02 !0.0001 0 0]2, ! "1, [0 0 b2]#*[0 b2 0] "[0 !0.0001 !0.02 0 0]2, b "[!0.02 0 0]2, ! "!1, [0 0 b2 ]#*[0 b2 0]"[0 0.0001 0 0]2, b "0 E singularity II, J"1, b "[!0.04 0]2, ! "1, [0 0 b 2 ]#*[0 b 2 0]"0,
Example 3.
b "0 E singularity II again,
f (z)"(z!1)(z!1.01) (z!0.99) "z!4z#5.9999z!3.9998z#0.9999,
These three equations describe the same system but in three di!erent domains. From Eq. (5.15) it is clear that +n\, n, n>,"+1, 2, 1, which is also true for Eqs. (5.16) and (5.17). Now perform the SCT, d-SCT, and p-SCT with "nite word-length arithmetic. SC¹:
(5.15)
F (c)"c!c, *"0.01, @
(5.16)
F (s)"0.999975s#0.01s!s. K
(5.17)
b "!0.04. The inertia is counted as follows. From b , b , and b it is clear that s\"1 and s>"1. From b and b we have n> b "0. Thus +n\, n, n>,"+1, 2, 1, according to Eq. (5.12).
801
H. Fan/Automatica 35 (1999) 791}807
Now perform the p-SCT with only two digits after the decimal point. p-SC¹: n"[1 0.01 !1 0 0]2, n "[2 0 !2 0 0]2, n "[0.02 0 0 0]2, n "[!0.02 0 0]2, n "0 E singularity II, J"1, n "[!0.04 0]2, n "0 E singularity II again, n "!0.04. The inertia count, and hence the conclusion, is the same as the d-SCT.
Then showing h "h is equivalent to showing (\R( (\R( "0. The recursion for h can be obtained by subh H (\R( tracting Eq. (2.5) from Eq. (2.13), noting that *"1 and RH g "R( g ∀j, RH H R( ( 0
R( g 0 RH "h #qJ R( ( !qJ h H H 0 H R( g H\ H ( ( R R 0
t #q ! H H 1
t H ; t !1 H
6. Conclusion Singular cases for the recently developed d-SCT and p-SCT are considered in this paper. &&Numerical singularity'' due to fast sampling is distinguished from &&algorithmic singularity''. Numerical singularity of the SCT is overcome by the d-SCT. Thus only algorithmic singularity of types I and II are considered in this paper. For type I singularity, a recent order recursive approach of Pal and Kailath is scaled and modi"ed. It is then transformed into the d-domain. Its limit as *P0 yields an s-domain version. Numerical examples are given to show the numerical advantage of the d-version against the original approach of Pal and Kailath for fast sampling. Finally, singularity type II is being tackled. A d-version is obtained from the traditional z-domain method, and an s-domain version results when taking the limit *P0. The resulting s-domain version coincides with the traditional approach, and is thus reconciled and uni"ed with the z-domain approach. A numerical example is also given for singularity type II algorithms.
0 R( g #2#(!1)RH\ R( ( 0 RH\
0 RH\ R( g R( ( 0 RH>
#2#(!1)RH\
t !q ! H H 1
t H t !1 H
0 RH> R( g R( ( 0 RH\
0 RH\ R( g R( ( 0
,
(A.2)
hI [1 0 2 0]h H , qJ "q !q "! H "! H H H [1 0 2 0]R( g g ( H RH
(A.3)
where
and the last terms in the brackets of Eq. (2.13) cancel. It is clear by Eq. (2.10) that h "0, which also implies that ( qJ "0. Thus recursion (A.2) can be evaluated for j"J to ( obtain
0 t h "q ! ( ( (\ 1 0
0 #2#(!1)R(\
R( g R( (
0 R(\ 0 R(\ R( g R( ( 0 R(>
Acknowledgements
t ( ; t !1 (
t !q ! ( ( 1
This work began when the author visited Systems and Control Group, Uppsala University, Uppsala, Sweden in 1995. The author wishes to thank Prof. P. Stoica for bringing up the importance of the singular cases and for providing a reprint of the reference, Stoica and Moses (1992), during his visit.
t ( #2#(!1)R(\ t !1 (
0 R(> R( g R( ( 0 R(\
0 R(\ R( g R( ( 0
.
(A.4)
Using Eq. (A.3) one obtains from Eq. (A.4) Appendix A. Proof of Lemma 1 When *"1, g "g ∀j is obvious comparing Eq. (2.12) H H with Eq. (2.4) and observing Eq. (2.10). De"ne h "h !h , H H H
j"J, J!2, 2, J!2t . (
(A.1)
[1 0 2 0]h [1 0 2 0]R( g (\"! R( ( qJ "! (\ [1 0 2 0]R( g [1 0 2 0]R( g ( (
t t ; ! ( q " ( q . 1 ( 1 (
(A.5)
802
H. Fan/Automatica 35 (1999) 791}807
Use Eqs. (A.4) and (A.5) to carry out recursion (A.2) further. One obtains 0
0
t h "q ! ( ( (\ 1 0 0
#(!1)R(\
0
t ( #2#(!1)R(\ t !1 (
R( g R( (
0 R(\
#(!1)R(\\
#(!1)R(\\
t ( t !1 (
0 R(> R( g R( ( 0 R(\
(\ t !1 (\
#q (\
0 R(\ R( g R( ( 0 R(\>
t (\ t !1 (\
t
0
0 R(\> R( g R( ( 0 R(\
0 R(\>
t t ( q ! (\ q (\ 2 ( 1
0
.
#(!1)
0 R(\ R( g R( ( #2 0
0 R( g R( ( 0 R(\
#2
0
(A.6)
On the other hand, recursion (A.2) can be executed backwards from its last stage, i.e., 0 h "qJ (\R( (\R(> 0
!qJ
R( g R( ( 0 0
(A.8)
t !l#1 ( q , (\J\ 1
where 14l4t !1 and Eq. (2.6) is used. (
(A.9)
(\R(>
0 R( g R( (
, #h (\R(>
(A.10)
where t "t !(t !1)"1 is used. Appending (\R(> ( ( a zero at each end and using Eq. (A.2) again on h , (\R(> Eq. (A.10) becomes 0
0
0
R( g 0 R( ( "q h 0 !q 0 (\R( (\R(> (\R(> R( g 0 0 R( ( 0 0 0 0
R( g R( ( 0
t t !1 ( q #(!1)J ( q #2 ( l l!1 (\
0
.
Continue the recursion, one sees that the terms with qJ in H Eq. (A.2) always cancel the corresponding terms in h of H Eq. (A.2) which are most shifted vectors of R( g for each j. R( ( The next most shifted vectors of R( g contribute to the R( ( "rst non-zero component of h . By Eq. (A.3) we can H\ deduce that
0 R(\> t R( g !q #q ! (\ (\ R( ( (\ 1 0
0
[1 0 2 0]h t (\"(!1) ( q qJ "! (\ [1 0 2 0]R( g 2 ( ( t #(!1) (\ q . (\ 1
qJ "(!1)J> (\J
0 R(\ R( g #2 R( ( 0
(A.7)
R( g R( ( 0 R(\ Using Eq. (A.7) in Eq. (A.3), one obtains
t !q ! ( ( 1
t !q ! (\ (\ 1
Note that in Eq. (A.6) the "rst terms in the "rst two brackets cancel the terms with qJ due to Eq. (A.5). (\ Therefore, the "rst (from the top) non-zero entry of h is due to the following terms: (\ 0 (!1)
R( g R( (
0 R(\ R( g R( ( 0 R(>
#q (\R(>
0 0 0
0 0 !q (\R(>
0 0 R( g R( (
803
H. Fan/Automatica 35 (1999) 791}807
we get
0
R( g R( ( 0
2 #q (!1) (\R(> 1
0
RH g 0 RH , t 51. "2a #q RH H !q 2a H H 0 H RH g H H\ RH RH H 0
0
(B.1)
0
The zeroth and ( j!2)th entry of a can be obtained H\ from Eq. (B.1), using Eq. (2.15), as
0 0
2 !q (!1) (\R(> 1
0
Rg R ( 0
#h . (\R(>
(A.11)
( (
Continue this process and group vectors of the same kind. It can be deduced that the coe$cients of these vectors are given by
0 0 R(\ R(> R( g ! R( g R( ( R( ( 0 0 R(> R(\ q
(\R(>
#(!1)
By the de"nition in Eq. (2.15), ! can be obtained from H\ Eqs. (B.2) and (B.3),
2 t ( q #2#(!1)R(\ qN ( (\R > 1 t !1 ( (
0
$ 0
!
t !l#1 t ( q #2#(!1)J ( qN (\J> 1 l (
$
0 R(\
(B.3)
terms:
q #(!1) (\J
q a "a ! H (a #! a ). H\ H\ H H\ 2 H H\RH\ H H RH>
$
0 0 J R(\J R( g ! R( g R( ( R( ( 0 0 R(\J J
R( g R( (
(B.2)
terms:
$
q a "a # H (a #! a ) H\ H 2 H RH> H H H\RH\
R(\ R( g R( ( 0
terms: q
#(!1) (\
t ( qN , 1 (
where 14l4t !1. However by Eqs. (A.5), (A.8), and ( (A.9) all these coe$cients equal to zero. Since Eq. (A.12) , 02 ]2, it is concluded contains all the terms in [02 , h 2 R( (\R( R( that h "0. (\R(
(A.12)
a !OH (a #! a ) a H H RH> ! "! H\ H\"! H H\ H H\RH\ H\ a #OH (a #! a ) a H H RH> H H H\RH\ H\
#! a ) a !OH (a H H\ H H\RH\ H H RH> "! . H !! a !OH (!a #! a ) H H H H H\RH\ H H RH> (B.4)
Appendix B. Proof of Lemma 2 Eq. (2.16) is easily solved from Eq. (2.15). Thus what needs to be proven is Eqs. (2.17) and (2.18). To this end we need to express Eqs. (2.4) and (2.5) in terms of a . H Thus adding Eq. (2.4) onto Eq. (2.5) and using Eq. (2.16),
For j"J, J!2, 2, J!2t #4, t '1 ( H g "0, i.e., by Eq. (2.15), H H\ a "!! a . H H\ H H
and
hence
(B.5)
804
H. Fan/Automatica 35 (1999) 791}807
Also, since g "0 in this range of j it is easily seen from H Eq. (2.15) that "! ""1, or !"1. Using these facts in H H Eq. (B.4) one easily concludes ! "! , j"J, J!2, 2, J!2t #4 H\ H (
(B.6)
which is exactly Eq. (2.17). To show Eq. (2.18), observe that for j"J!2t #2, t "1 and thus !"1 is still ( P H H P (!1)H> and a (!1)H valid. In addition, ! & & H H H\ a (Vijayan et al., 1991). Thus, although Eq. (B.5) does H not hold for a "nite *, but in the limit it still holds true. Therefore from Eq. (B.4), ! P! . (\R( (\R(>
(B.7)
Eqs. (B.7) and (2.17) imply Eq. (2.18).
Appendix C. Proof of Theorem 1 Left multiplying *\H¹\ to Eq. (2.12) and using the H identities (14) of Vijayan et al. (1991), 02
¹\ ¹ u " H H\ H\ 02
02
02
02 #*
I
I
02
0
0 #* u
0
" u
H\
H\ 0
"u . H
(C.1)
Left multiplication of *\H¹\ to Eq. (2.13) is slightly H more involved due to the shifted RH g terms. Thus let us RH H "rst investigate how those shifted RH g terms can be transRH H formed. Observe from Eq. (2.7), g 2"[02 , RH g 2, 02 ]"[02 , RH g 2 ] H RH RH H RH RH H
(C.2)
RH g 2"[RH g 2, 02 ]. H RH H RH
(C.3)
Thus, the "rst brackets in Eq. (2.13) contain only upward shifted RH g , or RH g . For k4t , one can calculate RH H H H 0 RH\I *\H¹\ RH g " ¹\ H H H 0 I
+*\H¹\ g ,"M u , H H H H
0H\I>;I I H\I> ¹ 0I;I 0I;H\I> H
where M is the product of the three matrices in the "rst H brackets and can be calculated as
0I\;H\I> I
H\I> 02 H\I>
k #2#*I k
* 2 * $ $
I H\I>
* 2 *
0I;H\I>
, (C.5)
where * denotes a &&do-not-care'' entry. Since u has H t zeros at the top, when Eq. (C.5) is multiplied by u the H H &&do-not-care'' entries would have no e!ect on the outcome. Then substitute Eqs. (C.5) into (C.4) to obtain, observing the de"nition of Eq. (3.7), 0 RH\I k *\H¹\ RH g "u #* H H H 1 0 I
u H #2 0
Iu H . 0 I
k k
(C.6)
The "rst brackets in Eq. (2.13) would then be transformed into the following, using Eq. (C.6) for k"t , H t !1, 2, 1, 0, H *\H¹\ H
RH g t H ! H 0 1 RH
0
;
t H t H
RHu H 0 RH
t !1 H t !1 H
t ; H m
u #* H
t H 1
#(!1)
t !m H l
RH\Ku H 0 RH\K
RH\u H 0 RH\
t !m H 1
0 RH RH g H
u t !1 H #2#*J H 0 l
#*RH\
;
t RH g #2#(!1)RH H H t H 0 RH\
u t H #2#*J H l 0
t " u #* H H 1
#*J (C.4)
k 1
I H\I>
* 2 *
#*RH
where
#*
* 2 * ; $ $
#*I
u H\
0I;H\I>
* 2 * M" $ $ H * 2 *
Ju H #2 0 J
u #* H
t !1 H 1
Ju H #2 0 J
#2#(!1)K
u H #2 0
Ju t !m H #2#*RH\K H 0 t !m J H
#2#(!1)RH
t H u. H t H
(C.7)
805
H. Fan/Automatica 35 (1999) 791}807
The terms [Ju2, 02 ]2 can be combined together with H J coe$cients given by, for 04l4t !1, H
R \J t *J(!1)K H m K H
*\H¹\ H
t !m H l
RH\J t! (t !m)! H H "*J (!1)K m!(t !m)! l!(t !m!l)! H H K
*J RH\J t !l " t (t !1)2(t !l#1) (!1)K H "0. H l! H H m K (C.8) That is, all terms on the right hand-side (RHS) of Eq. (C.7), except *RH (RH ) [RHu2, 02 ]2, add up to zero. Thus RH H RH
RHu H . 0 RH
(C.9)
Now turn attention to transforming the second brackets of Eq. (2.13) which contain only downward shifted RH g . RH H De"ne g "[02 , RH g 2 ]2. RH H RH RH H
(C.10)
Then the downward shifts of RH g in the second brackets RH H of Eq. (2.13) become upward shifts of g , i.e., RH H 0 RH\ t RH g #2#(!1)RH H RH H t H 0
0 t RH ! H RH g 1 RH H
0 t " RH ! H g 1 RH H
"(!1)RH
0
#(!1)RH
t H t H
0 RH g RH H
g RH H 0 RH
0
g t RH H ! H 0 1 RH
g #2 RH H 0 RH\ ,
0 RH . g RH H
RH\ RH g #2 RH H 0
0 RH RHg H . H g R #(!1)R "(!1)RH *RH (C.13) RH H 0 H R 0 RH Now g needs to be expressed in terms of u . To do so H H write u as H g u "*\H¹\ g "*\H¹\ RH H H H H H 0 RH 0 I 0 ; H\RH> "¹\ H\RH> RH ¹ *\H¹\ RH H H H g 0RH;RH 0RH;H\RH> RH H g RHg t t H #2#*RH H H , "g #* H (C.14) H 0 1 t 0 H H R where Eq. (C.5) and the fact that g has t zeros at the top H H are used. Actually, since g has t zeros at the top, from RH H H Eq. (C.12) it is seen that g has 2t zeros at the top. The H H non-zero part RHg can be solved from Eq. (C.14) starting H from the last element.
t H t H
u "g , H H H H
t u "g #* H g , H H\ H H\ 1 H H $ t t H g , u "g #* H g #2#*H\J H J H J 1 H J> j!l H H (C.15)
(C.11)
(C.12)
where 2t 4l4j. Note that in Eq. (C.15) we have de"ned H (RH)"0 for k't . RHg can then be solved from Eq. (C.15) I H H as RH>u H RH>u H #*j RHg "RHu #*j 0 #2 H H H H 0 0
Hu H , (C.16) H\RH where Hu "u and j is given by Eq. (3.9). H HH HJ Finally, left multiplying Eq. (2.13) by *\H¹\ and H using Eqs. (C.1), (C.9), (C.13), and (C.16), one obtains Eq. (3.6), where #*H\RH j
where (RH)"( RH ) is used. It is seen that the RHS of RH\J J Eq. (C.11) has the same form as the terms in the brackets of the LHS of Eq. (C.7). Thus de"ne g "*\H¹\ H H
0
0 RH RH g RH H 0 RH
RH\ t g #2#(!1)RH H H R H t H 0
0 t RH ! H RH g 1 RH H
H
*J RH\J (t !l)! H " t (t !1)2(t !l#1) (!1)K H H H l! m!(t !l!m)! H K
LHS of Eq. (C.7)"*RH
One can then use the results of Eqs. (C.4)}(C.9) to obtain
H H\RH 0
*RH [1 0 2 0] h H q "*RHq "! H H [02 , 1, 02 ] g RH H\RH H [1 0 2 0] *H¹ w t H H"! H . "! [02, 1, 02 ] *H\RH ¹ u u RH H\RH H H H RH
(C.17)
806
H. Fan/Automatica 35 (1999) 791}807
The last equality of Eq. (C.17) is true since ¹ is lowerH triangular and u has t zeros at the top. This proves the H H recursions (3.5)}(3.8). Now study how these recursions are entered. When J(N, i.e., singularity is not encountered at the very beginning, scaling by 1/* on g is needed. Left multiply( ing *\(¹\ to Eq. (2.10) and using Eq. (2.1), one obtains ( 1 u " +b #! bC,, w "b !! bC. ( * ( ( ( ( ( ( (
(C.18)
To eliminate bC from Eq. (C.18), note that the stages ( before Eq. (C.18) are regular and that u in Eq. (C.18) has ( the same form as the RHS of the scaled version of Eq. (1.12) as in Eq. (5.5). Thus, Eq. (1.14) applies, i.e., one obtains Eq. (3.3). Similar manipulation on w yields ( Eq. (3.4). On the other hand, if J"N, then g "g and left ( ( multiplication of *\(¹\ to Eq. (2.1) gives Eq. (3.2). ( To exit the singular type I recursion, one also needs to consider the J(N and J"N cases separately. First assume J(N. Then transforming Eq. (2.14) gives Eq. (3.12) and 1 bC " +*u !w ,. (\R( (\R( (\R( 2! (
(C.19)
Since we are now back in the regular case, the transformed version of Eq. (1.5), i.e., the scaled version of Eq. (1.12) as in Eq. (5.5), is in e!ect. Thus using Eq. (C.19) and the scaled version of Eq. (1.12) one has
0
b
#*
(\R(\
b 1 (\R(\ " +b #! bC , (\R( (\R( 0 * (\R(
1#! ! 1!! ! ( (\R( w ( (\R( u # , " (\R( (\R( 2 2*
(C.20)
where ! is given by Eq. (1.15) and the fact that (\R( ! "$1 is used. ( In case J"N, Eq. (C.19) becomes Eq. (3.10) and 1 bC " +u !w ,. (\R( 2! (\R( (\R( (
(C.21)
Accordingly, the factor 1/* on the RHS of Eq. (C.20) needs to be removed, i.e.,
0
b
#*
(\R(\
"b
b (\R(\ 0
Appendix D. Proof of Theorem 2 The recursions (4.7)}(4.9) are immediately obtained by taking the limit (*P0) of Eqs. (3.5), (3.6), and (C.17). Thus we only need to prove how to enter and exit the recursions. To enter the recursions of Eqs. (4.7) and (4.8), "rst consider J(N. Taking the limits of Eqs. (3.3) and (3.4) and using the facts that ! P(!1)H> and H 1/* (1#! / ! )P0 (Fan, 1997a), one obtains Eqs. H H> (4.4) and (4.5). Comparing Eq. (4.4) with Eq. (1.23) it is seen that k "l "0 for i even and that t "1 is not H G H G> ( a singular case. Now consider J"N. Taking the limit of Eq. (3.2) and considering Eqs. (1.12), (1.21) and (1.22), one obtains Eqs. (4.2) and (4.3). To exit the recursions (4.7) and (4.8), consider replacing t by (t !1)/2, i.e., J!2t by J!t #1 at the last stage. ( ( ( ( Note that although when t is reduced to being 1 (l H (\R(> has one zero at the top) it becomes consistent with the regular case, singularity is not exhausted at this step since RHl still cuts out the last zero at the top and Eq. (4.8) can H be executed one more time to produce m . Therefore (\R(\ for J(N, take the limit of Eq. (3.12) at J!t !1, ( Eq. (4.11) is obtained. On the other hand, observe from Lemma 2, Eq. (2.18) that in the limit ! has the same (\R( sign as ! . In other words, for * small ! ! ""! " ( ( (\R( (\R( in Eq. (3.13). Using Eq. (1.15) and the fact that m "0 for L P 0 in Eq. (3.13). Thus replacall n, (1!! ! )/2* & ( (\R( ing t by (t !1)/2, the limit of Eq. (3.13) is given by ( ( Eq. (4.12). Now for J"N the limit of Eq. (3.11) gives again Eq. (4.12). To obtain n , observe that since (\R(\ is still in the singular block, Eq. (3.10) is although n (\R(> valid but has no meaningful limit for J!t #1. On the ( other hand, since l does not exist, Eq. (3.10) does (\R(\ not have a limit for J!t !1 either. Thus n can( (\R(\ not be obtained from Eq. (3.10). This dilemma can be resolved by comparing Eqs. (4.2)}(4.3) with (1.21)}(1.22). Obviously n "m and [0, n2 ]"l2 in this case. If we , ( ,\ ( de"ne n " m and [0, n2 ]"l2 for j"J!2, J!4, H H H\ H 2, J!t(#1 and n(\R(\" m(\R(\, the singular type I recursions (4.7) and (4.8) can then be expressed in terms of n , j"J, J!1, 2, J!t !1. The desired n is H ( (\R(\ then naturally given again by Eq. (4.11). Thus the case J"N is consistent with J(N. The zero distribution rules of Eqs. (4.13) and (4.14) are derived by Delsarte et al. (1985) and used by Pal and Kailath (1992).
References
C
#! b (\R( (\R( (\R(
1#! ! 1!! ! ( (\R( u ( (\R( w " # . (\R( (\R( 2 2
(C.22)
Delsarte, P., Genin, Y. V., & Kamp, Y. G. (1985). Pseudolossless functions with application to the problem of locating the zeros of a polynomial. IEEE ¹rans. Circuits & Systems, 32(4), 373}381.
H. Fan/Automatica 35 (1999) 791}807 Fan, H. (1997a). E$cient zero location tests for delta-operator based polynomials. IEEE ¹rans. Automat. Control, 42(5), 722}727. Fan, H. (1997b). A normalized Schur-Cohn stability test for the delta operator based polynomials. IEEE ¹rans. Automat. Control, 42(11), 1606}1612. Goodwin, G. C., Middleton, R. H., & Poor, H. V. (1992). Highspeed digital signal processing and control. Proc. IEEE, 80(2), 240}259. Jury, E. I. (1982). Inners and stability of dynamic systems (2nd ed). Malabar, FL: Robert E. Krieger Publishing Co. Pal, D., & Kailath, T. (1992). Displacement structure approach to singular root distribution problems: The imaginary axis case. IEEE ¹rans. Circuits Systems I, 41(2), 138}148. Pal, D., & Kailath, T. (1994). Displacement structure approach to singular root distribution problems: The unit circle case. IEEE ¹rans. Automat. Control, 39(1), 238}245. Stoica, P., & Moses, R. L. (1992). On the unit circle problem: The Schur}Cohn procedure revisited, Signal processing, 26(1), 95}118. Vijayan, R., Poor, H. V., Moore, J. B., & Goodwin, G. C. (1991). A Levinson-type algorithm for modeling fast-sampled data. IEEE ¹rans. Automat. Control, AC-36(3), 314}321.
807
H. (Howard) Fan received the B.S. degree from Guizhou University, China, in 1976, and the M.S. and Ph.D. degrees from the University of Illinois, Champaign-Urbana, in 1982 and 1985, respectively, all in electrical engineering. From 1977 to 1978, he worked as a Research Engineer at the Provincial Standard Laboratory and Bureau of Guizhou Province in China. In 1978, he entered the Graduate School of the University of Science and Technology of China and then transferred to the University of Illinois where he was a Teaching Assistant and a Research Assistant from 1982 to 1985. He joined the University of Cincinnati in Electrical and Computer Engineering in 1985 as an Assistant Professor, where he is now a Professor. His research interests are in the general "elds of systems and signal processing, in particular signal representation and reconstruction, adaptive signal processing, signal processing for communications, and system identi"cation. He spent six months in 1994 with the Systems and Control Group, Uppsala University, Sweden, as a visiting researcher. He received the "rst ECE Departmental Distinguished Progress in Teaching Excellence Award at the University of Cincinnati in 1987, and served as an associate editor for IEEE ¹ransactions on Signal Processing from 1991 to 1994. He is a member of Tau Beta Pi and Phi Kappa Phi.
Singular root distribution problem for delta-operator based real polynomials H. Fan* Department of Electrical & Computer Engineering and Computer Science, University of Cincinnati, Cincinnati, OH 45221-0030, USA Received 14 November 1996; revised 17 October 1997; received in "nal form 5 November 1998
Delta operator based Schur}Cohn tests for root distribution of a polynomial in the singular cases I and II are developed. These singular case tests connect and unify the traditional discrete time Schur}Cohn tests and their continuous-time counterparts in the singular cases.
Abstract E$cient zero location tests for the delta-operator-based polynomials have been a recent topic of interest. These tests not only provide the needed e$cient tests for the delta-operator formulated systems, control, and signal processing problems, which are numerically better behaved than the traditional shift-operator based tests for fast sampling, they also provide missing links between classical tests in the z-domain and the s-domain such as the Schur}Cohn test (SCT) and the Routh test (RT). These will also be accomplished in this paper, but for the singular cases. Speci"cally, the recently developed singular case I Schur}Cohn test by Pal and Kailath (1994) is modi"ed and transformed into the delta-domain. As the sampling interval vanishes this new algorithm converges to an s-domain singular case I test associated with the recently developed s-domain SCT. Singular case II is also considered and the di!erence between the traditional approaches in the z-domain and the s-domain reconciled. Numerical examples are given to illustrate better numerical properties of the new test over Pal-Kailath recursion for fast sampling. 1999 Elsevier Science Ltd. All rights reserved. Keywords: Stability tests; Inertia tests; Delta operators; Singularities
1. Introduction The root distribution problem, or in other words the zero location tests, has been studied for over 100 years. Classical e$cient tests include the Schur}Cohn test (SCT) in the discrete-time shift operator case, denoted the z-domain, and the Routh test (RT) in the continuoustime case, denoted the s-domain (Jury, 1982). Recently, the d operator, i.e., d"(q!1)/* where qx(t)"x(t#1) 22222 * Corresponding author. Tel.: (513) 556-4765; fax: (513) 556-7326; e-mail: [email protected]. This paper was not presented at any IFAC meeting. This paper was recommended for publication in revised form by Associate Editor C. V. Hollot under the direction of Editor R. Tempo. This work is supported by the O$ce of Naval Research under Grant N00014-96-1-0241, and by the National Science Foundation under Grant INT-9321821.
and * is the sampling interval, has been used to develop such tests for the d-operator based discrete-time systems (Fan, 1997a, b). These new tests shall be referred to as in the d-domain or in terms of its frequency variable the c-domain. The new c-domain tests are computationally e$cient with O(N) computations where N is the system order. They are numerically better behaved than the traditional shift-operator based tests such as the SCT for fast sampling. In addition, they provide smooth transitions between the z-domain tests and the s-domain ones. For example, the limit (as *P0) of the d-SCT gives an s-domain SCT, called p-SCT (Fan, 1997a). However, these tests are all developed for the strongly regular case. They all fail in the singular cases. First let us review the strongly regular case of SCT, d-SCT, and p-SCT (Fan, 1997a) and show why they would fail in the singular cases. The formulation of this section will also be useful for the subsequent
0005-1098/99/$*see front matter 1999 Elsevier Science Ltd. All rights reserved PII: S 0 0 0 5 - 1 0 9 8 ( 9 8 ) 0 0 2 1 8 - 0
792
H. Fan/Automatica 35 (1999) 791}807
development. Consider an Nth-order real polynomial in the z-domain (assuming a O0): f (z)Oa z,#a z,\#2#a z#a . (1.1) ,\ , Using the vectors z O[z,, z,\, 2, z, 1]2 ,
Using the relationship c"(z!1)/* (Goodwin et al., 1992), it has been shown (Vijayan et al., 1991) that c "¹2 z where ¹ is a lower-triangular matrix whose , , , , l, kth entry is given by
(!1)J\I N!k (¹ ) " , 04l, k4N, l5k. , J I *,\I l!k
(1.10)
(1.2) aO[a , a , 2, a , a ]2 ,\ , Eq. (1.1) can be written as f (z)"z2 a. Here and through, out this paper, boldface letters will denote vectors whereas unboldfaced letters denote scalars. Furthermore, de"ne
To keep the identity f (z)"z2 a"*,c2 bO*,F(c), one , , needs to de"ne a"*,¹ b as in Vijayan et al. (1991). , C De"ne b "[b b 2b ]2 and b "[bC bC 2bC ]2, L L L LL L L L LL and
a OJ a, , J Oantidiag.+121,, (1.3) , and the lower-order parameter vectors a "[a a 2 L L L a ]2 and a "[a 2a a ]2 for all n"0, 1, 2, N!1. LL L LL L L Then the SCT is given by
Left multiply *\,¹\ and *\L\¹\ to the SCT, and , L> after some algebra, one obtains
a
b "*\L¹\ a , bC"*\L¹\ a . L L L L L L
0
b
(1.4)
a L "a #! a , n"N!2, N!3, 2, 0, L> L> L> 0 (1.5)
#*
,\ 0 0
,\ "a#! a , , 0
(1.11)
b ,\
b ,\ "b#! bC, , 0 0
(1.12)
#* b "! ,\ 0
b ,\ 0
(1#! ! ) (! #! ) C ,\ , b# ,\ , b # * *
(1.13)
and for n4N!3, a a ! "! LL , ! "! , , L , a a L
(1.6)
n>"var 1, 1!! , (1!! ) (1!! ), 2, , , ,\ , (1!!)G G G
(1.7)
or n\"C of entries in +a , ,a , ,\ 2 having the same sign as a , (1.8) where n> is the number of outside unit circle roots and n\ that of inside ones, and var+x, denotes the variation of signs in the sequence +x,. The singular cases occur when a "0 for some n. From Eq. (1.6) it is clear that L when this happens ! cannot be computed, therefore the L SCT cannot be continued. The same situations occur with the d-SCT of Fan (1997a). We "rst brie#y review the d-SCT. It turns out that the d-SCT can be easily obtained by a linear transformation of a scaled version of the SCT in Eqs. (1.4)}(1.8). First de"ne the d-operator quantities. c O[c,, c,\, 2, c, 1]2, , bO[b , b , 2, b , b ]2. ,\ ,
(1.9)
0
0
1 ! 0 #* b " 1# L> L * ! L> b 0 L
0
b L>
b ! ! b # L> L> ! L> L> b , L> ! 0 ! b L> L> L> where
(1.14)
1 !" [(!1)L>, (!1)L*, (!1)L\*, 2, !*L]b , L b L L (1.15) bC"*\,¹\ a "¹\J ¹ b. , , , , See Fan (1997a) for details. The inertia test is
(1.16)
b b b ,\ , ,\ ,\ , 2, n> b "var 1, b b b ,\ , b L L\ , (1.17) b L L n\"C of entries in +b , ,b , ,\ 2 having the same sign as b . (1.18) The singular cases of the d-SCT occur when b "0 for L some n. Similar to the SCT, from Eq. (1.15), it is seen that
H. Fan/Automatica 35 (1999) 791}807
in the singular cases ! cannot be computed, thus the L d-SCT cannot be continued. The following example illustrates one singular case for both the SCT and the d-SCT. Example 1. f (z)"z!4.999z#10.001z!10.002z#5.001z!1, (1.19) F(c)"c#0.004c#0.08c#0.448c#1.024c #1.024, *"0.25.
(1.20)
These two characteristic functions correspond to the same system/process but formulated using the shiftoperator and the d-operator, respectively. Now perform the SCT on f (z) and the d-SCT on F(c). SC¹: a"[1 !4.999 10.001 !10.002 5.001 !1]2, ! "1, a "[0 0.002 !0.001 !0.001 0.002]2, ! "?
a "0, d-SC¹:
b"[1 0.004 0.08 0.448 1.024 1.024]2, bC"[!1 0.004 0.032 0.064 0.256 1.024]2, ! "1, b "[0 0.032 0.32 0.768]2, b "0,
! "?.
The same situation may occur to the p-SCT in the s-domain as well. The p-SCT in the s-domain can be obtained by simply taking the limit *P0 of the d-SCT. By doing so ! P(!1)L> and Eq. (1.12) becomes L b Pn "[2m , 0, 2m , 0, 2]2, ,\ ,\
(1.21)
where n"[m , m , 2, m ]2 is the continuous-time sys , tem parameter vector which is the limit of b as *P0. Furthermore, if one de"nes n "[2m , 0, 2m , 0, 2]2, ,
(1.22)
then Eqs. (1.13) and (1.14) give the same limit 0
0
m n 0 P 0 " L> n ! L> L> m 0 L> b n L L
(1.23)
793
for n4N!2. Therefore, Eqs. (1.21)}(1.23) constitute the p-SCT. The inertia test is Eqs. (1.17) and (1.18) with b's replaced by m's. Note that m "0 for j odd and for all L H n. It is also clear from Eq. (1.23) that singularity occurs when m "0 for some n, for which case Eq. (1.23) cannot L be continued at stage n!2. Note that the above kind of singularities is not to be confused with the numerical singularity of the SCT due to fast sampling. To be speci"c, assume that the discretetime characteristic polynomial is obtained by sampling a continuous-time system or process whose characteristic polynomial is strongly regular in the sense of the p-SCT. When sampling is fast (*P0) the SCT gives "! "P1 and L a P0 for all n (Fan, 1997a). This is the numerical L singularity of the SCT due to fast sampling, which is now overcome by the d-SCT of Fan (1997a). There are two essential steps in the development of the d-SCT: (1) Scaling by 1/* and 1/* to eliminate such singularity; and (2) transformation to provide a meaningful limit as *P0. In other words, numerical singularity of the SCT due to fast sampling can be eliminated by scaling, i.e., if scaled properly as in Eqs. (21)}(25) of Fan (1997a), a as L obtained by Eqs. (1.4) and (1.5) do not vanish as *P0 when the underlying s-domain characteristic polynomial is strongly regular in the sense of p-SCT. Therefore, the kind of singularities that will be dealt within this paper is such that a "0 regardless of the L speed of sampling (without taking the limit *P0), and can be called &&algorithmic'', as opposed to &&numerical''. Speci"cally, two types of algorithmic singularities will be considered. In type I (or case I) a "a "2"a "0, L L LH b "b "2"b "0, and m "m "2"m "0 L L LH L L LH for some n and some j such that 0(j(n. Clearly Example 1 yields type I singularity. Type II (or case II) is type I with j"n, i.e., a "0, b "0, and n "0 for some L L L n. The strongly regular SCT, d-SCT, and p-SCT of Fan (1997a) as summarized above all fail in case of either Type I or Type II algorithmic singularities. In this paper the algorithmic singular cases will be investigated in the framework of the SCT, d-SCT, and p-SCT. The next three sections will be devoted to algorithmic singularity type I. In Section 2 an order recursive solution for the SCT due to Pal and Kailath (1994) will be scaled and modi"ed. It will then be transformed into the d-domain in Section 3, resulting in a new singular case I d-SCT. Its limit as *P0 naturally yields a new singular case I p-SCT in Section 4. For the singular case II, the traditional approach of taking derivative in the z-domain will be used to yield a similar method in the c-domain in Section 5. The traditional approach in the s-domain, however, is to use the original polynomial plus its derivative. This seeming discrepancy with the z-domain approach is easily reconciled by taking the limit of the new c-domain method. Numerical examples will be given to demonstrate better numerical behavior of the new tests over the traditional singular case SCT that
794
H. Fan/Automatica 35 (1999) 791}807
is based on the shift-operator, for fast sampling. Real polynomials (those with real coe$cients) will be assumed throughout this paper. Generalization to the complex case should not be too di$cult.
singularity. At the last stage of Eqs. (2.4) and (2.5) when there are no zeros at either end of g ,a and (\R( (\R( a are recovered by (\R(
2. Modi5ed singularity type I SCT
1 1 a " +g #h ,, a " +g !h ,. (\R( 2 (\R( (\R( (\R( 2! (\R( (\R( ( (2.8)
The singularity type I approach of Pal and Kailath (1994) for the (unscaled) SCT is "rst presented here. Assume that in either Eqs. (1.4) or (1.5) singularity type I is encountered at n#1"J, i.e.,
Then the regular SCT can resume. The zero distribution for a is counted by (
g Oa #! a , h Oa !! a , ( ( ( ( ( ( ( (
where n is the number of roots that are on the unit circle. If J(N, i.e., a is not at the beginning stage, then ( Theorem 2 of Stoica and Moses (1992) should be used together with Eq. (2.9) in determining zero distribution of a. To eliminate numerical singularity due to fast sampling and to facilitate later development, the above Pal}Kailath recursion needs to be scaled and modi"ed. The new version is as follows:
(2.1)
where the (J#1)-dimensional vector g would have ( t zeros at the top, i.e., g "0 for j"0, 1, 2, t !1. ( ( H ( When J"N, a is used in place of a per Eq. (1.4). By the ( de"nition of ! in Eq. (1.6), zero entries on top of g are ( ( possible only when ! "$1. It then follows that ( g "$g . In other words, g must also have t zeros at ( ( ( ( the bottom, i.e., g "0 for j"J!t #1, J!t # ( H ( ( 2, 2, J. Since only singularity type I is considered here, one then must have J J even: t 4 , ( 2
J!1 J odd: t 4 . ( 2
(2.2)
The Pal}Kailath recursion (Pal and Kailath, 1994) can then be written in the following vector form. For j"J, J!2, 2, J!2t #2: ( h h q " HH "! H , H g g H H\RH H RH
(2.3)
0 "g , g H H\ 0 0
(2.4)
+n\, n, n>,"+t , 0, t ,#+those of a ,, ( ( (\R(
1 g " g if J(N, g "g if J"N, ( * ( ( ( h "h . ( (
(2.10)
For j"J, J!2, 2, J!2t #2, ( hM q "! H , H gN H RH
t "t !1, H\ H
(2.5)
(2.6)
where 0 denotes a 2t -dimensional zero vector, and RH H RH g is a ( j!2t #1)-dimensional vector de"ned as RH H H RH g "[g , g , ,g ]2. RH H H RH H RH> 2 H H\RH
(2.7)
Obviously RH g "R( g for all j, due to Eqs. (2.4) and RH H R( ( (2.6). The second equality of Eq. (2.3) holds since ! "$1 and hence g "$g and h "Gh . The reH H H H H cursions of Eqs. (2.4) and (2.5) eliminate all the top zero entries of g , one at a time, thus overcoming the (
(2.11)
0
1 g " g , H\ * H 0 0
RH g 0 RH , h "h #q RH H !q H\ H H 0 H RH g H H H R R 0
(2.9)
(2.12)
1 q h " h # H H\ * H * 0
RH g t RH H ! H 0 1 RH
#2#(!1)RH
q ! H *
t H g H t H
0 t RH ! H RH g 1 RH H
#2#(!1)RH
0 0
RH g RH H RH\
0
RH\ RH g RH H 0
t H g . H t H
(2.13)
795
H. Fan/Automatica 35 (1999) 791}807
At the last stage a
(\R(
and a are recovered by (\R(
1 a " +*g #h ,, (\R( (\R( (\R( 2 1 a " +*g !h ,, if J(N, (\R( 2! (\R( (\R( (
1 h 2" [8 !1 1 8], *
q "8,
1 h 2" [7 7] *
E out of singularity,
1 " +g !h ,, a (\R( (\R( (\R( 2! (
if J"N.
(2.14)
Otherwise the new version is the same as the Pal}Kailath recursion. Similar to the scaled SCT of Fan (1997a), it is easily seen that the scaling factors 1/* and 1/* in Eqs. (2.10)}(2.14) only scale g by 1/*J> or 1/*J, (\J h by 1/*J, and q by * for J(N from their unscaled (\J H versions. Thus there is no essential change from the Pal}Kailath recursion in this kind of scaling except that the numerical singularity due to fast sampling is eliminated. Note that Eqs. (2.10)}(2.14) have no meaningful limit form as *P0. Thus, a linear transformation will be performed in Section 3. The additional terms of the g vector in Eq. (2.13) would facilitate such a transformaH tion. However, it is due to these additional terms in Eq. (2.13) that the equivalence of the new version (2.10)}(2.14) for *"1 with the Pal}Kailath recursion (2.1)}(2.8) is not obvious. The following lemma establishes such equivalence. Lemma 1. ¸et *"1 in Eqs. (2.10)}(2.14). ¹hen in the modi,ed Pal}Kailath recursions g "g for all j"J, H H J!2, 2, J!2t and h "h where g and h are ( (\R( (\R( H H given in Eq. (2.1). See Appendix A for a proof. This lemma establishes that, other than the scaling, t , q , and a are all ( ( (\R( unchanged from the Pal}Kailath recursion. Therefore the newly modi"ed version is equivalent to Pal}Kailath recursion. The following example is taken from Pal and Kailath (1994). Example 2. f (z)"z#z#z#2z#z#1, ! "!1,
g 2"[0 0 !1 1 0 0], t "2, E singularity type I, h 2"[2 2 3 3 2 2],
t "1,
1 g 2" [!1 1], *
1 a " +g #h ,, (\R( 2 (\R( (\R(
a2"[1 1 1 2 1 1],
1 g 2" [0 !1 1 0], *
q "2,
1 a2" [3 4], *
1 a 2" [4 3]. *
Since t "2 and "! "'1, n\"2 and n>"3. This is identical (other than the scaling) to Pal and Kailath (1994). Before leaving this section, we introduce the following interesting properties of the Pal}Kailath recursion which justi"es the use of ! in Eqs. (2.8) and (2.14) and will also ( be useful in Section 4. Lemma 2. In the Pal}Kailath recursions (2.1)}(2.7), de,ne a g "a #! a , h "a !! a , ! "! HH H H H H H H H H H a H
(2.15)
for all j. ¹hen 1 a " +g #h ,, H 2 H H
1 a " +g !h ,, H 2! H H H
! "! , j"J!2, 2, J!2t #2, H ( ( *P0
! && ! . (\R( (
(2.16) (2.17) (2.18)
See Appendix B for a proof. Note that Eqs. (2.15) and (2.16) are valid for any j including both those in the singularity I block and in the regular block.
3. Singularity type I d-SCT The linear transformation described in Section 1 can be applied to the modi"ed Pal}Kailath recursion of Section 2 to obtain a singularity type I recursion for d-SCT. First de"ne u "*\H¹\ g , w "*\H¹\ h H H H H H H
(3.1)
which are consistent with Eq. (1.11). Note that there are t zeros at the top of u for each j, but not at the bottom of H H u since ¹\ is lower-triangular. Left multiplying H H *\H¹\ to Eq. (2.12) and (2.13) one obtains the following H result.
796
H. Fan/Automatica 35 (1999) 791}807
Theorem 1 (Singular I d-SC¹). ¸et J be the order index at which the d-SC¹ becomes singular (¹ype I), i.e., b "0. If J"N, de,ne (\ u "b#! bC, w "b!! bC. ( ( ( (
(3.2)
If J(N, de,ne
1 ! 0 " 1# ( u * ! ( (>
*! b 0 ( ! ( b ! 0 (> (
*! b ( b # ( , (> ! b (> (>
(3.4)
0 #* u
0 u H\ 0
H\ 0
"u , H
#* w
w H\
0
0
b
#*
(\R(\
(3.10)
b 1#! ! (\R(\ " ( (\R( u (\R( 2 0
1!! ! ( (\R( w # , (\R( 2
RH Hu H 0 H\RH
0 RH RHu H 0 H R
,
0
b
(3.11)
(\R(\
#*
(3.12)
b 1#! ! (\R(\ " ( (\R( u (\R( 0 2
1!! ! ( (\R( w # . (\R( 2*
(3.13)
¹hen the regular d-SC¹ resumes. ¹he zero distribution of b is given by ( +n\, n, n>,"+t , 0, t ,#+those of b ,. ( ( (\R(
(3.6)
where lu is u cutting out the ,rst l entries, i.e., H H (3.7)
and t q "! H , H u H RH
b " +u #w ,, (\R( (\R( (\R(
0 0 RH RH #*j RH>u #* j RH>u H H H H 0 0 H H R > R >
u2"[u l, u l , 2, u ] H H H > H H
(3.9)
=hen u and w are obtained at the end of the (\R( (\R( recursions, b and b are recovered by, for (\R( (\R(\ J"N,
b " +*u #w ,, (\R( (\R( (\R(
RHu H !(!1)RH q "w #q H H 0 H RH
l
and for J(N,
H\ 0
#2#*H\RH j H H\RH
t J j "(!1)J H # (!1) G KG m G HJ 1 G GKGIGJ O KGIG J t K t K t H ; H 2 ! H . k k l
(3.5)
0
0
t t H ! H , 2 3
where u has t zeros at its top, and Eq. (2.2) is satis,ed ( ( for t . ¹he singular I d-SC¹ recursions are given by, for ( j"J, J!2, 2, J!2t #2, ( 0
t t j " H ! H , H 1 2
$ (3.3)
! 0 " 1! ( w ! (> (
t j "! H , H 1
t t j "! H #2 H H 1 1
! 0 b ( # ( b ! 0 ( (>
! b ( b ! ( , (> ! b (> (>
and the j l's are pre-calculated constants and are given by H
(3.8)
(3.14)
See Appendix C for a proof. When J(N, Theorem 2 of Stoica and Moses (1992) with a replaced by b and a replaced by b should L L be used together with Eq. (3.14) to determine the zero distribution of b. It is also important to recognize that Eq. (3.3) has exactly the same form as Eq. (1.14). In other words, in executing Eq. (1.14) for the regular case if the RHS gives more than one zero at the top, then singular type I is naturally entered as in Eq. (3.3). The total amount of computation for the singular block is roughly O((J!2t ) t #2Jt ) #ops. The following example illus( ( ( trates the numerical advantage gained using the new
797
H. Fan/Automatica 35 (1999) 791}807
d-version against the original Pal}Kailath recursion for fast sampling. Example 1 (Continued).
"ve decimal places would the correct conclusion be reached. This kind of improvement in numerical behavior is similar to that of Fan (1997a) and is expected to become more dramatic for a smaller *.
f (z)"z!4.999z#10.001z!10.002z#5.001z!1, (3.15) F(c)"c#0.004c#0.08c#0.448c#1.024c #1.024, *"0.25.
(3.16)
Firstly the new test in the c-domain is illustrated. Only four decimal places after the decimal point are saved in the calculation. b2"[1 0.004 0.08 0.448 1.024 1.024], bC2"[!1 0.004 0.032 0.064 0.256 1.024], ! "1E singularity, J"5"N, u2"b2#! bC2 "[0 0.008 0.112 0.512 1.28 2.048], t "1 w2"[2 0 0.048 0.384 0.768 0], q "!250 u2"[0.032 0.32 0.768 2.048], [0 0 w2 ]#*[0 w2 0] "[0 !36 !208 !512 !1024 0], w2"[!144 !256 !1024 0], E out of singularity, b2"[!71.984 !127.84 !511.616 1.024], ! "1.0004, [0 b2]#*[b2 0] "[0.0608 0.3713 0.9730 2.0484], b2"[0.2561 0.4609 2.0484], ! "!1.05, [0 0 b2]#*[0 b2 0] "[!0.0025 !0.0615 !0.3484 !0.41], b2"[!0.2459 !0.41], ! "0.5832, [0 0 b2]#*[0 b2 0] "[0 !0.4554 !1.8215], b2"!1.8215. Applying Eq. (1.18) to b , it is seen that in the sequence +b , b , b , there are two entries having the same sign as b . Thus b has n\"2 and n>"1. Finally, Eq. (3.14) yields that F(c) or b have n\"3 and n>"2. Now if one uses the original (unscaled) Pal}Kailath recursion on f (z) of Eq. (3.15) and still saves four decimal places after the decimal point, an erroneous conclusion will be reached (n\"1, n>"4). Only by saving at least
4. Singularity type I p-SCT Similar to the regular case, an s-domain version of singularity type I can be obtained by taking the limit of the singularity type I d-SCT as *P0. This limit shall be called the singularity type I p-SCT. It needs to be pointed out that in the limit as *P0 the rule of counting top zeros of Eq. (2.6) is to be changed. This can be seen by taking the limit of (C.1). The new rule is t "t !2 and t "t !2l. (4.1) H\ H H\J H Taking into consideration that t will always be odd (Pal ( and Kailath, 1992), it follows that t will be odd for all j. H Thus assuming b Pn , u Pl , and w Pm as *P0, H H H H H H taking the limit of Eqs. (3.5) and (3.6) yields the following result. Theorem 2 (Singular I p-SCT). ¸et J be the order index at which the p-SC¹ becomes singular (¹ype I), i.e., n "0. If J"N, i.e., if m "0, de,ne (\ l "[0, 2m , 0, 2m , 0, 2]2, (4.2) ( m "[2m , 0, 2m , 0, 2]2. (4.3) ( If J(N, de"ne
0 n m " ( n ! ( , (4.4) (> l m 0 ( (> 0 0 "2 , (4.5) m n ( ( where l has t zeros at its top, and t satis,es ( ( ( J even: t !1(J, J odd: t (J. (4.6) ( ( ¹he singular I p-SC¹ recursions are given by, for j"J, J!2, J!4, 2, J!t #3, ( 0 0 l H\ 0
"l , H
(4.7)
0 RH RHl H #q RHl , 0 "m #q H H 0 H H RH m 0 H\ RH l q "! H , H k H RH t "t !2, H\ H
(4.8)
(4.9) (4.10)
798
H. Fan/Automatica 35 (1999) 791}807
where RHl is l cutting out the ,rst 2t entries. At the end H H H of the above recursions, execute Eqs. (4.8) and (4.9) one more time for j"J!t #1 to obtain m . ¹hen ( (\R(\ n and n are recovered by (\R( (\R(\ n " m , (\R(\ (\R(\
(4.11)
(4.12)
0
n
(\R(
"l . (\R(>
¹hen the regular p-SC¹ can resume. ¹he zero distribution of n is given by: ( for q (0 and (t !1)/2 even, or q '0 and (t !1)/2 ( ( ( ( odd,
t !1 t !1 #1, 0, ( #+those of n ,; +n\, n, n>," ( (\R( 2 2 (4.13) for q '0 and (t !1)/2 even, or q (0 and (t !1)/2 ( ( ( ( odd,
t !1 t !1 +n\, n, n>," ( , 0, ( #1 #+those of n ,, (\R( 2 2 (4.14)
where n\, n, and n> are the numbers of the left-side plane zeros, the jw-axis zeros, and the right-side plane zeros, respectively. Note that there are t zeros at the top of l with a H H reduction of two zeros for each iteration. Thus (t !1)/2 ( replaces many functions of t in the previous sections. ( Also note that in case J(N, Theorem 2 of Stoica and Moses (1992) with a replaced by m and a replaced by L L m should again be used to count the zero distribution of n. Furthermore, t "1 is not a singular case, i.e., singu( larity is declared only for t "3, 5, 7, etc. The following ( example illustrates the utility of this new algorithm. Example 4. F(s)"s#s#s#s#s#s#s#2s#1, (4.15) n2"[1 0 1 0 1 0 1 0 1], n2"[1 0 1 0 1 0 2 0], m n2![n2, 0]"[0 0 0 0 0 0 !1 0 1 ], m GFFFHFFFI l2 singularity type I, t "5: m2"[2 0 2 0 2 0 4 0], l2"[0 0 0!1 0 1], m2"[4 0 2 0 4 0],
q "2,
t "3, q "4,
l2"[0 !1 0 1], t "1, m2"[6 0 4 0], q "6, m2"[16 0], E out of singularity, n2"[!1 0 1], n2"[8 0], regular: n "!8. Since there is no sign change from m to m , n\"1 and n>"0 for that part. For the singular block, since q '0 and (t !1)/2 even, Eq. (4.14) gives n\"2 and n>"3. Finally, the sequence +m , m , m , yields n\"1 and n>"1. Thus the root distribution of n is n\"4 and n>"4.
5. Singularity type II When a "b "n "0, singularity type II is encoun( ( ( tered. All regular as well as singular type I methods fail in this case. The traditional solution in the z-domain is to replace f (z) by its derivative f (z) since n> is the (> (> same for both f (z) and f (z) (Pal and Kailath, 1994; (> (> Stoica and Moses, 1992). On the other hand, the traditional solution in the s-domain is to replace F (s) by (> F (s)#F (s) since they have the same n> (Pal and (> (> Kailath, 1992). This seeming discrepancy can be reconciled by the d-domain version as a bridge. First de"ne the parameters of f (z) by a and note the relationship (> ( a "D a , (5.1) ( (> (> where
0
J#1
$
J
D " (>
\
0
0 $
1
.
(5.2)
0
Next transform Eq. (5.1) into the d-domain by left multiplying *\(¹\ to both sides and de"ne b "*\(¹\a : ( ( ( ( b "+*¹\D ¹ , +*\(\¹\ a ,"D b , ( ( (> (> (> (> (> (> (5.3) where the identity *¹\D ¹ "D is easily veri( (> (> (> "ed by straight forward multiplication using the de"nition (1.10). Identity (5.3) implies that b is precisely the ( coe$cients of F (c), the derivative of F (c). There(> (> fore for both the z-domain and the c-domain, if singularity type II is encountered one would then work with a or b , respectively, and continue the regular case ( ( recursion and count only n> a , the number of the outside (Y unit circle roots of a , or n> b , respectively. The root (Y ( distribution of a is then (Stoica and Moses, 1992) n\ a "n> a #s\, (Y
n> a "n> a #s>, (Y
na "J#1!2n> a , (Y (5.4)
799
H. Fan/Automatica 35 (1999) 791}807
where s\(s>) is the number of entries in +a , ,\ a , 2, a , having the same sign as a (!a ). ,\ (> The same holds true for b by replacing all quantities of a by those of b. However, since the regular case d-SCT recursion depends on two previous stages, this aspect needs further examination. Let the parameter vectors of the subsequent stages after b be b , b , etc. To calculate b one cannot ( (\ (\ (\ use (1.14) since the relationship between b and b (> ( is not due to the regular d-SCT. Instead the d-SCT has to be restarted based on b using a scaled version of ( Eq. (1.12) since J(N (see Fan, 1997a for details), i.e.,
0
#*
1 b (\ " +b #! b C,. ( ( 0 * (
b (\ To calculate b C, de"ne ( 0 1
$ DC " (> $ 0
0
2
.
\
0
(5.5)
(5.6)
J#1
Then b C can be expressed by ( b C"*\(¹\a "*\(¹\DC a ( ( ( ( (> (> "+*¹\DC ¹ , +*\(\¹\ a , ( (> (> (> (> (5.7) "+!D #*DC ,bC , (> (> (> where the equality of the "rst brackets in the last step is again easily veri"ed by straight forward multiplication. Substitute Eq. (5.7) into Eq. (5.5) and observe that, since b #! bC "0, bC "(!1/! )b , (> (> (> (> (> (> 0 b #* (\ b 0 (\ 1 ! ( +D !*DC , b " b # ( (> (> (> * *! (> 1 ! ! , (5.8) " 1# ( b ! ( DC b ( ! (> (> * ! (> (> where Eq. (5.3) is used. Furthermore, using Eq. (5.3) and writing the vectors componentwise, it can be veri"ed that
b b 0 ( b ! ( " C . (5.9) (> b 0 D b (> (> (> Substitute Eq. (5.9) into Eq. (5.8) to "nally obtain 0 0 b (\
0
1 ! #* b " 1# ( (\ * ! (> 0
! b ! b ( ! ( ( b # ( . (> ! 0 ! b (> (> (>
0 b (
(5.10)
It is interesting to see that Eq. (5.10) has exactly the same form as Eq. (1.14). Thus if one uses b to replace b "0, ( ( the regular case d-SCT can continue without any interruption. Eqs. (5.3) and (5.4) with a replaced by b, Eq. (5.10), and subsequently Eq. (1.14) comprise the singular type II d-SCT. The partial inertia count n> b is given ( by Eq. (1.17) with b replacing b and b replacing b. G G ( Therefore we have proven the following. Theorem 3 (Singular II d-SCT). In the regular d-SC¹ if b "0 at some stage J, then take derivative of the ( c-polynomial associated with b ,F (c)"b c(> (> (> (> #b c(#2#b c#b , and call the co(> (> ( (> (> e.cient vector of dF (c)/dc as b . ¹he regular d-SC¹ (> ( can then continue with b replacing b . ¹he zero distribu( ( tion of b is given by Eq. (5.4) with all a's replaced by b's, in which n> b is given by Eq. (1.17) with b replaced by b and ( G G b replaced by b . ( To make the singular type II d-SCT a bridge between the SCT and the p-SCT, observe that since b and (> ,a and a "(J#1)a b "(J#1) b ( (> (> ( (> have the same signs, the value of the RHS of Eq. (1.17) with b replaced by b is unchanged by the following G G modi"cation:
n> b "var 1, (
b b b ( , ( (\ , 2, b b b (> (> (
b (> ( b L ( L\ b b (> L L
(5.11)
"n> b
(>
such that Eq. (5.4) becomes n\ b "n> b #s\, (>
n> b "n> b #s>, (>
nb"J#1!2n> . b (>
(5.12) A corresponding version for a is by simply substituting b with a. Note that although n> is as expected, the b "n> b ( (> particular form of n> in Eq. (5.11) is not trivial. This b (> form states that in counting the zero distributions one can treat b , b , 2, etc., as subsequent stages of b , ( (\ (> similar to Eq. (5.10). This point of view is essential in connecting to the p-SCT. The singular type II p-SCT is easily obtained by taking the limit of the d-version as *P0. Since b Pn , it (> (> follows that b Pn , the derivative of n . Eq. (5.3) ( ( (> becomes n "D n , and Eq. (5.10) becomes ( (> (> 0 0 m n 0 P 0 " ( n ! ( , (5.13) (> m 0 (> b n (\ (\
800
H. Fan/Automatica 35 (1999) 791}807
which is the same as Eq. (1.23). In other words, if one replaces n "0 by n , then the regular case p-SCT recur( ( sion continues without any interruption. The zero distribution count would be the same as Eqs. (5.11) and (5.12) with n replacing b. Thus again we have proven the following result. Theorem 4 (Singular II p-SCT). In the regular p-SC¹ if n "0 at some stage J, then take derivative of the ( s-polynomial associated with n , F (s)"m s(> (> (> (> #m s(#2#m s#m , and call the coef(> (> ( (> (> ,cient vector of dF (s)/ds as n . ¹he regular p-SC¹ can (> ( then continue with n replacing n . ¹he zero distribution of ( ( n is given by Eq. (5.4) with all a's replaced by n's, in which n> n is given by Eq. (5.11) with all b's replaced by m's. (
It is important to point out that, since m "0 and L J m "0 for l odd, at any stage j the p-SCT can be viewed K J as a recursive procedure to test zero distribution of an s-polynomial with the coe$cient vector
1 0 n# 2 H n H\
1 " [m , m ,m ,m , ]2. 2 H H\ H H\ 2 (5.14)
This is certainly true for j"N in light of Eqs. (1.21) and (1.22). For j(N, if one applies the rules of Eqs. (1.21) and (1.22) to test the zero distribution of the vector on the RHS of Eq. (5.14), this vector will be separated into n and H n . Thus, the p-SCT is unchanged whether n and H\ H n are obtained as intermediate stages from the p-SCT H\ or as the beginning stages with Eq. (5.14). In other words, the statement leading to Eq. (5.14) is true for all j. In particular, for j"J#1 the singular-type II p-SCT developed above, due to the continuity in Eqs. (5.10) and (5.11), is then to test the zero distribution of an s-polynomial with the parameter vector +n #[n ], and the ( (> zero distribution count (5.12) with s\"s>"0. Clearly, this polynomial is +F (s)#F (s),, which is identi (> (> cal (except a factor of 2) to the traditional solution to singularity-type II (Pal and Kailath, 1992). Thus, the seemingly di!erent approaches in the z-domain and the s-domain for singularity-type II are reconciled. The following example demonstrates the singularitytype II algorithms in all three domains, thus showing the numerical advantages of the d-SCT and p-SCT.
a"[1 !4 5.9999 !3.9998 0.9999]2, ! "!0.9999, a " [0.00019999 !0.00059998 0.00059999 !0.0002]2. Even without going into the Pal/Kailath recursion, it is clear from a that unless one keeps 8 digits after the decimal point, a would be rounded to a "[0.0002 !0.0006 0.0006 !0.0002]2 which corresponds to a triple root at 1. This erroneous result is due to the sensitivity of the SCT to fast sampling. Now perform the d-SCT and keep only four digits after the decimal point. d-SC¹: b"[1 0 !1 0 0]2, bC"[0.9999 !0.02 !1 0 0]2, ! "!0.9999, [0 b2]#*[b2 0] "[0.0002 0.02 !0.0001 0 0]2, b "[0.02 !0.0001 0 0]2, ! "1, [0 0 b2]#*[0 b2 0] "[0 !0.0001 !0.02 0 0]2, b "[!0.02 0 0]2, ! "!1, [0 0 b2 ]#*[0 b2 0]"[0 0.0001 0 0]2, b "0 E singularity II, J"1, b "[!0.04 0]2, ! "1, [0 0 b 2 ]#*[0 b 2 0]"0,
Example 3.
b "0 E singularity II again,
f (z)"(z!1)(z!1.01) (z!0.99) "z!4z#5.9999z!3.9998z#0.9999,
These three equations describe the same system but in three di!erent domains. From Eq. (5.15) it is clear that +n\, n, n>,"+1, 2, 1, which is also true for Eqs. (5.16) and (5.17). Now perform the SCT, d-SCT, and p-SCT with "nite word-length arithmetic. SC¹:
(5.15)
F (c)"c!c, *"0.01, @
(5.16)
F (s)"0.999975s#0.01s!s. K
(5.17)
b "!0.04. The inertia is counted as follows. From b , b , and b it is clear that s\"1 and s>"1. From b and b we have n> b "0. Thus +n\, n, n>,"+1, 2, 1, according to Eq. (5.12).
801
H. Fan/Automatica 35 (1999) 791}807
Now perform the p-SCT with only two digits after the decimal point. p-SC¹: n"[1 0.01 !1 0 0]2, n "[2 0 !2 0 0]2, n "[0.02 0 0 0]2, n "[!0.02 0 0]2, n "0 E singularity II, J"1, n "[!0.04 0]2, n "0 E singularity II again, n "!0.04. The inertia count, and hence the conclusion, is the same as the d-SCT.
Then showing h "h is equivalent to showing (\R( (\R( "0. The recursion for h can be obtained by subh H (\R( tracting Eq. (2.5) from Eq. (2.13), noting that *"1 and RH g "R( g ∀j, RH H R( ( 0
R( g 0 RH "h #qJ R( ( !qJ h H H 0 H R( g H\ H ( ( R R 0
t #q ! H H 1
t H ; t !1 H
6. Conclusion Singular cases for the recently developed d-SCT and p-SCT are considered in this paper. &&Numerical singularity'' due to fast sampling is distinguished from &&algorithmic singularity''. Numerical singularity of the SCT is overcome by the d-SCT. Thus only algorithmic singularity of types I and II are considered in this paper. For type I singularity, a recent order recursive approach of Pal and Kailath is scaled and modi"ed. It is then transformed into the d-domain. Its limit as *P0 yields an s-domain version. Numerical examples are given to show the numerical advantage of the d-version against the original approach of Pal and Kailath for fast sampling. Finally, singularity type II is being tackled. A d-version is obtained from the traditional z-domain method, and an s-domain version results when taking the limit *P0. The resulting s-domain version coincides with the traditional approach, and is thus reconciled and uni"ed with the z-domain approach. A numerical example is also given for singularity type II algorithms.
0 R( g #2#(!1)RH\ R( ( 0 RH\
0 RH\ R( g R( ( 0 RH>
#2#(!1)RH\
t !q ! H H 1
t H t !1 H
0 RH> R( g R( ( 0 RH\
0 RH\ R( g R( ( 0
,
(A.2)
hI [1 0 2 0]h H , qJ "q !q "! H "! H H H [1 0 2 0]R( g g ( H RH
(A.3)
where
and the last terms in the brackets of Eq. (2.13) cancel. It is clear by Eq. (2.10) that h "0, which also implies that ( qJ "0. Thus recursion (A.2) can be evaluated for j"J to ( obtain
0 t h "q ! ( ( (\ 1 0
0 #2#(!1)R(\
R( g R( (
0 R(\ 0 R(\ R( g R( ( 0 R(>
Acknowledgements
t ( ; t !1 (
t !q ! ( ( 1
This work began when the author visited Systems and Control Group, Uppsala University, Uppsala, Sweden in 1995. The author wishes to thank Prof. P. Stoica for bringing up the importance of the singular cases and for providing a reprint of the reference, Stoica and Moses (1992), during his visit.
t ( #2#(!1)R(\ t !1 (
0 R(> R( g R( ( 0 R(\
0 R(\ R( g R( ( 0
.
(A.4)
Using Eq. (A.3) one obtains from Eq. (A.4) Appendix A. Proof of Lemma 1 When *"1, g "g ∀j is obvious comparing Eq. (2.12) H H with Eq. (2.4) and observing Eq. (2.10). De"ne h "h !h , H H H
j"J, J!2, 2, J!2t . (
(A.1)
[1 0 2 0]h [1 0 2 0]R( g (\"! R( ( qJ "! (\ [1 0 2 0]R( g [1 0 2 0]R( g ( (
t t ; ! ( q " ( q . 1 ( 1 (
(A.5)
802
H. Fan/Automatica 35 (1999) 791}807
Use Eqs. (A.4) and (A.5) to carry out recursion (A.2) further. One obtains 0
0
t h "q ! ( ( (\ 1 0 0
#(!1)R(\
0
t ( #2#(!1)R(\ t !1 (
R( g R( (
0 R(\
#(!1)R(\\
#(!1)R(\\
t ( t !1 (
0 R(> R( g R( ( 0 R(\
(\ t !1 (\
#q (\
0 R(\ R( g R( ( 0 R(\>
t (\ t !1 (\
t
0
0 R(\> R( g R( ( 0 R(\
0 R(\>
t t ( q ! (\ q (\ 2 ( 1
0
.
#(!1)
0 R(\ R( g R( ( #2 0
0 R( g R( ( 0 R(\
#2
0
(A.6)
On the other hand, recursion (A.2) can be executed backwards from its last stage, i.e., 0 h "qJ (\R( (\R(> 0
!qJ
R( g R( ( 0 0
(A.8)
t !l#1 ( q , (\J\ 1
where 14l4t !1 and Eq. (2.6) is used. (
(A.9)
(\R(>
0 R( g R( (
, #h (\R(>
(A.10)
where t "t !(t !1)"1 is used. Appending (\R(> ( ( a zero at each end and using Eq. (A.2) again on h , (\R(> Eq. (A.10) becomes 0
0
0
R( g 0 R( ( "q h 0 !q 0 (\R( (\R(> (\R(> R( g 0 0 R( ( 0 0 0 0
R( g R( ( 0
t t !1 ( q #(!1)J ( q #2 ( l l!1 (\
0
.
Continue the recursion, one sees that the terms with qJ in H Eq. (A.2) always cancel the corresponding terms in h of H Eq. (A.2) which are most shifted vectors of R( g for each j. R( ( The next most shifted vectors of R( g contribute to the R( ( "rst non-zero component of h . By Eq. (A.3) we can H\ deduce that
0 R(\> t R( g !q #q ! (\ (\ R( ( (\ 1 0
0
[1 0 2 0]h t (\"(!1) ( q qJ "! (\ [1 0 2 0]R( g 2 ( ( t #(!1) (\ q . (\ 1
qJ "(!1)J> (\J
0 R(\ R( g #2 R( ( 0
(A.7)
R( g R( ( 0 R(\ Using Eq. (A.7) in Eq. (A.3), one obtains
t !q ! ( ( 1
t !q ! (\ (\ 1
Note that in Eq. (A.6) the "rst terms in the "rst two brackets cancel the terms with qJ due to Eq. (A.5). (\ Therefore, the "rst (from the top) non-zero entry of h is due to the following terms: (\ 0 (!1)
R( g R( (
0 R(\ R( g R( ( 0 R(>
#q (\R(>
0 0 0
0 0 !q (\R(>
0 0 R( g R( (
803
H. Fan/Automatica 35 (1999) 791}807
we get
0
R( g R( ( 0
2 #q (!1) (\R(> 1
0
RH g 0 RH , t 51. "2a #q RH H !q 2a H H 0 H RH g H H\ RH RH H 0
0
(B.1)
0
The zeroth and ( j!2)th entry of a can be obtained H\ from Eq. (B.1), using Eq. (2.15), as
0 0
2 !q (!1) (\R(> 1
0
Rg R ( 0
#h . (\R(>
(A.11)
( (
Continue this process and group vectors of the same kind. It can be deduced that the coe$cients of these vectors are given by
0 0 R(\ R(> R( g ! R( g R( ( R( ( 0 0 R(> R(\ q
(\R(>
#(!1)
By the de"nition in Eq. (2.15), ! can be obtained from H\ Eqs. (B.2) and (B.3),
2 t ( q #2#(!1)R(\ qN ( (\R > 1 t !1 ( (
0
$ 0
!
t !l#1 t ( q #2#(!1)J ( qN (\J> 1 l (
$
0 R(\
(B.3)
terms:
q #(!1) (\J
q a "a ! H (a #! a ). H\ H\ H H\ 2 H H\RH\ H H RH>
$
0 0 J R(\J R( g ! R( g R( ( R( ( 0 0 R(\J J
R( g R( (
(B.2)
terms:
$
q a "a # H (a #! a ) H\ H 2 H RH> H H H\RH\
R(\ R( g R( ( 0
terms: q
#(!1) (\
t ( qN , 1 (
where 14l4t !1. However by Eqs. (A.5), (A.8), and ( (A.9) all these coe$cients equal to zero. Since Eq. (A.12) , 02 ]2, it is concluded contains all the terms in [02 , h 2 R( (\R( R( that h "0. (\R(
(A.12)
a !OH (a #! a ) a H H RH> ! "! H\ H\"! H H\ H H\RH\ H\ a #OH (a #! a ) a H H RH> H H H\RH\ H\
#! a ) a !OH (a H H\ H H\RH\ H H RH> "! . H !! a !OH (!a #! a ) H H H H H\RH\ H H RH> (B.4)
Appendix B. Proof of Lemma 2 Eq. (2.16) is easily solved from Eq. (2.15). Thus what needs to be proven is Eqs. (2.17) and (2.18). To this end we need to express Eqs. (2.4) and (2.5) in terms of a . H Thus adding Eq. (2.4) onto Eq. (2.5) and using Eq. (2.16),
For j"J, J!2, 2, J!2t #4, t '1 ( H g "0, i.e., by Eq. (2.15), H H\ a "!! a . H H\ H H
and
hence
(B.5)
804
H. Fan/Automatica 35 (1999) 791}807
Also, since g "0 in this range of j it is easily seen from H Eq. (2.15) that "! ""1, or !"1. Using these facts in H H Eq. (B.4) one easily concludes ! "! , j"J, J!2, 2, J!2t #4 H\ H (
(B.6)
which is exactly Eq. (2.17). To show Eq. (2.18), observe that for j"J!2t #2, t "1 and thus !"1 is still ( P H H P (!1)H> and a (!1)H valid. In addition, ! & & H H H\ a (Vijayan et al., 1991). Thus, although Eq. (B.5) does H not hold for a "nite *, but in the limit it still holds true. Therefore from Eq. (B.4), ! P! . (\R( (\R(>
(B.7)
Eqs. (B.7) and (2.17) imply Eq. (2.18).
Appendix C. Proof of Theorem 1 Left multiplying *\H¹\ to Eq. (2.12) and using the H identities (14) of Vijayan et al. (1991), 02
¹\ ¹ u " H H\ H\ 02
02
02
02 #*
I
I
02
0
0 #* u
0
" u
H\
H\ 0
"u . H
(C.1)
Left multiplication of *\H¹\ to Eq. (2.13) is slightly H more involved due to the shifted RH g terms. Thus let us RH H "rst investigate how those shifted RH g terms can be transRH H formed. Observe from Eq. (2.7), g 2"[02 , RH g 2, 02 ]"[02 , RH g 2 ] H RH RH H RH RH H
(C.2)
RH g 2"[RH g 2, 02 ]. H RH H RH
(C.3)
Thus, the "rst brackets in Eq. (2.13) contain only upward shifted RH g , or RH g . For k4t , one can calculate RH H H H 0 RH\I *\H¹\ RH g " ¹\ H H H 0 I
+*\H¹\ g ,"M u , H H H H
0H\I>;I I H\I> ¹ 0I;I 0I;H\I> H
where M is the product of the three matrices in the "rst H brackets and can be calculated as
0I\;H\I> I
H\I> 02 H\I>
k #2#*I k
* 2 * $ $
I H\I>
* 2 *
0I;H\I>
, (C.5)
where * denotes a &&do-not-care'' entry. Since u has H t zeros at the top, when Eq. (C.5) is multiplied by u the H H &&do-not-care'' entries would have no e!ect on the outcome. Then substitute Eqs. (C.5) into (C.4) to obtain, observing the de"nition of Eq. (3.7), 0 RH\I k *\H¹\ RH g "u #* H H H 1 0 I
u H #2 0
Iu H . 0 I
k k
(C.6)
The "rst brackets in Eq. (2.13) would then be transformed into the following, using Eq. (C.6) for k"t , H t !1, 2, 1, 0, H *\H¹\ H
RH g t H ! H 0 1 RH
0
;
t H t H
RHu H 0 RH
t !1 H t !1 H
t ; H m
u #* H
t H 1
#(!1)
t !m H l
RH\Ku H 0 RH\K
RH\u H 0 RH\
t !m H 1
0 RH RH g H
u t !1 H #2#*J H 0 l
#*RH\
;
t RH g #2#(!1)RH H H t H 0 RH\
u t H #2#*J H l 0
t " u #* H H 1
#*J (C.4)
k 1
I H\I>
* 2 *
#*RH
where
#*
* 2 * ; $ $
#*I
u H\
0I;H\I>
* 2 * M" $ $ H * 2 *
Ju H #2 0 J
u #* H
t !1 H 1
Ju H #2 0 J
#2#(!1)K
u H #2 0
Ju t !m H #2#*RH\K H 0 t !m J H
#2#(!1)RH
t H u. H t H
(C.7)
805
H. Fan/Automatica 35 (1999) 791}807
The terms [Ju2, 02 ]2 can be combined together with H J coe$cients given by, for 04l4t !1, H
R \J t *J(!1)K H m K H
*\H¹\ H
t !m H l
RH\J t! (t !m)! H H "*J (!1)K m!(t !m)! l!(t !m!l)! H H K
*J RH\J t !l " t (t !1)2(t !l#1) (!1)K H "0. H l! H H m K (C.8) That is, all terms on the right hand-side (RHS) of Eq. (C.7), except *RH (RH ) [RHu2, 02 ]2, add up to zero. Thus RH H RH
RHu H . 0 RH
(C.9)
Now turn attention to transforming the second brackets of Eq. (2.13) which contain only downward shifted RH g . RH H De"ne g "[02 , RH g 2 ]2. RH H RH RH H
(C.10)
Then the downward shifts of RH g in the second brackets RH H of Eq. (2.13) become upward shifts of g , i.e., RH H 0 RH\ t RH g #2#(!1)RH H RH H t H 0
0 t RH ! H RH g 1 RH H
0 t " RH ! H g 1 RH H
"(!1)RH
0
#(!1)RH
t H t H
0 RH g RH H
g RH H 0 RH
0
g t RH H ! H 0 1 RH
g #2 RH H 0 RH\ ,
0 RH . g RH H
RH\ RH g #2 RH H 0
0 RH RHg H . H g R #(!1)R "(!1)RH *RH (C.13) RH H 0 H R 0 RH Now g needs to be expressed in terms of u . To do so H H write u as H g u "*\H¹\ g "*\H¹\ RH H H H H H 0 RH 0 I 0 ; H\RH> "¹\ H\RH> RH ¹ *\H¹\ RH H H H g 0RH;RH 0RH;H\RH> RH H g RHg t t H #2#*RH H H , "g #* H (C.14) H 0 1 t 0 H H R where Eq. (C.5) and the fact that g has t zeros at the top H H are used. Actually, since g has t zeros at the top, from RH H H Eq. (C.12) it is seen that g has 2t zeros at the top. The H H non-zero part RHg can be solved from Eq. (C.14) starting H from the last element.
t H t H
u "g , H H H H
t u "g #* H g , H H\ H H\ 1 H H $ t t H g , u "g #* H g #2#*H\J H J H J 1 H J> j!l H H (C.15)
(C.11)
(C.12)
where 2t 4l4j. Note that in Eq. (C.15) we have de"ned H (RH)"0 for k't . RHg can then be solved from Eq. (C.15) I H H as RH>u H RH>u H #*j RHg "RHu #*j 0 #2 H H H H 0 0
Hu H , (C.16) H\RH where Hu "u and j is given by Eq. (3.9). H HH HJ Finally, left multiplying Eq. (2.13) by *\H¹\ and H using Eqs. (C.1), (C.9), (C.13), and (C.16), one obtains Eq. (3.6), where #*H\RH j
where (RH)"( RH ) is used. It is seen that the RHS of RH\J J Eq. (C.11) has the same form as the terms in the brackets of the LHS of Eq. (C.7). Thus de"ne g "*\H¹\ H H
0
0 RH RH g RH H 0 RH
RH\ t g #2#(!1)RH H H R H t H 0
0 t RH ! H RH g 1 RH H
H
*J RH\J (t !l)! H " t (t !1)2(t !l#1) (!1)K H H H l! m!(t !l!m)! H K
LHS of Eq. (C.7)"*RH
One can then use the results of Eqs. (C.4)}(C.9) to obtain
H H\RH 0
*RH [1 0 2 0] h H q "*RHq "! H H [02 , 1, 02 ] g RH H\RH H [1 0 2 0] *H¹ w t H H"! H . "! [02, 1, 02 ] *H\RH ¹ u u RH H\RH H H H RH
(C.17)
806
H. Fan/Automatica 35 (1999) 791}807
The last equality of Eq. (C.17) is true since ¹ is lowerH triangular and u has t zeros at the top. This proves the H H recursions (3.5)}(3.8). Now study how these recursions are entered. When J(N, i.e., singularity is not encountered at the very beginning, scaling by 1/* on g is needed. Left multiply( ing *\(¹\ to Eq. (2.10) and using Eq. (2.1), one obtains ( 1 u " +b #! bC,, w "b !! bC. ( * ( ( ( ( ( ( (
(C.18)
To eliminate bC from Eq. (C.18), note that the stages ( before Eq. (C.18) are regular and that u in Eq. (C.18) has ( the same form as the RHS of the scaled version of Eq. (1.12) as in Eq. (5.5). Thus, Eq. (1.14) applies, i.e., one obtains Eq. (3.3). Similar manipulation on w yields ( Eq. (3.4). On the other hand, if J"N, then g "g and left ( ( multiplication of *\(¹\ to Eq. (2.1) gives Eq. (3.2). ( To exit the singular type I recursion, one also needs to consider the J(N and J"N cases separately. First assume J(N. Then transforming Eq. (2.14) gives Eq. (3.12) and 1 bC " +*u !w ,. (\R( (\R( (\R( 2! (
(C.19)
Since we are now back in the regular case, the transformed version of Eq. (1.5), i.e., the scaled version of Eq. (1.12) as in Eq. (5.5), is in e!ect. Thus using Eq. (C.19) and the scaled version of Eq. (1.12) one has
0
b
#*
(\R(\
b 1 (\R(\ " +b #! bC , (\R( (\R( 0 * (\R(
1#! ! 1!! ! ( (\R( w ( (\R( u # , " (\R( (\R( 2 2*
(C.20)
where ! is given by Eq. (1.15) and the fact that (\R( ! "$1 is used. ( In case J"N, Eq. (C.19) becomes Eq. (3.10) and 1 bC " +u !w ,. (\R( 2! (\R( (\R( (
(C.21)
Accordingly, the factor 1/* on the RHS of Eq. (C.20) needs to be removed, i.e.,
0
b
#*
(\R(\
"b
b (\R(\ 0
Appendix D. Proof of Theorem 2 The recursions (4.7)}(4.9) are immediately obtained by taking the limit (*P0) of Eqs. (3.5), (3.6), and (C.17). Thus we only need to prove how to enter and exit the recursions. To enter the recursions of Eqs. (4.7) and (4.8), "rst consider J(N. Taking the limits of Eqs. (3.3) and (3.4) and using the facts that ! P(!1)H> and H 1/* (1#! / ! )P0 (Fan, 1997a), one obtains Eqs. H H> (4.4) and (4.5). Comparing Eq. (4.4) with Eq. (1.23) it is seen that k "l "0 for i even and that t "1 is not H G H G> ( a singular case. Now consider J"N. Taking the limit of Eq. (3.2) and considering Eqs. (1.12), (1.21) and (1.22), one obtains Eqs. (4.2) and (4.3). To exit the recursions (4.7) and (4.8), consider replacing t by (t !1)/2, i.e., J!2t by J!t #1 at the last stage. ( ( ( ( Note that although when t is reduced to being 1 (l H (\R(> has one zero at the top) it becomes consistent with the regular case, singularity is not exhausted at this step since RHl still cuts out the last zero at the top and Eq. (4.8) can H be executed one more time to produce m . Therefore (\R(\ for J(N, take the limit of Eq. (3.12) at J!t !1, ( Eq. (4.11) is obtained. On the other hand, observe from Lemma 2, Eq. (2.18) that in the limit ! has the same (\R( sign as ! . In other words, for * small ! ! ""! " ( ( (\R( (\R( in Eq. (3.13). Using Eq. (1.15) and the fact that m "0 for L P 0 in Eq. (3.13). Thus replacall n, (1!! ! )/2* & ( (\R( ing t by (t !1)/2, the limit of Eq. (3.13) is given by ( ( Eq. (4.12). Now for J"N the limit of Eq. (3.11) gives again Eq. (4.12). To obtain n , observe that since (\R(\ is still in the singular block, Eq. (3.10) is although n (\R(> valid but has no meaningful limit for J!t #1. On the ( other hand, since l does not exist, Eq. (3.10) does (\R(\ not have a limit for J!t !1 either. Thus n can( (\R(\ not be obtained from Eq. (3.10). This dilemma can be resolved by comparing Eqs. (4.2)}(4.3) with (1.21)}(1.22). Obviously n "m and [0, n2 ]"l2 in this case. If we , ( ,\ ( de"ne n " m and [0, n2 ]"l2 for j"J!2, J!4, H H H\ H 2, J!t(#1 and n(\R(\" m(\R(\, the singular type I recursions (4.7) and (4.8) can then be expressed in terms of n , j"J, J!1, 2, J!t !1. The desired n is H ( (\R(\ then naturally given again by Eq. (4.11). Thus the case J"N is consistent with J(N. The zero distribution rules of Eqs. (4.13) and (4.14) are derived by Delsarte et al. (1985) and used by Pal and Kailath (1992).
References
C
#! b (\R( (\R( (\R(
1#! ! 1!! ! ( (\R( u ( (\R( w " # . (\R( (\R( 2 2
(C.22)
Delsarte, P., Genin, Y. V., & Kamp, Y. G. (1985). Pseudolossless functions with application to the problem of locating the zeros of a polynomial. IEEE ¹rans. Circuits & Systems, 32(4), 373}381.
H. Fan/Automatica 35 (1999) 791}807 Fan, H. (1997a). E$cient zero location tests for delta-operator based polynomials. IEEE ¹rans. Automat. Control, 42(5), 722}727. Fan, H. (1997b). A normalized Schur-Cohn stability test for the delta operator based polynomials. IEEE ¹rans. Automat. Control, 42(11), 1606}1612. Goodwin, G. C., Middleton, R. H., & Poor, H. V. (1992). Highspeed digital signal processing and control. Proc. IEEE, 80(2), 240}259. Jury, E. I. (1982). Inners and stability of dynamic systems (2nd ed). Malabar, FL: Robert E. Krieger Publishing Co. Pal, D., & Kailath, T. (1992). Displacement structure approach to singular root distribution problems: The imaginary axis case. IEEE ¹rans. Circuits Systems I, 41(2), 138}148. Pal, D., & Kailath, T. (1994). Displacement structure approach to singular root distribution problems: The unit circle case. IEEE ¹rans. Automat. Control, 39(1), 238}245. Stoica, P., & Moses, R. L. (1992). On the unit circle problem: The Schur}Cohn procedure revisited, Signal processing, 26(1), 95}118. Vijayan, R., Poor, H. V., Moore, J. B., & Goodwin, G. C. (1991). A Levinson-type algorithm for modeling fast-sampled data. IEEE ¹rans. Automat. Control, AC-36(3), 314}321.
807
H. (Howard) Fan received the B.S. degree from Guizhou University, China, in 1976, and the M.S. and Ph.D. degrees from the University of Illinois, Champaign-Urbana, in 1982 and 1985, respectively, all in electrical engineering. From 1977 to 1978, he worked as a Research Engineer at the Provincial Standard Laboratory and Bureau of Guizhou Province in China. In 1978, he entered the Graduate School of the University of Science and Technology of China and then transferred to the University of Illinois where he was a Teaching Assistant and a Research Assistant from 1982 to 1985. He joined the University of Cincinnati in Electrical and Computer Engineering in 1985 as an Assistant Professor, where he is now a Professor. His research interests are in the general "elds of systems and signal processing, in particular signal representation and reconstruction, adaptive signal processing, signal processing for communications, and system identi"cation. He spent six months in 1994 with the Systems and Control Group, Uppsala University, Sweden, as a visiting researcher. He received the "rst ECE Departmental Distinguished Progress in Teaching Excellence Award at the University of Cincinnati in 1987, and served as an associate editor for IEEE ¹ransactions on Signal Processing from 1991 to 1994. He is a member of Tau Beta Pi and Phi Kappa Phi.