Statistical Analysis of the Overdetermined Yule-walker Estimator of Frequencies of Multiple Sinusoids

Copyright © IFAC Identification and System Parameter Estimation, Beijing, PRC 1988

STATISTICAL ANALYSIS OF THE OVERDETERMINED YULE-WALKER ESTIMATOR OF FREQUENCIES OF MULTIPLE SINUSOIDS P. Stoica*, T. SOderstrom**, F-N. Ti*** and B. Friedlandert *Department of Automatic Control, Polytechnic Institute of Bucharest, Splaiul Irulependentei 313, R-77206 Bucharest, Romania **Automatic Control and Systems, Analysis Group, Department of T echnology, Uppsala University, p.a. Box 534, S-751 21 Uppsala, Sweden ***Automatic Control and Systems Analysis Group, Department of Technology , Uppsala University, p.a. Box 534, S-751 21 Uppsala, Sweden , 01/ leave fmm Department of Mechanical and Electrical Engineering, Shandong College of Light Industry , Jinan , PRC tSaxfrY Computer Corporation, 255 San Geronimo Way, SUlznyvale, CA 94086, USA

Abstract . The estiaation of sinusoidal frequencies by the overdetermined Yule -Walker (OTV) .ethod is considered. An explicit expression is derived for the asymptotic covariance aatrix of the esti.ation errors . The effect of increasing the number of YuleWalker equations on estiaation accuracy is analyzed. The asymptotic estimation accuracy of the OTV .ethod is co.pared to the best achievable accuracy corresponding to the Cra.er-Rao bound . Keywords. system identification; para.eter esti.ation; Walker .ethods; asymptotic distribution

sinusoidal frequencies;

Yule-

Finally, we note that related studies of accuracy of the OTV esti.ates of ARMA system para.eters have been presented recently by the authors in Stoica et al (1985, 1987b) . As is well known, the sinusoidsin-noise process can be .odelled as a "degenerate" ARMA process whose poles (and zeros) are located on the unit circle . The theory derived in Stoica et al (1985, 1987b) does not appply directly to such degenerate ARMA processes, and a separate treat.ent is required .

INTRODUCTION The problem of esti.ating the frequencies of .ultiple sinusoids fro. noisy .easure.ents has received considerable attention recently, Kay and Marple (1981), Haykin and Cadzow (1982), Cadzow (1980, 1982), Tufts and Ku.aresan (1982), Chan and Langford (1982), Rife and Boorstyn (1976), Friedlander (1984), Friedlander and Porat (1984), Stoica et al (1986a, 1986b, 1986) SOderstro. and Stoica (1988). The overdeter.ined Yule-Walker (OTV) .ethod, Kay and Marple (1981), Cadzow (1982), Tufts and Ku.aresan (1982), Chan and Longford (1982), Friedlander (1984), SOderstro. and Stoica (1988), is often used for this purpose due to its co.putational si.plicity and satisfactory statistical accuracy.

THE OYW METHOD

2

Consider the sinusoids x(t)

The accuracy properties of the OTV .ethod have been studied so far only by si.ulation. A for.al accuracy theory for the OTV estiaates of sinusoid frequencies appears to be lacking . Our aim here is to fill this gap . We derive an explicit expression for the asymptotic covariance aatrix of the estiaation errors. The effect of increasing the nu.ber of TV equations on estiaation accuracy is analyzed in detail. This analysis provides valuable verification of the e.pirically established fact that estiaation accuracy i.proves when the nuaber of TV equations is increased, cadzow (1980, 1982), Tufts and Ku.aresan (1982), Chan and Longford (1982), Friedlander (1984), Stoica et al (1986a) SOderstro. and Stoica (1988). Coaparison of the asymptotic accuracy of the OTV estiaates with the best achievable accuracy corresponding to the Cra.er-Rao Lower Bound, Rife and Boorstyn (1976), Stoica et al (1986) is also included.

