A Unified Performance Analysis of Subspace-Based DOA Estimation Algorithms in Array Signal Processing

In the last decade, the subspace approach has found prominence in the problem of estimating directions of arrival using an array of sensors. Many subspace methods have been proposed and improved; the most attractive ones among these are MUSIC, Min-Norm, State-Space Realization (TAM) and ESPRIT. However, performance analyses are required for justifying and comparing these methods before applying them. Early performance justifications and comparisons were based on simulations. In recent years, many excellent analytical studies have been reported, but these studies have one or more of the following restrictions: (i) assume asymptotic measurements, (ii) analyze some specific parameter perturbation directly instead of through the perturbation of the appropriate subspace, (iii) evaluate individual algorithms using different approximations (so it is hard to compare the analyses of different methods), (iv) involve complicated mathematics and statistics which result in difficult expressions. In our attempt to obtain a unified, nonasymptotic analysis to subspace processing algorithms in a greatly simplified and self-contained fashion, we 1. classify these algorithms into category by the subspace they use orthogonalsubspace processing and signal-subspace processing. We then derive expressions for the first-order perturbation of the signal and orthogonal subspaces using a matrix approximation technique. These formulas provides a common foundation for our analysis of all the DOA estimation algorithms mentioned above. 2. define three approaches by the numerical procedure these algorithms exploit extrema-searching, polynomial-rooting approach, matrix-shifting approach. We establish a common model for each approach and analyze these common models (instead of individual algorithms), and specialize the results for each algorithm. 3. provide a first-order relationship between subspace perturbations and direction-of-arrival perturbations. 4. use the perturbation formulas to derive variance expressions for DOA estimates for all the algorithms. We make the comparisons and discussions among these algorithms and approaches with our theoretical prediction and numerical simulations. The tractable formulas derived in this analysis provide insight into the performance of the algorithms. Simulations verify the analysis.


Abstract
In the last decade, the subspace approach has found prominence in the problem of estimating directions of arrival using an array of sensors. Many subspace methods have been proposed and improved; the most attractive ones among these are MUSIC, Min-Norm, State-Space Realization (TAM) and ESPRIT. However, performance analyses are required for justifying and comparing these methods before applying them. Early performance justifications and comparisons were based on simulations. In recent years, many excellent analytical studies have been reported, but these studies have one or more of the following restrictions: (i) assume asymptotic measurements, (ii) analyze some specific parameter perturbation directly instead of through the perturbation of the appropriate subspace, (iii) evaluate individual algorithms using different approximations (so it is hard to compare the analyses of different methods), (iv) involve complicated mathematics and statistics which result in difficult expressions. In our attempt to obtain a unified, nonasymptotic analysis to subspace processing algorithms in a greatly simplified and self-contained fashion, we 1. classify these algorithms into category by the subspace they use -orthogonalsubspace processing and signal-subspace processing. We then derive expressions for the first-order perturbation of the signal and orthogonal subspaces using a matrix approximation technique. These formulas provides a common foundation for our analysis of all the DOA estimation algorithms mentioned above.
2. define three approaches by the numerical procedure these algorithms exploit -extrema-searching, polynomial-rooting approach, matrix-shifting approach.
We establish a common model for each approach and analyze these common models (instead of individual algorithms), and specialize the results for each algorithm.
3. provide a first-order relationship between subspace perturbations and direction-of-arrival perturbations. 4. use the perturbation formulas to derive variance expressions for DOA estimates for all the algorithms. We make the comparisons and discussions among these algorithms and approaches with our theoretical prediction and numerical simulations.
The tractable formulas derived in this analysis provide insight into the performance of the algorithms. Simulations verify the analysis.

Statement of the Problem
This dissertation analyzes the statistical performance of subspace processing algorithms for estimating directions-of-arrival (DOA) from sensor array signals.
Given an array of L sensors, the problem is to estimate the directions of P uncorrelated plane-wave signals which simultaneously impinge on the array. By taking measurements simultaneously at all sensors, at M different instants in time, an L x M (sensor versus snapshot) data matrix can be formed. The first step in the subspace approach to the estimation of directions-of-arrival is to perform an eigenvalue decomposition on the covariance estimates of the data matrix or, equivalently, perform a singular value decomposition directly on the data matrix.
The second step is to obtain a signal-subspace spanned by the eigen (singular) vectors associated with P largest eigen (singular) values and an orthogonal-subspace spanned by the eigen (singular) vectors associated with L -P smallest eigen (singular) values. Since these two subspaces are orthogonal complements of each other, either of them can be used to estimate the directions of arrival. After acquiring the appropriate subspace, each method estimates the interesting parameters in a different fashion of exploiting the underlying signal propagation model.
In the presence of the noise introduced by observation, quantization, etc., the direction-of-arrival estimators cannot be error-free. Different subspace methods 1 CHAPTER 1. INTRODUCTION 2 have different sensitivities to perturbations. An analytical justification and comparison of their performances is therefore important. Based on a comprehensive study of their statistical error, we develop a unified approach to the performance analyses of these subspace methods.

Scope of the Problem
Array signal processing has been a very active research field for several decades motivated by the increasing demand for radar, sonar, radio, seismic, speech and other signal processing in applications of electronic surveillance, telecommunications, astronomy, geographies, oceanographies, and bio-medical science. In the last ten years, the subspace approach has found its prominence in the problem of estimating directions of arrival, which is one of the most important topics in array processing. Many subspace methods have been proposed and improved. Numerous The first topic addresses the perturbation of the subspaces, which provides the common foundation for tl~e analyses of various algorithms. We then classify those various algorithms into three categories based on their estimation procedures: extrema-searching, polynomial-rooting and matrix-shifting. The analyses of these types of algorithms form the topics of latter three chapters.

Summary of the Other Works
Motivated by the increasing popularity of subspace processing algorithms in DOA estimation, many researchers turned to study the performance of these algorithms.
Early performance justification and comparisons are based on simulations [1], which cannot avoid dependence on the selection of simulation examples. In recent years, analytical performance justification and comparisons have been actively developed.
Among many other excellent analyses, the following works have specific contributions: • Kaveh and Barabell [2] used the first and second moments of the perturbed signal eigenvalues and eigenvectors to derive a first order approximation to the mean and variance of the estimated spectrum for MUSIC, and the mean for Min-Norm based on first order perturbation of the null-spectrum in the signal neighborhoods.
• Porat and Friedlander [3,4] later developed a similar approach for Min-Norm and MUSIC using first-order Taylor series approximations of estimated parameters around their true values.
• Clergeot, Tressens and Ouamri [5] analyzed the performances of MUSIC and Min-Norm based on a first-order linear perturbation method to calculate the perturbed signal eigenvectors.
• Hua and Sarkar [6] presented an analysis of Min-Norm by calculating the perturbation of the pseudo-inverse of the truncated data matrix.
• Stoica and Nehorai [11,12,13] derived an expression for the covariance of MUSIC estimate of the direction of arrival and used it for performance comparison of MUSIC with the Maximum-Likelihood Method and the Cramer-Rao Bound.
• Farrier, Jeffries and Mardani [14,15] presented a second-order analysis to the expected value of MUSIC spectrum.
• Tufts, Vaccaro and Kot [16] proposed a backward error analysis to Min-Norm Linear Prediction of time series using matrix approximation.
Many of these excellent methodologies and results, especially the simplicity and efficiency of the backward error analysis, inspired this research. However, despite of their great successes, most of these performance analyses have one or more of the following restrictions: • assume asymptotic measurements which may not be realistic in practice, • analyze some specific parameter perturbation directly instead of the perturbation of the appropriate subspace, • evaluate individual algorithms using different approximations, and so it lS hard to compare the analyses of the different algorithms.
• involve complicated mathematics and statistics which result in complicated expressions.
Therefore, the goal of this research was to provide a non-asymptotic, unified and simplified performance analysis.

