Digital Signal Processing 18 (2008) 835–843 Contents lists available at ScienceDirect Digital Signal Processing www.elsevier.com/locate/dsp Computation of the exact Cramer–Rao lower bound for the parameters of a nonsymmetric half-plane 2-D ARMA model Aydin Kizilkaya ∗ Electrical and Electronics Engineering Department, Pamukkale University, 20070 Kinikli, Denizli, Turkey a r t i c l e i n f o a b s t r a c t Article history: Available online 15 May 2008 Keywords: Cramer–Rao lower bound (CRLB) Fisher information matrix (FIM) Nonsymmetric half-plane (NSHP) Two-dimensional (2-D) Autoregressive (AR) model Moving average (MA) model Autoregressive moving average (ARMA) model Parameter estimation Homogeneous random fields The Cramer–Rao lower bound (CRLB) that gives the minimal achievable variance/standard deviation for any unbiased estimator offers a useful tool for an assessment of the consistency of parameter estimation techniques. In this paper, a closed-form expression for the computation of the exact CRLB on unbiased estimates of the parameters of a twodimensional (2-D) autoregressive moving average (ARMA) model with a nonsymmetric half-plane (NSHP) region of support is developed. The proposed formulation is mainly based on a matrix representation of 2-D real-valued discrete and homogeneous random field characterized by the NSHP ARMA model. Assuming that the random field is Gaussian, the covariance matrix of the NSHP ARMA random field is first expressed in terms of the model parameters. Then, using this matrix structure, a closed-form expression of the exact Fisher information matrix required for the CRLB computation of the NSHP ARMA model parameters is developed. Finally, the main formulas derived for the NSHP ARMA model are rearranged for its autoregressive and moving average counterparts, separately. Numerical simulations are included to demonstrate the behavior of the derived CRLB formulas.  2008 Elsevier Inc. All rights reserved. 1. Introduction Modeling and representation of two-dimensional (2-D) homogeneous random fields based on parametric 2-D models have gathered considerable interest in a wide range of image and signal processing areas. The most popular models used for the characterization and modeling purposes of 2-D random fields are autoregressive (AR), moving average (MA), and autoregressive moving average (ARMA) models. These parametric models perform the estimation processes over some regions of support such as quarter-plane (QP), nonsymmetric half-plane (NSHP), and noncausal plane. As addressed thoroughly in [1], models with NSHP region of support (ROS) yield good spectral estimates and correct power spectrum densities resulting from its efficient spectral factorization and modeling abilities. Thus, in this work we deal with the parametric models having a NSHP ROS. In general, there are eight classes of NSHP ROS [1]. All our discussions are related to the following form:     ℜ+⊕ = (i , j )|i  1, ∀ j ∪ (i , j )|i = 0, j  0 , (1) where (i , j) stands for a pair of integers located on a regular 2-D discrete lattice depicted in Fig. 1. Note that the ROS defined by (1) has infinite dimension. In practice, finite-dimensional case of this support is used with parametric models. 2-D finite-order AR, MA, and ARMA models with NSHP ROS have found numerous applications in texture processing, image coding, image modeling, high-resolution spectral estimation, and filtering, etc. [1–13]. From the parametric modeling * Fax: +90 258 296 3262. E-mail address: akizilkaya@pau.edu.tr. 1051-2004/$ – see front matter doi:10.1016/j.dsp.2008.05.002  2008 Elsevier Inc. All rights reserved. 836 A. Kizilkaya / Digital Signal Processing 18 (2008) 835–843 Fig. 1. NSHP ROS (shaded region). point of view, an important issue that arises in estimation processes is to evaluate the accuracy of algorithms used for estimating the model parameters. A prominent tool to evaluate the performance of any parameter estimation methods from a theoretical manner is the Cramer–Rao lower bound (CRLB). The CRLB sets a lower bound on the variance/standard deviation of the parameter estimates achieved by an unbiased estimation procedure for a given data set [14]. According to the number of data samples, two bounds related to the CRLB are defined. These are asymptotic and exact CRLBs that are convenient for large and small amounts of data samples, respectively. Comparing the variance of parameters obtained by any unbiased estimator to the CRLB provides a benchmark for the achievable estimation accuracy of that estimator. The objective of this paper is to provide a closed-form expression for computing the exact CRLB on the achievable estimation accuracy of the NSHP ARMA model parameters, given a finite dimensional observed realization. To the author’s best knowledge, this problem has not been delved into the literature so far, although its MA counterpart [9,15] and QP ARMA case [16] have been explored. In [9], the problem of estimating the parameters of 2-D real-valued homogeneous MA random fields is investigated and a CRLB computation procedure is also presented for the consistency analysis of the proposed estimator associated with the NSHP MA parameter estimation. An extension of the CRLB procedure in [9] is introduced in [15] for the case of complex-valued MA fields. In [16], a matrix-based formula is derived to compute the exact CRLB on the parameters of a QP ARMA model. Currently, the formulation presented with this work appears as the generalized version of [16]. First, we express the covariance matrix of the observed random field in terms of the NSHP ARMA model parameters, based on a matrix representation of this observed random field. Then, assuming the ARMA field is Gaussian, we derive a closed-form expression for the exact Fisher information matrix (FIM) so that the exact CRLB is obtained by its inverse. In addition, we offer some computation procedures so as to calculate the exact CRLBs on unbiased parameter estimates of the NSHP AR and MA models, by using the main formulas derived for the NSHP ARMA model. The rest of this paper is organized as follows. In Section 2, parametric representation of the NSHP ARMA field and its covariance matrix description by means of the ARMA model parameters are considered. A closed-form expression for computing the exact FIM as well as the exact CRLB on the variance of the unbiased parameter estimates of a NSHP ARMA model is derived in Section 3. In Section 4, the exact CRLB computation procedure presented for the NSHP ARMA model parameters are respectively described for the NSHP AR and MA models. In Section 5, we illustrate the use of derived formulas by some numerical examples. Finally, concluding remarks are drawn in Section 6. 2. Parametric representation of the NSHP ARMA field and its covariance matrix computation Let {x(n1 , n2 ), (n1 , n2 ) ∈ Z } be a real-valued discrete and homogeneous random field represented by a finite-order innovations-driven NSHP-ROS ARMA model. In that case, the (n1 , n2 )th sample of this field satisfies x(n1 , n2 ) = −  (a) (i , j )∈ℜ+⊕ −(0,0) ai , j x(n1 − i , n2 − j ) +  bk,l w (n1 − k, n2 − l), (2) (b) (k,l)∈ℜ+⊕ where a0,0 = b0,0 = 1, and {w (n1 , n2 )} is the innovations field of {x(n1 , n2 )}. The innovations field is also called the driving noise of the ARMA model, and is assumed as a zero-mean, real-valued white Gaussian noise field with variance σ 2 . The (a) (b) respective coefficient sets {ai , j |(i , j ) ∈ ℜ+⊕ } and {bk,l |(k, l) ∈ ℜ+⊕ } correspond to the AR and MA parameters of the NSHP (a) (b) ARMA model. ℜ+⊕ and ℜ+⊕ represent the AR and MA ROSs, that are respectively defined by     (a) ℜ+⊕ = (i , j )|i = 0, 0  j  p 2 ∪ 1  i  p 1 , − p 2  j  p 2 ,     (b) ℜ+⊕ = (k, l)|k = 0, 0  l  q2 ∪ 1  k  q1 , −q2  l  q2 . (3) (4) 837 A. Kizilkaya / Digital Signal Processing 18 (2008) 835–843 (a) (b) From (1), (3), and (4), it is ℜ+⊕ ⊂ ℜ+⊕ and ℜ+⊕ ⊂ ℜ+⊕ . In addition, (p 1 , p 2 ) and (q1 , q2 ) stand for the orders of the AR and MA parts of the NSHP ARMA model, respectively. It should be noted that the samples of the output field (2) must be recursively computable [17]. Taking this fact into account and considering the AR ROS defined by (3), we determine the recursively computable (N 1 × N 2 )-point output field in a vector form        x = x N 1 − 1, N 2 − 1 − ( N 1 − 1) p 2 , . . . , x N 1 − 1, −( N 1 − 1) p 2 , x N 1 − 2, N 2 − 1 − ( N 1 − 2) p 2 ,   T . . . , x N 1 − 2, −( N 1 − 2) p 2 , . . . , x(1, N 2 − 1 − p 2 ), . . . , x(1, − p 2 ), x(0, N 2 − 1), . . . , x(0, 0) , (5) where N 1 and N 2 denote the number of samples in the vertical and horizontal directions of the observed 2-D random field, respectively. In order to compute the sample values {x(n1 , n2 )} of the output field (5), the boundary conditions must be chosen as x(n1 , n2 ) = 0, 0, {(n1 , n2 )|n1 < 0, ∀n2 } ∪ {n1 = 0, n2  −1}, 1  n1  N 1 − 1, n2 < −n1 p 2 . (6) It is worthwhile to note that the boundary conditions (6) are defined by taking into consideration of the recursive computability of the (N 1 × N 2 )-point observed field (5), which are an acceptable set of initial conditions [17]. The driving noise field required to generate the output field (5) with the boundary conditions determined by (6) is therefore defined in a vector form    w = w ( N 1 − 1, N 2 + q2 − 1), . . . , w N 1 − 1, −( N 1 − 1) p 2 − q2 , w ( N 1 − 2, N 2 + q2 − 1),    T  . . . , w N 1 − 2, −( N 1 − 1) p 2 − q2 , . . . , w (−q1 , N 2 + q2 − 1), . . . , w −q1 , −( N 1 − 1) p 2 − q2 . (7) In the light of above statements, the parameter vector of the observed field {x(n1 , n2 )} is then given by θ = [a0,1 , . . . , a0, p 2 , a1,− p 2 , . . . , a1, p 2 , . . . , a p 1 ,− p 2 , . . . , a p 1 , p 2 , σ 2 , b0,1 , . . . , b0,q2 , b1,−q2 , . . . , b1,q2 , . . . , bq1 ,−q2 , . . . , bq1 ,q2 ] T . (8) T In (5), (7), and (8), [·] denotes the transpose operator, and the respective vectors x, w, and θ are N 1 N 2 −, [ N 2 + 2q2 + ( N 1 − 1) p 2 ]( N 1 + q1 )-, and [2p 1 p 2 + p 1 + p 2 + 2q1 q2 + q1 + q2 + 1]-dimensional column vectors. Indeed, the observations (2) can also be rewritten as  ai , j x(n1 − i , n2 − j ) = (a)  bk,l w (n1 − k, n2 − l). (9) (b) (i , j )∈ℜ+⊕ (k,l)∈ℜ+⊕ Using (5)–(7), the difference equation (9) is rearranged in a matrix-vector form as (10) Ax = B w , where A is a N 1 N 2 × N 1 N 2 symmetric matrix and B is a N 1 N 2 × [ N 2 + 2q2 + ( N 1 − 1) p 2 ]( N 1 + q1 ) rectangular matrix, whose (p , r ) entries are respectively defined by [ A ] p ,r = ⎧ ai , j , ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ a (1 − δ p −k,r −l ), ⎪ ⎨ i, j ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ for (i , j ) = (0, 0) and (i , j ) = (1, − p 2 ), p = i ( N 2 + p 2 ) + j + 1, . . . , N 1 N 2 , r = N 1 N 2 − p + i ( N 2 + p 2 ) + j + 1, (a) for (i , j ) ∈ ℜ+⊕ − {(0, 0), (1, − p 2 )}, p = i ( N 2 + p 2 ) + j + 1, . . . , N 1 N 2 , r = N 1 N 2 − p + i ( N 2 + p 2 ) + j + 1, k = N 1 N 2 − l + i ( N 2 + p 2 ) + j + 1, l = N 1 N 2 − cN 2 + 1, . . . , N 1 N 2 − cN 2 + j + p 2 i , c = 1, 2, . . . , N 1 − i − 1, otherwise, (11) 0, ⎧ (b) bk,l , for (k, l) ∈ ℜ+⊕ , ⎪ ⎪ ⎪ ⎨ p = sN 2 + 1, . . . , (s + 1) N 2 , [ B ] p ,r = r = [ N 2 + 2q2 + ( N 1 − 1) p 2 ]( N 1 + k) − p 2 [( N 1 − 1) + s( N 1 − 2)] − q2 (2s + 1) − p + l + 1, ⎪ ⎪ ⎪ s = 0, 1, 2, . . . , N 1 − 1, ⎩ 0, otherwise. (12) x = A −1 B w . (13) From (10), we can write 838 A. Kizilkaya / Digital Signal Processing 18 (2008) 835–843 It should be noted that the matrix A defined by (11) is in the form of an anti-lower triangular matrix and the entries along its anti-diagonal are all a0,0 = 1. Its determinant is always different from zero and has a value of det( A ) = (−1) N1 N2 2 ( N 1 N 2 +3) (14) . −1 As a result of (14), the existence of the inverse of matrix A, i.e., A , is always ensured. Applying the multiplication and expectation operations on (13), as in [16], the covariance matrix of the observed ARMA field is expressed via the model parameters by Γ x = σ 2 A −1 B ( A −1 B )T . (15) It is clear that the observation covariance matrix (15) is symmetric because of Γ x = Γ T x. 3. CRLB on the NSHP ARMA model parameters In this section we derive an exact form of the FIM to compute the exact CRLB on the variance of the unbiased parameter estimates of the NSHP ARMA model. Assuming the driving noise of the ARMA model as a zero-mean white Gaussian noise field with variance σ 2 implies that the observed field {x(n1 , n2 )} will also be zero-mean Gaussian field. Thus, for a real-valued zero-mean Gaussian random field, the (m, n)th entry of the FIM is calculated by [18]    1 −1 ∂Γ x −1 ∂Γ x , Γ J (θ) m,n = tr Γ x 2 ∂θm x ∂θn (16) where tr{·} denotes the trace operator and θ m is the mth element of θ . For the parameter vector θ in (8), the elements of the FIM J (θ) are obtained by 1  (m, n)  2p 1 p 2 + p 1 + p 2 + 2q1 q2 + q1 + q2 + 1. It follows from (16) that we need to compute the partial derivatives of observation covariance matrix (15) with respect to the parameters in (8). 3.1. Partial derivatives of the observation covariance matrix w.r.t. AR parameters Taking the fact Γ x = Γ xT into consideration, the partial derivatives of Γ x w.r.t. the AR parameters {ai , j } are achieved by [16]       T  ∂Γ x ∂A ∂A Γ + A −1 Γ = −σ 2 A −1 , ∂ ai , j ∂ ai , j ∂ ai , j (17) where Γ is the normalized version of Γ x , and is given by Γ = A −1 B ( A −1 B )T . (18) (a) Using the form of (11), the derivatives of the matrix A w.r.t. the AR parameters {ai , j |(i , j ) ∈ ℜ+⊕ − (0, 0)} are expressed by (19) as follows: For (i , j ) = (1, − p 2 ) ∂A [C y ] p ,r = = ∂ ai , j  1, p = i ( N 2 + p 2 ) + j + 1, . . . , N 1 N 2 , r = N 1 N 2 − p + i ( N 2 + p 2 ) + j + 1, 0, otherwise. (19a) (a) For (i , j ) ∈ ℜ+⊕ − {(0, 0), (1, − p 2 )} ⎧ 1 − δ p −k,r −l , r = N 1 N 2 − p + i ( N 2 + p 2 ) + j + 1, ⎪ ⎪ ⎪ ⎪ p = i ( N 2 + p 2 ) + j + 1, . . . , N 1 N 2 , ⎪ ⎨ ∂A k = N 1 N 2 − l + i ( N 2 + p 2 ) + j + 1, [C y ] p ,r = = l = N 1 N 2 − cN 2 + 1, . . . , N 1 N 2 − cN 2 + j + p 2 i , ∂ ai , j ⎪ ⎪ ⎪ ⎪ c = 1, 2, . . . , N 1 − i − 1, ⎪ ⎩ 0, otherwise. (19b) (a) The matrix index y in (19) takes its values to i and j, and is calculated by y = (2p 2 + 1)i + j. Thus, for (i , j ) ∈ ℜ+⊕ − (0, 0), the range of y is to be 1  y  2p 1 p 2 + p 1 + p 2 . Rewriting (17) using the expressions (18) and (19), we have ∂Γ x ∂Γ x = = σ 2 D y, ∂θ y ∂ ai , j (20) where D y = − A −1 C y Γ + ( A −1 C y Γ ) T .   (21) A. Kizilkaya / Digital Signal Processing 18 (2008) 835–843 839 3.2. Partial derivative of the observation covariance matrix w.r.t. the driving noise variance Differentiating (15) w.r.t. the driving noise variance σ 2 , we get ∂Γ x ∂Γ x = = σ 2 D 2p 1 p 2 + p 1 + p 2 +1 , ∂θ 2p 1 p 2 + p 1 + p 2 +1 ∂σ 2 (22) where D 2p 1 p 2 + p 1 + p 2 +1 = C 2p 1 p 2 + p 1 + p 2 +1 = Γ σ2 (23) . 3.3. Partial derivatives of the observation covariance matrix w.r.t. MA parameters The partial derivatives of the observation covariance matrix Γ x upon the MA parameters {bk,l } are given by [16]       T  ∂Γ x ∂B ∂B −1 2 −1 T −1 −1 T =σ A ( A B) + A ( A B) . ∂ bk,l ∂ bk,l ∂ bk,l (24) (b) Considering (12), the derivatives ∂ B /∂ bk,l for (k, l) ∈ ℜ+⊕ − (0, 0) are obtained as follows: ⎧ 1, s = 0, 1, 2, . . . , N 1 − 1, ⎪ ⎨ ∂B p = sN 2 + 1, . . . , (s + 1) N 2 , = (25) [C z ] p ,r = r = [ N 2 + 2q2 + ( N 1 − 1) p 2 ]( N 1 + k) − p 2 [( N 1 − 1) + s( N 1 − 2)] − q2 (2s + 1) − p + l + 1, ⎪ ∂ bk,l ⎩ 0, otherwise. The matrix index z in (25) depends on the values of k and l, and is obtained by z = (2q2 + 1)k + l + 2p 1 p 2 + p 1 + p 2 + 1. Thus, the range of z is 2p 1 p 2 + p 1 + p 2 + 2  z  2p 1 p 2 + p 1 + p 2 + 2q1 q2 + q1 + q2 + 1. Using (25), we can rewrite (24) in a more compact form ∂Γ x ∂Γ x = = σ 2 Dz, ∂θ z ∂ bk,l (26) where D z = A −1 C z B T + (C z B T ) T ( A −1 ) T .   (27) It should be noted that the matrices D i defined by (21), (23), and (27) are all N 1 N 2 × N 1 N 2 symmetric matrices. Moreover, [C y ] p ,r (19) and [C z ] p ,r (25) imply the (p , r ) entries of the matrices C y and C z , respectively. Moreover, δi , j in (11) and (19b) refers to the dirac-delta function defined by δ0,0 = 1 and δi , j = 0 for (i , j ) = (0, 0). Thus, substituting (15), (18), (20), (22), and (26) into (16), we arrive at the following exact formulation:   J (θ ) m,n = 1 2 tr Γ −1 D m Γ −1 D n   (28) which corresponds to a closed-form expression for the FIM of the real-valued, discrete, and homogeneous Gaussian ARMA random field (2). Note that the elements of FIM in (28) are found for 1  (m, n)  2p 1 p 2 + p 1 + p 2 + 2q1 q2 + q1 + q2 + 1. The FIM J (θ) obtained by (28) will be symmetric, since the matrices defined by (18), (21), (23), and (27) are symmetric. As mentioned earlier in [16], the number of computations required to obtain the elements of FIM in (28) can be reduced almost 50% by utilizing the symmetry of FIM. In this respect, (28) is initially calculated for 1  m  2p 1 p 2 + p 1 + p 2 + 2q1 q2 + q1 + q2 + 1 and m  n  2p 1 p 2 + p 1 + p 2 + 2q1 q2 + q1 + q2 + 1. Then, employing the computed results, the other elements of FIM are acquired by [ J (θ)]m,n = [ J (θ )]n,m , for 2  m  2p 1 p 2 + p 1 + p 2 + 2q1 q2 + q1 + q2 + 1 and 1  n  m − 1. We have thus completed the derivation of the closed-form expression for the exact FIM of a homogeneous Gaussian random field characterized by the parameters of a NSHP ARMA model. The exact CRLB on the variance of unbiased estimates of the NSHP ARMA model parameters can now be calculated by CRLB{θ} = diag  −1  J (θ ) , (29) where diag{·} refers to the diagonal elements of a matrix. It is noted that the CRLB on the variance of parameter vector θ works with the diagonal elements of the inverse of the FIM. Remark 1. Some equations introduced up to now are in the same notations of [16] since we pursue the same approach in derivations. On the other hand, the derived exact CRLB computation procedure with this paper is the generalized version of [16] to the NSHP ARMA case. Remark 2. It is also worth noting here that the boundary conditions stated by (6) have admittedly a crucial role on deriving the coefficient matrices and so on the observation covariance matrix. They affect the forms of mentioned matrices as well as the results to be attained by above formulas. In other words, if the boundary conditions are not chosen accurately then the results associated with the CRLBs of the model parameters may not be trusty. 840 A. Kizilkaya / Digital Signal Processing 18 (2008) 835–843 4. Computation of the exact CRLBs for the parameters of NSHP AR and MA models The main formulas presented in the last section can simply be rearranged for the NSHP AR and MA models, separately. At this point, maintaining the same notations and assumptions in Sections 2 and 3, the formulas derived up to now can be reorganized for the NSHP AR model by choosing (q1 , q2 ) = (0, 0) and for the NSHP MA model by taking ( p 1 , p 2 ) = (0, 0), separately. The rearranged formulas for the NSHP AR and MA models are briefly presented as follows. 4.1. NSHP AR model case Choosing (q1 , q2 ) = (0, 0) in (2) we obtain the difference equation related to a NSHP AR random field. In this case, parameters modeling this random field are given as θ (a) = [a0,1 , . . . , a0, p 2 , a1,− p 2 , . . . , a1, p 2 , . . . , a p 1 ,− p 2 , . . . , a p 1 , p 2 , σ 2 ] T . (30) Using (11), the covariance matrix of the NSHP AR field is expressed in terms of the AR model parameters by (a) Γ x = σ 2 A −1 ( A −1 )T . (31) Normalizing (31) with the driving noise variance 2 σ gives Γ (a) = A −1 ( A −1 )T . (32) Next, replacing the expressions (8), (15), and (18) with (30), (31), and (32), respectively, and using them in (17), (20)–(23), and (28), the elements of FIM (28) for the parameters of the NSHP AR field are achieved by 1  (m, n)  2p 1 p 2 + p 1 + p 2 + 1. At last, the exact CRLB on the variance of NSHP AR field parameters is found by (29). 4.2. NSHP MA model case The formulas presented before sections can also be employed for a NSHP MA random field, by ( p 1 , p 2 ) = (0, 0). In that case, the difference equation (2) turns into the difference equation of a NSHP MA random field. Hence, the parameter vector modeling this random field is in the form of θ (b) = [σ 2 , b0,1 , . . . , b0,q2 , b1,−q2 , . . . , b1,q2 , . . . , bq1 ,−q2 , . . . , bq1 ,q2 ] T . (33) For ( p 1 , p 2 ) = (0, 0), the AR coefficient matrix A in (15) is vanished and thereby the covariance matrix of the NSHP MA field is stated by means of the MA model parameters as follows: (b) Γ x = σ 2 B BT (34) where B is the N 1 N 2 × ( N 2 + 2q2 )( N 1 + q1 ) rectangular matrix whose (p , r)th entry ⎧ (b) bk,l , for (k, l) ∈ ℜ+⊕ , ⎪ ⎪ ⎪ ⎨ s = 0, 1, 2, . . . , N 1 − 1, [ B ] p ,r = p = sN 2 + 1, . . . , (s + 1) N 2 , ⎪ ⎪ ⎪ r = ( N 2 + 2q2 )( N 1 + k) − q2 (2s + 1) − p + l + 1, ⎩ 0, otherwise (35) is obtained by taking ( p 1 , p 2 ) = (0, 0) in (12). Thus, using (35), the partial derivatives of the observation covariance matrix (34) w.r.t. the parameters given by (33) are found as follows: (b) (b) ∂Γ x ∂Γ x = = σ 2 D1, ∂θ 1 ∂σ 2 (36) (b) (b) ∂Γ x ∂Γ x = = σ2Dv, ∂θ v ∂ bk,l (b) (k, l) ∈ ℜ+⊕ − (0, 0), (37) where D1 = C 1 = Γ (b) σ2 (38) , Γ (b) = B B T , T (39) T T D v = C v B + (C v B ) ⎧ ⎪ 1, ⎨ ∂B = [C v ] p ,r = ∂ bk,l ⎪ ⎩ 0, , (40) s = 0, 1, 2, . . . , N 1 − 1, p = sN 2 + 1, . . . , (s + 1) N 2 , r = ( N 2 + 2q2 )( N 1 + k) − q2 (2s + 1) − p + l + 1, otherwise. (41) 841 A. Kizilkaya / Digital Signal Processing 18 (2008) 835–843 Table 1 CRLB on the parameters of ARMA field Parameters Data sizes → (N 1 , N 2 ) σ2 = 1 a0,1 = 0.09 a1,−1 = −0.3 a1,0 = 0.05 a1,1 = 0.1 b 0, 1 = − 0 . 3 b1,−1 = 0.2 b1,0 = −0.5 b1,1 = 0.4 (10, 10) (20, 20) (30, 30) (40, 40) 0.02490 0.04709 0.03506 0.05392 0.03585 0.04925 0.04275 0.04764 0.03991 0.00521 0.01149 0.00735 0.01205 0.00823 0.01066 0.00765 0.00863 0.00675 0.00226 0.00512 0.00309 0.00521 0.00358 0.00458 0.00306 0.00357 0.00269 0.00126 0.00289 0.00170 0.00289 0.00200 0.00254 0.00163 0.00195 0.00144 (10, 10) (20, 20) (30, 30) (40, 40) 0.14142 0.10227 0.09084 0.11337 0.10477 0.07071 0.04965 0.04344 0.05362 0.04740 0.04714 0.03279 0.02855 0.03512 0.03063 0.03536 0.02448 0.02126 0.02611 0.02263 Table 2 CRLB on the parameters of AR field Parameters Data sizes → (N 1 , N 2 ) σ2 = 1 a0,1 = 0.25 a1,−1 = −0.48 a1,0 = 0.1 a1,1 = 0 In (41), [C v ] p ,r denotes the (p , r)th entry of the matrix C v , i.e., C v = ([C v ] p ,r ) N 1 N 2 ×( N 2 +2q2 )( N 1 +q1 ) that corresponds to the (b) derivative of matrix B with respect to the MA parameters {bk,l |(k, l) ∈ ℜ+⊕ − (0, 0)}. Also, the matrix index v in (40) and (41) takes its values upon v = (2q2 + 1)k + l + 1, i.e., the range of v is 2  v  2q1 q2 + q1 + q2 + 1. Replacing the respective expressions (8) and (15) with (33) and (34), (16) turns into  J (θ (b) )  = m,n 1 2 tr   (b) −1 ∂Γ (xb)  (b) −1 ∂Γ (xb) Γx . Γ x (b) ∂θ m ∂θ n(b) (42) Thus, substituting (34), (36), and (37) into (42), the closed-form expression for computing the FIM elements of a real-valued zero-mean Gaussian MA field is given by  1  (b) −1 tr (Γ ) D m (Γ (b) )−1 D n . (43) 2 Note that calculation of (43) are realized over 1  (m, n)  2q1 q2 + q1 + q2 + 1, by exploiting the expressions defined by (38)–(40). Finally, using (33) and (43) instead of (8) and (28) in (29), the exact CRLB on the variance of the parameters of a NSHP MA field is accomplished. Remark 3. It is interesting to note that the CRLB computation procedure originated for the parameters of a NSHP MA field is similar in nature with previously proposed computation procedure in [9] for the NSHP case. It is also possible to argue this fact from the expressions (5)–(7), which are of central importance on defining the coefficient matrices required for the related calculations. In this respect, assigning ( p 1 , p 2 ) = (0, 0) in (5) and (7) will provide the same notations offered in [9]. Moreover, in this case, the boundary conditions defined by (6) is to be convenient for the representations in [9]. This statement is also confirmed in the following section by numerical examples. Consequently, the rearranged formulas will be valuable to make a performance analysis of some parameter estimation techniques such as the ones given in [10].  J (θ (b) )  m,n = 5. Numerical examples We studied on some parametric models for computing the exact CRLBs on its parameters, as a function of the size of data fields by choosing ( N 1 , N 2 ) = {(10, 10), (20, 20), (30, 30), (40, 40)}. In Tables 1–3, respectively, the results of CRLBs on the standard deviations of the parameters of NSHP ARMA, AR, and MA models are presented by computing (CRLB)0.5 . Note that the MA model parameters considered in Table 3 are taken from [9,10]. It is observed from Table 2 that numerical results achieved by (N 1 , N 2 ) = (30, 30) are the same with [9] supporting the assertion given by Remark 3. 6. Conclusions In this paper, we have developed a closed-form expression for computing the exact CRLB on the variance of unbiased parameter estimates of a NSHP ARMA model. In this framework, we firstly expressed the covariance matrix of a NSHP ARMA 842 A. Kizilkaya / Digital Signal Processing 18 (2008) 835–843 Table 3 CRLBs on the parameters of the MA fields addressed in [9,10] Parameters [9] [10] Data sizes → (N 1 , N 2 ) (10, 10) (20, 20) (30, 30) (40, 40) b 0, 1 = − 0 .9 b1,−1 = 0.1 b1,0 = −0.5 b1,1 = 0.4 0.19721 0.19814 0.11110 0.13239 0.14464 0.07413 0.04725 0.04916 0.06210 0.05427 0.04778 0.02318 0.03171 0.04063 0.03343 0.03556 0.01495 0.02344 0.03021 0.02426 σ2 = 1 b 0, 1 = 0 .7 b0,2 = −0.14 b0,3 = −0.048 b1,−3 = 0 b1,−2 = −0.48 b1,−1 = −0.4496 b1,0 = 0.268 b1,1 = 0.1 b1,2 = 0 b1,3 = 0 0.21211 0.19096 0.14727 0.13670 0.15034 0.15917 0.17770 0.14053 0.13631 0.13412 0.13478 0.07443 0.05654 0.06356 0.05591 0.05173 0.06154 0.06414 0.06347 0.06317 0.05811 0.05166 0.04791 0.03515 0.04145 0.03516 0.03100 0.03806 0.04022 0.04132 0.04103 0.03706 0.03204 0.03562 0.02581 0.03083 0.02564 0.02214 0.02756 0.02942 0.03067 0.03037 0.02720 0.02324 σ2 = 1 random field through the model parameters. Then, this derivation is employed in order to derive an exact CRLB computation procedure for the parameters of a Gaussian ARMA random field represented by the NSHP ARMA model. It is shown that the main formulas developed for the NSHP ARMA model can be exploited to compute the exact CRLBs for the parameters of NSHP AR and MA models, individually, by making simple rearrangements on them. It should be noted that the proposed CRLB computation procedures are easy to use especially for the estimation processes with small-size data fields, and it is hoped that the proposed formulas will be valuable for assessing the performance of algorithms dealing with the parameter estimation of NSHP AR, MA, and ARMA models. Acknowledgment The author is thankful to the anonymous referees for their constructive comments and suggestions on the presentation of this paper. References [1] M.P. Ekstrom, J.W. Woods, Two-dimensional spectral factorization with applications in recursive digital filtering, IEEE Trans. Acoust. Speech Signal Process. 24 (1976) 115–128. [2] T.L. Marzetta, Two-dimensional linear prediction: autocorrelation arrays, minimum-phase prediction error filters, and reflection coefficient arrays, IEEE Trans. Acoust. Speech Signal Process. 28 (1980) 725–733. [3] J.M. Francos, A.Z. Meiri, B. Porat, A unified texture model based on a 2-D Wold-like decomposition, IEEE Trans. Signal Process. 41 (1993) 2665–2678. [4] J.M. Francos, A. Narasimhan, J.W. Woods, Maximum likelihood parameter estimation of textures using a Wold-decomposition based model, IEEE Trans. Image Process. 4 (1995) 1655–1666. [5] J.M. Francos, A. Narasimhan, J.W. Woods, Maximum-likelihood parameter estimation of discrete homogeneous random fields with mixed spectral distributions, IEEE Trans. Signal Process. 44 (1996) 1242–1255. [6] A.H. Kayran, Two-dimensional orthogonal lattice structures for autoregressive modeling of random fields, IEEE Trans. Signal Process. 44 (1996) 963–978. [7] T. Nakachi, K. Yamashita, N. Hamada, Asymmetric half-plane lattice modeling based on 2-D Levinson algorithm, IEEE Trans. Circuits Syst. II Analog Digit. Signal Process. 44 (1997) 865–868. [8] A.H. Kayran, I. Erer, Optimum asymmetric half-plane autoregressive lattice parameter modeling of 2-D fields, IEEE Trans. Signal Process. 52 (2004) 807–819. [9] J.M. Francos, B. Friedlander, Parameter estimation of two-dimensional moving average random fields, IEEE Trans. Signal Process. 46 (1998) 2157–2165. [10] A. Kizilkaya, On the parameter estimation of two-dimensional moving average random fields, IEEE Trans. Circuits Syst. II Exp. Briefs. 54 (2007) 989–993. [11] R.L. Kashyap, Characterization and estimation of two-dimensional ARMA models, IEEE Trans. Inform. Theory 30 (1984) 736–745. [12] L. Xu, M.R. Azimi-Sadjadi, Two-dimensional modeling of image random field using artificial neural network, in: Proc. IEEE Int. Conf. Acoust., Speech, Signal Process. (ICASSP ’93), Minneapolis, USA, 1993, pp. 581–584. [13] Y.-S. Chung, M. Kanefsky, On 2-D recursive LMS algorithms using ARMA prediction for ADPCM encoding of images, IEEE Trans. Image Process. 1 (1992) 416–422. [14] S.J. Orfanidis, Optimum Signal Processing, second ed., McMillan, New York, 1988. [15] J.M. Francos, Cramer–Rao bound on the estimation accuracy of complex-valued homogeneous Gaussian random fields, IEEE Trans. Signal Process. 50 (2002) 710–724. [16] A. Kizilkaya, A.H. Kayran, Computation of the exact Cramer–Rao lower bound for 2-D ARMA parameter estimation—I: The quarter-plane case, IEEE Trans. Circuits Syst. II Exp. Briefs. 53 (2006) 23–27. [17] D.E. Dudgeon, R.M. Mersereau, Multidimensional Digital Signal Processing, Prentice Hall, Englewood Cliffs, NJ, 1984. [18] B. Porat, B. Friedlander, Computation of the exact information matrix of Gaussian time series with stationary random components, IEEE Trans. Acoust. Speech Signal Process. 34 (1986) 118–130. A. Kizilkaya / Digital Signal Processing 18 (2008) 835–843 843 Aydin Kizilkaya received the B.Sc. degree from the Black Sea Technical University, Turkey, in 1994, and the M.Sc. degree from the Pamukkale University, Turkey, in 1997, both in electrical and electronics engineering, and the Ph.D. degree in electronics and communication engineering from the Istanbul Technical University, Turkey, in 2006. In 2002, he joined the Istanbul Technical University, and he was appointed as a Research and Teaching Assistant from 2002 to 2006. He is now in Pamukkale University as an Assistant Professor of electrical and electronics engineering. His research interests include statistical digital signal processing, multidimensional signal processing, Cramer–Rao lower bounds for parameter estimation, spectrum estimation, system identification, and circuits and systems.