•[ a

following

sin(~kt+'

k=l k

signal

k

consisting

)

of

•

(2 . la)

where a , ,

€ R, ~ € (0,11), ~ . _~ . for i_j (2 . 1b) k k k 1) Let y(t) denote the noise-corrupted .easureaent of x(t)

y(t) = x(t)+e(t)

t=l, 2, . .. ,

(2 . 2)

where (e(t)) is assuaed to be a sequence of independent and identically distribut,d randoa variables with zero .ean and variance ~ . We also assu.e that x(t) and e(s) are uncorrelated for any t and s . As is well known, x(t) obeys the following autoregressive difference equation, Kay and Marple (1981), Tufts and Ku.aresan (1982), Stoica and SOderstro. (1983): \035

P. Stoica et al.

1036

A(q -1 )x(t)

=0

(2.3a)

-1

whefe q denotes [q- x(t)=x(t-1)] and degree n=2m given by

the_ unit A(q 1 ) is

-1" A(q) 1+a 1q -1 + . . . +a q -n n

m

rr (1-2coswkq

-1

+q

-2

a

delay operator polynomial of

~_

(2.3b)

)

k=l It follows from (2.2) and (2 . 3) that yet) obeys the following autoregressive moving-average (ARNA) difference equation

n complex conjugate zeros of A(z) are chosen as the estimates of the sinuseidal frequencies {wi} ' The estiaation errors of {wi} will evidently depend on the estiaation errors of {at') ' A foraula .relating t~e covariance matrix of {w · to that of (a.) , when (wi) is determined from (at, t as described a~ove, is presented in Stocia et al 1986b) . The availability of this formula makes it possible to restrict attention to the problem of analyzing the estiaation errors of (ail and this is what we will do in this paper . It is worth noting that since A(z) has conjugate unit-modulus roots we must have

The poles of (2.4), i . e . the f£gts of A(z), are located on the unit circle at e- k, k=l, .. . ,m . Next we briefly review the weighted OYW method for estimating {w · }. The essential step of this method is the estimaEion of the {a . } parameters by solving the following minimization problem, Kay and Marple (1981), Cadzow (1980, 1982), Friedlander (1984), Soderstrom and Stoica (1988): •• • 2

= an-l..

(2 . 7) i=O, . .. ,m (a =1) o If these constraints on the paraaeters {ai' are taken into account, then the OYW systea of equations (2.5a) has to be properly modified, Stoica et al (1986a, 1986b) . The constrained estimates are often, but not always, more accurate than the unconstrained estimates; see Stoica et al (1986a) for details. In this paper we will concentrate on the study of the accuracy properties of the unconstrained OYW estimate (2 . 5a) . The constrained estimate can be analyzed in a similar way, Stoica et al (1986a) . a.

1

(2.4)

(2 . 5a)

min{IRe+rN Q)

e

ASYMPTOTIC DISTRIBUTION

3

where

r

R

"=

,

M

>n

(2.5b)

First we establish the asymptotic (for N-) distribution of the OYW estimate (2 . 5a) in the case M<- . Then we treat the case where M-, which appears somewhat more complicated.

(2 . 5c)

Let e denote the true parameter.vector . Similarly, let Rand r denote the matrix.R and the vector r, with the sample covariance Irk) reRlaced by the theoretical covariances Irk}' rk=E{y(t)y(t+k»). Since (cf (2 . 2), (2 . 3a) )

r n+M-l

(3 . 1)

R9+r = 0 rk ~

N=X