Significance and Contributions
The novelty of the performance analysis presented in this dissertation is that it encompasses a variety of problems and provides a unified framework. A theoretically and computationally simple analysis is common to all the subspace methods, which are formulated in a general framework. In the next few paragraphs, the main results of the present work are listed: The expressions for the perturbation of the signal-subspace and of orthogonalsubspace are derived, which provide a common foundation for our analysis of all the DOA estimation algorithms mentioned above. This approach is also applicable to several other methods beyond subspace processing methods.
2. The performance is analyzed in a greatly simplified and self-contained fashion, which provides the tractable formulae and insight for the algorithms.
Unnecessarily complicated mathematics and statistics have been successfully avoided.
3. The performance is analyzed for the case in which only a finite number of measurements are available, which is the common restriction for many array signal processing applications owing to some practical limitation, for instance, nonstationarity. This work makes the assumption of high signal-to-noise ratio which is reasonable in some practical situations, and is also of theoretical interest. Furthermore, the analyses are shown to agree with simulations down to a fairly low signal-to-noise ratio near threshold, which indicates the high SNR assumption is only sufficient but may not be necessary.
4. To achieve the goal, we have also categorized three different types of subspace methods by their procedure of solving the problem: extrema-searching, polynomial-rooting and matrix shifting.

5.
We have derived a unified framework for different algorithms in each category, which is a necessary step toward unifying all the subspace methods.
6. As a verification, the predicted performance from our analysis is shown to be in good agreement with measured performance from simulations over a wide range of signal-to-noise ratios, over various array geometries, over different observation lengths, over various angular separations.
During the course of this work, the application of the Min-Norm Linear Prediction technique to a sensor array of arbitrary geometry has been documented.
Previous papers describing the Min-Norm technique apply it to linear arrays, and use a polynomial rooting procedure to obtain the estimates. It was shown in [ 17] that the Min-Norm technique is applicable to arrays of arbitrary geometry by using an extrema-searching procedure similar to the MUSIC algorithm.

Organization of the Dissertation
This dissertation is arranged as outlined below: Chapter !!'. Provides a comprehensive review of array processing tracing back to its origin in time series analysis in the 18th century when Prony initiated the linear prediction of time series, up to and including the state-of-art with emphasis on an historical perspective of major trends. The classifications based on subspaces and numerical procedures of the algorithms are described.
Chapter 9: The concept of subspace is introduced using the principal component decomposition. An equivalence of the covariance based approach (which uses the eigendecomposition) and the direct data based approach (which uses the singular value decomposition) is discussed. The perturbation analyses of subspaces developed in this chapter will play a central role through the whole dissertation.
The perturbation of the orthogonal-subspace is presented through a first-order perturbation expansion. However, it is interesting to see that the perturbation of the signal-subspace is not derived in the same way, but rather from the orthogonality between the orthogonal subspace and signal subspace. The results in this chapter provide a common tool with which to analyze all the subspace methods.
Chapter .+: The extrema-searching approach of the orthogonal-subspace based methods is defined, reviewed and analyzed. x, y -sensor position coordinates for an array of arbitrary geometry.

Assumption and Signal Model
Given an array of L sensors and M snapshots of measurements taken at each sensor, the objective is to estimate the directions of P uncorrelated plane-wave signals which simultaneously impinge on the array.
The following assumptions are made throughout the whole dissertation unless being otherwise stated.
• The L array sensors are identical to each other. The number of sensors is greater than the number of signal sources (L > P).
• The P signal wavefronts s1c are narrow-band (with the center wavelength Ac) plane-waves (far-field), uncorrelated from source to source.
• The observation noise ni at each sensor is additive complex white-Gaussian with zero-mean and variance a! (ta! for independent real and imaginary parts), independent from sensor to sensor and from the signals.
Under the above generic assumptions, the signal arriving at the i-th sensor In the presence of additive observation noise, the data matrix is Y=Y+N.
The covariance matrix of the array data vector y is (2.4) where R, = E[ssH ] (P x P) is the covariance matrix of the signal vector. In practice, the covariance matrix is estimated from (2.5)

Pre-Subspace Methods
The theory of array processing has its origin in time-series analysis and spectral estimation, so it can be traced back to the 18th century when Gaspard Riche, the Baron de Prony published his work on fitting a sum of sinusoids to data in 1795 [18]. But most of the Pre-Subspace Methods are restricted to the case of a uniform line array.

