Maximum Likelihood Parameter Estimation of the ... - CiteSeerX

3 downloads 4540 Views 573KB Size Report
all i 6= j, and the spectral density function of the process fs(1;0) i. (n)g is k(1;0) i. (!). More generally, a model for the evanescent eld which corresponds to the ...
Maximum Likelihood Parameter Estimation of the Harmonic, Evanescent and Purely Indeterministic Components of Discrete Homogeneous Random Fields  Joseph M. Francosy, Anand Narasimhanzand John W. Woodsx

Abstract This paper presents a maximum-likelihood solution to the general problem of tting a parametric model to observations from a single realization of a 2-D homogeneous random eld with mixed spectral distribution. On the basis of a 2-D Wold-like decomposition, the eld is represented as a sum of mutually orthogonal components of three types: purelyindeterministic, harmonic, and evanescent. The suggested algorithm involves a two-stage procedure. In the rst stage we obtain a suboptimal initial estimate for the parameters of the spectral support of the evanescent and harmonic components. In the second stage we re ne these initial estimates by iterative maximization of the conditional likelihood of the observed data, which is expressed as a function of only the parameters of the spectral supports of the evanescent and harmonic components. The solution for the unknown spectral supports of the harmonic and evanescent components reduces the problem of solving for the other unknown parameters of the eld to linear least squares. The CramerRao lower bound on the accuracy of jointly estimating the parameters of the di erent components is derived, and it is shown that the bounds on the purely-indeterministic and deterministic components are decoupled. Numerical evaluation of the bounds provides some insight into the e ects of various parameters on the achievable estimation accuracy. The performance of the Maximum-Likelihood algorithm is illustrated by Monte-Carlo simulations and is compared with the Cramer-Rao bound. Index terms: ML estimation of 2-D random elds, 2-D Wold decomposition, 2-D mixed spectral distributions, purely-indeterministic elds, harmonic elds, evanescent elds, Cramer-Rao bound. This work was partially supported by NSF grant MIP-9120377. ECE Department, Ben-Gurion University, Beer-Sheva 84105, Israel. z I.B.M Thomas J. Watson Research Center, NY 10598. x ECSE Department, Rensselaer Polytechnic Institute, Troy NY 12180

 y

1

1 Introduction In this paper, we consider the problem of tting a parametric model to observations from a single realization of a two-dimensional (2-D) complex valued discrete and homogeneous random eld with mixed spectral distribution. This fundamental problem is of great theoretical and practical importance. It arises in several areas of radar and sonar processing, and the special case of real valued 2-D random elds arises quite naturally in terms of the texture estimation of images [9]. The general problem of random elds' parameter estimation has received considerable attention. Most approaches reported to date fall into one of two categories. They either try to t noise driven linear models (2-D auto-regressive (AR), moving average (MA), or auto-regressive moving average (ARMA)) to the observed eld, or they treat the special case of estimation of the parameters of sinusoidal signals in white noise. Noise driven linear models have absolutely continuous spectral distribution functions, and hence, are inappropriate for the general problem considered here. Parameter estimation techniques of sinusoidal signals in additive white noise include the periodogram based approximation (applicable for widely spaced sinusoids) to the maximum likelihood (ML) solution [2], extensions to the Pisarenko harmonic decomposition [3], or the singular value decomposition [5]. These methods rely heavily on the white noise assumption, and are therefore not applicable here, since in our more general setting, the noise is colored, and a-priori unknown. Note that covariance based estimation procedures must assume knowledge of the true covariances. If these are unknown, substituting them with the sample covariances is incorrect, since it is well known [6] that even under the Gaussian assumption, the sample covariances are not consistent estimates of the covariance function if the spectral distribution function has discontinuities. The 2-D Wold-like decomposition [1] implies that any 2-D regular discrete and homogeneous random eld can be represented as a sum of two mutually orthogonal components: A purelyindeterministic eld and a deterministic one. The deterministic component is further orthogonally decomposed into a half-plane deterministic eld and a countable number of mutually orthogonal evanescent elds. This decomposition results in a corresponding decomposition of the spectral measure of the regular random eld into a countable sum of mutually singular spectral measures. The spectral distribution function of the purely-indeterministic component is absolutely continuous, while the spectral measure of the deterministic component is singular with respect to the Lebesgue measure, and therefore it is concentrated on a set of Lebesgue measure zero in the frequency plane. For practical applications, the \spectral density function" of the regular eld's deterministic component can be assumed to have the form of a countable sum of 1-D and 2-D delta functions. The 1-D delta functions are singular functions supported on curves in the 2-D spectral domain. The 2-D delta functions are singular functions supported 2

on discrete points in the spectral domain. In this paper, we consider the problem of estimating the parameters of the di erent components of the decomposition from a single realization of the eld. In general, an unbiased estimator of the eld parameters will require joint estimation of the parameters of the harmonic, evanescent and purely-indeterministic components. We present a conditional maximumlikelihood solution to this simultaneous parameter estimation problem, for the case in which the purely-indeterministic component is a complex valued Gaussian random eld. The algorithm is a two-stage procedure. In the rst stage we obtain suboptimal initial estimates for the parameters of the spectral support of the evanescent and harmonic components. The initial estimates are obtained by solving the set of overdetermined 2-D normal equations for the parameters of a high-order linear predictor of the observed data. In the second stage we re ne these initial estimates by iterative maximization of the conditional likelihood of the observed data. This maximization requires the solution of a highly nonlinear least-squares problem. By introducing appropriate parameter transformations the nonlinear least-squares problem is transformed into a separable least-squares problem [11], [12]. In this new problem, the solution for the unknown spectral supports of the harmonic and evanescent components reduces the problem of solving for the other unknown parameters of the eld to linear least squares. Hence, the solution of the original least-squares problem becomes much simpler. The proposed method is useful even when the separation between the spectral supports of any two deterministic components is less than 1=N in each dimension (for an N  N observed eld). We also present the Cramer-Rao lower bound (CRB) for this estimation problem. We show that the bounds on the deterministic and purely-indeterministic components are decoupled, and derive closed form CR bounds on the accuracy of estimating the parameters of the harmonic and evanescent components of the eld. An early discussion on the problem of analyzing 2-D homogeneous random elds with discontinuous spectral distribution functions can be found in [7]. There, harmonic analysis is employed to analyze the long-lag sample covariances, since for such lags the contribution of the purely-indeterministic component is assumed to be insigni cant. In this framework, the detection problem for a special case of evanescent elds is also discussed. The idea in [7] is to rst test for the existence of the deterministic components. If such components are detected, their parameters are estimated and their contribution to the sample covariances is removed. Next, the spectral density function of the purely-indeterministic component can be estimated from the \corrected" sample covariances. In [9], a similar periodogram based approach was used. The paper is organized as follows. In section 2 we brie y summarize the results of the 2-D Wold-like decomposition, which establish the theoretical basis for the suggested solution. Then, in section 3, assuming that the purely-indeterministic component is a complex valued Gaussian AR eld, we formalize the parameter estimation problem and derive the conditional maximum 3

likelihood estimator in the presence of a single evanescent eld. In section 4 we elaborate on the problem of estimating the parameters of the evanescent random eld. Section 5 describes an iterative solution for the parameters of the spectral support of the evanescent and harmonic components and its initialization algorithm. In section 6 the Cramer-Rao bound is derived, and in section 7 we present some numerical examples to illustrate the performance of the suggested algorithm, and the behavior of the derived bounds.

