4574
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 56, NO. 10, OCTOBER 2008
Multidimensional Multiple-Order Complex Parametric Model Identification Denis Kouamé, Member, IEEE, and Jean-Marc Girault, Member, IEEE
Abstract—This paper presents a way to access both the multiple-order and parameters of a multidimensional complex number autoregressive (AR) model through matrix factorization. The principle of this technique consists of the transformation of the multidimensional model to a pseudo simple-input simple-output AR model, then performing factorization of the covariance matrix of the data. This factorization then leads to a recursive form of the parameter and order estimation. This paper makes two principal contributions. The first is a generalization of one dimensional factored form algorithm, and the second is that it makes it possible to access all the possible different orders and parameters of a multidimensional complex number AR model of any dimension, whereas classical approaches are limited to at most four-dimensional models. Computer simulation results are provided to illustrate the behavior of this method.
mension. The problem of multidimensional order estimation is of great importance in several areas, and many works have been devoted to it—e.g., [7] and [9]. Based on matrix factorization, the technique presented here can yield a straightforward recursive-like form of the algorithm for both dimensions (time, space, ) and orders. It allows easy AR modeling of, for example, three, four (or more) dimension signals. The properties of the algorithm are introduced and its behavior is illustrated and using a numerical example. with II. PRELIMINARY The following notations will be used. Bold type will be used for vectors. Thus, when no confusion is possible
Index Terms—Autoregressive, model, multidimensional, order, parameter.
I. INTRODUCTION
T
HIS paper is a generalization of the works recently reported in [1], where the investigations were restricted to a model with the same order in each dimension. Multidimensional parametric modeling has achieved a high degree of popularity in the literature due to its success in the areas of spectral analysis and image processing. Over 20 years in the field of multidiaumensional parametric modeling, two-dimensional toregressive (2D-AR) modeling has received much attention in many areas, particularly in digital signal- and image-processing areas. Many works have therefore merged, e.g., [2]–[5] and relevant references therein. Although three-dimensional AR (3D-AR) has received some attention, studies concerning more than three dimensions are rare, particularly in contrast to their potential applications [6]–[8]. This is mainly due to the high complexity of the generalization of 2D-AR cases. The aim of this paper is to propose a general framework for estimation of general -dimensional (ND) AR model parameters and orders. The key point here is that, unlike the study reported in [1], the model orders may be different in each diManuscript received March 23, 2007; revised December 12, 2007. First published July 9, 2008; current version published September 17, 2008. The associate editor coordinating the review of this paper and approving it for publication was Dr. Peter Handel. D. Kouamé is with the University Paul Sabatier Toulouse 3, IRIT-UMR 5505, F-31062 Toulouse, France (e-mail:
[email protected]). J.-M. Girault is with the University of Tours, LUSSI, FRE CNRS 2448, F-37032 Tours, France (e-mail:
[email protected]). Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org. Digital Object Identifier 10.1109/TSP.2008.928088
and generally any set of index variables will . Consequently, the value of be denoted will be denoted an -dimensional field at . Also will be denoted and will be denoted . Finally, i.e., ( and ) will be denoted . III. MULTIPLE ORDER MULTIDIMENSIONAL COMPLEX AR IDENTIFICATION A. Methodology Presentation The multiple-order multidimensional model problem may be formulated using multidimensional AR models. A complex-number multidimensional or ND-AR represents , the components of the complex number signal at location , as a linear combination of the complex and an number components additive noise , where and is a set of neighbors excluding (0,0, ,0). This paper is , primarily concerned with developing an approach for . We consider a although the results are also valid when second-order stationary multidimensional complex AR process (ND-AR) defined by (1) is a field of zero-mean constant variance-indepenwhere are complex numdent random noise and the parameters bers and provide a stable system. In the following, attention is focused on the first hyperplane model, without loss of gen-
1053-587X/$25.00 © 2008 IEEE
KOUAMÉ AND GIRAULT: MULTIDIMENSIONAL MULTIPLE-ORDER COMPLEX PARAMETRIC MODEL IDENTIFICATION
erality (the methodology may be applied to the other hyperplanes), i.e., the set of neighbors is . Here, one can have
vectors with dimension 1 to
4575
defined as follows:
(9) (2) i.e., the model orders are not necessarily identical in all direcand the tions. Our aim here is to estimate both the using a recursive-like algorithm, which generalizes matrix decomposition [10], [1]. First define such that
Also define the following vectors in which elements of and are stacked:
(3) Note that Let also
and 1 are part of these vectors.
(4)
is explained below. in (8) is a diagonal matrix. Remark 2: is a column vector. Its dimension is • •
. (10)
th column. is the in (8) contains loss functions • The diagonal matrix for order 0 to , as seen below. The components of its ditypes of periodicity, each type being with agonal have 1. This actually means by definition (see (4),(6), period 1)-periodicity in each dimenand (7)) that there is an ( sion of such a process, as explained in Remark 1. • Accordingly there is the same periodicity in the columns of and in the components of each column. • Due to the structure of the model defined in (3), consists of part or all of the true parameters of the model, depending on the model order and the dimension of , as explained below. Equation (8) is achieved from successive decompositions. Indeed, using (6) in (7)
Thus, (1) is equivalent to
(11) (5)
and finally
From classical least squares estimation theory (see [11]), it is known that
(6) where denotes the Hermitian matrix transpose. Remark 1: From these definitions, it is clear that the defined vectors (and also the matrices defined from these vectors) are all 1 and . For , this is due to the fact pseudoperiodic over the matrix elements have been stacked dimension by dimension. 1) can be seen explicitly in (3). The periodicity over ( It behaves as if the actual model parameters have been filled 1) elements, and thus with zeros to constitute sets of ( and all resulting matrices and vectors have the same type of periodicity. The data covariance matrix can then be defined as
(12) or similarly (13) Thus (11) may also be written
(14)
(7) where Assuming that , the size of this matrix is can then be written in factored form as follows:
.
(15)
(8) where denotes the Hermitian matrix transpose and is an upper triangular matrix with all diagonal elements equal to unity. The elements of this upper triangular matrix are column
(16)
4576
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 56, NO. 10, OCTOBER 2008
is equivalent to a loss function. Indeed, (12) is also equivalent to (17) Defining (18) yields (see Appendix A) (19) and generally (20) Thus, the are positive. Now, returning to (14), , which is in the form of (7), is then iteratively decomposed in the same way. Ultimately, after taking the inverse of the matrices involved in the decomposition, (8) is obtained. and, namely, the inverse of its elements is equivalent Thus to the loss function (see Appendix A) and, according to its periodicity, its minimum in each direction will exhibit the relevant order as in classical estimation theory. is an upper triangular complex As mentioned above, number matrix, the elements of which are . is a diagonal matrix, the elements of which are the in. It is verse of the loss function, as mentioned above, i.e., important to note that the similarity with the case when the order is identical in all the dimensions is only apparent. Here, every thing has to be analyzed dimension by dimension. Since the compohave types of periodicity, nents of the covariance matrix 1 (see Remark 1), our model has the each type with period following characteristics (see Appendixes A and D). • The true parameters as defined in (3) are not directly the but are included in these columns. This is columns of due to the definitions of the model in (3). The true model parameter vector is filled with sets of zeros to construct . , Indeed, assuming the order of the model is of the model are arranged in then the actual parameters , where is the number of the the vector denoted column containing the estimates of the true model parameter vector. That is
(21) “useful” components, This vector has which are actual parameters. For convenience let us assume that
(22)
is without the last set of zeros and provides the i.e., estimates of the true model parameter vector. This is illustrated in Fig. 1, where we chose for convenience , and . Still in the 3-D case, for given (see Fig. 2), the “super blocks” of parameter consists of “blocks” of “subblocks” of useful elements. All these super blocks, blocks, and subblocks are completed by 1). This can be easily sets of zeros so that their size is ( generalized for any , where the number of “subblocks” will be . This parameter is partially or totally included in each column depending on the number of elements column. of , i.e., , as defined The columns of in (9) and (10), consist of part (or all) of the elements of the parameter vector depending on the fact that the dimension of is lower (or greater respectively) than . Clearly, if the dimension of that of the parameter vector is lower than that of the column of , then the model parameters are . This holds provided all included in the model order is found. This is the , as explained below. task played by the elements of Assuming the model order is found, the column of containing the estimated parameters is the first column of such that , where is defined in (22). Due to ( )-periodicity, the “useful” parameters are also present at . A general expression may of course be given for , but it depends on the arrangement of the elements in (3). For instance, if we consider the parameter organization shown in Fig. 2, it can be seen that the number of “useful” components in a 1 1 . Since the total number of desired block is components is , to access them totally we need to choose close to
where Int represents the integer part of the fraction. This expression is modulo . We illustrate this in the example given below (Section IV). • Returning to the model order, the elements of are 1 in each direction (dipseudoperiodic, with a period has types of periodicity, as explained mension); i.e., in Remark 1. This is a direct consequence of arrangement defined in (6). The minimum within of the element of each ( ) type of periodicity constitutes the actual order . Thus, the minimum of each periodicity has to be found to find the model order. It is now clear that the order recursion is achieved from the above. The dimension (time, space, )-recursion may be decomposition. This last recursion is obtained by using obtained by defining and writing the relevant recursions. From (7), it can be written as follows: (23)
KOUAMÉ AND GIRAULT: MULTIDIMENSIONAL MULTIPLE-ORDER COMPLEX PARAMETRIC MODEL IDENTIFICATION
4577
B. Summary of the Proposed Algorithms
Fig. 1. Parameter vector as defined in (3) in the case ; p ; p ) = (2; 3; 4).
N
=
3; m
=
6; (p
Fig. 2. Parameter vector as defined in (3) in the case m; (p ; p ; p ).
N
= 3, for any
Defining the variables as in the one-dimensional case [12], , and , where the asterisk denotes the complex conjugate, at the rank one can be expressed by update expression, the matrix
(24) From these recursions, only elements of with physical meanings are retained. To speed up computations, the recursions may be applied only to the elements with physical meanings. The overall complexity of the parameter estimation algorithm, in terms of number of operations , for data of size , is [10]
Two algorithms can be obtained from the above. The first, referred to as the batch algorithm, deals with parameter estimation in the case of a stationary model, i.e., when model parameters do not vary. The second, referred to as the recursive algorithm, deals with recursive parameters estimation. In both cases, the order is estimated simultaneously. The proposed algorithms can be summarized as follows, given -dimensional complex data with representing an AR model (1). 1) Batch Algorithm: 1. Choose an arbitrary number (maximum possible order) 2. Stack elements of according to (4) to obtain vector according to (6) according to (7) 3. Compute according to (8) and take the inverse of 4. Decompose . Any decomposition may be used. See an example in Appendix B 5. Find the types of periodicity in the (diagonal) elements of 6. Find the minimum in each type of periodicity , con7. Find the column of , the number of which is taining the useful parameters (i.e., the parameters of the model); the number may be chosen as (see Fig. 2)
where denotes the integer part of the fraction, is the direction in which elements of the data are stacked. Due 1)-periodicity, the columns 1 to ( 2 1 also contain useful parameters. Note that this 1 , the number number should not be greater than of columns of . 2) Recursive Algorithm: This algorithm is based on the and Bierman technique [10]. For readability, we set . 1. Same as step 1 of batch algorithm
complex additions complex multiplications complex divisions
2. Set 3. Initialize: and is (25)
where “complex operations” (additions, multiplications, divisions) means operations on complex numbers. To summarize, the computational load of the algorithm is nearly [1.5 1) . ( Remark 3: The methodology developed here may be straightforwardly applied to any general one-dimensional complex models (e.g., ARX, ARMAX, etc.) to estimate their parameters and orders. Finally, the properties in terms of uniqueness, exactness, and convergence of the parameter estimates developed above are proven in Appendixes C and D.
;
;
and ; where identity matrix.
for
is a large number
for for
for
; 4. Construct a vector 5. 6. 7. 8.
by stacking elements of
4578
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 56, NO. 10, OCTOBER 2008
TABLE I COMPARISON BETWEEN THE PROPOSED ALGORITHM AND YULE WALKER ALGORITHM IN 2-D CASE
9. 10. for
11. for
skip to step 11 skip to step 10 are intermediate variables 12. Same as step 6 of batch algorithm 13. Same as step 7 of batch algorithm IV. NUMERICAL SIMULATIONS
TABLE II COMPARISON BETWEEN THE PROPOSED ALGORITHM AND YULE WALKER ALGORITHM IN 2-D CASE, CONTINUATION
For readability, we performed simulations on two-, three-, and four-dimensional cases. A. 2-D Case To illustrate the exactness of the parameter estimation, we compare here the proposed method and well-known Yule Walker algorithm in the complex version; see [13]. The model used was
where is a 256 256 complex field driven by a complex . We took different Gaussian random field with variance values of in 0.01; 0.1; 1; 10 and . The theoretical parameter values and results are shown in Tables I and II. Tables III and IV show the extracted minima in each direction is the versus the standard deviation of the additive noise. minima in the first direction allowing access to , and is the same in the second direction and allows access to . The first lines of the Tables are order 0, the second lines are order 1, etc. As can be seen, for , there is no significant change from line 3, i.e., and similarly . The behavior of is illustrated in the 3-D example below. From these results, it can be seen that the proposed approach has a very good numerical stability and the effect of noise is limited. B. 3-D Case For simplicity in the presentation of the parameters of the model, we set . We
considered the following 3-D complex AR model as defined in (3) by:
where is a 256 256 256 complex field driven by a complex . The parameter Gaussian random field with variance vector is thus arranged as in Table V.
KOUAMÉ AND GIRAULT: MULTIDIMENSIONAL MULTIPLE-ORDER COMPLEX PARAMETRIC MODEL IDENTIFICATION
4579
TABLE III LOSS FUNCTION VERSUS IN THE 2-D CASE, WITH m = 10; p = 2; p = 3
TABLE IV LOSS FUNCTION VERSUS IN THE 2-D CASE, WITH m = 10; p = 2; p = 3
Fig. 3. Loss function (elements of diagonal matrix D m = 6). Thin curve shows all the elements of D . Periodicity in each dimension (here three dimensions) provides the “true” order here (2,3,4). The thick curves illustrate the behavior of D in each dimension.
TABLE V THEORETICAL PARAMETERS IN THE 3-D CASE
Fig. 4. Loss function (elements of diagonal matrix D m = 6) zoomed from Fig. 3. The curves illustrate the behavior of D in each dimension. The minima provide the “true” order here (2,3,4).
and following Fig. 4, the first column containing all “useful” . parameters is Table VI shows the estimated parameters extracted from at column . As can be seen, the results show that the estimates are close to theoretical values. C. 4-D Case For readability, we show a low-order model for this four-dimensional (4-D) case. The used model was
The results are shown in Fig. 3 and Table VI. is shown in . As can be seen, exhibits three types of Fig. 3 with periodicity here. Finding the minimum in each type provides the , as shown “true” order of the model, here in Fig. 4. The number of useful parameters is ,
where is a 32 32 32 32 complex field driven by a complex Gaussian random field with variance . The theoret-
4580
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 56, NO. 10, OCTOBER 2008
TABLE VI ESTIMATED PARAMETERS IN THE 3-D CASE,
tion error. The search for the “true” order of the model consists of looking for the “knee” (see [11]) of the minimum value of or its variants (Akaike the loss function, i.e., information criterion, final prediction error, etc.). From (26)
m=6
(27) and therefore (28)
(29) (30)
ical and estimated parameter vectors are shown in Table VII. The order found was (1,1,1,2). As can be seen, here again the estimates are close to theoretical values of the parameters and the estimated order was correct.
(31) and thus, using (27) (32)
V. CONCLUSION AND GENERAL REMARKS We have presented here a technique for accessing both the multiple-order and parameters of a multidimensional complex number AR model through matrix factorization. This approach is a generalization of the case where the model order is the same in each dimension. The proposed approach is an order recursive algorithm, which allows simultaneous access to the parameters of the model with any order from zero to a given order in each dimension of the model. The dimension (time, space, etc.) recursive form of the algorithm is straightforward. The complexity of the algorithm increases with the number of dimensions of the system. This is the price that has to be paid to obtain a simultaneous estimate of the multiple order and parameters of a general model. The algorithm may be speeded up by using other types of fast computation decomposition. The algorithm is finally illustrated using simulation data in two-, three-, and four-dimensional cases.
This is a possible procedure to decompose a matrix -dimension A in the form . 1. for 2. for
4. for
Let us first analyze the model in a given dimension. Let us and also denote the signal in the th-dimension . This behaves as a one-dimensional system. Let us therefore consider such a system, i.e., a 1-D -order AR model in the form of (26) denotes the regression vector. Let and
APPENDIX B BATCH DECOMPOSITION OF
3. for skip 3 skip 2
APPENDIX A LOSS FUNCTION
where of
Returning to our multidimensional model, the regression vector is -pseudoperiodic, and this periodicity may be seen in each dimension. When considering the th dimension, , (32) holds, as it holds in all the other dimensions. will both decrease periodiThis means that the elements of related to each dimension. cally with period
be an estimate the predic-
skip 4 skip 1 lowertri ; end;
diag
; for
where lowertri(A), diag(A) and denote lower triangular part of A, diagonal of A and complex conjugate symbol, respectively.
KOUAMÉ AND GIRAULT: MULTIDIMENSIONAL MULTIPLE-ORDER COMPLEX PARAMETRIC MODEL IDENTIFICATION
4581
TABLE VII 4-D SIMULATION: THEORETICAL PARAMETERS (WITH p = 1; p = 1; p = 1; p = 2), m = 4; = 1, AND ESTIMATED PARAMETERS. THE ESTIMATED ORDER WAS ALSO fp ; p ; p ; p g = f1; 1; 1; 4g
APPENDIX C UNIQUENESS OF THE DECOMPOSITION
By considering the elements of the matrices, (38) is
Define (33) Assume
(34) , and hence the Therefore, from (7) and (8), , where demodel parameter matrix . notes The loss function matrix (see Appendix A) is thus . From (34), we have (35) Now assume that there exists another parameter matrix , which has exactly the same unit upper triangular structure as but which gives smaller loss function matrix, i.e.,
(39) where denotes the complex conjugate symbol. , Since elements of matrices are loss functions, then as shown in Appendix A. Thus the only possibility for (39) is . This yields , where is the identity matrix. Therefore and thus provides uniqueness. Similar proof in a different context and in a real case may be found in [14]. APPENDIX D EXACTNESS AND ASYMPTOTIC PROPERTIES
(36) or (37) For readability, we ignore the argument in the matrices, and the elements of the matrices will be denoted with the relevant . lower case letter, i.e., . By definition, is a lower triNow define angular matrix with all its diagonal elements . Thus (37) is (38)
The theory presented in Section III-A shows that the parame, are estimates of the ters obtained, namely, the elements of model parameters in the least squares sense. Provided the correct order is found, the model parameter can be extracted from , as shown in Section III-A. the relevant column of Now assume that the data are actually generated by the model of . Let us (5) via a true parameter embedded in column extract the parameters from this column and call them . From (5) and for ease, we write (40) where is a field of zero-mean constant variance-independent random noise. We know from (10) and (13) that the estimated parameter is given by . For readability, define .
4582
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 56, NO. 10, OCTOBER 2008
Thus, from (13)
(41) Defining
and inserting (40) in (41) gives
(42) Remembering that means ( and and, , and ), yields . Since is a field of zero-mean constant variance-independent random noise, when , the term will converge to its expected value, that is, zero, and thus will converge to , the true parameter. ACKNOWLEDGMENT The author would like to thank the anonymous reviewers for their helpful suggestions.
[5] A. H. Kayran and S. R. Parker, “Optimum quarter-plane autoregressive modeling of 2-D fields using four-field lattice approach,” IEEE Trans. Signal Process., vol. 45, pp. 2363–2373, 1997. [6] H. K. Kwan and Y. C. Lui, “Lattice predictive modeling of 3-D random fields with application to interframe predictive coding of picture sequences,” Int. J. Electron., vol. 66, pp. 489–505, 1989. [7] B. Choi, “An order-recursive algorithm to solve the 3-D Yule Walker equations of causal 3-D AR models,” IEEE Trans. Signal Process., vol. 47, no. 9, pp. 2491–2502, 1999. [8] C. L. Nikias and R. R. Raghuveer, “Multi-dimensional parametric spectral estimation,” Signal Process., vol. 9, pp. 191–205, 1985. [9] B. Aksasse, Y. Stitou, Y. Berthoumieu, and M. Najim, “3-D AR model order selection via rank test procedure,” IEEE Trans. Signal Process., vol. 54, no. 7, pp. 2672–2677, Jul. 2006. [10] G. J. Bierman, Factorization Methods for Discrete Sequential Estimation. New York: Academic, 1977. [11] L. Ljung, System Identification. Theory for Users. Englewood Cliffs, NJ: Prentice-Hall, 1987. [12] J. P. Remenieras, D. Kouame, J. M. Gregoire, and F. Patat, “Multiphase pipe flow velocity profile measurements by Doppler ultrasound containing a high level of colored noise,” in Proc. IEEE Int. Ultr. Symp., 1996. [13] S. L. Marple, Digital Spectral Analysis With Application. Englewood Cliffs, NJ: Prentice-Hall, 1987. [14] S. Niu, L. Ljung, and A. Bjorck, “Decomposition methods for solving least-squares parameter estimation,” IEEE Trans. Signal Process., vol. 44, no. 11, pp. 2847–2852, Nov. 1996. Denis Kouamé (M’97) graduated in computer integrated in manufacturing engineering from the Ecole d’Ingénieurs de Tours, Tours, France, in 1992. He received the M.Sc. degree and the Ph.D. degree in signal processing and medical ultrasound imaging from the University of Tours, Tours, France, in 1993 and 1996, respectively. After having been Associate Professor at the University of Tours, he is currently a Professor at the University of Toulouse 3, France. His research interests include mono and multidimensional signal and image processing, ultrasound and medical imaging techniques with relevant applications, in biomedical signal and image analysis.He is regularly involved in the Technical Program committees of major congress in these fields.
REFERENCES [1] D. Kouamé, C. Garnier, J. M. Grégoire, and J. M. Girault, “Multidimensional complex number parametric model order and parameters estimation,” in Proc. IEEE ICASSP, 2006, vol. 3, p. III. [2] B. F. McGruffin and B. Liu, “An efficient algorithm for two dimensional autoregressive spectrum estimation,” IEEE Trans. Acoust., Speech, Signal Process., vol. 37, pp. 106–117, 1989. [3] J. H. McClellan, “Mulitimensional spectral estimation,” Proc. IEEE, vol. 29, pp. 454–456, 1981. [4] C. W. Therrien, “Relation between 2-d and multichannel linear prediction,” IEEE Trans. Acoust., Speech, Signal Process., vol. ASSP-70, pp. 1029–1039, 1982.
Jean-Marc Girault (S’98–M’00) received the master’s degree in signal processing and biological and medical imaging from the University of Angers, France, in 1996 and the Ph.D. degree in “sciences de l’ingénieur” from the University of Tours, Tours, France. He is Lecturer at the University of Tours in signal processing in the I.U.T. of Blois. His research area concerns ultrasound signal processing in biomedical applications and particularly in emboli detection.