Classical Fourier Analysis -Periodogram
In 1898, Schuster introduced the idea of the periodogram for determining the periodicities of meterological phenomenon. At that time the calculation of the periodogram was computationally a very expensive procedure. But with the advent of digital computers and after the discovery of fast Fourier transform (FFT) algorithm [19], the periodogram has become the standard method for frequency estimation and is probably the most frequently used method even in recent days. When the noise in the observed data is Gaussian distributed, the periodogram has been shown to produce a maximum likelihood estimate of the frequency when only a single sinusoid is present. But the periodogram fails to distinguish between two closely spaced frequencies (separated by less than ;~, where N is the number of observations and T is the observation interval) and provides only a single frequency estimate instead of two. The inability to separate the two closely spaced sinusoids is due to the length of the data record ( N < 2; I sin (J" -sin 6; I), and is not related to the signal-to-noise ratio. When enough data is available, the periodogram of sufficiently zero padded data sequence can provide reasonably good estimates.

Conventional Beamforming
The conventional method of mapping the monochromatic field power as a function of angle from the array axis is the phase-to-phase beamforming operation, which consists of phasing the narrow band sensor outputs to a set of plane wavefronts whose directions span those of interest. It employs a procedure known as delayand-sum processing in an attempt to steer a beam in a particular direction. The concept is based on the physics of wave propagation phenomenon whereby sensors at different locations with respect to the incoming plane wave receive (nearly) the same signal, but with different time-delays with respect to some reference due to the propagation path lengths. If delays are inserted at each of the sensors to exactly compensate for the propagation and differential receiver delays, and outputs of the delay elements summed to form a scalar output, the energy in the signal waveforms will add voltage-wise while the measurement noise in each sensor, being uncorrelated from sensor to sensor, adds power-wise. The steering vector can be solved from a nonlinear constrained maximization problem of finding a unit-norm vector w such that the array out. put power is maximum. Such a steering vector has a form of (2.6) The idea behind conventional beamforming as a DOA estimation technique can be seen in the weight vector construction. Given measurements from the sensor array, and the set of known (calculated or calibrated) weight vectors, scan the region of 6-space in which the source may present (i.e., V(J E 0) for relative maximal output power (2.7) The values of fJ which maximizes PBF(O) are DOA estimates.
The beam.forming technique is conceptually and computationally simple, but in the presence of multiple sources, the method breaks down completely because maxima are no longer guaranteed to be at the true source DOAs.

Maximum Entropy Method
In 1967, Burg [20] pointed out that the conventional beam.forming has the problem of utilizing a truncated correlation estimate in lag space which results in a smoothing of true spectral function in frequency space. The extrapolated correlation lags should commit least with respect to observed data. Motivated by the results of information theory, he proposed a constrained maximization problem, i. e., find a P(f) which maximizes subject to the constraints max Jw log P(f)df, P(f) -W r(kLlr ) = 1_: P(f)ei 2 rflctl.r df, for k = -L, ... , L, (2.8) (2.9) where the 21+1 correlation lags are assumed known and Llr = 1/2W. The resulting variational problem obtained by using standard Lagrange multiplier techniques can be written as follows : (2.10) The solution is straightforward U1 P(f) = jjufR-lajj2' (2.11) where R is the Toeplitz correlation matrix of known lags, and def [ jT U1 = 1, 0, ... , 0 , (2.12) (2.13) For uniform linear arrays and m known correlation lags, the maximmientropy method (MEM) is equivalent to order m auto-regressive (AR(m)) mo fitting and Linear Prediction where the extremal problem (2.14) has as its solution the optimal w given by U1, (2.15) which leads to a power spectrum identical to (2.11) . To complete the MM formulation, the DOA estimates are associated with the minima of the null-spMum (or maxima of inverse null-spectrum) (2.16)

Minimum Variance Method
In 1969, Capon [21] argued that, in a filtering interpretation, the cc, iventional beamforming at a given angle gives equal weight to all angles, so that•he undesirable contribution from the other sources is not rejected. He proposd a linear estimator which minimizes the interference from angles outside the regi: >. of interest.
min E[llwHy(t) 11 2 ] subject to wH a(O) = 1, w (2.17) The optimal solution is easily found using the standard Lagrange mult1>lier techniques. Minimizing (2.18) with respect to w yields (2.19) where (2.20) and the DOA estimates are associated with the minima of the null-spectrum (or maxima of inverse null-spectrum) (2.21)

Subspace Based Algorithms
The subspace approach to direction-of-arrival (DOA) estimation has found prominence because of its higher resolution and lower computation compared with other classical methods. It was pioneered by the work of Pisarenko [22].

Subspace Decomposition
The subspace decomposition can be performed on either the direct data matrix by a singular value decomposition (SVD), or on the covariance matrix by an eigenvalue decomposition. These two approaches are essentially equivalent are discussed in Chapter 3. For simplicity, we will use direct data matrix approach.
The singular value decomposition of a data matrix Y (2.22) where the diagonal elements of ~. are arranged as singular values in decreasing order. The largest P singular values correspond to the signal powers (plus noise for noisy case), and the L -P smallest singular values (zero for noise-free case, close to a" for noisy case) are corresponding to the noise power. Therefore the eigenspace spanned by singular vectors is divided into two orthonormal subspaces: the signal subspace spanned by singular vectors associated with largest P singular values and the orthogonal subspace spanned by singular vectors associated with L-P smallest singular values.

MUSIC
The MUSIC [23,24,25] algorithm first estimates the signal subspace from array measurements, then estimates the parameters from the intersections between the array manifold and the estimated signal subspace. This search is typically carried out by computing a weighted Hermitian norm using the direction vectors for each angle of interest and kernel obtained from orthogonal singular vectors of the data matrix. It employs the fact that for signal directions 61i:, · k = 1, · · ·, P, Thus the DOA estimates of the MUSIC algorithm are given by the P zeros of the null-spectrum (2.23) In this and future equations, the symbol (J without a subscript IS a scalar variable which represents a possible direction of arrival, while the subscripted symbol 61i:, k = 1, · · ·, P refers to the actual directions of arrival in the noise-free data. The advantage of MUSIC is that no constraint on the geometry of the sensor array is required. The array pattern can be absolutely arbitrary as long as it is known or calibrated. Though performance advantages are substantial, they are achieved at a considerable cost in computation (search over parameter space) and storage (of array calibration data). Up to now MUSIC is the most popular and often cited work. It has been widely studied and new versions have been developed such as Weighted-MUSIC [26], and Root-MUSIC [1].

Minimum-Norm
The Min-Norm method [27,28] is to apply the linear prediction method to uniform line array data. The first step is to identify a single vector d, with first element equal to unity and minimum Euclidian norm, in the span of the orthogonal subspace.
The polynomial D(z) with coefficient vector d has L-1 roots. The P roots on the unit circle contain the DOA information and the L -P -1 extraneous roots tend to be uniformly distributed within unit circle.
If we partition the signal and orthogonal subspace vectors as follows where gH and cH are the first rows of U, and U 0 , respectively, then from the orthogonality of U, and U 0 , To obtain DOA estimates, a polynomial rooting is then performed, i.e., i=l (2.25) In the case of noisy data, the P roots which are closest to the unit circle are chosen as the signal-roots and the rest are regarded as noise roots. There are other more elaborate signal-root selection procedures. The directions of arrival are found from the angles of signal roots.
The Min-Norm method is also applicable to arbitrary array geometry by searching for the P zeros of the null-spectrum over fJ [17] (2.26)

Root-MUSIC Algorithm
As an improvement of MUSIC algorithms, Root-MUSIC, forms and roots the following null-spectrum polynomial [1] L-1 The polynomial PRM(z) has 2(1-1) roots. Unlike Min-Norm, Root-MUSIC always chooses the P roots with largest amplitudes inside unit-circle. This choice results in a bias in the radial direction of the estimated roots since when white Gaussian observation noise is present the signal-roots will be perturbed inside and outside the unit-circle. However, DOA estimates are only functions of the angles of the roots, not the radii. Thus the radial bias does not affect the DOA estimates obtained by Root-MUSIC.

State-Space Realization
The state space approach includes a covariance approach and a direct-data approach [29,30,31,32]. When autCH:orrelation estimates are used, the covariance matrix is Toeplitz (yielding the TAM algorithm) but this matrix does not have the low-rank property even with noise-free data; for the case of using covariance estimates (R = J.t Y Hy), the covariance approach is equivalent to direct-data approach (sometimes called DDA).
The State-Space approach can be used with an ESPRIT-type array (see next subsection), but what follows is based on a uniform line array for the purpose of analysis.

28)
A state-space model for a plane-wave signal propagating in a sensor array can be derived as where the xis the state-vector with the initial value x 1 (t) = s(t) and h = [1 , · · · , 1] is row-vector of sensor gains, and The diagonal entries of F in a diagonal realization (in general, eigenvalues of the system matrix in an arbitrary realization) contain the information of arrival direction. The data matrix is formed and factored as follows hFL-1 0 and C are respectively the observability matrix and controllability matrix. Let QT be 0 except the first row and oi be 0 except the last row. From (2.31) it is clear that the following shift-invariance is true (2.32) Therefore a transition matrix F is obtained from (2.32) (2.33) In general, the factors of the data matrix are obtained from a singular value decomposition (2.22). By taking the principal components and using the concept of state-variable balancing I o = u.:r:1 (2.34) and then F can be solved for as (2.35) The eigenvalues of F provide the estimates of arrival direction.

ESPRIT
Estimation of Signal Parameters via Rotational Invariant Techniques (ESPRIT) [33,34,35,36,37,38,39,40,41] mitigates the computational and storage requirements of MUSIC by exploiting an underlying displacement invariance structure of an array of sensors. In order to do this, another array which is a displaced version of the first array is needed. This could be implemented by one array with pairwise matched sensor doublets. To describe mathematically the effect of the shift-invariance of the sensor array, it is convenient to describe the array as being comprised of two subarrays identical in every respect although physically displaced from each other by a known displacement quantity~-The signals received at each subarray are The matrix F is diagonal P x P matrix of the phase delays between the subarrays for the P wavefronts and is given by Note that F is a unitary matrix (operator) that relates the measurements from subarray z and those from x. Now define x(t)

(2.37)
A matrix-pencil is formed as (2.38) where Uu (Uaz) are the principal components from U.i (U:z:)· If we pre-multiply (2.38) by a full-rank (P x L) TH, then the generalized eigenvalues of the matrixpencil in (2.38) are equal to the eigenvalues of (2.39) Since we use a uniform line array here to simplify analysis, the best selection of U a:z: and U"' are the first L -1 and the last L -1 rows of U,, respectively. Then Realization [42,43]. Here, we choose T = U az as in [39] to yield (2.40) ESPRIT has certain advantages over MUSIC: • It does not require knowledge of the array geometry and element characteristics.
• It is computationally much less complex.
However, there are still some strong constraints of array geometry: the corresponding sensors in each subarray must have identical characteristics, and they must be equally displaced, in parallel and in the same plane.

7 Matrix-Pencil Method
The ideas that the Matrix-Pencil Method [44,45,46] and ESPRIT method exploit are identical despite the fact that each of them has many different versions.
ESPRIT was proposed from an array geometrical design point of view while the Matrix-Pencil method was developed in applying results from linear algebra to the problem of direction-of-arrival estimation. In this sense, we can say that the

Matrix-Pencil Method is a generalized version of State-Space Realization and ES-
PRIT. The Matrix-Pencil Method forms two matrices with certain relationship in between, then solves the generalized eigenvalue problem of matrix-pencil by employing that relation. Three versions of the Matrix-Pencil Method summarized in [47] are as follows: where .\ can be solved from a generalized eigenvalue problem (2.42)

Y=(~)
and take the principal components from the singular value decomposition to the same data (2.38) using ESPRIT notation. A matrix pencil is then formed as (2.43) and ..\ is solved as an eigenvalue of F in (2.38). This solution is from signal subspace invariance as in ESPRIT (for this reason, we combine the analysis of the Matrix-Pencil method into the analysis ESPRIT).
3. A total least-squares approach [47] (which is also applicable to State-Space Realization and ESPRIT).

Remark
Many researchers have been making efforts to unify these subspace processing al- properties. The polynomial rooting approach is only applicable to the case of a linear equispaced sensor array. Comparative studies of some of above methods can be found in [48,42,43,49,50].

Classification by Subspace
The subpace-based algorithms can be classified into two categories by the subspace which an individual algorithm utilizes: • In the light of the above classification, the analysis therefore starts with the perturbation of each subspace in Chapter 3.

Classification by Numerical Procedure
The subspace-based algorithms can also be classified into three approaches by the numerical procedure which an individual algorithm exploits: • Extrema-Searching Approach: MUSIC and Min-Norm searching algorithm.
We will establish a common model for each approach, analyze the common model and specialize the results to each algorithm in Chapters 4, 5, and 6.

Chapter 3
Perturbation of Subspaces

Introduction
All the subspace processing algorithms utilize the properties of subspaces obtained by subspace decomposition. To analyze the different subspace processing algorithms, it is important to analyze the subspaces. It not only makes the analysis better founded, but also provides a common basis for performance comparison of vector [17], and then examining the perturbations of the orthogonal subspace [51] and· of the signal-subspace [52].
The effects of perturbations on an estimation algorithm can be classified as "above threshold" or "below threshold". An algorithm is said to be above threshold when the mean-squared error (MSE) of the parameters is close to the Cramer-Rao Bound. On the other hand, when the MSE is much larger than the C-R bound, the algorithm is said to be below threshold. The threshold phenomenon is a characteristic of nonlinear estimation algorithms. A threshold analysis of Min-Norm algorithm is given in [53]. In this work, we assume that the SNR is high enough so that the algorithm is not in threshold.

Subspace Review
Assume there are P uncorrelated plane waves simultaneously incident on an L sensor array. Take M snapshot simultaneously at each sensor. The data matrix is formed as where the subscript is spatial index of the sensors and the index in parenthesis is temporal for snapshots. The covariance estimate of the data matrix and its eigenvalue decomposition are where E. are the eigenvectors associated with P non-zero eigenvalues, which span the signal subspace, while E 0 are the eigenvectors associated with the zero eigenvalues, which define the orthogonal subspace. The term "noise-subspace" has been previously used to describe what we call the orthogonal subspace. Since this subspace is the orthogonal complement of the signal-subspace, and since it is well defined whether or not any noise is present in the problem, we believe the term "orthogonal-subspace" is more appropriate.
There are two important properties of the subspaces which the subspace based algorithms exploit: and • The vectors in E, span the signal subspace which is identical to the array manifold (column space of A(O)).
• The vectors in E 0 span the orthogonal subspace which is orthogonal to the array manifold. This orthogonality can be expressed as The above result can be seen from the fact that Even with presence of additive observation noise, under the ergodic assumption, the orthogonal-subspace is still (asymptotically) orthogonal to the signal manifold.
A simple proof can be stated as follows: where e, is the eigenvector ofR associated with eigenvalue .A,. Since A(8)R. (Lx P) has rank of P, there is a P x L matrix D of rank P which satisfies The subspace decomposition can also be performed on the direct-data matrix Y shown in (3.1) by a singular value decomposition. Let the SVD of Y be denoted as y = m:v" = ( u, u,) ( ~· ~ )( ~n.
Then we get the same results as before: • The vectors in U, span the signal subspace which is identical to array manifold (column space of A (8)).
• The vectors in U 0 span the orthogonal subspace which is orthogonal to the array manifold. This orthogonality can be expressed as a(8~JHU 0 = 0, fork= 1, ... ,P.
In a noisy environment, Obviously, we can see and -2 -:E =AM, U=E. (3.8) We assume that the singular vectors are normalized so that UHU =I and iJHiJ = I.
As shown above the subspaces obtained from the eigenvalue decomposition of a covariance matrix and the subspaces from a singular value decomposition of the direct data matrix are the same. In this work, all the algorithms used will be based on the direct-data. They will have the same performance as the covariance-based algorithm described above, with perhaps slightly better numerical properties. One reason that the direct-data algorithm might not be often used in practice is that the dimensions of the matrices grow wit · the data length (number of observations), while the covariance matrix has fixed dimensions (number of sensors). However, for the purpose of analysis, the direct-data algorithm is much easier to deal with.
Because we do not restrict our analysis on uniform line array geometry, only the forward data formulation will be considered in this chapter.

Perturbation of Orthogonal-Subspace
Recall that the columns of U 0 are an arbitrary orthonormal basis for the orthogonal-subspace of the noise-free signal matrix Y, and the columns of U 0 are the estimated orthogonal-subspace vectors associated with smallest singular values :E of the noisy data matrix Y. We can write where 6. U 0 is the perturbation in the estimated orthogonal-subspace vectors. The additive noise in the data matrix Y = Y + N is transformed in a highly nonlinear way by the SVD to produce the noisy singular vectors U 0 [54,55]. It is noted that the perturbations in estimated subspaces which produce errors in DOA estimates are those orthogonal to the original noise-free subspace under consideration.
For orthogonal subpace-based algorithms, only the signal-subspace component of the perturbation in estimated orthogonal subspace is relevant; for signal-subspace based algorithms, only the orthogonal-subspace component of the perturbation in the estimated signal-subspace is relevant. At high SNR, we seek an expression for AU 0 which is linear in the noise matrix N and which is the projection of U 0 onto U •. Such an expression can be found by a forward error analysis using a Taylor series expansion for the SVD. This approach has been used for a related problem and was found that the resulting expression was extremely cumbersome [56]. In our previous derivation, we started with a minimization problem for which U 0 is the solution and by approximating this problem, we obtained a simple expression for AU 0 • This so-called backward error analysis was first used in [16] to develop a performance analysis for linear prediction algorithms suggested by Tufts. It has been found that forward and backward error analyses give about the same performance predictions, but the expressions resulting from the backward analysis are much simpler than those resulting from the forward analysis .
In this section, we give a simple derivation of AU 0 based on perturbation expansion theory suggested by Dr. G. W. Stewart. This derivation relies on the assumption of high SNR (which all other perturbation analyses assume either implicitly or explicitly) .
Let us seek AU 0 in the form of U .P, where P is of order of the noise N (similarly define AVo = v.P) . Then pre-multiply y = ui:vH by U!1 to obtain (3.9) Using the fact that i: 0 = A'E 0 , (3.9) can be written as (3.10) Substituting 6U 0 = U.P and 6V 0 = V.P into (3.10), we get (3.11) By neglecting second order terms and using the fact that U!1Y = 0, we get (3.12) Finally, we post-multiply (3.12) by V. and simplify to obtain (3.13) We can solve for P from (3.13) as follows (3.14) and then we have ~Uo = -U,I:; 1 V,HNHU 0 • (3 .15) Note that (3.15) expresses ~U 0 as a linear function of the noise matrix.
We emphasize the fact that (3.15) is a general expression for the perturbation of the orthogonal-subspace due to additive noise in the data matrix. This expression can be used to analyze the performance of any algorithm which estimates the orthogonal-subspace from data.

Perturbation of Signal-Subspace
We start the investigation of the signal-subspace perturbation by assuming U, = U, + ~ U,. We now seek ~ U 1 in the form of U 0 Q. Using the orthogonality between the the orthogonal-and signal-subspaces, we have -H-U0 U, = 0. (3.16) Equivalently, we have (3.17) Substitute ~U 0 = U,P and ~U, = U 0 Q into (3.17), (3.18) or (3.19) Notice the noise-free signal-and orthogonal-subspaces are orthogonal, so the first term and the last term of above equation are zero. Since U !1U 0 = I and U !1'U, = I, we then have (3.20) and this gives us (3.21) Notice that in deriving the perturbation of signal-subspace, no further approximation was made, so the perturbation expansion for the signal-subspace is as good as the perturbation expansion for the orthogonal-subspace.

Summary
This chapter presents the perturbation analysis of subspaces obtained from a singular value decomposition of a data matrix formed using noisy data. The expressions for the perturbed subspaces are derived using perturbation theory. We first derived the orthogonal-subspace perturbation through a first-order perturbation expansion, then derived the signal-subspace perturbation using the orthogonality between the two subspaces. The columns of the perturbation matrix for the orthogonal-subspace are in signal-subspace while the columns of the perturbation matrix for signal-subspace are in the orthogonal-subspace. To first-order, the perturbation of subspaces is linear in the observation noise. The analysis of the subspace perturbation will provided a common ground for the comparison of various subspace processing algorithms.

Performance of Extrema
Searching Algorithms

Introduction
People had used the classical Fourier analysis and conventional beamforming to process sensor array signals and to estimate the directions of arrival until the late 1960's when Burg's Maximum-Entropy Methods [20] and Capon's Minimum-Variance Method [21] were proposed in order to increase the resolution in array signal processing. In the 1970's, a high-resolution approach -subspace processingemerged, pioneered by Pisarenko [22]. In the late 1970's, Schmidt [23,24] and, independently, Bienvenu and Kopp [25] suggested the MUSIC algorithm which became a landmark of subspace processing. In the early 80's, another important method -Minimum-Norm Linear Prediction was developed by Tufts and Kumaresan [27,28].
Pisarenko's method, MUSIC and Min-Norm were the first three subspace processing methods. They all utilize the orthogonality between the signal manifold and orthogonal-subspace to find the directions of arrival. MUSIC and Min-Norm give much better performance than Pisarenko's method. MUSIC and Min-Norm have gained popularity because of their excellent performance and have been widely used and extensively cited. Also they have attracted most of performance analyses, [2,57,4,58,14,15,11,12,13,59,51] for MUSIC and [17,60,8] for Min-Norm searching algorithm.
Originally, Pisarenko's method and the Min-Norm method used the polynomial rooting procedure based on uniform line array geometry; while MUSIC used the extrema searching procedure based on arbitrary (calibrated) array geometry. Later, a new version of MUSIC -Root-MUSIC was developed by Barabell [1] and Min-Norm Linear Prediction has been applied to arbitrary (calibrated) array geometry by Li,Vaccaro and Tufts [17]. A similar application of Pisarenko's method to arbitrary array geometry is easy to show [see Appendix B].
The major advant age of the searching algorithm is that array geometry is arbitrary as long as the sensor locations are calibrated, i.e. the sensor locations ( In this chapter, we analyze the searching algorithms of MUSIC and Min-Norm.
Two versions of them -a MUSIC algorithm based on direct-data and a Min-Norm algorithm for arbitrary array geometry will be used. A description of these algorithms can be found, respectively, in [17], [51] and [61]. Pisarenko's Method, because of using different data formulation, is analyzed in Appendix B.

Common Model
Now consider a common problem, which has been referred as Weighted MUSIC [26].
The algorithm consists of searching for the P zeros over 0 of the null-spectrum (4.1) where W is the weighting matrix and 0 is a scalar variable which represents a possible direction of arrival. Then by (3.3) Both MUSIC and Min-Norm can be considered as special cases of this algorithm by choosing the weights appropriately.
MUSIC: the weighting matrix is the identity matrix W=I.
Min-Norm: weighting matrix is rank 1 and has the form of Equation (4.1) provides a common basis for our analysis of MUSIC and Min-Norm searching algorithms.

Perturbation of the Angle Estimates
In a noisy environment, the estimated angles of arrival are denoted as perturbations from the true directions of arrival as (4.2) where the D..81c is the perturbation of the k-th direction of arrival. For the purpose of analysis, we adopt the following two steps in analyzing D..81c: 1. Derive the perturbation of the DOA's from P(fJ, U 0 ) by one Newton step initialized at true arrival angle 81c.

Approximate the Newton
Step by a linear function of the orthogonal-subspace perturbation D..U 0 • Step 1. Perturbation of the DOA's from one Newton step of P(81c).
We assume at high SNR that the estimated angles can be found by one step of the Newton algorithm for finding the minima (no longer zeros) of the estimated nullspectrum function P(O, U 0 ) when the algorithm is initialized at a true direction, 81c. The Newton method attempts to find a direction B1c which causes the first derivative of the null-spectrum function P(O, U 0 ) to be zero. The resulting Ok will then minimize P ( 8, U 0 ). This gives the following formula for the estimated angles aP(8t,Uo) ae' (4.3) The first-and second-order partial derivatives of P(O, U 0 ) with respect to f) are computed as follows, and (4 .5) where the superscripts (l) and ( 2 ) stand for first-and second-order derivatives, respectively, with respect to the scalar variable 8. Note that W is not a function of 0, but it could be a function of U 0 , as it is in the case of Min-Norm.
From (    where the nominal term T 0 = N(OA:, V 0 ) = 0, T 1 represents terms in which !::..Vo appears once, T 2 represents terms in which t::..V 0 appears twice, and T3 represents terms in which AU 0 appears three times. For a first-order perturbation analysis, we ignore the terms in T 2 and T 3 and only consider the terms in T 1 • Most of these terms are zero when calculated at()= 61e since they contain the factor a(61e) 8 U 0 • In particular, all terms in T 1 containing AW are zero at () = 61e. Since the nominal term in ( 4.11) is zero, terms in T 1 are the first-order perturbation AN we are looking for. There are only two non-zero terms in T 1 when() = 61e, and the result IS AN -a( 1 l(fJ1e) 8 U 0 W AU!1 a(61e) -a(01e) 8 AU 0 WU!1 a( 1 l(01e) ( 4.12) Finally, using (4.8) (4.10) and (4.12)

Statistical Performance
Since we have used a first-order perturbation analysis, the predicted bias of an    (4.20) In above analysis, only the first-order partial derivative of a( 0) with respect to 0 is needed. When this derivative is evaluated at a true direction 01c, the result is aP) ( 01c) = ei1f-(z; ain 9 1r.+ll; ccl811r.) j 2; (x, cos 01c -Yi sin 01c).
where d is the spacing between adjacent sensors.

Analytical Comparison of MUSIC and Min-N orm
If we divide (  In general, the mean-squared error of the estimated directions-of-arrivals from the Min-Norm algorithm are lower-bounded by the corresponding mean-squared error from MUSIC. In the above equations we are making the correspondence

Numerical Examples
In this section, we give a simulation example to demonstrate the algorithms as well as to verify the performance analysis. A twenty-element uniform circular array (shown in Figure 4.1) is used where the distance between adjacent sensors is one half the signal wavelength. The first sensor of the array is located at the origin of the X-Y plane, and angles are measured with respect to the Y axis. Two sources are considered at 0.2 and 0.5 radians. Twenty snapshots of array data were used to estimate the directions of arrival for a given trial, data matrices of dimension 20 x 20 were formed, and 100 trials were run at each value of SNR. Here SNR is defined as 10 log(::\-) where the er~ is the variance of the complex, additive noise. u,.
The data is generated according to the following formula 2 Yi(n) = L e1¥.(z; sin81+11; coe8 1 )Hkn + ni(n) A:=l where (x,, y,) are the sensor locations, Ac is the center wavelength of the narrowband signal, and <Pi m are independent random phase angles uniformly distributed in the interval (-?r,?r). Figure 4.2 shows a null-spectrum function used in the extrema search approach of subspace based algorithms.
The theoretical mean-squared error for each estimated direction was computed using (4.18) and (4.19) for each trial. The theoretical mean-squared error is a random variable because the signal is random. The mean and standard deviation of the theoretical mean-squared error were computed for each set of 500 trials, and it was found that the standard deviation was less the 5% of the mean. Thus in the figures, we only plot the mean of theoretical mean-squared error.

Summary
In this chapter, a non-asymptotic statistical performance analysis using matrix ap-

Introduction
I_>olynomial rooting is another orthogonal-subspace processing approach to direction of arrival estimation for a uniform line array, which has better performance than extrema searching if both approaches are applied to a uniform line array.
But the searching approach is also valid for arbitrary array geometry, which is an important property in many applications. The rooting approach is similar to the searching approach in many aspects except that the DOA's are determined from the roots of a polynomial formed from the intersection of the signal manifold and the orthogonal-subspace. MUSIC was originally proposed as a searching algorithm; the uniform line array version -Root-Music was developed by Barabell [1]. Pisarenko's Method (see Appendix A) and Min-Norm were proposed as polynomial rooting approaches, and later an extrema search approach was developed for these methods [17,60,8]. Search algorithms look for minima or maxima of a function in one-dimension so that finite sample bias can be induced; however, this is easily seen not to be a limitation of orthogonal-subspace approach per se, but rather a consequence of using a computationally attractive one-dimensional search.
The rooting algorithms are better suited to the geometric nature of the problem in the sense that they actually look for a global minimum or maximum over a multi-dimension space. In other words, all the roots of the polynomial are simultaneously sought. This is important because distinct roots do not imply distinct extrema in the spectrum function. H there are two closely spaced roots, the spectrum may have only one extremum resulting in an apparent loss in resolution using a search procedure. However, the rooting approach can obtain (and resolve) true roots. Therefore, locating the roots, and using their angular location to obtain the DOA's is preferable in the case where the estimation is performed near threshold.
In general, for a polynomial with the order greater than 3, iterative root finding techniques are required, which means the improvement of resolution is gained at a high expense of costly computation. Besides, as long as the roots are distinct, i.e., above threshold, the extrema-search algorithm has the same performance as polynomial-rooting algorithm, as we shall show in this chapter.
We have developed a unified analysis of the polynomial-rooting approach [63] which we will elaborate on in this chapter.
In this chapter, we will first develop a common form of these algorithms for analysis (Section 5.2); then we provide a first-order relation between subspace perturbations and signal-root perturbations (Section 5.3). This is done with some care because the signal roots of the spectral polynomial always have multiplicity two. The statistical performance of the algorithms and the relation to the extremasearch approach, will be derived in section 5.4 and 5.5, respectively. Finally, in Section 5.6, the perturbation of signal roots of the spectral polynomial wil be discussed.

Common Model
Min where W is a weighting matrix and A is a scalar factor. For future reference we define the signal roots of P(z) to be {r 1 }, i = 1, · · ·, P. These roots correspond to the P different sources, and are on the unit circle. The other roots {r 1 }, i = P+l, ···,Lare not on the unit circle. Min-Norm and Root-MUSIC can be obtained as special cases of (5.1) and (5.2) by choosing W and A as follows.

Min-Norm
Note that the choice of Wand A yields the spectral polynomial P(z) corresponding to Min-Norm. Half of the roots of P(z) are identical to the roots of Min-Norm polynomial in (2.25), and the other half are at locations reflected through the unit circle.

(L-P)
W =I and A= llhll2 (5.3) where h is the coefficient vector (with first element of unity) for the polynomial

H(z). H(z) is the causal spectral factor of P(z), namely, P(z) = H(z)H(z-1 )•.
The vector h can be easily obtained from the roots {ri}·

Perturbation of the Angle Estimates
In Chapter 3, a formula was given which showed, to first-order, the perturbation of the derivative of P(z). We are only interested in the roots at z1c = e'A;"" 9 '" k corresponding to noise-free DOA's.
These roots will be on the unit-circle if the null-spectrum has distinct minima, and off the unit circle otherwise. Since it is only the angle of the roots which contains DOA information, it does not matter if the roots are on the unit circle or not, and so the analysis given in the this section involving the derivative of the spectral polynomial is appropriate whether or not the null spectrum has distinct minima.
The discussion above considers the case when the perturbed null-spectrum has distinct minima for each DOA. However, we However, we have pointed out previous that this is not necessary for polynomial rooting to work. In the case when the perturbed null-spectrum polynomial does not possess distinct minima for each DOA, it is still appropriate to look at the roots of the derivative of the spectral polynomial.
In what follows, we calculate the perturbations in the roots of the derivative of the spectral polynomial as a function of the perturbation in the orthogonal subspace. The following steps are used in the calculation: 1. Derive a first-order expression of a~~z) as a function of the perturbation of the orthogonal subspace.
2. Derive a first-order expression for the perturbation of a~~z) as a function of the perturbation in the signal-roots.
3. Equate the expressions derived in 1 and 2 above to get the desired relationship between subspace and signal-root perturbations.
Step 1. Perturbation of the derivative of the spectral polynomial as a function of the orthogonal-subspace perturbations.
We start with The first derivative with respect to z is where a( 1 l(z) = ~~) = ( 0 1 2z · · · (Ll)zL-2 t . (5.5) Substitute U 0 = U 0 + AU 0 and W = W +AW into (5.4) to obtain aP(z, Uo) az (5 .6) where T 0 represents terms in which AU 0 does not appear, T 1 represents terms in which AU 0 appears once, T 2 represents terms in which AU 0 appears twice, and T 3 represents terms in which AU 0 appears three times. For a first-order perturbation analysis, we ignore the terms in T 2 and T 3 • T 0 is the nominal term which equals zero, and so we only need to consider the terms in T 1 • Most of these terms are zero when calculated at z = r1c since they contain the factor a(z-1 )TU 0 • In particular, all terms in T 1 containing AW are zero at z = r1c. There are only two non-zero terms in T 1 when z = r1c, and the result is that to the first order in AU 0 [-z-2 a( 1 l(z-1 )TU 0 WA U!'" a(z) + a(z-1 f A Uo WU!'" a( 1 l(z)]lz=r 1 (5.7) Step 2. Perturbation of the derivative of the spectral polynomial as a function of root perturbations.

• 5
Ar; tl.S z,r1c G z,r1c (5.18) Step 3. Equating the perturbation of the derivative of the spectral polynomial as a function of signal-roots and orthogonal-subspace perturbations.

Relation to Extrema-Searching Algorithms
As shown above, using L'Hospital's rule, we have (5.28) and (5.30). Then the meansquared error expressions of MUSIC and Min-Norm (5.29) (5.31) can be further written as (5.32) and {5.33) These expressions are identical to the mean-squared error expressions derived for extrema-searching algorithms when they are applied to the uniform line array.
Therefore the mean-squared error ratio holds for polynomial-rooting approach

Mean-Squared Error of Root Perturbation
From the previous two sections, we can get

{5.34)
This is from a derivation where ~r.c is linear in the perturbation of orthogonalsubspace, and thus linear in the zero-mean observation noise. Under such an assumption, {5.34) gives us the mean-squared error of estimated roots E(~r.c) 2 •

{5.35)
Min-Norm (RMN): (5.36) Root-MUSIC (RMS): (5.37) However, when the mean-squared error of root estimates is of interest, we need to consider bias. The zero-mean assumption of root estimation error should be avoided. We now present an analysis of the mean-squared error of root estimation error (including the implied contribution of bias squared) of the spectral polynomial (not its derivative). We use similar steps to those used in analyzing the meansquared error of the root estimation. The major difference is we are now dealing with the spectral polynomial itself, instead of dealing with the derivative of the spectral polynomial as we did in the mean-squared error analysis of the angle of the roots.

(5.39)
Here we keep the second-order perturbation terms since we are dealing with the "spectral" function, and we neglect the higher-order perturbation term (for Root-MUSIC, this term is zero anyway). The result is (5.40) Step 2. Perturbation of signal-roots of spectral polynomial Evaluate this function at the signal-roots rA: and neglect the higher order terms.

AG(rA:)
Substitute (3.15) in (5.43) and take the expectation Comparing the (5.37} and (5.45), we find that E(jAr1cl} 2 /Cf is L-P time larger than E(A81c} 2 for MUSIC which indicates a large bias exists in the radial direction as verified in simulations. The radial nature of the error makes the extrema in the spatial spectrum less distinct. This renders procedures that examine extrema using a search procedure less attractive. Notice that (5.36} and (5.46} are identical, which indicates the Min-Norm root estimates are unbiased.

Numerical Examples
A twenty-element uniform line array (with d = 'Ac/2} is used as shown in agrees well with predicted mean-squared error (shown by the lines) . Again, our analysis starts from a high-SNR assumption, but the simulation results show that the analytical expressions give accurate results over a wide range of SNR. We also tried the Min-Norm algorithm using a forward-backward data matrix (instead of (3.1)} and then its performance was improved.

Summary
In this chapter, a non-asymptotic statistical performance analysis using matrix ap- State-Space Realization (sometimes called Toeplitz Approximation Method and Direct Data Approximation) was first suggested by Kung, Arun and Rao in 1983 [29, 30], and further studied by Rao [64,65,10], separately by Foka [31,32] and by Le Cadre [66]. This approach applies the stochastic realization technique to direction of arrival estimation. The property of minimum roundoff noise and coefficient sensitivity of state-space structures was summarized by Jackson for filter design in [67]. Arun proved that a newly developed state-variable balancing technique has the minimum sensitivity with respect to parameter quantization [ 68].

Estimation of Signal Parameters via Rotational Invariant Techniques (ESPRIT)
was proposed by Paulraj, Roy and Kailath [33,34,35,36,37,38,39,40,41], and also studied by Zoltowski [49]. This approach suggests a shift-invariant array which can be considered as two subarrays. The corresponding sensors in each subarray are identical, spaced equally in the same direction so that the second subarray can viewed as a shifted version of the first one. The shift-invariance property is therefore utilized.
The Matrix-Pencil Method was proposed by Ouibrahim, Weiner and Sarkar [44,45,46], and summarized by Hua and Sarkar [47]. It forms two data matrices with a certain relationship between them, and solves the generalized eigenvalue problem for the matrix-pencil by exploiting that relationship. Different matrix pencil methods are obtained by using the delay relation (which results in the same algorithm of ESPRIT), the moving-window relation, or the summation relation [46].
The statistical performance of individual signal subspace algorithms have been analyzed: State-Space Realization [65,69,70,10] and Matrix-Pencil [4 7]. Bhaskar Rao and Hari [50,71] recently presented an interesting analysis of ESPRIT and TAM but their result is based on an asymptotic argument valid only for large data records.
In general, signal-subspace based algorithms utilize a shift-invariant structure of the signal-subspace. So our analysis of performance is based on the signal-subspace perturbation we derived in Chapter 3. Recall that the derivation assumes finite measurements. We begin our analysis of signal subspace methods by establishing a conunon model of shift invariance.

Common Model
Define a transition equation AF=B. (6.1) The transition matrix is solved as (6.2) State-Space Realization A= or and B = o! (6.3) ESPRIT (6.4) If the transition matrix is defined as FA=B (6.5) the solution is where the superscript ~ stands for right pseudo-inverse.

Perturbation of the Angle Estimates
In analyzing the perturbation of estimated DOA, we adopt the following three steps: We can cancel AF and B, and neglect the second-order perturbation term 6A6F to obtain A6F + 6AF = 6B.
Finally, we can solve for 6F as where the superscript t denotes left pseudo-inverse.

If the transition equation is
FA=B a similar result can be derived: where the superscript ~ denotes right pseudo-inverse.

State-Space Realization
The transition equation in the noise-free case for State-Space Realization is

{6.7)
and the transition matrix is {6.8) Therefore the perturbation of the transition matrix is (6.9) where from {2. 35 In addition with ort = E;iu!t and oit = E;iu;t.

ESPRIT
The transition matrix of ESPRIT in the noise-free case is obtained from and in noisy case, from We can expand the above equation to obtain Expand both sides of the above equation yields We can cancel the first terms on both sides, keep the first-order perturbation terms, and drop the higher-order terms to obtain which can be re-arranged as The last term on right hand-side is zero because it satisfies the noise free transition equation, and so we obtain (6.11) Step 2. Perturbation of the eigenvalues of the transition matrix.

State-Space Realization
The eigenvalue perturbation of (6. Ct~[ a{1~f3t ] (6.17) where and where the Doz and Doz are defined from the rows in D 0 corresponding to subarray X and Z.

Statistical Performance
Under the first-order approximation, the bias of an estimated direction of arrival lS = 0.

Numerical Examples
In general, the array geometry of the Matrix Shifting approach requires some type of shift-invariance. For instance, the ESPRIT array shown in Figure 6.1 has a displacement invariance (sensors can be grouped in pairs with identical displacements). In the example in this section, we use the uniform line array with same configuration as in Chapter 5 for the purpose of comparison.  ESPRIT had virtually identical performance (the middle lines on the graphs). We also tried the Min-Norm algorithm using a forward-backward data matrix (instead of (3.1)) and then its performance was the same as SSR and ESPRIT.

Summary
In this chapter, a non-asymptotic statistical performance using matrix approximation to a common model for signal-subspace based matrix-shifting algorithms

Comparative Examples
The configuration of the experiment is: a twenty-element uniform line array (with Under high SNR, all the algorithms appear to be unbiased in DOA estimation, so the important statistical characteristic is the mean-squared error of the estimation. With mean-squared error as the performance measure the subspace based algorithms get ranked as MUSIC, State-Space Realization and ESPRIT, Min-Norm (for the purpose of comparison, we used forward-only data). However, as indicated in a preliminary study of threshold (low SNR) performance, that at low SNR, the estimates' bias, rather than mean-squared error, assumes more importance in determining resolution. In this this case, with the bias as performance measure the algorithms can be ranked in the order of Min-Norm, ESPRIT and State-Space Realization, MUSIC.

Concluding Remarks
In summary, we have • categorized DOA estimation algorithms into orthogonal subspace based algorithms and signal subspace based algorithms according to the subspace being utilized.
• derived a first-order expression of the subspace perturbations induced by observation noise.
• classified the algorithms into extrema searching, polynomial rooting and matrix shifting approaches by the numerical procedures they exploit.
• obtained a linear relationship between perturbation of DOA and perturbation of a given subspace.
• established a general expression of mean-squared error of DOA estimation which can be specialized to all the algorithms by substituting the appropriate parameters.
• demonstrated that the mean-squared error of DOA estimation by extrema searching and by polynomial rooting are equivalent.
• proved that the mean-squared error of DOA estimation by MUSIC is lower than that by Min-Norm.
The following areas should be pursued: • sensor error analysis; • low SNR analysis {threshold analysis); • coherent interference (or multipath time-delay) analysis; • wideband signal performance. Appendix B

Pisarenko's Method
The analysis presented in this research can be easily applied to Pisarenko's method [22]. Pisarenko's Method is the first subspace based estimation algorithm, originately designed for the harmonic· retrieval problem. Pisarenko's Method applied principal component decomposition on a covariance matrix of dimension ' . (P + 1) x (P + 1)]. Then the P + 1-th eigenvector Up+i is in the orthogonalsubspace which is perpendicular to the signal manifold at a source direction. The covariance matrix, limited to be dimension of [(P + 1) x (P + 1)], can be estimated from array data in a spatial smoothing manner of shifted subarray data [73]. Pisarenko's method was originally designed as a polynomial-rooting algorithm for uniform line array, but it can be easily applied to extrema-searching algorithm for the array of arbitrary geometry, like we did for the Min-Norm method.

(B.1)
For a uniform line array the above equation can viewed as a polynomial evaluated · lrrd . e at the signal-roots with z = e'-Te" •m t. The searching algorithm is to find these minima over 8 where 8 is in The statistical performance can be analyzed in the similar procedure as in the analyses for MUSIC and Min-Norm.

Perturbation of Angle and Root
The relation between an arrival angle and a signal root of the characteristic polynomial (or eigenvalue of the signal transition matrix) can be derived as follows: The noise-free signal roots are

Cramer-Rao Bound Calculation
The formula for calculating the Cramer-Rao Lower Bound used in this research is suggested by Clergeot, Tressens and Ouamri. The interested reader can find the detailed derivations in [5]. For uncorrelated signals, the lower bound of meansquared error can be expressed as (D.l) where Pk is the power of k-th signal calculated from k-th diagonal element of Rs = E(ssH). Il 0 is the projector on orthogonal-subspace, defined as (D.2) where U 0 are the left-singular vectors associated with L-P smallest singular values from the noise-free data matrix. The approach we used, independent from data, is