2 The Homogeneous Random Field Model The presented random eld model is derived based on the results of the Wold-type decomposition of 2-D regular and homogeneous random elds, [1]. In this section we brie y summarize the results of [1]. Let fy(n; m); (n; m) 2 Z 2g, be a complex valued homogeneous random eld. Let y^(n; m) be the projection of y(n; m) on the Hilbert space spanned by those samples of the eld that are in the \past" of the (n; m)-th sample, where the \past" is de ned with respect to the totally ordered, non-symmetrical-half-plane support, i.e.,   [  (i; j )  (s; t) i (i; j ) 2 (k; `)jk = s; ` < t (k; `)jk < s; ?1 < ` < 1 : (1) The innovation with respect to the de ned support and total order is given by u(n; m) = y(n; m) ? y^(n; m) and its variance is denoted by 2. If E jy(n; m) ? y^(n; m)j2 = 2 > 0; the eld fy(n; m)g is called regular. The eld is called deterministic if E jy(n; m) ? y^(n; m)j2 = 0: A regular eld fy(n; m)g is called purely-indeterministic if y(n; m) 2 Spfu(s; t)j(s; t)  (n; m)g, where Spfg denotes the closure of the span. These de nitions result in the following decomposition theorem: Theorem 1 [8]: Let fy(n; m); (n; m) 2 Z 2g be a 2-D regular and homogeneous random eld. Then there exists a deterministic random eld fv(n; m)g and an innovations eld fu(n; m)g such that fy(n; m)g can be uniquely represented by the orthogonal decomposition

y(n; m) = w(n; m) + v(n; m) where

w(n; m) =

X (0;0)(k;`)

a(k; `)u(n ? k; m ? `)

(2) (3)

and P(0;0)(k;`) ja(k; `)j2 < 1; a(0; 0) = 1. The eld fw(n; m)g is purely-indeterministic and regular. The elds fw(n; m)g and fv(s; t)g are mutually orthogonal for all (n; m) and (s; t). 4

De nition 1: A 2-D deterministic random eld feo(n; m)g is called evanescent w.r.t. the

NSHP total-order o if it spans a Hilbert space identical to the one spanned by its column-tocolumn innovations at each coordinate (n; m) (w.r.t. the total order o).

The concept of column-to-column innovations of deterministic elds is best illustrated using the following example. Let f(i)j ? 1 < i < 1g be an in nite two sided sequence of i.i.d. Gaussian random variables with zero mean and unit variance. Let also the 2-D random eld fy(k; l)g be de ned by y(k; l) = (k). It is clear that y^(k; l) = y(k; l ? 1) = (k) = y(k; l). Therefore the eld fy(k; l)g is deterministic. On the other hand it is obvious that y(k; l) is not a vector in the Hilbert space spanned by the eld samples y(s; t) for s < k since this Hilbert space is spanned by f(i)j?1 < i < k ? 1g and contains no information about (k). The innovation of y(k; l) with respect to the Hilbert space spanned by the eld samples y(s; t) for s < k is what we call the column-to-column innovation of the deterministic eld y at the coordinate (k; l). Hence, in this example the eld fy(k; l); (k; l) 2 Z 2g is an evanescent eld. It can be shown that it is possible to de ne a family of NSHP total-order de nitions such that the boundary line of the NSHP has rational slope. Let ; be two coprime integers, such that 6= 0. The angle  of the slope is given by tan  = = . (See, for example, Fig. 1.) Each of these supports is called rational non-symmetrical half-plane (RNSHP). We denote by O the set of all possible RNSHP de nitions on the 2-D lattice, (i.e. the set of all NSHP de nitions in which the boundary line of the NSHP has rational slope). For the special case in which  = =2 the transformation is obtained by interchanging the roles of columns and rows. The introduction of the family of RNSHP total-ordering de nitions results in a corresponding decomposition of the deterministic component of the random eld: Theorem 2 [1]: Let fv(n; m)g be the deterministic component of a 2-D regular and homogeneous random eld. Then fv(n; m)g can be uniquely represented by the following countably in nite orthogonal decomposition X v(n; m) = p(n; m) + e( ; )(n; m) : (4) ( ; )2O

The random eld fp(n; m)g is half-plane deterministic, i.e., it has no column-to-column innovations w.r.t. any RNSHP total-ordering de nition. The eld fe( ; )(n; m)g is the evanescent component which generates the column to column innovations of the deterministic eld w.r.t. the RNSHP total-ordering de nition ( ; ) 2 O. Hence, if fy(n; m)g is a 2-D regular and homogeneous random eld, then y(n; m) can be

5

uniquely represented by the orthogonal decomposition

y(n; m) = w(n; m) + p(n; m) +

X

( ; )2O

e( ; )(n; m) :

(5)

Place Figure 1 here In the following, all spectral measures are de ned on the square region K = [?1=2; 1=2]  R [?1=2; 1=2]. The spectral representation of y(n; m) is given by y(n; m) = exp[2j (n! + K m )]dZ (!;  ), where Z (!;  ) is a doubly orthogonal increments process, such that dFy (!;  ) = E [dZ (!;  )dZ  (!;  )]. Fy (!;  ) is the spectral distribution function of fy(n; m)g. Let f (!;  ) be the corresponding spectral density function, which is the Lebesgue 2-D derivative of Fy (!;  ). F s(!;  ) denotes the singular part in the Lebesgue decomposition of Fy (!;  ). Let L be a set of Lebesgue measure zero in K , such that the measure de ned by F s(!;  ) is concentrated on L. Theorem 3 [1]: The spectral measures of the decomposition components in (5) are mutually singular. The spectral distribution function of the purely-indeterministic component is absolutely continuous, while the spectral measure of the deterministic component is concentrated on the set L of Lebesgue measure zero in K . Moreover, since both the half-plane deterministic eld and all the evanescent elds in the decomposition (5) are components of the deterministic component of the regular eld, their spectral measures are concentrated on subsets of the set L. The de nition of the evanescent eld and Theorem 3 imply that the spectral measure of the evanescent component that generates the column-to-column innovations of the deterministic component for ( ; ) = (1; 0) is a linear combination of spectral measures of the form dFe ; (!;  ) = k(!)d!dF s ( ) ; where F s( ) is a one-dimensional singular spectral distribution function and k(!) is a 1-D spectral density function. In other words, the spectral distribution function of each evanescent component is separable: it is absolutely continuous in one dimension and singular in the orthogonal one, (or a linear combination of such separable distribution functions). From Theorem 3 we have that the spectral measure of each evanescent component of the regular eld is concentrated on a set of Lebesgue measure zero. For practical applications we can exclude singular-continuous spectral distribution functions from the framework of our treatment. Hence, the \spectral density function" of the evanescent eld e(1;0) has the countable sum form, fe ; (!;  ) = Pi ki(1;0)(!)( ? i(1;0)); where () is a Dirac delta function. A model for this (1 0)

(1 0)

6

evanescent eld is given by, [1],

e(1;0)(n; m) =

(1;0) IX

i=1

s(1i ;0)(n)ej2mi

;

(1 0)

(6)

where the 1-D purely-indeterministic processes fs(1i ;0)(n)g, fs(1j ;0)(n)g are mutually orthogonal for all i 6= j , and the spectral density function of the process fs(1i ;0)(n)g is ki(1;0)(!). More generally, a model for the evanescent eld which corresponds to the RNSHP de ned by ( ; ) 2 O is given by ; IX  ; j 2 i (n +m ) ( ; ) (7) si (n ? m )e e( ; )(n; m) = (

(

)

)

2+ 2

i=1

where the 1-D purely-indeterministic processes fsi( ; )(n ?m )g, fsj( ; )(n ?m )g are mutually orthogonal for all i 6= j , and I ( ; ) is in nite in general. Hence, the \spectral density function" of each evanescent eld has the form of a countable sum of 1-D delta functions which are supported on lines of rational slope in the 2-D spectral domain. In the following we assume that each of the 1-D purely-indeterministic processes si( ; ) obeys a nite order autoregressive (AR) model. Thus, for example, the purely-indeterministic modulating process of e(1;0) is given by

s(1;0)(n) = ? i

(1;0) VX i

t=1

a(1i ;0)(t)s(1i ;0)(n ? t) + ?(1i ;0)(n) ;

(8)

where ?(1i ;0)(n) is a 1-D white innovations process. One of the half-plane-deterministic eld components, which is often found in physical problems is the harmonic random eld P X (9) h(n; m) = Cpej2(n!p+mp ) p=1

where the Cp 's are mutually orthogonal random variables, E jCpj2 = p2, and (!p; p) are the spatial frequencies of the pth harmonic. In general P is in nite. The parametric modeling of deterministic random elds whose spectral measures are concentrated on curves other than lines of rational slope, or discrete points in the frequency plane, is still an open question to the best of our knowledge. Theorem 1 implies that the most general model for the purely-indeterministic component of a regular homogeneous random eld is the MA model (3). However, if its spectral density function is strictly positive on the unit bicircle and analytic in some neighborhood of it, a 2-D AR representation for the purely-indeterministic eld exists as well [10]. In the following, we 7

assume that the above requirements are satis ed. Hence the purely-indeterministic component autoregressive model is given by X w(n; m) = ? b(k; `)w(n ? k; m ? `) + u(n; m) (10) (0;0)(k;`)

where fu(n; m)g is the 2-D white innovations eld whose variance is 2.

3 The Conditional Maximum Likelihood Estimator 3.1 Problem De nition and Assumptions The orthogonal decompositions of the previous section imply that if we exclude from the framework of our model those 2-D random elds whose spectral measures are singular continuous, or are concentrated on curves other than lines of rational slope, y(n; m) is uniquely represented by y(n; m) = w(n; m)+ h(n; m)+ P( ; )2O e( ; )(n; m). The problem of estimating the ( ; ) pairs of the di erent evanescent components is beyond the scope of the present paper. In order to keep the notations as simple as possible, we restrict our attention to the case in which it is a-priori known that ( ; ) = (1; 0) for all the evanescent components. The more general problem of estimating the eld parameters in the presence of evanescent elds which are characterized by unknown ( ; ) parameters, will be discussed in a forthcoming paper. Hence, the problem faced here is the parameter estimation of the harmonic and evanescent components (those of ( ; ) = (1; 0)) of the eld in the presence of an unknown colored noise generated by the purely-indeterministic component, jointly with estimating the purely-indeterministic component parameters. When expressed in the general form (9), the coecients fCpg of the harmonic component are complex valued, mutually orthogonal random variables. However, since in general, only a single realization of the random eld is observed we cannot infer anything about the variation of these coecients over di erent realizations. The best we can do is to estimate the particular values which the Cp 's take for the given realization; in other words we might just as well treat the Cp's as unknown constants. Finally, we note that a maximum likelihood solution to our parameter estimation problem involves maximization of the exact likelihood function. However, this is a formidable task due to the complexity in representing the covariance matrices in terms of the model parameters of the di erent components. For large enough data records the exact likelihood function can be approximated by the conditional likelihood function. Since this approach results in a more tractable solution, we have chosen it for the above parameter estimation problem. 8

We next state our assumptions and introduce some necessary notations. Let fy(n; m)g, (n; m) 2 D where D = f(i; j )j0  i  S ? 1; 0  j  T ? 1g be the observed random eld. Note however that the observed eld just as well could have any arbitrary shape. Assumption 1: The purely-indeterministic component is a complex valued Gaussian AR eld, whose model is given by (10) with (k; `) 2 SN;M nf(0; 0)g, where SN;M = f(i; j )ji = 0; 0  j  M g Sf(i; j )j1  i  N; ?M  j  M g, and N; M are a-priori known. The driving noise of the AR model is a complex valued Gaussian eld such that its real and imaginary components are independent real Gaussian white noise elds each with zero mean and variance 2 . Assumption 2: The number P of harmonic components in (9), as well as the number I (1;0) of evanescent components in (6), are a-priori known. Assumption 3: The 1-D purely-indeterministic processes fs(1i ;0)g, of the type (8), are all assumed to be complex valued Gaussian AR processes of known orders, Vi(1;0). The driving noise of each of the AR models is an independent, zero mean, complex valued Gaussian process, such that its real and imaginary components are independent real Gaussian white noise processes each ( i ; ) with zero mean and variance 2 . In the proposed algorithm we take the approach of rst estimating a non-parametric representation of the 1-D purely-indeterministic processes fs(1i ;0)g, and only in a second stage the AR models of these processes are estimated. Hence, in the rst stage we estimate the particular ?1 , i = 1 : : : I (1;0) take for the given realization, i.e., we values which the processes fs(1i ;0)(n)gSn=0 treat these as unknown constants. To simplify the presentation of this section, we shall describe the solution for I (1;0) = 1, i.e., in (6), i  1 . Hence in the following we omit all the subindices i. Further, since we shall only deal with the case where ( ; ) = (1; 0), we shall omit the notation (1; 0) as superscripts and subscripts from our derivation, up to section 7. Thus, the parameters to be estimated are fCp; !p; pgPp=1, ?1 , fb(k; `)g  , fs(n)gSn=0 (k;`)2SN;M ,  2. We denote this vector of unknown parameters by  . De ne 2

(1 0) 2

u = [u(N; M ); : : :; u(N; T ?1?M ); u(N +1; M ); : : : : : :; u(N +1; T ?1?M ); : : :; u(S ?1; T ?1?M )]T (11)

9

The vector y is similarly de ned. 2 ::: y(N; 0) y(N ? 1; 2M ) 66 y(N; M ? 1) y(N; M ) ::: y(N; 1) y(N ? 1; 2M + 1) 66 ... 66 6 Y = 666 y(N; T ? M ? 2) : : : y(N; T ? 1 ? 2M ) y(N ? 1; T ? 1) 66 y(N + 1; M ? 1) : : : y(N + 1; 0) y(N; 2M ) 66 . .. 64 y(S ? 1; T ? M ? 2) : : : y(S ? 1; T ? 1 ? 2M )

::: :::

y(0; 2M )

::: y(0; 0) ::: y(0; 1) ... ::: y(0; T ? 1) ::: y(0; T ? 1 ? 2M ) ::: y(1; 2M ) ::: y(1; 0) ... : : : y(S ? 1 ? N; T ? 1) : : : y(S ? 1 ? N; T ? 1 ? 2M )

::: y(N ? 1; 0) ::: ... : : : y(N ? 1; T ? 1 ? 2M ) ::: y(N; 0) ...

3 77 77 77 77 77 : 77 77 75

Also, we set 2 ej2[N! +M ] ej2[N! +M ] : : : ej2[N!P +MP ] 66 ... ... ... ... 66 66 ej2[N! +(T ?1?M ) ] ej2[N! +(T ?1?M ) ] : : : 66  ej2[(N +1)! +M ] : : : Eh = 666 ej2[(N +1)! +M ] .. .. .. .. 66 66 .. .. .. .. 64 ej2[(S?1)! +(T ?1?M ) ] ::: : : : ej2[(S?1)!P +(T ?1?M )P ] 2 j2M 3 66 je2(M +1) 77 e 77 W = 6666 77 ; ... 4 5 ej2(T ?1?M ) 1

1

2

1

1

1

1

2

2

1

2

2

1

10

2

(12)

3 77 77 77 77 77 ; (13) 77 77 77 5

(14)

2 66 W 66 W 0  6 Ee = 66 W 66 0 ... 4 and

W

3 77 77 77 ; 77 75

(15)

b = ?[b(0; 1); : : :; b(0; M ); b(1; ?M ); : : :; b(1; M ); : : : : : : ; b(N; ?M ); : : : ; b(N; M )]T : (16)

3.2 Conditional ML Estimation in the Presence of a Single Evanescent Component Since u(n; m) is assumed to be Gaussian,

 1 SX  ?1 T ?X 1?M 1 2 p(Y ; ; D n D1) = (2)jD j exp ? 2 ju(n; m)j : n=N m=M 1

(17)

The conditional MLE of  is found by maximizing (17), or equivalently by minimizing J () = P (n;m)D ju(n; m)j2, where D1 = f(i; j )jN  i  S ? 1; M  j  T ? 1 ? M g, and D n D1 is the set of required initial conditions. Thus, only actually occurring values of the observed eld are used in the estimation procedure. Using this method we sum the squares of only jD1j values of u(n; m), but this slight lost of information will be unimportant if the size of the observed eld, jDj, is large enough. Using (10), u(n; m) is given by u(n; m) = P(k;`)2SN;M b(k; `)w(n ? k; m ? `) with b(0; 0) = 1. Since w = y ? h ? e, we have   X u(n; m) = b(k; `) y(n ? k; m ? `) ? h(n ? k; m ? `) ? e(n ? k; m ? `) (k;`)2SN;M  P X X = b(k; `) y(n ? k; m ? `) ? Cpej2[(n?k)wp +(m?`)p ] p=1 (k;`)2SN;M  ?s(n ? k)ej2(m?`) X = b(k; `)y(n ? k; m ? `) 1

(k;`)2SN;M

?

P X

p=1

Cp

X (k;`)2SN;M

b(k; `)ej2[(n?k)wp+(m?`)p ]

11

X

?

b(k; `)s(n ? k)ej2(m?`)

X (k;`)2SN;M = b(k; `)y(n ? k; m ? `) (k;`)2SN;M 0 1 P X X b(k; `)e?j2(k!p+`p )A ej2(n!p+mp ) ? Cp @ 0p=1 (k;`)2SN;M 1 X ? j 2 ` A ej2m ?@ b(k; `)s(n ? k)e (k;`)2SN;M

De ne now the following transformations: X p = Cp b(k; `)e?j2(k!p +`p ) (k;`)2SN;M

(n) =

X (k;`)2SN;M

b(k; `)s(n ? k)e?j2` :

(18)

(19) (20)

Let B (ej2! ; ej2 ) = P(k;`)2SN;M b(k; `)e?j2(!k+`) . The assumptions made in section 2 as to the properties of the spectral density function of the purely-indeterministic eld imply that the eld AR model is such that B (z1; z2) is minimum phase. (We implicitly assume here that the nite support B (z1; z2) de ned above retains this property of the in nite support lter which corresponds to the AR model (10)). Since B (ej2!; ej2 ) is non-zero on the unit bicircle, and in particular at the frequencies of the harmonic components, the transformation (19), of the Cp's to p's is one-to-one. The transformation (20) is also one-to-one since given N initial values of the process fs(n)g, each newly introduced s(n) results in a unique (n). The idea of using a transformation of the type (19) was developed in one-dimension for estimating the parameters of harmonic signals in colored noise by Chatterjee et al. [14], as well as by Kay and Nagesha [15]. Let S 0N;M = SN;M nf(0; 0)g. We can therefore rewrite (18) in the following form

u(n; m) = y(n; m)+

X (k;`)2S N;M 0

Let and

P X

b(k; `)y(n?k; m?`)? pej2(n!p+mp )?(n)ej2m (n; m) 2 D1 : p=1

(21)

iT h ;  = 1; 2 : : : ; p

(22)

h iT :  = (N ); (N + 1) : : : (S ? 1)

(23)

12

Since J () = uHu, we obtain by writing (21) for all (n; m) 2 D1, the following matrix representation for J (): J () = jjy ? Yb ? Eh ? Ee jj2 : (24) Thus the transformations (19) and (20) allow us to minimize the objective function J () with respect to b, ,  , and the deterministic component spectral support parameters, f!p; pgPp=1;  , instead of minimizing it with respect to the original problem parameters. The properties of the above transformations guarantee that both minimizations will result in the same minima for J . De ne D  = [YEhEe] and 1 = [bTT T]T. Then we can rewrite (24) as

J () = jjy ? D1jj2:

(25)

Because of the fact that the objective function is a quadratic function of 1, the minimization over 1 can be carried out analytically for any given value of D. Using the well known solution to the least-squares problem we have that  ?1 (26) ^1 = DHD DHy will minimize J () over 1 . By inserting (26) into (25) we nd that the minimum value of J () is given by (27) Jmin (f!p; pgPp=1;  ) = yH(I ? D(DHD)?1 DH)y: Here D is assumed to be full rank so that (DHD)?1 exists. Thus, maximization of the likelihood function is achieved by minimizing the new objective function Jmin (f!p; pgPp=1;  ), which is a function only of the deterministic component spectral support parameters. We have thus shown that the minimization problem (24) which is obtained after taking the transformations (19), (20) is separable since its solution can be reduced to a minimization problem in the nonlinear deterministic component's spectral support parameters, f!p; pgPp=1;  , only, while b;; can then be determined by solving a linear least-squares problem. This new minimization problem is of a considerably lower complexity. A broad discussion on the subject of separable, can be found in [11]  nonlinear,  least-squares minimization problems P P and [12]. Since Jmin f!p; pgp=1;  is a nonlinear function of f!p; pgp=1;  this optimization problem cannot be solved analytically and we must resort to numerical methods. In order to avoid the enormous computational burden of an exhaustive search, we use the two-step procedure which is described in section 5. In the discussion above we assumed that the noise variance 2 is known. If it is not known, it can be estimated. The maximum likelihood estimate of 2 is derived by maximizing (17) with 13

respect to 2. Using the estimated frequencies, and (27) we have that ^ ^2 = (S ? NJ)(min T ? 2M ) :

(28)

Thus (26) and (28) establish the estimate for the autoregressive model of the purely-indeterministic component of the eld. Using the estimated frequencies of the harmonic component and the transformation (19), a complete estimate for the parameters of the harmonic component is obtained. The solution for the parameters of the evanescent eld is more involved and it is given in the next section.

4 Estimating the Parameters of the Evanescent Component De ne

8 PM ?j 2` k=0 `=0 b(0; `) e P M ? j 2 ` : `=?M b(k; `) e 1kN :

< B (k) =

Let us rewrite (20)

(n) =

X M `=0



b(0; `)s(n)e?j2`

= s(n)B (0) + =

N X k=0

N X k=1

+

N X M X k=1 `=?M

B (k)s(n ? k)

B (k)s(n ? k)



b(k; `)s(n ? k)e?j2`

n = N; : : : ; S ? 1 :

(29)



(30)

Note, that as a result of the preceding stages of the algorithm, the estimated values of the (n)'s and B (k)'s are available, rather than the true values. We next present two di erent approaches for estimating the parameters of the 1-D purely-indeterministic process fs(n)g. Following (8) and Assumption 3, we have that

s(n) = ?

V X t=1

a(t)s(n ? t) + ?(n) :

14

(31)

Substituting (31) into (30), we obtain # " X N V X (n) = B (k) ? a(t)s(n ? k ? t) + ?(n ? k) k=0 V X

= ?

= ?

t=1 V X t=1

t=1

a(t)

N X

k=0

B (k)s(n ? k ? t) +

a(t)(n ? t) +

N X k=0

N X k=0

B (k)?(n ? k)

B (k)?(n ? k)

n = N; : : : ; S ? 1 :

(32)

Hence, (32) implies that solving the problem of estimating the unknown parameters of the 1-D purely-indeterministic process associated with the evanescent component, is equivalent to solving the above one-dimensional ARMA equation. In the equivalent problem, the MA parameters fB (k)gNk=0 have previously been estimated, the driving Gaussian noise source is of an unknown variance ((1;0))2, and the \observations" are the f(n)g, for n = N; : : : ; S ? 1. In [17](pp. 205 { 208) it is shown how the exact likelihood function of an ARMA process can be computed, from the ARMA model parameters, by using the Kalman lter. Hence, maximization of the likelihood function with respect to the unknown parameters, fa(t)gVt=1; (1;0), will result in a ML estimate of the ARMA parameters. The alternative approach is simple to implement, but suboptimal. Moreover, this solution can be used to initialize the above ML algorithm. De ne the 1-D function B(z1) = PNk=0 B (k)z1?k . Note that the parameters of this function were previously estimated. Hence, if this estimated function is minimum phase then instead of solving (30), we can solve the equivalent problem obtained by deconvolving the sequence fB (k)gNk=0 out of the sequence f(n)gSn=?N1 . The result of ?1 . However, the system (30) is clearly an underdeterthis deconvolution is the sequence fs(n)gSn=0 mined system having S unknowns with only S ? N equations. Hence, we obtain a solution for ?1 by ltering the sequence f (n)gS ?1 through the IIR the required unknown sequence fs(n)gSn=0 n=N ? 1 lter B (z1) with zero initial conditions. Since fs(n)g is an AR process, a standard AR model parameter estimation algorithm can now be applied to estimate its parameters. (In the present paper we use the conditional maximum likelihood estimator [18]).

5 The Solution for the Spectral Support Parameters of the Deterministic Component In section 3 we concluded that the minimization problem (24) which is obtained after taking the transformations (19), (20) is separable since its solution can be reduced to a minimization 15

problem in the nonlinear deterministic component's spectral support parameters, f!p; pgPp=1;  , only, while b;; can then be determined by solving a linear least squares problem. Hence, the rst step in solving the presented estimation problem is the minimization of Jmin with respect to the unknown spectral support parameters of the deterministic component. Since Jmin is a nonlinear function of the deterministic component's spectral support parameters this optimization problem cannot be solved analytically and we must resort to numerical methods. In general, Jmin has a complicated multi-modal shape. Hence, in order to avoid the enormous computational burden of an exhaustive search, we used the following two-step procedure. In the rst stage we obtain a suboptimal initial estimate for the parameters of the spectral support of the deterministic component. This stage is implemented by solving the system of overdetermined 2-D normal equations for the parameters of a high-order linear predictor of the observed data. In the second stage we re ne these initial estimates by an iterative numerical minimization of the objective function Jmin . In our experiments we used the conjugate gradient method of Fletcher and Reeves [19] (p. 253). Note that only for the case of a quadratic objective function, the conjugate gradient procedure is guaranteed to converge in at most N steps. For our problem, we simply restart the algorithm using new gradients, until the objective function becomes appreciably small. As is well known, this type of iterative optimization procedure converges to a local minimum, and does not guarantee global optimality, unless the initial estimate is suciently close to the global optimum. As we show in section 7, the initial estimates provided by the solution of the overdetermined high-order normal equations appear to provide a good initial starting point (i.e., one which leads to convergence to the global minimum) as long as the local signal to noise ratio is suciently high, and the frequencies of the di erent deterministic components are not too close. We next describe the initialization algorithm. Formalizing the 2-D linear prediction problem for some NSHP predictor with support S = f(k; `)jk = 0; 0 < `  Lg Sf(k; `)j1  k  K; ?L  `  Lg results in the following linear system of equations: X a(k; `)r(i ? k; j ? `) = ?r(i; j ) (i; j )  (0; 0) (33) (k;`)2S

where fa(k; `)g(k;`)2S are the linear predictor coecients. Rewriting the system (33) in a matrix form for all (i; j ) 2 S results in the well known 2-D Yule-Walker equations. Including in this system additional equations for (i; j )  (0; 0) such that (i; j ) 62 S , results in an overdetermined system. The overdetermined Yule-Walker method is a modi cation of the basic Yule-Walker method, which was reported [4] to lead to a considerable increase in the estimation accuracy of the frequencies of harmonic signals in white noise for 1-D signals. It is further concluded in [4], that the asymptotic accuracy of the estimates will increase with the number of Yule-Walker equations 16

used and with the model support. Intuitively, it can be expected that increasing the predictor support (i.e., increasing K and L), will improve the accuracy of the estimates of the deterministic component's spectral support, since the covariances for large lags contain \useful information" about the deterministic component. In order to solve (33) for the linear predictor parameters, the covariances of the observed eld must be available. Since the covariance functions of the observed eld are unknown they must be estimated from the data itself. Hence, in (33) we replace the true covariances by their estimates. However, as noted earlier, the ergodic property of the sample covariances does not generally hold in the present case. This is due to the result [6], that the sample covariances of a Gaussian process are consistent estimates of the covariance function if and only if the spectral distribution function of the process has no discontinuities. Clearly, this requirement does not hold in the present case. Nevertheless, since in practice the above method produces accurate estimates of the spectral support of the deterministic component, and since these estimates are used only to initialize the iterative minimization of the cost function Jmin, the above theoretical problem is avoided. Since (33) is an overdetermined system, it is solved in the least-squares sense for the linear predictor coecients. We then look for the peaks of 1=jA^(ej2! ; ej2 )j2, to obtain the required initial estimates of the deterministic component spectral support parameters.

6 The Cramer-Rao Lower Bound A conditional ML algorithm for estimating the parameters of the harmonic and evanescent components in the presence of a 2-D circular Gaussian AR noise, jointly with estimating the AR model parameters, was suggested in the previous sections. In this section we derive the performance bound for this algorithm. As is well known the Cramer-Rao Bound (CRB), provides a lower bound for the covariance matrix of any unbiased estimator of the model parameters. The bound is given by the inverse of the so-called Fisher Information Matrix (FIM) [16]. In the following derivation, the number of components of the evanescent eld is I (1;0), the number of harmonic components is P , and the 2-D AR model support is SN;M . The conditional CRB is derived from the conditional probability density function of the observed process. Rewriting (17) we obtain 2  X X ; b ( k; ` )[ y ( n ? k; m ? ` ) ? d ( n ? k; m ? ` )] p(Y ; ; DnD1) = (21)jD j exp ? 12 (n;m)2D (k;`)2SN;M (34) 1

1

17

were we de ne

d(n; m) = h(n; m) + e(n; m) : (35) Note that in the present framework fd(n; m)g is the mean component of fy(n; m)g. Collecting the parameters of the harmonic component into vector representations, we have

c = [C1; : : :; CP ]T ;

(36)

fi = [!i; i]T ; ! h = [f1T ; : : :; fPT ]T :

(37) (38)

The parameter vectors of the evanescent components are de ned by

 e = [1; 2; : : :; I ; ]T ;

(39)

s0i = [si(0); si(1); : : : ; si(N ? 1)]T ; si = [si(N ); si(N + 1); : : : ; si(S ? 1)]T ; i h xi = s0i T ; siT T ;

(40) (41)

(1 0)

h

s = s1T ; s2T ; : : :; sI

;

(1 0)

i T T

:

(42) (43)

Also let

d = [d(0; 0); : : : ; d(0; T ? 1); d(1; 0); : : : ; d(1; T ? 1); : : :; : : : ; d(S ? 1; 0); : : : ; d(S ? 1; T ? 1)]T ;

(44) h = [h(0; 0); : : : ; h(0; T ? 1); h(1; 0); : : : ; h(1; T ? 1); : : : ; : : :; h(S ? 1; 0); : : : ; h(S ? 1; T ? 1)]T ; (45) ei = [ei(0; 0); : : :; ei(0; T ? 1); ei(1; 0); : : : ; ei(1; T ? 1); : : : ; : : :; ei(S ? 1; 0); : : : ; ei(S ? 1; T ? 1)]T : (46) Rewriting (35) in vector form we have,

d=h+

(1;0) IX

i=1

ei :

(47)

Let dR = Refdg, dI = Imfdg, d~ = [dRT dI T ]T . In a similar way we de ne the vectors hR, hI , h~ , eiR, eiI , e~i, xiR, xiI , x~ i, cR, cI , ~c, sR, sI , ~s. Also let bR(k; `) = Refb(k; `)g, bI (k; `) = Imfb(k; `)g. Let us denote by  the mean component parameter vector, i.e.,  = [~cT ! Th ~sT  Te ]T : 18

Note that the initial conditions vectors si0; i = 1; : : : ; I (1;0), are not parameters to be estimated, as they are assumed known. Taking now the partial derivatives of the conditional log likelihood function w.r.t. the elements of , we have 8 0 1 < X X @d(n ? k; m ? `) A @ ln P = 1 @ b ( k; ` ) 2  :  (n;m)2D (k;`)2SN;M @(i) @(i) 0 1 X @ b(k; `)[y(n ? k; m ? `) ? d(n ? k; m ? `)]A 0(k;`)2SN;M 1 X @ X b(k; `)[y(n ? k; m ? `) ? d(n ? k; m ? `)]A + (n;m)2D (k;`)2SN;M 19 0  (n ? k; m ? `) = X @d A; : @ (48) b(k; `)  @  ( i ) (k;`)2SN;M 1

1

Taking the partial derivative w.r.t. the AR process driving noise variance parameter yields 8 1  @ 2 ln P  < X 0 X @d ( n ? k; m ? ` ) 1 @ A b(k; `) ?E @2@(i) = 4 E :  @  ( i ) (n;m)2D (k;`)2SN;M 0 1 X @ b(k; `)[y(n ? k; m ? `) ? d(n ? k; m ? `)]A 0(k;`)2SN;M 1 X @ X + b(k; `)[y(n ? k; m ? `) ? d(n ? k; m ? `)]A (n;m)2D (k;`)2SN;M 19 0 (n ? k; m ? `) = X @d A; @ b(k; `)  @  ( i ) (k;`)2SN;M = 0: (49) 1

1

Similarly, taking the partial derivatives w.r.t. the real and imaginary parts of the AR process parameters we nd that



2

?E @ (i@)@blnRP(k; `)



 X  @d(n ? k; m ? `)  = ? 12 E @ (i) (n;m)2D  X  b(k; `)[y(n ? k; m ? `) ? d(n ? k; m ? `)] (k;`)2SN;M   X  X ? 12 E b(k; `) @d(n ?k; m ? `) @ (i) (n;m)2D (k;`)2SN;M 1

1

19

[y(n ? k; m ? `) ? d(n ? k; m ? `)]



 X 1 ? 2 E [y(n ? k; m ? `) ? d(n ? k; m ? `)] (n;m)2D  X  (n ? k; m ? `)  @d  b (k; `) @(i) (k;`)2SN;M  X  X  1 ? 2 E b(k; `)[y(n ? k; m ? `) ? d(n ? k; m ? `)] (n;m)2D (k;`)2SN;M  @d(n ? k; m ? `)  @(i) = 0: (50) 1

1

Similar derivation w.r.t. bI (k; `), yields





2

? E @ (@i)@blnI P(k; `) = 0 :

(51)

Thus the conditional FIM is block diagonal. Hence, the conditional CRB on the harmonic and evanescent components parameters is decoupled from the bound on AR process parameters. Therefore, the conditional CRB's on the deterministic component, and on the AR component parameters are obtained by inverting the FIM blocks which correspond to the deterministic and the AR parameters, respectively. Asymptotic CRB on the parameters of an AR eld with an NSHP support is given in [13]. In the following we concentrate on deriving the conditional CRB on the parameters of the harmonic and evanescent components. Taking the partial derivatives w.r.t. the mean component parameters we nd that   2 ? E @ @(i)ln@ P(j ) 10 1 0  X X X @ = 12 b(k; `) @d (n ?k; m ? `) A b(k; `) @d(n ?k; m ? `) A @ @(i) @ (j ) (n;m)2D (k;`)2SN;M (k;`)2SN;M 10 1 0  (n ? k; m ? `) X X @ X @d 1 @d ( n ? k; m ? ` ) A@ A b(k; `) + 2 b(k; `)   @  ( j ) @  ( i ) (k;`)2SN;M (n;m)2D (k;`)2SN;M 0 1 0 1   (n ? k; m ? `) X @ X  X 2 @d @d ( n ? k; m ? ` ) A@ A : = 2 Re b (k; `) b(k; `)   @  ( i ) @  ( j ) (n;m)2D (k;`)2SN;M (k;`)2SN;M (52) 1

1

1

Thus, the above derivation of the conditional FIM reveals that the bounds on both the amplitude and the frequency parameters of the harmonic components, as well as the bounds on the 20

parameters of the evanescent components, are functions of the frequency response of the colored noise model at the frequencies of the spectral support of the deterministic components, and of the derivative of the frequency response at these frequencies. Let d be the \ ipped around" version of d, i.e.,

d = [d(S ? 1; T ? 1); : : : ; d(S ? 1; 0); : : : ; : : :; d(1; T ? 1); : : : ; d(1; 0); d(0; T ? 1); : : : ; d(0; 0)]T : Note that d = Kd, where

3 2 0 : : : 0 1 7 66 66 0 : : : 1 0 777 K = 66 .. ... 77 ; . 5 4 1 ::: 0 0 is the exchange matrix. Let 0k denote a k-dimensional column vector of zeros. Let also, h i b0 = 0TM ; 1; b(0; 1); : : : ; b(0; M ); 0TT ?(2M +1) T ; h i b1 = b(1; ?M ); : : :; b(1; 0); : : : ; b(1; M ); 0TT ?(2M +1) T ; ... ... h i bN ?1 = b(N ? 1; ?M ); : : : ; b(N ? 1; 0); : : : ; b(N ? 1; M ); 0TT ?(2M +1) T ; bN = [b(N; ?M ); : : :; b(N; 0); : : : ; b(N; M )]T ; and

h

(53)

(54)

(55)

i

b = bT0 ; bT1 ; : : : ; bTN T : (56) Let b  denote the conjugate of b . De ne the following S  T  (S ? N )  (T ? 2M ) matrix 2  66 b 66 0 66 66 0 6 B = 666 0. 66 .. 66 66 66 4 0

0 ::: b  0 0 ... 0

0 ... 0 b  0 ... ... 0

j j j j j j j j

0 T 0T : : : b  0 0 b  ...

0

0

... 0

0T j 0 j ... j 0 j b  j 0 j ... j 0 j

3

:::

j : : : 0T(S?N?1) 7 77 j ::: 0 77 77 j 77 77 j 77 ; ... j 77 77 j 77 75 j 0 j b 

where each of the S ? N blocks is an S  T  (T ? 2M ) matrix. 21

(57)

Using these notations (52) can be written in the following matrix form   H  2   ?E @@(i)ln@P(j ) = 22 Re @@d(i) B B H @@(dj ) :

(58)

We next evaluate @@d(i) for each of the parameters of the deterministic component. We begin with the harmonic components. Let 2 3 j 2[0! +0 ] j 2[0! +0 ] j 2[0!P +0P ] e e : : : e 66 7 ej2[0! +1 ] : : : ej2[0!P +1P ] 777 66 ej2[0! +1 ] ... ... ... ... 66 77 66 j2[0! +(T ?1) ] j2[0! +(T ?1) ] 7 j 2[0!P +(T ?1)P ] 7 e e : : : e 6 77 Eh = 666 ej2[1! +0 ] (59) 77 ; j 2  [1 ! +0  ] j 2  [1 ! +0  ] P P e : : : e 66 77 ... ... ... ... 66 77 66 77 ... ... ... ... 4 5 ej2[(S?1)! +(T ?1) ] ::: : : : ej2[(S?1)!P +(T ?1)P ] 1

1

2

2

1

1

2

2

1

1

1

1

2

1

2

2

2

1

i.e., the ith column of Eh consists of the values of the ith harmonic component evaluated for all (s; t) 2 D. We therefore have h~ = c~ ; (60) where 2 R 3 h ?EhI E =4  I  R 5 ; (61)

Eh

Eh

and EhR = RefEhg, EhI = ImfEhg. Taking the partial derivatives of d~ with respect to the harmonic component amplitude parameters we nd that

@ d~ =  ; @ c~(`) `

(62)

where ~c(`) is the `th element of ~c, and ` is the `th column of . Hence,

@ d = @ dR + j @ dI @ c~(`) @ ~c(`) @ ~c(`) = R` + jI` ; where

R` = [`(0); : : : ;`(ST ? 1)]T 22

(63)

I` = [`(ST ); : : :;`(2ST ? 1)]T :

(64)

 1 = [0; 1; : : : ; (S ? 1)]T 1T  2 = 1S [0; 1; : : : ; (T ? 1)]T ;

(65)

Let

where 1T and 1S are T -dimensional and S -dimensional column vectors of ones, respectively, and

is the Kronecker product. In other words,  1 is the vector of the rst indices of the elements of d in (44), and  2 is the vector of the second indices of the elements of d. Taking now the partial derivatives w.r.t. the harmonic frequencies yields @ hR = ?2diag( ) cR(p)E I + cI(p)E R 1 hp hp @!p @ hR = ?2diag( ) cR(p)E I + cI(p)E R 2 hp hp @p @ hI = 2diag( ) ?cI (p)E I + cR(p)E R 1 hp hp @!p @ hI = 2diag( ) ?cI (p)E I + cR(p)E R ; (66) 2 hp hp @p where diag( 1), (diag( 2) ), is a ST  ST matrix whose diagonal is the vector  1, ( 2). cR(p), (cI (p)); is the pth element of cR, (cI ), and EhRp, (EhIp), is the pth column of EhR, (EhI ). Hence,

@ d = @ hR + j @ hI @!p @!p @!p @ d = @ hR + j @ hI : @p @p @p Similarly for the parameters of the ith evanescent component, we rst de ne 2 j20i 3 66 ej21i 77 77 66 e 6 j 2  2   i = 6 e i 77 ; W 77 66 ... 75 64 ej2(T ?1)i i; Eei = ISS W

where ISS is an S  S identity matrix.

23

(67)

(68)

(69)

We therefore have

~ei = i x~ i ;

(70)

where

2 R 3 ei ?EeIi E 5 (71) i=4  I Eei EeRi ; and EeRi = RefEeig, EeIi = ImfEeig. Taking now the partial derivatives w.r.t. the ith evanescent component frequency yields @ eiR = ?2diag( ) E I x R + E Rx I  2 ei i ei i @ei = ?2diag( 2)eiI @ eiI = 2diag( ) E Rx R ? E I x I  2 ei i ei i @ei = 2diag( 2)eiR : (72) Hence,

@ d = @ ei R + j @ ei I @ei @ei @ei = j 2diag( 2)ei :

(73)

Because the initial conditions vectors si0; i = 1; : : :; I (1;0), are not parameters to be estimated (as they are assumed known), we nd that taking the partial derivatives of d~ with respect to elements of ~si yields @ d~ = @ d~ (74) @ sR(`) @ x~ (` + N ) = i;`+N ; i

i

@ d~ @ d~ = (75) @ sIi (`) @ x~i(S + ` + N ) = i;S+`+N ; where i;` is the `th column of i, sRi(`) denotes the `th element of Refsig, and sIi (`) denotes the `th element of Imfsig. Hence, @ d = @ dR + j @ dI @ sRi(`) @ sRi(`) @ sRi(`) = Ri;`+N + j Ii;`+N ;

24

(76)

and

@d = @ sIi (`)

where R i;` I i;`

h = h =

R i;S +`+N

+j

I i;S +`+N

;

i

? 1) T iT i;` (ST ); : : :; i;` (2ST ? 1) : i;` (0); : : : ; i;` (ST

(77)

(78)

Note that the bound on the variance of estimating the parameters of the 1-D modulating purely-indeterministic processes is given in terms of their non-parametric representation. Substituting (63), (67), (73), (76), and (77) into (58), we obtain the FIM block which corresponds to the parameters of the deterministic component. Since the conditional FIM is block diagonal, the lower bound on the accuracy of estimating the deterministic component parameters is obtained by inverting (58), after the above substitutions have been made.

7 Numerical Examples In this section, we investigate the behavior of the conditional CRB and the performance of the suggested ML algorithm using some speci c examples. First, we investigate the CRB as a function of the spectral support parameters of the harmonic and evanescent components, and the shape of the purely-indeterministic component spectral density. In the second part of this section we illustrate the performance of the ML algorithm by Monte-Carlo simulations, and compare the variance of the ML algorithm estimation errors with the lower bound given by computing the CRB for these examples.

7.1 The Bounds as Functions of the Harmonic and Evanescent Frequencies In this section we rst investigate the bound on the harmonic component parameters, as a function of frequency, for a xed data size of 1616 samples. The harmonic component comprises a single exponential of unit amplitude, and no evanescent components are present. For each of the two di erent AR models of the purely-indeterministic component listed in Table I, the frequency of the exponential is varied in the square region K = [?1=2; 1=2]  [?1=2; 1=2], and the bound on 25

the estimation error variance of each of the harmonic component parameters is computed for each spatial frequency the exponential assumes. The results are illustrated in gures 2 and 3. Note that both for the narrow-band AR eld, and the wide-band AR eld, the shape of the bound as a function of frequency matches the shape of the spectral density of the AR eld. In other words, the lower bound on the estimation error variance of any of the exponential parameters becomes higher, and hence the estimation more dicult, as the local SNR given by

jCk j2 SNR(!k ; k ) = S (ej2! k ; ej 2k )

(79)

decreases. Here, S (ej2!k ; ej2k ) = jB(ej !;ej  )j denotes the spectral density function of the 2-D AR eld. 2

2

2

2

Place Table I here

Place Figure 2 here

Place Figure 3 here In the next example, we investigate the bound on the frequency parameter of a single evanescent component embedded in the above wide-band, and narrow-band AR modeled purelyindeterministic components, as a function of the evanescent component frequency parameter, for a xed data size of 16  16 samples. The evanescent component has spectral support parameters ( ; ) = (1; 0), and its modulating one-dimensional process is a zero mean, unit variance Gaussian white noise process. No harmonic component is present. For each of the two di erent AR models of the purely-indeterministic component listed in Table I, the frequency parameter of the evanescent component,  (1;0), is varied in the interval [?1=2; 1=2], and the bound on its estimation error variance is computed for each value  (1;0) assumes. The results are depicted in Figure 4. Note that both for the narrow-band AR eld, and the wide-band AR eld, the shape of the bound as a function of frequency matches the spectral density of the AR eld. In other words, the lower bound on the estimation error variance of  (1;0) becomes higher as the power 26

of the purely-indeterministic component increases relative to the power of the evanescent eld embedded in it.

Place Figure 4 here

7.2 Performance Examples of the ML Algorithm In this section we illustrate the performance of the ML algorithm by Monte Carlo simulations, and by comparing the variance of its estimation errors with the CRB. The experimental results are based on 100 independent realizations of the purely-indeterministic component, and of the modulating 1-D purely-indeterministic process of each evanescent eld. Similar results were obtained when a single realization of the modulating 1-D purely-indeterministic process of each evanescent eld was used in all 100 experiments. We consider two sets of test data, represented as 32  32 realizations of the elds. Example 1: Consider a eld which consists of the sum of a purely-indeterministic component modeled by the narrow bandwidth 2-D AR model with support S1;1 whose parameters are listed in Table I, two exponentials of equal amplitudes, and two evanescent components. The frequencies of the two harmonic components are (!1; 1) = (0:15; 0:25) and (!2; 2) = (0:16; 0:26), which are far away from the peak of the spectral density of the purely-indeterministic component. The frequencies of the two evanescent components are  (0;1) = 0:1 and  (1;0) = ?0:4. The modulating 1-D purely-indeterministic process s(1;0)(n) is a rst-order AR process with parameter 0:9  exp(j 20:4), whose input is a unit variance Gaussian white noise process. The process s(0;1)(m) is a unit variance Gaussian white noise. Example 2: Consider the following case. The purely-indeterministic component is the wideband AR eld with support S1;1 whose parameters are listed in Table I, there is a single evanescent component and two exponentials of unequal amplitudes. Here, the harmonic component frequencies are (!1; 1) = (0:25; 0:25) and (!2; 2) = (0:26; 0:26), which are close to the peak of the spectral density of the purely-indeterministic component. The evanescent component has spectral support parameters ( ; ) = (1; 0), with 1(1;0) = ?0:4. The modulating 1-D purelyindeterministic process s(1;0)(n) is a rst-order AR process with parameter 0:9  exp(j 20:4), whose input is a unit variance Gaussian white noise process.

Place Table II here 27

The experimental results show that the estimates obtained by the proposed conditional ML algorithm are essentially unbiased as the experimental bias is much smaller than the standard deviation of the experimental results. Furthermore, the results show that the CRB on the error variance in estimating the spectral support parameters is within the 0.95 con dence interval of the experimental variance. Note however that the experimental error variances of the harmonic components' amplitudes are slightly higher than the CRB, probably due to the small data sizes used in these experiments.

8 Conclusions In this paper we have presented a conditional maximum-likelihood solution to the general problem of tting a parametric model to observations from a single realization of a 2-D homogeneous random eld with mixed spectral distribution. The Cramer-Rao lower bound on the accuracy of jointly estimating the parameters of the di erent components was derived, and it was shown that the estimation of the purely-indeterministic and deterministic components are decoupled. The results presented in this paper provide a useful tool for estimating the parameters of 2-D random elds with mixed spectral distributions, and a method for assessing the achievable accuracy of this estimation procedure.

References [1] J. M. Francos, A. Z. Meiri and B. Porat, \A Wold-Like Decomposition of 2-D Discrete Homogeneous Random Fields," Annals Appl. Prob., vol. 5, pp. 248-260, 1995. [2] A. M. Walker, \On the Estimation of a Harmonic Component in a Time Series with Stationary Independent Residuals," Biometrika, vol. 58, pp. 21-36, 1971. [3] S. W. Lang and J. H. McClellan, \The Extension of Pisarenko's Method to Multiple Dimensions," Proc. Int. Conf. Acoust. , Speech, Signal Processing, 1982, pp. 125-128. [4] P. Stoica, R. L. Moses, B. Friedlander and T. Soderstrom, \Maximum Likelihood Estimation of the Parameters of Multiple Sinusoids from Noisy Measurements," IEEE Trans. Acoust., Speech and Signal Process., vol. 37, pp. 378-391, 1989. [5] R. Kumaresan and D. W. Tufts, \A Two-Dimensional Technique for Frequency-Wavenumber Estimation," Proc. IEEE , vol. 69, pp. 1515-1517, 1981. 28

[6] J. L. Doob, Stochastic Processes, John Wiley & Sons Inc., 1953. [7] M. B. Priestley, \The Analysis of Two-Dimensional Stationary Processes with Discontinuous Spectra," Biometrika, vol. 51, pp. 195-217, 1964. [8] H. Helson and D. Lowdenslager, \Prediction Theory and Fourier Series in Several Variables," Acta Math. vol. 106, pp. 175-213, 1962. [9] J. M. Francos, A. Z. Meiri and B. Porat, \A Uni ed Texture Model Based on a 2-D Wold Like Decomposition," IEEE Trans. Signal Proc., vol. 41, pp. 2665-2678. [10] M. P. Ekstrom and J. W. Woods, \Two-Dimensional Spectral Factorization with Applications in Recursive Digital Filtering," IEEE Trans. Acoust., Speech and Signal Processing, vol. 24, pp. 115-128, 1976. [11] G. H. Golub and V. Pereyra, \The Di erentiation of Pseudo-Inverse and Nonlinear Least Square Problems whose Variables Separate," SIAM J. Numer. Anal., vol. 10, pp. 413{432, April 1973. [12] M. R. Osborne, \Some Special Nonlinear Least Square Problems," SIAM J. Numer. Anal., vol. 12, pp. 571{592, September 1975. [13] P. Whittle, \On Stationary Processes in the Plane," Biometrika, vol. 41, pp. 434-449, 1954. [14] C. Chatterjee, R. L. Kashyap and G. Boray, \Estimation of Close Sinusoids in Colored Noise and Model Discrimination," IEEE Trans. Acoust., Speech and Signal Process., vol. 35, pp. 328-337, 1987. [15] S. M. Kay and V. Nagesha, \Spectral Analysis Based on a Canonical Autoregressive Decomposition," Proc. Int. Conf. Acoust. , Speech, Signal Processing, Toronto, 1991, pp. 3137-3140. [16] S. Zacks, The Theory of Statistical Inference, John Wiley and Sons, 1971. [17] B. Porat, Digital Processing of Random Signals, Prantice-Hall, 1994. [18] S. M. Kay, Modern Spectral Estimation, Prantice-Hall, 1988. [19] D. G. Luenberger, Linear & Nonlinear Programming, 2nd ed., Addison-Wesley, 1984. [20] T. W. Anderson, The Statistical Analysis of Time Series, John Wiley & Sons, 1971.

29

FUTURE

θ PRESENT n

m PAST

Figure 1: RNSHP support.

30

Table I: Purely-indeterministic components parameters.

2 b(0,1) b(1,-1) b(1,0) b(1,1) Narrowband AR 1 0:9  exp(j 20:4) 0 0:9  exp(j 20:3) 0:81  exp(?j 20:3) Wideband AR 1 ?0:4  exp(j 20:4) 0 ?0:3  exp(j 20:3) 0:12  exp(?j 20:3)

31

Wide Band AR Spectral Density

-5

x 10 CRB (omega_0)

6 4 2

0 Frequency (omega)

-0.5 -0.5

0.5 0 Frequency (nu)

-5

CRB (nu_0)

x 10

2

0 0.5

0.5

0 omega_0

-0.5 -0.5

0 nu_0

1

0 0.5

0.5

0 omega_0

CRB (Amplitude real comp.)

0 0.5

2

-0.5 -0.5

0 nu_0

-0.5 -0.5

0 nu_0

0.015 0.01 0.005 0 0.5

0.5

0 omega_0

Figure 2: The CRB on the parameters of a single exponential in a wide-band 2-D AR eld.

32

Narrow Band AR Spectral Density

CRB (omega_0)

10000 5000 0 0.5

0.5

Frequency (omega)

-0.5 -0.5

CRB (nu_0)

0.02 0 0.5 omega_0

0.5 -0.5 -0.5

0 0.5

0 nu_0

0.5

0

0 Frequency (nu)

0.04

0

0.02

omega_0

CRB (Amplitude real comp.)

0

0.04

-0.5 -0.5

0 nu_0

-0.5 -0.5

0 nu_0

30 20 10

0 0.5

0.5

0 omega_0

Figure 3: The CRB on the parameters of a single exponential in a narrow band 2-D AR eld.

33

Narrow Band AR PI Component 1.4

0.035

1.2

0.03

1

CRB (Evanescent Freq.)

CRB (Evanescent Freq.)

Wide Band AR PI Component 0.04

0.025

0.02

0.8

0.6

0.015

0.4

0.01

0.2

0.005 -0.5

0 Frequency

0 -0.5

0.5

0 Frequency

0.5

Figure 4: The CRB on the frequency parameter of a single evanescent eld in wide-band and narrow-band AR elds.

34

Table II: Deterministic components estimation results. Parameters First !1 harmonic 1 component Real Imag Second !2 harmonic 2 component Real Imag First Evanescent component  ( ; ) Second Evanescent component  ( ; )

Orig.

0.15 0.25 1 0 0.16 0.26 1 0 0 1 0.1 1 0 -0.4

CRB

Example 1 Bias

6.9865e-08 7.4476e-08 6.2165e-04 3.6350e-03 6.9287e-08 6.7335e-08 6.4339e-04 3.4100e-03 4.5160e-08 5.9676e-08

5.408e-05 5.477e-05 5.412e-03 9.799e-03 5.478e-05 5.463e-05 5.632e-03 1.052e-02 4.292e-05 4.418e-05

Var

7.6399e-08 7.8123e-08 1.0022e-03 5.4130e-03 7.8378e-08 7.7931e-08 1.0662e-03 5.7620e-03 4.8076e-08 6.1000e-08

35

Orig.

0.25 0.35 0.3 0 0.26 0.36 0.5 0 1 0 -0.4 -

CRB

Example 2 Bias

3.3129e-05 3.3636e-05 3.3353e-02 1.6959e-01 1.2598e-05 1.2853e-05 3.2483e-02 1.7993e-01 2.3338e-07 -

1.210e-03 1.234e-03 1.960e-03 1.001e-01 7.828e-04 7.998e-04 1.915e-04 1.002e-01 6.217e-05 -

Var.

3.9733e-05 3.9882e-05 4.3187e-02 3.5453e-01 1.6616e-05 1.7350e-05 4.5418e-02 3.0118e-01 3.0158e-07 -

Captions Figure 1: RNSHP support. Table I: Purely-indeterministic components parameters. Figure 2: The CRB on the parameters of a single exponential in a wide-band 2-D AR eld. Figure 3: The CRB on the parameters of a single exponential in a narrow band 2-D AR eld. Figure 4: The CRB on the frequency parameter of a single evanescent eld in wide-band and narrow-band AR elds. Table II: Deterministic components estimation results.

36

Joseph M. Francos was born on November 6th, 1959 in Tel-Aviv, Israel. He received the B.Sc. degree in Computer Engineering in 1982, and the D.Sc. degree in Electrical Engineering in 1990, both from the Technion-Israel Institute of Technology, Haifa.

From 1982 to 1987, he was with the Signal Corps Research Laboratories, Israeli Defense Forces. From 1991 to 1992 he was with the Department of Electrical Computer and Systems Engineering, at the Rensselaer Polytechnic Institute, as a Visiting Assistant Professor. During 1993, he was with Signal Processing Technology, Palo Alto, CA. In 1993 he joined the Department of Electrical and Computer Engineering, Ben-Gurion University, Beer-Sheva, Israel, where he is now a Lecturer. He also held visiting positions at the MIT Media Laboratory, and at the Electrical and Computer Engineering Department, University of California, Davis. His current research interests are in parametric modeling and estimation of 2-D random elds; random elds theory; parametric modeling and estimation of nonstationary signals; image modeling; and texture analysis and synthesis.

37

Anand Narasimhan was born in Madras, India. He received a B.E. degree in Electrical and

Electronics Engineering from the College of Engineering, Madras, India in 1986, a M.S. degree in Electrical Engineering from the University of Rhode Island in 1988, a M.S. degree in Mathematics and a Ph.D. degree in Electrical Engineering from the Rensselaer Polytechnic Institute in 1991 and 1992 respectively.

He has been with the Personal Systems group at the I.B.M Thomas J. Watson Research Center, New York, since July 1992. His current research interests include wireless communication systems architecture, problems in spectrum estimation, and computer telephony. He was a recipient of the I.B.M. Graduate Fellowship (1990-1992), and is a member of Sigma Xi and the IEEE.

38

John W. Woods received the Ph.D. degree from the Massachusetts Institute of Technology in

1970. Since 1976 he has been with the ECSE Department at Rensselaer Polytechnic Institute, where he is currently Professor and Acting Director of the Center for Image Processing Research. Dr. Woods was co-recipient of 1976 and 1987 Senior Paper Awards of the now IEEE Signal Processing (SP) Society. He served on the editorial board of Graphical Models and Image Processing and was chairman of the Seventh Workshop on Multidimensional Signal Processing in 1991. He served as Technical Program Co-chairman for the First IEEE International Conference on Image Processing (ICIP) held in Austin, Texas in November 1994. He received the Signal Processing Society Meritorious Service Award in 1990. He is currently on the editorial board of the IEEE Transactions on Video Technology. He received a Technical Achievement Award from the IEEE Signal Processing Society in 1994. He is a member of Sigma Xi, Tau Beta Pi, Eta Kappa Nu, and the AAAS. Dr. Woods is a Fellow of the IEEE.

39