N-k [y(t)y(t+k) t=l

(2.5d)

and where Q is a positive definite weighting mat!Ax T, N denotes the number of data points and IVIQ=V Qv . There are several numerically robust algorlthms which can be used to solve this problem, see for example Kayand Marple (1981), Cadzow (1982), Tufts and Kumaresan (1982), Chan and For the Langford (1982), Stoica et al (1986b) . purpose of this paper it suffices to note that the solution e of (2 . 5a) is given by the leastl/~quares solution of the weighted OYW equations (Q is a square root of Q) 1/2 "

Q

Re = -Q

coaplex-

1/2 '

r

or more explicitly 'T • -l'T e = -(R QR) R Qr

(2 . 6a)

and

the estimation errors {rk-r} are Bartlett (1946), Anderson and k Walker (1964), Stoica et al (1987), we can write for M
o

since

(11 fN),

T

-(R QR)

Once the {a . } parameters are estimated, one can in principle obtain estimates of the frequencies {w . } from (2.3b) . Indeed, (2.3b) induces a map from {all tQ {w . } which could be used in estimating {w.l f!qa (ai).l However,. since the polynomial lA(q ) corresponding to e is not guaranteed to have all its zeros on the unit circle, it is not possible to use that map exactly . A simple idea to circumvent this difficulty is to use the mfollowing approximation: the angular positions (w i )i=1 of the

•

It follows from (3 . 2) and standard results on convergence in distribution, Cramer ( 1951) , (see alsQ Theorem A. l of Stoica et al (1987)), that fN(e-e) has asymptotically a normal distribution, with zero mean and covariance matrix T

The inverse in (2 . 6) exists, at least for large N, since asymptotically the matrix R has full rank, Chan and Langford (1982), Stocia and Soderstrom (1983) .

•

R Q{fN[(r-r)+(R-R)e]}+O(l/fN) (3 . 2)

P = (R QR) (2 . 6b)

-1 T

-1 T

T

R QSQR(R QR)

-1

(3 . 3)

where. S is the asymptotic covariance matrix of {fN[ (r-r)+(R-R)e]) . An explicit expression for S can be obtained as follows . The (k,l)-th element of S is given by n S1.l = lim E{N[ [a . (rn+k_·-r . )] ~ N.i=O 1. l. n+k-l. n

.[[ a.(rn+l_j-r 1_ ' )]) j=O J n+ ) n

n

N [ .[ aiajcov(rn+k_i , rn+l_j) i=O )=0

1
1037

Statistical Anal YS is o f th e O verd ete rmin ed Yule-walke r Estimator

To proceed we need the following result, proved in Stoica et al (1987) : The nor.alized esti.ation errors IH(~n-r) , n=O,l, ... M, are as~ptotically jointly nor.al!y distributed with zero .ean and covariance matrix V.

~TR is clearly bounded and nonsingular for all M, it follows that the sa.e is true for ~TQR,

o

~.in(Q) T 1 T IQI. T < --M- R R < rf. QR < ~ R < • for all M

(3.9) Further.ore, for any vector v=O{£), where £ is so.e s.all nu.ber , we have IQvl

4

4

2~ (Qi-j+Qi+))+~ 61,)+(~-2~ )6 i ,06 j ,0

Vij

(3 . 10)

Thus, for large N and M we get from (3 . 2) [using (3 . 10))

The i,j element of V is given by 2

Ivl = 0(£) . < IQI..

rnm1(9-9)

O
=

1 T 1 -1 1 T •. IM -[HR QR+O(lN)] [JRR QIN{(r-r)+(R-R)91+0(~)]

( 1 i=j

61, · ).

~

,0

1 T

-[~ QR]

i_j

-1

T

•

•

[~Q/N{(r-r)+(R-R)9I]+0(~)

and (3.11 ) Using this result in (x(t)x(t+k)l, we get S

kl

(3.4) ,

and

denoting

A k

n 2~ 2 r

n r a' a . (Q . . +r , .) + i=O j=O 1) k+)-1-1 2n+k+1-1-)

n

+~4r

where the inverse exists (cf (3 . 9)) .

Q

n

raa .6 i=O j=O 1 n+k-1,n+l-j

It follows from (3 . 11) that to obtain a result similar to (3 . 5), M should increase more slowly than N. To proceed with the analysis and be .ore exact about the rate at which M can increase with H, we need to consider the covafiance .atrix of the first term in (3 . 11), say P. Paralleling the calculations in (3.2), (3 . 3) we get (3 . 12)

MP where we used the convention a . =O if in . In establishing the last equality in (3.5) we used the fact that (cf (2 . 3a)) n (3 . 6) for all k r a . Qk ' = i=O 1 -1

°

where P and S are defined by (3 . 3) and (3.7), respectively . The matrix P is bounded since ~

lIax

sup ~4IA(eiw) 12 <.

(S)

for all M

w

Hote that (3 . 5) can be rewritten as ~

2

E{A(q

-1

)e(t-k)A(q

-1

)e(t-l)1

(3.13)

(3.7)

Thus, S is the MIM covariance matrix associated with a mov4ng-ay~ra~e process with spectral density funchon ~ IA(e ) I . Hext we consider the case where M increases without bound as H tends to infinity . In other words, we let M depend on H, M=M(H), and allow M(H) to tend to infinity as H+M. One may expect that the results derived above will still hold, provided that the ratio M(H)/H goes to zero at an 'appropriate' rate as H.·. This 'appropriate rate' is, however, problem dependent and , therefore, difficult to determine . For Q=I we will ~?~~6that if M(H) does not increase faster than H ,6>0, then the results above remain valid. However, no atte.pt will be made to derive a tight upper bound on M (H) for a general Q. When M+M the matrix (RTQR) (in (3 . 3)) beco.es unbounded, unless it is properly nor.alized . The nor.alizing factor depends on Q. We assu.e that for all M (including M·.) IQI • ~ .ax r IQ 1) " I <• " 1

(3 . 8)

)

(The use of these mild conditions on Q will si.plify the analysis considerably . ) Then , since

Grenander and szego (1958), Moses et Stoica et al (1987a), and thus

al

(1987),

~ (S) T 2 ~ (S)IQI 2 ~ RTQSQR < ~ R Q R < .a~ • RTR < • (3 . 14) However, P may tend to zero as M+M . For instance this happens for the couonly used choice Q=1. Indeed, for Q=I the central matrix in (3 . 12) becomes

r~n . ~TSR

~n+M- l1

l '

I

JI!" •

r

'r

l 1

a

n

al 1 0

°a

M

j

a

a

n

l

J

° 1 I

0 a

H~n

a' I . l 1 r n+M- l

1

r

r

l

M

n

.a

1

1! j

1038

P. Stoica et al.

o

(3.15)

where • denotes entries which may be non-zero . To establish the second equality we used (3.6) and the fact that r =0 for n#O. Since P=O( l/M) (cf (3.15», the ¥ir~t term in (3 . 11) is O(l/lM) . Thus, to conclude in the case of Q=I that the second term in (3 2 11) is negligible we must require that [M(N)] IN~O as N~- . For most other choices of Q (satisfying (3 . 8» we_expect that the first term in (3.11) is 0(1) (i . e. P is 0(1» which means that in such cases it is sufficient to require that M(N)/N+O as N+-. As shown above the appropriate rate at which M (N) can increase when N increases is 'problem dependent' and cannot be readily derived for a general Q matrix . In any case, if M does not increase 'too fast' compared to N, the second term in (3.11) can be neglected and then we can conclude from (3.11), (3.12) that for large N ill M the covariance matrix of IN(9-e) is still given by P. The results on the asymptotic properties of the OYW estimate (2 . 6), which were derived above are su..arized in the following theorem . Theorem 3 1 Consider the sinusoids-in-noise process (2 . 1), (2.2) and its ARMA representation (2.4) . Let 9 be the weighted OYW estimate of the (ai} parameters (2 . 6) . Assume that the matrix Q satisfies the condition (3.8) (for M<- this condition on Q is trivially satisfied). Assume also that either M is finite or M increases without bound as N increases but does so no~ 'too fast' (for Q=I not 'too fast' means that M IN~ as N+-; for many other Q matrices, it is sufficient that M/N~O as N~-). Then for large N, the estimation errors (9-e) are (approximately) Gaussian distributed with zero mean and covariance matrix PIN, where P is defined by (3 . 3), (3 . 5) .

4

ACCURACY ASPECTS

Theorem 3 . 1 can be used to evaluate various choices of M and Q by comparing the accuracies of the resulting estimates . In particular it is of interest to clarify the role of increasing M, the number of the Yule-Walker equations . The weighting matrix Q also influences the accuracy and we will discuss the feasibility of an optimally weighted OYW estimator . It follows from (3.3) tbat for small M and large N the estimation errors (e-e) are of order l/fN. For large N ~ M (3 . 11), (3 . 12) imply that, under the cQnditions of Theorem i~A' the estimation errors (9-9) are of order l/N ,for some small positive 6. Therefore, the accuracy improves considerably in the limit as ~- . The estimation errors do not necessarily decrease monotonically with increasing M. See the numerical examples in Section 5. However, they generally tend to decrease when M increases and in the limit (as M+-) they become much smaller than the estimation errors corresponding to a finite (small) M. The discussion above provides a theoretical verification of the empirical observations that estimation accuracy improves considerably as the number of YW equations increases, Cadzow (1980, 1982), Tufts and Kumaresan (1982), Chan and Langford (1982), Friedlander (1984), SOderstrom and Stoica (1988).

Next we consider the problem of assessing the statistical efficiency of the OYW method . The Cramer-Rao lower bound (CR~¥~ on the estimation errors (~-9) is of order l/N ,Rife and Boorstyn (1976), Stoica et al (1986). Thus it appears that the OYW estimator is asymptotically very inefficient (in the statistical sense) as, 1/2+6

19-910yw/le-9ICRLB * N

»1

(for large N) (4 . 1)

However, the validity of (4.1) depends on the tightness of the upper bound imposed on M by Theorem 3.1. To be more exact, consider the case where Q=I . Then it follows from the calculations in the previous section that for large N and M we have 1

19-91 0yw = 0(M7N)

(4.2) N1/2 - 6 ,

provided M does not increase faster than for some small 6>0. However, this upper bound on the rate at which M can tend to infinity may be too conservative . In practice one often uses M=~N for some constant ~<1. 11 Theorem 3.1 would remain valid for this choice of M~ then H9-91 0yW WOUld. be of t~12same order of magnItude as Ne-9I CRLS ' I . e . O(l/N ) (cf (4 . 2», and the ratio (4 . 1) would remain bounded as N~- . In Section 5 we evaluate the asymptotic covariance matrix P for some sinusoids-in-noise processes, with Q=I and increasing M. For large values of M, P approaches PCRLB but never reaches it. Thus, the OYW method is asymptotically inefficient even if its statistical accuracy improves considerably when M increases. Let us remark that in Stoica et al (1987) we give some considerations pertaining to the choice of the matrix Q. 5

NUMERICAL EXAMPLES

In the previous section we have derived the asymptotic covariance matrix P of the parameter vector 9. From P we can evaluate the variances of the corresponding frequency estimates (w.} as described in Stoica et al (1986b). We have ~one so for the cases shown in Table 5. 1. In all cases the signal consists of two sine ~aves and white measurement noise with variance A =1 . The number of data, N, was assumed to be 2000 . Case

a1

a

2 3

1012 10/2 10/2

4 5

"'1

"'2

10/2

0 . 23.

10/2

0 . 53" 0 . 33.

10/2

0.23" 0 . 23"

12 12

12 12

0 . 23" 0.23"

0 . 53" 0.33"

6 7

12 120

12 12

0 . 23" 0 . 23"

0 . 26" 0 . 53"

8 9

120 120

12 12

0 . 23" 0 . 23"

0 . 33" 0.26"

I~t!l!: ~

2

0.26"

1. Amplitudes and frequencies for the cases

considered . In Figure 5.1 we illustrate the variances of wl and "'2 as functions of M. The Cramer-Rao lower bounds, RIfe and Boorstyn (1976), Stoica et al (1986), are also shown . Case 7 is illustrated in Figure 5 . 2, where we display the variances for a large number of M values . Let the variances given in Figure 5. 1

Statistical Analysis of the Qverdetermined Yule-walker Estimator "1

be denoted by P~i' (i=I,2), where "1=2000 denotes the correspondlng nuaber of data points. S~ilarlY, the Craaer-Rao bounds are denoted by p 1. The variances and the Craaer-Rao bounds for "'i any other value of ", say "2' can be siaply obtained froa those of Figure 5.1, 5.2 as

provided "1

and

The nuaerical results presented in the figures confira the following properties of the OYW estiaates, which were predicted by the theory developed. The as~ptotic variances of the OYW frequency estiaates decrease significantly with increasing M. However, they do not always decrease aonotonically. The variances are saall, indicating accurate frequency estiaates, when the signal-to-noise ratio is large or the frequencies are well separated.

(i=I,2)

(these equations hold 'sufficiently large').

are

The variances obtained with the OYW aethod aay be considerably larger than the Craaer-Rao lower bound, cf for exaaple !ases 6 and 9 in Figure 5.1. In fact POYW~I/SHR for large M and" while the Craaer-Rao lower bound is ~I/SHR for large

".

• 04

r:

• . e·· ....... .04

.

••-.. L==========·~··-~·~··-~-~·-=·--~··-=-~-=--""-·

---------

L-==-"'-·===-=--oc""'--====·=··"',,·=·_=-=-o:-=-=-c:c-"'·==o

.I .....'.

,.4 __";'::"~'~-~'::::.-=I~""'~'::':"::.'.::"::.":..:'::."'::'''='::''':.':'-:::::,;_-,

I...... . • 04

.

.

... -_.------ - - - ,..... L:::::.;=-=-=-·=·-=·-.::---=-~-==:;.:.::=::;:---:;---;;---;

.. ..'

'04~==--=-=-~--== '~ - "'= -~=='=-~'=--=-=~-==-=-~~=--~

~_~'::':o.:"::': '-:;:..::'-:::::.;':...'-="~"'::"':'['::":":'-'-",,:,'" ;.:.'''.:I_ _ _ _,

.04

I

I ........•04 ....

....

•04 ,04•

•• __...:..••.::'.:..:..... ::;::.:....:...-:...:.l.:"::"'=..:-'c.'-'-"_"~'_'_'''-,I,--.-___--.

.•

.•

....

.... Figure 5.1 . The variances of the frequency estiaates as functions of M.

CIIF2-N-

1039

P. Stoica et al.

1040

U Grenander and M Rosenblatt (1957) Statistical Analysis of Stationary Almquist and Wiksell, Stockholm.

Time

U Grenander and G Szego (1958) Toeplitz Forms and Their Applications. California Press, Berkeley, California .

Univeristy

Series .

of

E J Hannan (1970) Multiple Time Series . John Wiley & Sons, New York . S. Haykin and J . A. Cadzow, eds . (1982) Proc IEEE, Special Issue on Spectral Estimation, 70, no 9.

.....!• .....::.'7.--'----::----::---::---::----::----::-"'-~

Figure 5.2 . The variances as a functions of M, with higher resolution, case 7 .

6

CONCLUSIONS

We presented an accuracy analysis for the problem of estimating the frequencies of sinusoidal signals corrupted by measurement noise, using the overdetermined Yule-Walker method. An explicit formula for evaluating the asymptotic estimation error covariance matrix was presented. This formula provides a useful reference point for Monte- Carlo studies of the OYW method . It verifies the experimentally observed fact that estimation accuracy will generally improve as the number of Yule-Walker equations is increased . Finally, the accuracy of the OYW method was compared to the best achievable accuracy corresponding to the Cramer-Rao bound. It was shown that the OYW method is statistically inefficient, even for a large number of Yule-Walker equations . ACKNOWLEDGEMENTS The work by B. Friedlander was supported by the Office of Naval Research under contract NOOOI4-84C-0408 . REFERENCES T W Anderson and A M Walker (1964) On the Asymptotic Distribution of the Autocorrelations of a Sample from a Linear Stochastic Process. Annals of Math Statist, vol 35, no 3, pp 1296 - 1303. M S Bartlett (1946) On the Theoretical Specification and Sampling Properties of Auto- correlated Time-Series . J Roy Statist Soc, Ser B, vol 8, pp 27-41. J A Cadzow (1980) High performance Estimation - A New ARMA Method. IEEE Trans Acoustics, Speech and Signal Processing, vol ASSP-28, no 5, pp 524-529 . J A C4dzow (1982) Spectral Estimation : An overdetermined Rational Model Equation Approach. Proc IEEE, vol 70, no 9, pp 907-939 .

S Kay and S L Marple (1981) Spectrum Analysis - A Modern Perspective . vol 69, pp 1380-1419 . L Ljung (1977) Some Limit Theorems for Functionals Processes. Report LiTH - ISY-I-0167, Electrical Engineering, Linkoping Linkoping .

Proc

vol

IEEE,

of Stochastic Department of University,

R Moses, P Stoica, B Friedlander and T Soderstrom (1987) An Efficient Linear Method for ARMA Spectral Estimation. Proc IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Dallas. D C Rife and R R Boorstyn (1976) Multiple Tone Parameter Estimation from Discrete-Time Observations . Bell Sys Tech Journal, vol 55, pp 13891410 . T SOderstrom and P Stoica (1988) System Identification. Prentice-Hall, Hemel Hempstead. P Stoica, B Friedlander and T Soderstrom (1986a) On Instrumental Variable Estimation of Sinusoid Frequencies and the Parsimony Principle . IEEE Trans on Automatic Control, vol AC-31, pp 793 - 795. P Stoica, B Friedlander and T Soderstrom (1986b) Asymptotic Properties of High-Order Yule-Walker Estimates of Frequencies of Multiple Sinusoids . Proc IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), Tokyo. P Stoica, B Friedlander and T Soderstrom (1987a) Approximate Maximum Likelihood Approach of ARMA Spectral Estimation . Int J of Control, vol 45, pp 1281-1310. P Stoica, B Friedlander and T Soderstrom (1987b) Optimal Instrumental Variable Multistep Algorithms for the Estimation of AR Parameters of an ARMA Process . Int J of Control, vol 45, pp 2083-2107. P Stoica, R Moses , B Friedlander and T Soderstrom (1986) Maximum Likelihood Estimation of the Parameters of Multiple Sinusoids from Noisy Measurements . Third ASSP Workshop on Spectrum Estimation and Modelling, Boston. P Stoica and T Soderstrom (1983) Optimization with Respect to a Covariance Sequence : Algorithmic Aspects and Some Applications . Report UPTEC 83112R, Department of Technology , Uppsala University, Uppsala, Sweden. An abbreviated version has appeared in Automatica, vol 21, pp 671-675, 1985 .

Y T Chan and R P Langford (1982) Spectral Estimation via the High-Order Equations. IEEE Trans Acoustics, Speech Processing, vol ASSP-30, no 5, pp 689-698 .

Yule-Walker and Signal

p Stoica , T Soderstrom and B Friedlander (1985) Optimal Instrumental Variable Estimates of the AR Parameters of an ARMA Process . IEEE Trans Automatic Control , vol AC-30, pp 1066-1074 .

B Craaer (1951) Mathematical Methods of Statistics. University Press, Princeton, NJ .

Prince ton

P Stoica, T Soderstrom and F N Ti (1987) Overdetermined Yule-Walker estimation of frequencies of multiple sinusoids : some accuracy aspects . Report UPTEC 87107R , Department of Technology , Uppsala University, Uppsala , Sweden .

B Friedlander (1984) The Overdetermined Recursive Instrumental Variable Method. IEEE Trans Automatic Control, vol AC-29, no 4, pp 353-356. B Friedlander and B Porat (1984) The Modified Yule-Walker Method for ARMA Estimation . IEEE Trans Aerospace Electronics AES-20, no 2, pp 158-173 .

Spectral Systems,

D W Tufts and R Kumaresan (1982) Estimation of Frequencies of Multiple Sinusoids : Making Linear Prediction Perform Like Maximum Likelihood . Proc IEEE , vol 70, no 9, pp 975-989 .