An Order-Recursive Algorithm for Estimating Pole- Zero ... - IEEE Xplore
154
IEEE TRANSACTIONS ON ACOUSTICS, SPEECH, AND SIGNAL PROCESSING. VOL. ASSP-35, NO. 2, FEBRUARY 1987
An Order-Recursive Algorithm for Estimating PoleZero Models Abstract-Thispaper dealswiththepole-zeroestimationofadiscrete-timelinearsystemfromameasuredinput-outputrecord. It is shown that the minimization of the squared equation error €or a recursive (n, n - 1) filter can be implemented by an order-recursivealgorithrn using scalar products of records derived from the data. The algorithm is based on the Gram-Schmidt orthogonalization of an intertwinedKrylovsequencewhichconsists of successivelytime-shifted or circularly time-shifted versions of the input and output records.
R
I.INTRODUCTION ECURSIVE digital (n, n - 1) filters may be characterized by an input-output relationship of the form n
In this paper, an order-recursive algorithm for the minimization of (2) is presented which, as its counterpart in
r31 for the no input does not construct the Gramian explicitly. The algorithm is based On general linear vector space considerations involving an intertwined Krylov [6], [7] sequence and a unitary operator. 11. MAINRESULT In what follows, we work in Z2(Z), the Hilbert space of biinfinite sequences with scalar product +cc
(x, Y> =
n- 1
r=
--m
(3)
4Y:'
and norm 1x1 = (x, x)"'. A is the biinfinite shift operator, defined by
where ut, u,, and e, are, respectively, the input, the output, and an error term which represents effects not ac(A),= x,- 1. (4) counted for. Conventionally, one normalizes the filter by It is clear that A is unitary, since taking a, = 1, but for reasons that will become clear in the next section the monic normalization an = 1 will be (A,,AY) = (x, Y ) (5) adopted. In linear system identification, the objective is for all x , y. to estimate the parameters a k and bk from the given input Defining the monic polynomial p ( z ) and the polynomial and output sequences, and a method commonly used is 4(2) as the least squares technique, i.e., the parameters ak and bk n- 1 n- 1 are found by minimizing the squared equation error n-1
under the constraint an = 1. This results in normal equations with a symmetric block-Toeplitz Gram matrix of dimension 2n. If there is no input, the Gramian is Toeplitz of dimension n and an algorithm of the Levinson [l] type can be used to generate the predictorcoefficients ak for all n. In the block-Toeplitz case, block-predictors can be generated by block-Levinson [2] recursive relationships. In two important papers [3], 141, Cybenko points out that the Levinson-type algorithmsoflinear prediction (here the no input case) make implicit use of the Szego [5] recurrences for polynomials orthogonal over the unit circle. Moreover, the recursive relationships described in [3] do not first construct the Gramian but work directly on the data, without matrix manipulations. Manuscript received October 4, 1985; revised June 3, 1986. The author is with the Laboratory of Electromagnetic~and Acoustics, University of Ghent,Sint-Pietersnieuwstraat 41, B9000 Ghent,Belgium. IEEE Log Number 8610888.
J , = [el2= Ip(~)v + q(A)uI2.
(7)
From elementary linear least squares approximation theory [SI, we know that the minimum of (7) with an = 1 is obtained when the following orthogonality conditions are satisfied, i.e.: (eopt,A k v)
=
(eopt,A ~ U=)
fork
=
o
0, 1,
*
*
,n
-
(8)
1.
-
Defining&, = A ~ aUn d h k c l = fork = 0 , I , * we now consider the intertwined Krylov [7] sequence
(fo,fi,
,
. , A"- Iu, Anv)
(u, u , Av, Au, =
*
*
*
,.hn-t,.hn).
(9)
If the numerator and denominator polynomials of the ( n , - n - 1) filter are relatively prime and theerror term e, vanishes, then it is Seen by inspection Of (l) that the tors (fo, ,f i ) are linearly independent for i up to 2n
0096-3518/87/0200-0154$01.00 O 1987 IEEE
155
FOR ESTIMATING POLE-ZERO MODELS
KNOCKAERT: ORDER-RECURSIVE ALGORITHM
- 1 and are linearly dependent in the limiting case i = Axj - , and A x j - 2. Defining the intermediary quantities 2n only.Hence,itispossibleto perform the Grampj = < A j- xj) (1 8) Schmidt orthogonalization procedure (GSO) [8] on the sequence (A), which yields the set of mutually orthogonal and vectors xi recursively defined by vj = ( h j , xj) ( 19) i-1 19
and keeping in mind that A is unitary, all xj are orthogonal and, since xj - A x j - is in Cj - (xj,A x j -2) = (xj, xj) = wj, it is not difficult to show that the unknown quantities with wk (xk,xk). It is clear that the GSO scheme implies that eopt = x Z n , are the solution of the following linear system: since is orthogonal t o 6 for i = 0, , 2 n - 1 which is exactly what is needed in (8). Note also that the related ( n , n) filter problem of minimizing
,
-
n
under the constraint b, = 1 is solved by , C = -x2,,+ The following theorem shows that the GSO procedure (10) can alternatively be implemented by an order-recursive procedure which takes advantage of the fact thatA is a unitary operator. Theorem: The vectors xi satisfy a six-term recursive relationship of the form xj+2
- Axj
= cY.x. 1 ]+I
- * .,
(12)
with initial values x-1 = x-2 = 0
(13)
x0 = u x1
The scalar
E
=u
(14)
+ EU.
(15)
= O*
(16)
is such that (xo7
uj+2- a. =
cy.
*+ 1 + Pjvjk.
II (A)U + rj(A)u, (22) where I I j (2) and rj(z)are polynomials of degree [j/2], xi =
[ ( j - 1)/2], respectively. Inserting (22) in the recurrence (12), and making explicit that this, in light of (1), should be valid for all choices of u and u , leads to the following corollary. The polynomials IIj(z) and rj(z)satisfy the recursive relationship
x,+, - f i j
=
cYjjxJ+1
+ P j x , + yjzxj-1 + S j f i j - 2 (23)
Proof: Introducing the linear spaces
cj = span {f,,
*
,$}
= span
(x,,
(21)
Looking at (lo), we see that x j can be written as
+ Pjxj + y j h j - 1 + a j h j - 2
j = 0,1,
Note that for j = 0 ( j = l), a 2 (3)-dimensional subset of (20) needs to be solved. Finally, making thescalarproduct of (12) with A x j yields the additional recurrence
* *
, xj> (17) with X - ]
it is readily seen that the vectors xj + , x j , A x j - , A x j - 2 , and^^+^ - A x j are all in E j + l . It is also clear that the vectors xj+ , x j , and xj + - A x j are orthogonal to the subspace AXj - 3 . Since xj- 1, xj are orthogonal to subspace C j - 3 , it follows from the unitary character of A that A x j - , A x j - are orthogonal to AC,-3. Hence, we have found five vectors, all belonging to Cj + l , which are orthogonal to the subspace AXj - 3 . In other words, a space of dimension j 2 contains 5 vectors which are orthogonal to a subspace of dimensionj 2. This is only possible if these vectors are linearly dependent [9]. Since the intertwined Krylov sequence is supposed to be linearly independent, this is only possible if there exists a recursive relationship of the form (12).
= X - 2 = 0, and
x, = 1, x1 = E
(24)
x, = 0 , x,= 1
(25)
for I I j ( z ) , for rj(z). Note that, for the( n , n a, = 1
- 1) filter with monic constraint
+
Q.E.D. Equations for the parameters aj, Pj, y j , and Sj are obtained by taking the scalar products of (12) with xj+ ],xj,
and for the (n, n) filter with monic constraint 6, = 1
P(z) = II 2n + 1 (z),
rzn+ I(z).
The transfer function is in both cases given by
(27)
156
ACOUSTICS, SPEECH, AND SIGNAL PROCESSING, VOL. ASSP-35, NO. 2, FEBRUARY 1987
IEEE TRANSACTIONS ON
as
TABLE I
(nu), = u,-1
Initialization X.
(nu), = uj)-
(v, u ) = - ( u , u)iwo XI = u EU wo =
E
a1
+
1,2,
,N
' * *
- 1
n
P I = (Axe, X I ) Solve wla0 = -pl
-vo
woo0 =
+
Ut =
+
x2 = A X , (YOXI Pox0 w2 = wo + " O P I + P0.0 V I = ('4x1, X I ) ~2 = ( A x , , ~ 2 ) Solve w z a l wzyl= - p z W l P l + PtYl = - V I
C ckzfc k= 1
+ PIPI + WOYl
= AI + "1x2 a 3 = w1 + " l P 2 v2 = (Axz, x,) ~3 = (Ax,, X , )
X3
0
+ PlVl
w,pi
p(z) =
x,+Z
V,+I
= ( A x t + ~ .X ~ + I )
(30)
(3 1)
IT
k= 1
(z
- 2;')
It is clear that the minimization of lei2 can still be implemented with the above algorithm with Q instead of A . Note that z "p(z still represents the pole polynomial independently of N , while the roots of qN(z-') are dependent of the data length. Therefore, formula (28) no longer holds in the cyclical shift case, although the coefficients ck can still be recovered.
AX,-,
-')
( A x ; + , ,X i + d
TABLE 11
1) - 4
1.
n
End Do
Adds 6N-5 15N+8 12(N + $)(n -
-
=0
~ l i ~ , + l
= ai
P;+Z =
Wi-lYi
+ qPi+I + P A
0,+2
N
and
+ w,-26; = 0 2 End Loop + + @,xi + ~ , A x ; - +I
= 2n = AX,
9
n
Doi=2,2n-2 Solve wi+lcli + w,+,yi = - p , + , w ; p i + piyi + w,& = -vi If i
*
where the polynomials p ( z ) and qN(z)are given by
+ PIXI + Y l h O
+ pipi +
0, 1,
e = p(Q)u + qN(n)u = 0, =
Loop
Wj+Icri
l =
If (30) is exact, it is easy to prove [lo] that the modified Kalman equation error
+
w2al
(29)
1.
The impulse response of an ( n , n - 1) filter restricted to a data length N is given by
= (XI, XI)
vo = ('4x0, Xo)
n = l n=2 n >2
t =
=u
Mults 6N, 15N + 19 12(N + 2 ) ( n
Divs 3
- 2) +
11
10 8n - 2
111. DISCUSSION A summary of the algorithm is given in Table I, and a detailed computation count is done in Table 11. It is seen that the computational complexity is O ( n N ) , and hence, for data sequences of fixed length N , the complexity is linear in the system order YE. By comparison, classical methods for orthogonalization [3] have # ( A h 2 ) computational complexities. As has been stated in the Introduction, the algorithm works directly on the data, without first constructing the block-Toeplitz Gramian which in E) Hence,thealitself already requires ~ ( N Y mult/adds. gorithm implements the GSO process, which is equivalent to the Cholesky decomposition of the inverseof the Gramian, with the same order of computational complexity as needed for the construction of the Gramian itself. The algorithm is based on scalar products in Z2(Z) of infinitely long data records. Therefore, as one of the reviewers pointed out, the algorithm's accuracy might be susceptible to truncation of therecords, especially for poles near the unit circle (large time constants). These difficulties can be overcome, at least in the impulse response case (u, = S,,), by replacing the infinite shift operator A with the circular shift operatorQ, defined
EXAMPLES IV. NUMERICAL A. ExampleOne Input: 3-' for t 2 0, 0 f o r t < 0. output: 4-' - 2' -,for t 2 0, 0 for t < 0. Input and output are truncated: u, = u, = 0 for t > 29. The true transfer function is - 8 + 2.6666662 H(z) = 8 - 6z-1 + z - 2
-'
A (2, 1) approximation yields
'
- 8.000004 + 2.6666672 H(z) = 8.000004 - 6.000002~z -2 ' B. Example Two Input: ut = 0 f o r t # 0, u, = 1. Output: 4-'- + 2~ sin (7rt/3)/J?;.for t 2 0, o for t < 0. The impulse response is truncated: u, = 0 for t > 29. is added to u, Weak zero-mean white noise (a = before processing. The true transfer function is
'+
A (2, 1) approximation yields
H(z)
=
0.904803 f 0.5297312-' 3.619209 - 2.405112z-' + z-*
157
KNOCKAERT: ORDER-RECURSIVE ALGORITHM FOR ESTIMATING POLE-ZERO MODELS
REFERENCES
A (3, 2) approximation gives
H(z) = -3.994757 - 1 . 9 9 6 2 3 4 ~ ~ +’0 . 0 0 0 6 8 6 ~ - ~ -15.979004 11.9889432-’ - 5 . 9 9 5 2 2 5 ~ -+ ~z-~’
+
\
C. Example Tu2ree Input: u, = 0 f o r t = 1, 2, * * , N - 1, u, = 1. Output: sin (0.20t) sin (0.21t) f o r t = 0, 1 , * - * , N - 1. Weak zero-mean white noise ( 0 = is added to U, before (circular shift) processing. The true pole polynomial is
+
z 4 - 3 . 9 1 6 1 9 5 ~+~ 5 . 8 3 4 1 4 2 ~-~ 3.9161952, + 1.0. A (4, 3) approximation with N = 30 yields
z 4 - 3 . 8 8 7 6 5 8 ~+~ 5 . 7 4 7 2 6 2 ~ ~
- 3.8269492
+ 0.968994.
[l] N. Levinson, “The Wiener (root mean square) error criterion in filter design and prediction,” J . Math. Phys., vol. 25, pp. 261-278, 1947. [2] P. Delsarte, Y . Genin, and Y. Kamp, “A polynomial approach to the generalized Levinson algorithm based on the Toeplitz distance,”IEEE Trans. Inform. Theory, vol. IT-29, pp, 268-278, 1983. [3] G. Cybenko, “A general orthogonalization method with applications to time series analysis and signal processing,” Math. Comput., vol. 40, pp, 323-336, 1983. “Restrictions of normal operators, Pad6 approximation and au[4] -, toregressive time series,” SIAM J . Math. Anal., vol. 15, July 1984. 151 G.Szego, OrthogonalPolynomials. Providence, RI:AMS Colloquium Publications, XXIII, American Mathematical Society, 1959. [6] A. S. Householder, The Theory of Matrices in Numerical Analysis. New York: Dover, 1975. [7] L. Knockaert, “A recursivealgorithmforlinearsystemidentification,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP34, pp. 492-498, June 1986. [8] P. J . Davis, Interpolation and Approximation. New York: Blaisdell, 1963. [9] P. R.Halmos, Finite Dimensional Vector Spaces. New York: Van Nostrand Reinhold, 1958. [lo] L. Knockaert,“Parametricmodeling of electromagneticwaveforms,” Electromagnetics, vol. 4, pp. 415-430, 1984.
A (4, 3) approximation with N = 300 yields
z4 - 3 . 9 1 6 0 0 1 ~+~ 5.8335562’ - 3.9155972
+ 0.99793.
V . CONCLUSION A O(Nn) order-recursive algorithm for estimating polezero models based on the minimization of the squared equation error has been presented. The derivation makes use of the unitary character of time-shift operators. As with all least squares algorithms, it works well in a lownoise environment. For poles on or close to the unit circle; truncation effects have been avoided by replacing the biinfinite shift operator with the circular shift operator. Applied to a sum of closely related sinusoids promising results are obtained, even for relatively short data lengths.
Luc Knockaert (“81) received the PhysicalEngineer degree and Telecommunications Engineer degree from the State University of Ghent, Belgium, in 1974 and 1977, respectively. From 1979 to 1984 he was with the UniversitC Nationale du Zaire. His current interests are in the signal identification and electromagnetic problem solving areas.
IEEE TRANSACTIONS ON ACOUSTICS, SPEECH, AND SIGNAL PROCESSING. VOL. ASSP-35, NO. 2, FEBRUARY 1987
An Order-Recursive Algorithm for Estimating PoleZero Models Abstract-Thispaper dealswiththepole-zeroestimationofadiscrete-timelinearsystemfromameasuredinput-outputrecord. It is shown that the minimization of the squared equation error €or a recursive (n, n - 1) filter can be implemented by an order-recursivealgorithrn using scalar products of records derived from the data. The algorithm is based on the Gram-Schmidt orthogonalization of an intertwinedKrylovsequencewhichconsists of successivelytime-shifted or circularly time-shifted versions of the input and output records.
R
I.INTRODUCTION ECURSIVE digital (n, n - 1) filters may be characterized by an input-output relationship of the form n
In this paper, an order-recursive algorithm for the minimization of (2) is presented which, as its counterpart in
r31 for the no input does not construct the Gramian explicitly. The algorithm is based On general linear vector space considerations involving an intertwined Krylov [6], [7] sequence and a unitary operator. 11. MAINRESULT In what follows, we work in Z2(Z), the Hilbert space of biinfinite sequences with scalar product +cc
(x, Y> =
n- 1
r=
--m
(3)
4Y:'
and norm 1x1 = (x, x)"'. A is the biinfinite shift operator, defined by
where ut, u,, and e, are, respectively, the input, the output, and an error term which represents effects not ac(A),= x,- 1. (4) counted for. Conventionally, one normalizes the filter by It is clear that A is unitary, since taking a, = 1, but for reasons that will become clear in the next section the monic normalization an = 1 will be (A,,AY) = (x, Y ) (5) adopted. In linear system identification, the objective is for all x , y. to estimate the parameters a k and bk from the given input Defining the monic polynomial p ( z ) and the polynomial and output sequences, and a method commonly used is 4(2) as the least squares technique, i.e., the parameters ak and bk n- 1 n- 1 are found by minimizing the squared equation error n-1
under the constraint an = 1. This results in normal equations with a symmetric block-Toeplitz Gram matrix of dimension 2n. If there is no input, the Gramian is Toeplitz of dimension n and an algorithm of the Levinson [l] type can be used to generate the predictorcoefficients ak for all n. In the block-Toeplitz case, block-predictors can be generated by block-Levinson [2] recursive relationships. In two important papers [3], 141, Cybenko points out that the Levinson-type algorithmsoflinear prediction (here the no input case) make implicit use of the Szego [5] recurrences for polynomials orthogonal over the unit circle. Moreover, the recursive relationships described in [3] do not first construct the Gramian but work directly on the data, without matrix manipulations. Manuscript received October 4, 1985; revised June 3, 1986. The author is with the Laboratory of Electromagnetic~and Acoustics, University of Ghent,Sint-Pietersnieuwstraat 41, B9000 Ghent,Belgium. IEEE Log Number 8610888.
J , = [el2= Ip(~)v + q(A)uI2.
(7)
From elementary linear least squares approximation theory [SI, we know that the minimum of (7) with an = 1 is obtained when the following orthogonality conditions are satisfied, i.e.: (eopt,A k v)
=
(eopt,A ~ U=)
fork
=
o
0, 1,
*
*
,n
-
(8)
1.
-
Defining&, = A ~ aUn d h k c l = fork = 0 , I , * we now consider the intertwined Krylov [7] sequence
(fo,fi,
,
. , A"- Iu, Anv)
(u, u , Av, Au, =
*
*
*
,.hn-t,.hn).
(9)
If the numerator and denominator polynomials of the ( n , - n - 1) filter are relatively prime and theerror term e, vanishes, then it is Seen by inspection Of (l) that the tors (fo, ,f i ) are linearly independent for i up to 2n
0096-3518/87/0200-0154$01.00 O 1987 IEEE
155
FOR ESTIMATING POLE-ZERO MODELS
KNOCKAERT: ORDER-RECURSIVE ALGORITHM
- 1 and are linearly dependent in the limiting case i = Axj - , and A x j - 2. Defining the intermediary quantities 2n only.Hence,itispossibleto perform the Grampj = < A j- xj) (1 8) Schmidt orthogonalization procedure (GSO) [8] on the sequence (A), which yields the set of mutually orthogonal and vectors xi recursively defined by vj = ( h j , xj) ( 19) i-1 19
and keeping in mind that A is unitary, all xj are orthogonal and, since xj - A x j - is in Cj - (xj,A x j -2) = (xj, xj) = wj, it is not difficult to show that the unknown quantities with wk (xk,xk). It is clear that the GSO scheme implies that eopt = x Z n , are the solution of the following linear system: since is orthogonal t o 6 for i = 0, , 2 n - 1 which is exactly what is needed in (8). Note also that the related ( n , n) filter problem of minimizing
,
-
n
under the constraint b, = 1 is solved by , C = -x2,,+ The following theorem shows that the GSO procedure (10) can alternatively be implemented by an order-recursive procedure which takes advantage of the fact thatA is a unitary operator. Theorem: The vectors xi satisfy a six-term recursive relationship of the form xj+2
- Axj
= cY.x. 1 ]+I
- * .,
(12)
with initial values x-1 = x-2 = 0
(13)
x0 = u x1
The scalar
E
=u
(14)
+ EU.
(15)
= O*
(16)
is such that (xo7
uj+2- a. =
cy.
*+ 1 + Pjvjk.
II (A)U + rj(A)u, (22) where I I j (2) and rj(z)are polynomials of degree [j/2], xi =
[ ( j - 1)/2], respectively. Inserting (22) in the recurrence (12), and making explicit that this, in light of (1), should be valid for all choices of u and u , leads to the following corollary. The polynomials IIj(z) and rj(z)satisfy the recursive relationship
x,+, - f i j
=
cYjjxJ+1
+ P j x , + yjzxj-1 + S j f i j - 2 (23)
Proof: Introducing the linear spaces
cj = span {f,,
*
,$}
= span
(x,,
(21)
Looking at (lo), we see that x j can be written as
+ Pjxj + y j h j - 1 + a j h j - 2
j = 0,1,
Note that for j = 0 ( j = l), a 2 (3)-dimensional subset of (20) needs to be solved. Finally, making thescalarproduct of (12) with A x j yields the additional recurrence
* *
, xj> (17) with X - ]
it is readily seen that the vectors xj + , x j , A x j - , A x j - 2 , and^^+^ - A x j are all in E j + l . It is also clear that the vectors xj+ , x j , and xj + - A x j are orthogonal to the subspace AXj - 3 . Since xj- 1, xj are orthogonal to subspace C j - 3 , it follows from the unitary character of A that A x j - , A x j - are orthogonal to AC,-3. Hence, we have found five vectors, all belonging to Cj + l , which are orthogonal to the subspace AXj - 3 . In other words, a space of dimension j 2 contains 5 vectors which are orthogonal to a subspace of dimensionj 2. This is only possible if these vectors are linearly dependent [9]. Since the intertwined Krylov sequence is supposed to be linearly independent, this is only possible if there exists a recursive relationship of the form (12).
= X - 2 = 0, and
x, = 1, x1 = E
(24)
x, = 0 , x,= 1
(25)
for I I j ( z ) , for rj(z). Note that, for the( n , n a, = 1
- 1) filter with monic constraint
+
Q.E.D. Equations for the parameters aj, Pj, y j , and Sj are obtained by taking the scalar products of (12) with xj+ ],xj,
and for the (n, n) filter with monic constraint 6, = 1
P(z) = II 2n + 1 (z),
rzn+ I(z).
The transfer function is in both cases given by
(27)
156
ACOUSTICS, SPEECH, AND SIGNAL PROCESSING, VOL. ASSP-35, NO. 2, FEBRUARY 1987
IEEE TRANSACTIONS ON
as
TABLE I
(nu), = u,-1
Initialization X.
(nu), = uj)-
(v, u ) = - ( u , u)iwo XI = u EU wo =
E
a1
+
1,2,
,N
' * *
- 1
n
P I = (Axe, X I ) Solve wla0 = -pl
-vo
woo0 =
+
Ut =
+
x2 = A X , (YOXI Pox0 w2 = wo + " O P I + P0.0 V I = ('4x1, X I ) ~2 = ( A x , , ~ 2 ) Solve w z a l wzyl= - p z W l P l + PtYl = - V I
C ckzfc k= 1
+ PIPI + WOYl
= AI + "1x2 a 3 = w1 + " l P 2 v2 = (Axz, x,) ~3 = (Ax,, X , )
X3
0
+ PlVl
w,pi
p(z) =
x,+Z
V,+I
= ( A x t + ~ .X ~ + I )
(30)
(3 1)
IT
k= 1
(z
- 2;')
It is clear that the minimization of lei2 can still be implemented with the above algorithm with Q instead of A . Note that z "p(z still represents the pole polynomial independently of N , while the roots of qN(z-') are dependent of the data length. Therefore, formula (28) no longer holds in the cyclical shift case, although the coefficients ck can still be recovered.
AX,-,
-')
( A x ; + , ,X i + d
TABLE 11
1) - 4
1.
n
End Do
Adds 6N-5 15N+8 12(N + $)(n -
-
=0
~ l i ~ , + l
= ai
P;+Z =
Wi-lYi
+ qPi+I + P A
0,+2
N
and
+ w,-26; = 0 2 End Loop + + @,xi + ~ , A x ; - +I
= 2n = AX,
9
n
Doi=2,2n-2 Solve wi+lcli + w,+,yi = - p , + , w ; p i + piyi + w,& = -vi If i
*
where the polynomials p ( z ) and qN(z)are given by
+ PIXI + Y l h O
+ pipi +
0, 1,
e = p(Q)u + qN(n)u = 0, =
Loop
Wj+Icri
l =
If (30) is exact, it is easy to prove [lo] that the modified Kalman equation error
+
w2al
(29)
1.
The impulse response of an ( n , n - 1) filter restricted to a data length N is given by
= (XI, XI)
vo = ('4x0, Xo)
n = l n=2 n >2
t =
=u
Mults 6N, 15N + 19 12(N + 2 ) ( n
Divs 3
- 2) +
11
10 8n - 2
111. DISCUSSION A summary of the algorithm is given in Table I, and a detailed computation count is done in Table 11. It is seen that the computational complexity is O ( n N ) , and hence, for data sequences of fixed length N , the complexity is linear in the system order YE. By comparison, classical methods for orthogonalization [3] have # ( A h 2 ) computational complexities. As has been stated in the Introduction, the algorithm works directly on the data, without first constructing the block-Toeplitz Gramian which in E) Hence,thealitself already requires ~ ( N Y mult/adds. gorithm implements the GSO process, which is equivalent to the Cholesky decomposition of the inverseof the Gramian, with the same order of computational complexity as needed for the construction of the Gramian itself. The algorithm is based on scalar products in Z2(Z) of infinitely long data records. Therefore, as one of the reviewers pointed out, the algorithm's accuracy might be susceptible to truncation of therecords, especially for poles near the unit circle (large time constants). These difficulties can be overcome, at least in the impulse response case (u, = S,,), by replacing the infinite shift operator A with the circular shift operatorQ, defined
EXAMPLES IV. NUMERICAL A. ExampleOne Input: 3-' for t 2 0, 0 f o r t < 0. output: 4-' - 2' -,for t 2 0, 0 for t < 0. Input and output are truncated: u, = u, = 0 for t > 29. The true transfer function is - 8 + 2.6666662 H(z) = 8 - 6z-1 + z - 2
-'
A (2, 1) approximation yields
'
- 8.000004 + 2.6666672 H(z) = 8.000004 - 6.000002~z -2 ' B. Example Two Input: ut = 0 f o r t # 0, u, = 1. Output: 4-'- + 2~ sin (7rt/3)/J?;.for t 2 0, o for t < 0. The impulse response is truncated: u, = 0 for t > 29. is added to u, Weak zero-mean white noise (a = before processing. The true transfer function is
'+
A (2, 1) approximation yields
H(z)
=
0.904803 f 0.5297312-' 3.619209 - 2.405112z-' + z-*
157
KNOCKAERT: ORDER-RECURSIVE ALGORITHM FOR ESTIMATING POLE-ZERO MODELS
REFERENCES
A (3, 2) approximation gives
H(z) = -3.994757 - 1 . 9 9 6 2 3 4 ~ ~ +’0 . 0 0 0 6 8 6 ~ - ~ -15.979004 11.9889432-’ - 5 . 9 9 5 2 2 5 ~ -+ ~z-~’
+
\
C. Example Tu2ree Input: u, = 0 f o r t = 1, 2, * * , N - 1, u, = 1. Output: sin (0.20t) sin (0.21t) f o r t = 0, 1 , * - * , N - 1. Weak zero-mean white noise ( 0 = is added to U, before (circular shift) processing. The true pole polynomial is
+
z 4 - 3 . 9 1 6 1 9 5 ~+~ 5 . 8 3 4 1 4 2 ~-~ 3.9161952, + 1.0. A (4, 3) approximation with N = 30 yields
z 4 - 3 . 8 8 7 6 5 8 ~+~ 5 . 7 4 7 2 6 2 ~ ~
- 3.8269492
+ 0.968994.
[l] N. Levinson, “The Wiener (root mean square) error criterion in filter design and prediction,” J . Math. Phys., vol. 25, pp. 261-278, 1947. [2] P. Delsarte, Y . Genin, and Y. Kamp, “A polynomial approach to the generalized Levinson algorithm based on the Toeplitz distance,”IEEE Trans. Inform. Theory, vol. IT-29, pp, 268-278, 1983. [3] G. Cybenko, “A general orthogonalization method with applications to time series analysis and signal processing,” Math. Comput., vol. 40, pp, 323-336, 1983. “Restrictions of normal operators, Pad6 approximation and au[4] -, toregressive time series,” SIAM J . Math. Anal., vol. 15, July 1984. 151 G.Szego, OrthogonalPolynomials. Providence, RI:AMS Colloquium Publications, XXIII, American Mathematical Society, 1959. [6] A. S. Householder, The Theory of Matrices in Numerical Analysis. New York: Dover, 1975. [7] L. Knockaert, “A recursivealgorithmforlinearsystemidentification,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP34, pp. 492-498, June 1986. [8] P. J . Davis, Interpolation and Approximation. New York: Blaisdell, 1963. [9] P. R.Halmos, Finite Dimensional Vector Spaces. New York: Van Nostrand Reinhold, 1958. [lo] L. Knockaert,“Parametricmodeling of electromagneticwaveforms,” Electromagnetics, vol. 4, pp. 415-430, 1984.
A (4, 3) approximation with N = 300 yields
z4 - 3 . 9 1 6 0 0 1 ~+~ 5.8335562’ - 3.9155972
+ 0.99793.
V . CONCLUSION A O(Nn) order-recursive algorithm for estimating polezero models based on the minimization of the squared equation error has been presented. The derivation makes use of the unitary character of time-shift operators. As with all least squares algorithms, it works well in a lownoise environment. For poles on or close to the unit circle; truncation effects have been avoided by replacing the biinfinite shift operator with the circular shift operator. Applied to a sum of closely related sinusoids promising results are obtained, even for relatively short data lengths.
Luc Knockaert (“81) received the PhysicalEngineer degree and Telecommunications Engineer degree from the State University of Ghent, Belgium, in 1974 and 1977, respectively. From 1979 to 1984 he was with the UniversitC Nationale du Zaire. His current interests are in the signal identification and electromagnetic problem solving areas.
