Eigenvalue Pairing for Direction Finding with Vector Sensor Arrays

This thesis introduces a novel sorting pairing method for the azimuth-elevation estimates obtained from the Estimation of Signal Parameters via Rotational Invariance Techniques (ESPRIT) based closed-form source localization algorithm. The ESPRIT algorithm estimates Direction of Arrival (DOA) angles from a given source using arbitrarily spaced three-dimensional arrays of vector sensors, whose locations need not be known. In estimating the DOAs there can be miss-pairings of DOA azimuth and elevation direction cosines in certain cases. In such cases, the incident angles calculated can be permuted which produces miss matching of angles or misses. The sorting pairing method presented exploits the order of eigenvalues and thereby the DOAs estimated from the ESPRIT algorithm. The sorting pairing method is compared to two other pairing methods: the traditional pairing method, a method that is simplistic in nature when dealing with misses; and the exhaustive pairing method, a method that exhaustively attempts to find the eigenvalue pairing that yields the least amount of error. Simulation results provides strong evidence that the sorting pairing method gives good results over a much larger range of signal-to-noise ratios than the traditional method, and has an accuracy that is nearly co-linear with the exhaustive method.

5. Chapter 5 concludes the thesis, and evaluates the research and its contribution. This chapter also discusses implications of the results and additional areas for potential research.

Hydrophones
Most detection systems in air or free space use electromagnetic waves, but these generally absorb quickly in salt water; however, sound waves, unlike electromagnetic waves, can travel great distances, sometimes thousands of kilometers, before being dissipated underwater. These underwater sound waves are reflected by objects and are produced by machinery, making sound useful for active and passive detection. Some undersea applications reject active sonar for two reasons: • Active sonar pulses travel far beyond the maximum detection range; thus, targets can intercept these pulses at great distance and avoid detection.
• Active sonar is not covert, and sources actively transmitting energy in the water may be located and classified.
In order to achieve sensing underwater sound waves, one would require a hydrophone, a common sensor employed for sonar, which is at its barest an underwater microphone. Like an in air microphone, a hydrophone measures pressure only and sound waves passing over a hydrophone introduce changes that are measured and used for detection. Hydrophones that have no directionality are called omni-directional hydrophones and they are quite common due to the fact that they are easy to build, maintain, and analyze. Omnidirectional hydrophones measure the acoustic pressure as a scalar quantity.

Vector Sensors
Performance of a configuration of sensors can be improved upon by increasing the amount of information that is measured by each sensor. For acoustic measurements, the particle velocity of the medium in which the sound is traveling through can provide additional information about the direction of sound arrival. This type of sensor that measures the particle velocity is known as a vector sensor.

Description of Vector Sensors
Acoustic vector sensors measure the pressure as well as non-scalar, vector, components of an acoustic field such as a particle velocity. The non-scalar components of an acoustic field cannot be obtained using a single acoustic pressure sensor which has no directionality, an omni-directional sensor, but can be by a single acoustic vector sensor which contains one omni-directional hydrophone measuring pressure and two or three orthogonal velocity hydrophones measuring the components of particle velocity [1]. In the past few decades extensive research has been conducted on the theory and design of vector sensors [2] [3], and vector sensors have been primarily used for underwater target localizations, acoustic communications, and SOund Navigation And Ranging SONAR applications.
The concept of vector sensors has been within the acoustic community since the 1930s with the seminal paper by H.F. Olson [4]. In this work Olson provides the set up and derivation of how a vector sensor would work. Olson describes a velocity sensor structure that is a tube that contains a magnetic mass suspended by springs where any vibration along the axis of the tube will cause the mass to move; thus, inducing a current in a wire coil. This induced current yields a measurement of velocity along the velocity sensor axis. Olson proved that the response of a high quality microphone with uniform sensitivity over a wide frequency range set up in such a system can measure the velocity component of a sound wave.
Vector sensors can be categorized into two general categories: inertial and gradient [5]. Inertial vector sensors measure the acceleration or velocity by responding to particle motion, while gradient sensors use a finite-difference approximation to estimate the gradients and thereby the direction of the acoustic field. There are obvious benefits and drawbacks to each category of sensor, but inertial sensors offer a broad dynamic range. However, proper supporting and packaging of the sensor without affecting its response to motion is an issue. This is due to the fact that inertial sensors do not distinguish between acoustic and non-acoustic sources such as package vibrations; therefore, they must be properly shielded from such disturbances. This shielding, and supporting and packaging issue makes it difficult to manufacture accurate and small inertial sensors for higher frequencies. Gradient sensors, unlike inertial sensors can be manufactured at significantly smaller sizes as they do not suffer the same drawbacks as their counterpart, inertial sensors. Gradient sensors smaller size makes them more suitable for higher frequencies; however, the drawback to gradient sensors is that their method of determining the 'vector' of an acoustic source, finite-difference approximation, limits their operational use, e.g. their dynamic range.
The proposed use of vector sensors takes advantage of vector sensor components at the receiver and are not limited to a particular sensor type, inertial or gradient.

Current Applications of Vector Sensors
Vector sensors, despite their history dating back to 1931 [4], took 38 years to have interest from the United States Navy where they have been used in acoustic sensor systems such as the Directional Frequency Analysis and Recording (DIFAR), the passive Vertical Line Array DIFAR (VLAD), and the DIrectional Command Activated Sonobuoy System (DICASS) [6] [7]. There are also other uses for vector sensors such as a receiver for underwater acoustic communication [8].
These systems that use vector sensors can do a wide variety of tasks such as can be detected as well as the localization accuracy without increasing the array aperture [9].
The vector sensor for the purpose of this thesis is assumed to be consisting of three sensors: two identical co-located but orthogonally oriented velocity hydrophone plus another pressure hydrophone. This is covered in greater detail in Chapter 2.

ESPRIT Algorithm overview 1.4.1 ESPRIT Algorithm and ESPRIT Algorithm for Vector Sensors
DOA estimation or source localization is important for a wide variety of reasons and there are two subspace methods that are typically used, MUltiple SIgnal Classification (MUSIC) [10] that was popularized by Schmidt [11], and the Estimation of Signal Parameters via Rotational Invariance Techniques (ESPRIT) Algorithm. The ESPRIT method, first described in [12] and later described in more detail in [13,14], can estimate source DOAs of narrowband sources that are highly correlated. The ESPRIT algorithm, an eigenstructure source localization algorithm, has been expanded over the years and adapted for a wide variety of different systems and purposes. One of the more recent purposes for vector sensors from Wong et. al. [15,1]

Results and Limitations
Results that compare and contrast three different pairing methods are explored in Chapter 3: classical, sorting, and exhaustive. The classical pairing method is the pairing method that is described in [1] and the sorting method is a new method presented in this thesis. Finally the exhaustive method runs through every possible permutation and selects the permutation matrix that presents the least amount of error.

Contributions of thesis
The contribution of this thesis is a new pairing method that relies upon sorting.
This method performs significantly better than the traditional pairing method.
Simulation results in Chapter 4 show that the sorting pairing method performs better than the traditional method and that the sorting pairing method performs almost as well as the exhaustive pairing method.

CHAPTER 2
Review of Existing ESPRIT Algorithm Developments and Eigenvalue Pairing for VSAs

Assumptions for the ESPRIT Algorithm
To simplify the derivation and the entire document, the following assumptions are made throughout the entire thesis.

Coordinate Systems
For this thesis, the 3-Dimensional Cartesian coordinate system consisting of x, y, and z axes will be used. This system uses θ to refer to the angle from the vertical z-axis, elevation angle, and φ to refer to the angle from the x-axis in the xy-plane, azimuth angle.

Sensor and environmental Model
• The entire signal model is in a free-space environment where sound waves travel in a homogeneous, isotropic fluid wholespace which implies direct-path propagation only. Constructive and deconstructive interference or multi-path propagation is not assumed. It is also assumed that the speed of the sound wave is uniform across all sensors.
• The impinging signal is a narrow-band acoustic source that operates on one frequency f where the wavelength λ = 1/f . In practice this means that the signal is sufficiently band-limited to allow narrowband processing in the frequency domain. Mathematically speaking this would mean that the bandwidth β is significantly less than the carrier frequency f c . Such a band-limited signal may be obtained by pre-filtering or computing the Discrete Fourier Transform of the time series measurements. This implies that the unit-vector from each sensor to a source or K sources is the same, regardless of the sensor location. In practice this would require sources to be in the far-field whose distances are much greater than the length of the array.
• The K source signals impinging upon the array of L vector sensors has the stipulation that L > K at possibly arbitrary and possibly unknown locations in a 3-D region with arbitrary geometry.

Vector Sensor Array
• In this thesis it is assumed that each vector sensor consists of three components: one pressure hydrophone and two velocity sensors. In each vector sensor, all three components are located at the same point; in other words all three sensors are co-located. In practice this requires that the component spacing is small compared to the minimum wavelength, which is set by the highest operating frequency.
• Each vector-sensor is modeled as a single point. In practice this requires the sensor dimensions to be small when compared with the minimum wavelength.
• The signal response of each velocity hydrophone is proportional to the cosine of the angle between the velocity axis and the source. Cosine velocity response results from measuring the velocity, and is performed only along one axis.
• The axes of the two velocity hydrophones are orthogonal. In practice this is true where each vector-sensor is a static unit.
• The vector sensor array consists of L vector sensors and while all locations of the vector sensors need not be known, all the velocity hydrophones, as previously stated, need to be identically oriented, i.e. the array of L vector sensors are located at arbitrary and possibly unknown locations in a 3-Dimensional region with completely arbitrary geometry, and all vector hydrophones are oriented in the same direction.

Derivation of noise-free VSA data 2.2.1 Signal Measurement Model
Each vector sensor's (possibly unknown) location can be defined as where p l denote the position of the l th vector sensor in 3-Dimensional space, and each position is defined for the purpose of this paper in units of wavelengths.
The three-component vector hydrophone would produce the following 3 × 1 manifold with regard to the k th source impinging from (θ k , φ k ): where θ denotes the kth source's elevation angle measured from the vertical z axis parametrized by θ k ∈ (0, 2π] and φ denotes the kth source's azimuth angle parametrized by φ k ∈ [−π/2, π/2]. Note that u(θ k , φ k ) and v (θ k , φ k ) or u and v respectively can be referred to as direction cosines. This nomenclature will be used later in the chapter and in Chapter 3.

Intervector hydrophone spatial phase factor
Each vector sensor is located at a particular 3-Dimensional (possibly unknown) location p l (1). Given the different sensor locations and the nature of the signal, each sensor will have a different phase. In order to simulate this in the signal model, an intervector hydrophone spatial phase factor is needed. The phase response of the l th vector sensor from a plane wave impinging from angle θ and φ is given by The phase responds for all L sensors may be combined into a vector This intervector hydrophone spatial phase factor as described in (3) can now be employed on a 3L × 1 array for the entire L-element vector hydrophone array in order to spatially relate the received signal on each individual sensor.
where ⊗ symbolizes the Kronecker-product, or tensor product, that is defined as: where A is an i × j matrix and B is a k × m matrix. This operation is possible because the phase vector a is simply the measurement vector for the corresponding pressure-sensor array. With the total of K ≤ L co-channel signals, the entire array would yield a 3L × 1 vector measurement z(t) at time t: where n(t) = 0, P k denotes the kth signal's power, σ k (t) represents a zero-mean unit-variance complex random process, λ refers to the signals' wavelength, c represents the propagation speed, and ψ k denotes the kth signal's uniformly distributed random carrier phase. The total of N snapshots (where N > K) taken at the distinct instants t n , n = 1, · · · , N , yields a 3L X N data matrix This matrix is a representation of the received acoustic energy from the vectorhydrophone sensor field. In it contains information that can be used to determine and estimate the DOA's through direction angles (θ k , φ k ) where k = 1, · · · , K.

ESPRIT Algorithm
This section will estimate the DOA's given the measurements from Z through the ESPRIT algorithm and makes no assumptions about array geometry, element characteristics, DOA's, or signal correlations.

Overview of ESPRIT Algorithm
DOA estimation is important for a wide variety of reasons and there are two methods that are typically used, the single subspace method, MUltiple SIgnal Classification (MUSIC) [1] that was popularized by Schmidt [2], and the

Estimation of Signal Parameters via Rotational Invariance Techniques
(ESPRIT) Algorithm [3]. The MUSIC algorithm yields asymptotically unbiased and efficient estimates. The MUSIC algorithm estimates the signal subspace from the array measurements and then estimates the parameters of interest from the intersections between the array manifold and the estimated signal subspace; however, in order to use the MUSIC algorithm, a priori knowledge of the array geometry and element characteristics is required [4]. Therefore it is impossible to do DOA estimation using MUSIC when the array geometry and element characteristics are unknown.
The ESPRIT approach first proposed in [3] is similar to MUSIC as it exploits the underlying data model and produces generally unbiased estimates. [3,4] in computationally efficient manner. The advantage of using ESPRIT over MUSIC is that • ESPRIT does not require a priori knowledge of the array geometry and element characteristics.
• It is computationally less complex since it does not need the search procedure that MUSIC uses.
These benefits explain why ESPRIT has had a significant amount of research over the past two and a half decades. The ESPRIT algorithm is not perfect as it assumes that the signal source is in the far field; in other words, the ESPRIT algorithm can only operate with plane waves.
Despite the above limitations, there are evident reasons why the ESPRIT algorithm would be favored in a wide variety of cases over other algorithms. These same limitations give rise to the special conditions under which the ESPRIT algorithm can offer a unique solution to the DOA estimation problem.

ESPRIT Algorithm Derivation
It is important to realize that the goal of the ESPRIT algorithm is to effectively remove the effects of q, shown in (3), the intervector phase factor between the sensors; thus, recovering the direction-cosines. By recovering the direction-cosines, the estimated DOA of the source can then be calculated.

Overview and Intuitive Analysis of the ESPRIT algorithm to estimate DOA's
Recall that each column the data matrix Z is a 3L × 1 vector. Each vector can be divided up into three L × 1 sub-vectors. These L × 1 sub-vectors are related to one another by invariant factors q that are dependent only on the positional locations p l of the vector hydrophones and not the direction cosines of the sources.
The 3L × 1 array manifold described as a from equation (5) can be rewritten as: Now that a is formed in a more explicit way one can see that the three sub-vectors are each related by the same q(θ k , φ k ) factor. Furthermore, one can form three L × K data blocks out of the matrix A, shown in (9), as follows where S j is an L × JL sub-array selection matrix where O m,n denotes an m × n zero matrix, and I m denotes an m × m matrix identity matrix. This matrix multiplication divides up A into 3 sub-blocks, which can be shown using (12) to be interrelated and described as:

ESPRIT derivation for Vector Sensors
An alternative derivation to the literature derivation is as follows. Recall that (15) where Φ u is a diagonal matrix containing the u direction cosines. Also, In the absence of noise, equation (7) shows that Z = AS, and the columnspace of Z is the same as the column space of A or where col(Z) denotes column-space for Z. If the singular value decomposition (SVD) is taken of Z, it can be written as follows where U 1 contains the left singular vectors of Z corresponding to the non-zero singular values. Recall that a property of the SVD is that col(Z) = col(U 1 ) and using (17) this gives Note that the relationship between the eigenvalue and singular value decomposition is that the left singular vectors of Z are the eigenvectors of ZZ H . Thus [U, E] = eig(ZZ H ) could be used in place of the SVD.

From (19) the following is true
where T is a non-singular change-of-basis matrix. Recall that given this fact, U 1 can be partitioned such as the following From (20), and through manipulation and substitution of the first and third equations into (15) one can get the resultant equation: Combining (24) and (25) yields U 11 = U 31 F u , which may be solved for F u to Using (25) and the fact that Φ u is diagonal, the eigenvalues of F u are the diagonal Similarly through substitution of the second and third equations from (23) into A 2 = A 3 Φ v (16) then one can get the resultant equation: where the expression T Φ v T −1 is given the name F v , Using a derivation similar to that for F u , it can be shown that Therefore, the eigenvalues of F v are the diagonal elements of Φ u .
From inspection (25) and (28) show that F u and F v have the same eigenvector matrix T . This is true when the u direction cosines and v direction cosines are listed in corresponding order. If the eigenvectors of F u and F v are computed separately, it is unlikely that they will be in corresponding order. In general, the eigenvectors and eigenvalues of F v will have to be permuted to get into corresponding order with the eigenvectors and eigenvalues of F u . Put differently, [T 1 , , where T 2 = T 1 P and P is a permutation matrix, and the diagonal elements of P Λ 1 P T are in corresponding order with the diagonal elements of Λ 2 . Note that this will be covered more in the next section, and in much greater detail in Chapter 3.

Pairing Problem
First recall known facts about permutations. Let π represent the permutation function that takes integers i = 1, . . . , N into some permuted order defined by π(i). Let e i be the ith standard basis vector in n . The permutation matrix corresponding to the permutation function π is Also recall some facts regarding permutations matrices • Permutation matrices are orthogonal such that • Given a matrix A=[a 1 · · · a n ], multiplication on the right by P permutes the columns of A as follows: • A permutation matrix and its inverse may be used to permute the diagonal elements of a diagonal matrix.
If P is defined by (30) then

Eigenvalue Pairing
To develop the concept of eigenvalue pairing consider two different matrices, F and G, that have the same eigenvectors but different eigenvalues.
where the columns of T 1 are eigenvectors of the matrices F and G, and the corresponding eigenvalues are the diagonal elements of Λ and D. Assume that the eigenvalues of F and G are real numbers and correspond to each other in the following order of pairs Suppose that the columns of T 1 are permuted into a new matrix T 2 and that the diagonal elements of D are permuted by a permutation matrix P into a new matrix D where Note that the diagonal elements of D are now a permutation of the diagonal elements of D. That is to say that d = d π(i) where π is the permutation function and defined in (33); thus, the columns of T 2 are eigenvectors of G and the corresponding eigenvalues are the diagonal elements of D. As a result of the preceding In the eigenvalue pairing problem the computed eigendecomposition of G is (38) and not (35). Therefore given the eigendecomposition G in (38) and the eigendecomposition of F in (35) the eigenvalue pairing problem is to find the permutation matrix P or its associated permutation function π.
From equation (37) the permutation matrix P and, therefore, the permutation function π can be obtained from the eigenvector matrices T 1 and T 2 The classical pairing method algorithm can be described by the following: given the matrices F and G with identical eigenvectors, suppose that the eigendecomposition of each matrix can be defined as where the permutation matrix that puts the columns of T 2 into the same order as the columns of T 1 is where the elements of P will be 0 or ±1, and where a −1 in the ijth element of P implies that the ith column of T 2 is the negative of the jth column of T 1 . Note that the negative signs cancel out in the computation of the diagonal matrix Now that D has been permuted by P , the diagonal elements of D correspond to the diagonal elements of Λ.
The implementation in MATLAB M-code for calculating equation (42)

Consider Noisy Data
Recall from equations (8) and (10) that the additive noise, n(t), to the data snapshot, z(t), was previously noiseless, n(t) = 0. Now suppose that n(t) = 0 and data matrix, which contains additive noise, is Z noise . When the DOAs are estimated using the ESPRIT algorithm from the steps described in (25) or (28) and the information provided about the eigenvalue permutation in section 2.4 are applied to Z noise then the estimated noisy DOAs are the eigenvalues of the following matrices:F where ∆ 1 and ∆ 2 represent perturbations due to noise. Note that one of the eigenvalue sets is ordered; the other has to be put into corresponding order. Let F = F u and G = F v . In this case the eigendecomposition will yield In this particular case when the matrix P is calculated it is no longer a permutation matrix. The matrix P will instead take on values other than 0 and 1, which are the only element values of a permutation matrix.
In order to form a permutation matrix from this P that has values other than 0 and 1 one has to modify the classical noise-free algorithm to handle noisy matrices.
The simplest approach, which is used by Wong et al. [5], is to compute the matrix P and then in each row put a 1 in place of the element with the largest magnitude and replace all other elements with zeros. In other words let (j i ) denote the row index of the matrix element with the largest absolute value in the ith column of P such that The following MATLAB code can be added to the previous function pairing methods % Create the permutation matrix such that it consists of all % zeros and only 1's at the maximum (j_k,k) locations for k=1:length(T1) [m,ind]=max(abs(P(:,k))); P(:,k)=zeros(length(T1),1); P(ind,k)=1; end This resulting matrix will likely be a permutation matrix. The only way that this process may fail is when a permutation matrix has the maximum elements from two different columns occur in the same row. This phenomena occurs when there are larger perturbations, or low SNR, in the matrices U 1 or U 2 . The terminology that will be used in this thesis to describe the above phenomena is a miss. The authors in [5] suggest that with a two source problem K = 2 a miss occurs with a near-zero probability if the computer represents numbers by more than a few bits and that if a miss should occur, however rare, that it shall be handled by performing the pairing arbitrarily, i.e. coin-toss.
For the case of noiseless or the case of Gaussian additive noise with an infinite number of snapshots the approximation of P becomes exact which is to say that the step where the maximum value in P is found is not needed. [3] A. Paulraj, R. Roy, and T. Kailath, "Estimation of signal parameters via rotational invariance techniques -ESPRIT," in Proc. Nineteenth Asilomar Conference on Circuits, Systems and Computers, Aslomar, CA, November 1985.

CHAPTER 3 Improved Eigenvalue Pairing Methods for DOA Estimation for VSAs
Recall from Chapter 2 the situation where we assume noisy data and form the resulting permutation matrices from the ESPRIT calculated results. Let P represent a 'pseudo-permutation' matrix where P is a matrix that has values other than 0 and 1, which are the only element values of a permutation matrix. In order to form a permutation matrix, P from the 'pseudo-permutation' matrix P has to be modified to take on values of 0 and 1. The easiest approach to replace elements within the matrix P to values of 0 or 1 is calculated by first computing the matrix P and then in each row put a 1 in place of the element with the largest magnitude and replace all other elements in that row with zeros. This resulting matrix will likely be a permutation matrix.
The only way that this process may fail is when a permutation matrix has the maximum elements from two different columns occur in the same row. This phenomenon, a miss, occurs when there are larger perturbations, low SNR in the data.
This chapter will present the method used by Wong et al.
[1] to randomly assign tied elements or what will be referred to as missed elements in this thesis, e.g., coin-toss, in the event of a miss. The following two sections present respectively a new eigenvalue pairing method based on sorting, and an exhaustive eigenvalue method, respectively.

Coin-toss Eigenvalue Pairing Method
A miss, or what is referred to as a miss in this work, according to the authors in [1] suggest that with a two source problem K = 2 a miss occurs with a near-zero probability if the computer represents numbers by more than a few bits and that if a miss should occur, however rare, that it shall be handled by performing the pairing arbitrarily, i.e. coin-toss. The methodology to implement a coin-toss can be best demonstrated with an example.
Let P be a matrix which is a K × K or 2 × 2 matrix P = 1 1 0 0 where there exists a miss in P in row 1 where there are two ones and in row 2 where there are two zeros. Once a row is identified to have a miss, an element is randomly selected with probability 1/T where T is the number of missed elements.
After the element is selected, it is moved to an empty row, a row that has all 0s and no 1s, and that has a probability of 1/T as well. This is done iteratively until there exists no missed rows. In order to ensure that the permutation matrix is well formed one can multiply P by P T which should by definition return the identity matrix I.
By the nature of the problem the maximum number of elements that can be in a missed row is T = K; therefore, if one extends K to a much larger number there exists a condition in which T becomes much larger as well. If this occurs, then the probability of the randomly distributed coin-toss matrix P matching the real or actual permutation matrix P real becomes smaller according to the following equation: where P r is the probability of the matrix P matching P real . Thus, if K increases there is a higher probability for more mismatches.

New Approach to Eigenvalue Pairing 3.2.1 Eigenvalue Pairing for Noise-Free Real-Valued Vectors via Sorting
Given the above issues with the coin-toss method and the likelihood of a mismatch of DOAs with the coin-toss method an alternative method is proposed.
Consider two vectors x and y in N whose elements are identical up to a permutation; where P is some permutation matrix with corresponding permutation function π.
For example, consider By looking at the corresponding vectors (51) and the definition (50) the permutation function π and the permutation matrix P for these two vectors are defined The permutation function for any two such vectors x and y may be calculated using their sorting index vectors. A sorting index vector s x corresponding to the given vector x shows how to arrange the elements of x in ascending order: x(s x (i)) ≤ x(s x (j)) for i < j.
Likewise the sorting index vector for y is denoted by s y . The permutation that relates the vectors x and y,y(i) = x(π(i)), may be written in terms of the sorting index vectors s x and s y π(s y (i)) = s x (i), i = 1, . . . , n.
Now consider the vectors given in (51). Their sorting vectors would be which, when s x and s y are put with (54), yields π(2) = 2, π(4) = 1, π(1) = 4, π(3) = 3, which gives the correct result, Equation (54) may be used to solve the eigenvalue pairing problem. Once again suppose that the matrices F and G have the following eigendecompositions where T 2 = T 1 P for some permutation matrix P with corresponding permutation function π (see (30)). The diagonal elements ofD must be permuted in order to have the elements ofD correspond to the diagonal elements of Λ. The permutation is as follows:d where the permutation function π corresponds to the permutation matrix P . Consider the following matrix M 1 where P represents the permutation matrix, and T 2 represents the associated eigenvectors from the eigendecomposition of the matrix G. The equation shows that M 1 is a diagonal matrix whose diagonal elements, m i , are a permuted version of the diagonal elements of Λ m i = λ π(i) , i = 1, . . . , n.
Let x and y be the diagonal elements of Λ and M 1 , respectively, and let s x and s y be the corresponding sorting index vectors. Once again (54) may be used to calculate the permutation function π that relates the eigenvalues of F and the eigenvalues of M 1 from s x and s y .
Consider another matrix M 2 where P represents the permutation matrix, and T 1 represents the associated eigenvectors from the eigendecomposition of the matrix F . The equation shows that M 2 is a diagonal matrix whose diagonal elements,m i , are a permuted version of the the diagonal elements ofD. Similarly (62) can be written as Let x and y be the diagonal elements ofD and M 2 , respectively, and let s x and s y be the corresponding sorting vector, and then (54) may be used to find the permutation π that relates the eigenvalues of M 2 and the eigenvalues of G.

Eigenvalue Pairing for Noisy Real-Valued Vectors via Sorting
Once again consider two vectors x and y whose elements are identical up to a permutation function π that is to say Now suppose though that only the perturbed versions of these vectors are known, where δx and δy are perturbations of x and y. That is to say that the two vectors x and y are perturbed by some amount δx or δy. Recall from (54) that the sorting vectors can be calculated from x and y. Now assume that the same methodology from (53) with the associated sorting matrix and apply the perturbation δx to x such that The sorting vector forx would be instead of the sorting vector for x This error or transposition in the sorting matrix occurs only when δx and δy are large enough to cause elements within the vector to swap places with one another.
That is to say that large perturbations in x and y will cause changes in the sorting index vectors, which will cause unavoidable errors in the estimated permutation function; thus, it is proposed to use (54) even in the case of noisy vectors.
Using the aforementioned realization, a novel sorting pairing method is explored. Once again consider matrices F and G whose eigenvectors are identical, and suppose that for both F and G only the perturbed matrices are known, F = F + ∆F andG = G + ∆G. Let the eigenvalue decomposition of each perturbed matrix be given byF where the diagonal elements ofΛ are x 1 , and the diagonal elements ofD are x 2 .
Similarly to (53) the sorting index vectors, s x1 and s x2 , are calculated from x 1 and Recall that a sorting index vector shows how to arrange elements of a given vector, x 1 or x 2 , in ascending order. Also recall that the columns ofT 1 andT 2 are the eigenvectors of F and G, and that the corresponding eigenvalues are the diagonal elements ofΛ andD respectively.
Calculate the matrices M 1 and M 2 that are formed by taking the eigenvectors T 1 from matrix F and replacing them for the eigenvectorsT 2 found in the matrix G and the eigenvectorsT 2 from matrix G and replacing them for the eigenvectors T 2 found in the matrix F .
Let y 1 be the diagonal elements of M 1 and y 2 be the diagonal elements of M 2 , and let s y1 and s y2 be sorting index vectors for y 1 and y 2 respectively.
Once the sorting index vectors are calculated the permutation functions, π 1 and π 2 can be calculated as follows: π 1 (s y1 (i)) =s x1 (i), i = 1, . . . , n π 2 (s y2 (i)) =s y2 (i), i = 1, . . . , n Now that the permutation functions, π 1 and π 2 , have been calculated, a test can be done to see which permutation function might be used. If the permutation functions are equal to one another, that is to say that then the permutation function π = π 1 . If π 1 is not equal to π 2 then one must calculate the error from each different permutation function in order to determine which permutation function to use, Once the error values α and β are calculated then the following permutation function assignment is made based on the following The advantage of the proposed permutation pairing method via sorting is that, unlike the classical approach, the proposed eigenvalue pairing method via sorting always returns a valid permutation function.

Exhaustive Pairing Method
The most exhaustive method is to permute through all possible pairing combinations. By utilizing a multitude of options and comparing the calculated error to one another, the best permutation option can be chosen. Consider a repository matrix that contains the order of eigenvalues for n different numbers. Then the following equation summarizes the number of distinct ways in which k eigenvalues can be picked in a distinct manner: where the equation assumes distinct numbers. In the case of the paper, it is assumed that n is finite; therefore, k is also finite. The variations of permutation listing, i.e. (1 2); (2 1), can be defined as π i where π is the exhaustive list of permutation matrices and i is the selection of a given permutation within the exhaustive permutation matrix.
The method described in this paper to determine the proper permutation matrix from an exhaustive matrix of permutation matrices can be found by calculating the error between one set of eigenvalues and another set of eigenvalues whose order is permuted through all possible combinations. That is to say where K is the number of eigenvalues, Error i is the error, F u and F v are labeled as x and y respectively, M 1 refers to equation (60) and M 2 refers to equation (62) and π i is the unique permutation as described above.

CHAPTER 4
Simulation Results

Sensor Locations
The simulation results presented in the figures below have the sensor setup that is described in

Number of Sources and Respective Parameters
In the results below there exist three separate scenarios consisting of 2, 3, or 4 sources, respectively, that impinge upon the sensor field. For each source there exists an angle pair, (θ, φ), and a direction cosine pair (u, v). Recall the relationship between the angle pair of direction cosine pair from Chapter 2: Note that the first two angle pairs and direction cosine pairs shown in Table 1 below are the same as those used in [1]. Following the structure laid out by Wong et al.
the additional pairs are designated accordingly such that The following in Table 1 are the parameters used in the results below. For all cases the signal power is set equal to unity, P = 1, and the sources impinge upon the previously outlined and described 13-element non-uniformly spaced 3-D array

Simulations
The simulation results in 5 through 12 illustrate the effectiveness of the pre- ing method gives away misses, its standard deviation and bias are not shown. The traditional method is considered useless in this case and the data is not presented.

Two Sources
The two source scenario is the same as the simulation presented in [1] as previously stated. The angle values can be found in Table 1 in the Two Source row. The results are as follows:   Table 1. Note that the SNR must be greater than 20 dB for the traditional method to give valid results.  Table 1.

Four Source
The Four source scenario includes the first two angle pair or direction cosine pair that the Wong et al. presents and two additional angle pairs which follows the described methodology for choosing direction cosines that are separated by 0.08 by Wong et al. The angle values can be found in Table 1 Table 1 ship, would be resolved and identified with high probability if both the estimation standard deviation and the bias are under approximately 0.02. The proposed sorting pairing method resolves these closely spaced sources for all SNR's at or above approximately 10 dB. Above the SNR resolution thresholds the estimation for both the standard deviation and bias both decrease for the both the exhaustive and the sorting pairing method fairly linearly with increasing SNR values. Since biases are typically one order of magnitude lower than the standard deviations the high biases are not the limiting factor for resolution of the estimated direction of arrival angles.  Figure 11. RMS standard deviation of {û 1 ,û 2 ,û 3 ,û 4v1 ,v 2 ,v 3 .v 4 } versus SNR: four closely spaced equal-powered uncorrelated narrow-band sources with angular values found in Table 1. There are three different pairing methods presented: classical, sorting, and exhaustive pairing method. The traditional method is not used below 30 dB (see Fig 10).

Overall Analysis
Given the fact that the traditional method was only able to avoid producing misses at very high SNR, it is quite evident that the sorting and exhaustive method performs better than the traditional method. Due to the fact that the lower limit where the algorithm can successfully resolve closely spaced sources depends on the composite rms standard deviation value and the composite rms bias being under 0.02 it is quite evident from inspection that the sorting pairing method presented as well as the exhaustive method performed as well with SNR = 30 or 40 dB, and better in all other SNR values.
In all cases the sorting and exhaustive pairing method were relatively colinear and would be able to resolve DOAs at approximately the same SNR. This is important to note as the exhaustive method is considered to be the best that can be achieved, which in turn means that in the presented scenarios, the sorting method would be the best pairing method out of the traditional and sorting pairing method. Finally, it is important to note that the exhaustive permutation method is computationally taxing as the number of possible permutations grows exponentially to the number of sources. Due to this fact, as the number of sources increase the number of computations must increase significantly more than the number of sources; thus, the sorting method, whose computational complexity increases approximately linearly to the number of sources, has a practical advantage.

Future Work
Future work could involve applying the same pairing method for the root-MUSIC algorithm for a similar case, but unlike the ESPRIT case where the vector sensor locations are unknown and the orientation is known, the root-MUSIC case has the sensor locations known and the orientation unknown. Despite the conclusion that the results would be likely similar in the root-MUSIC case as is found in the ESPRIT case, it would be of interest to research the future use of this sorting method in the root-MUSIC algorithm.
Another area that would be significantly interesting would be to formulate the CR-Bounds for this particular sensor setup. The lack of CR-bounds creates the ambiguity of how well the algorithm performed with the new pairing method.
It would be of interest to see in future work how the results in this paper compare to the CR-bound for this field of sensors and different presented scenarios.
It is quite evident that the current algorithm has its limitations where the direction cosines can be resolved only to difference of 0.08. Therefore, it would be of great interest to research how high-resolution, a resolution smaller than 0.08 for close direction cosines, ESPRIT algorithm would perform with this new pairing method.
There are other eigenvalue pairing problems that can be found in other fields including but not limited to chemistry, physics, and engineering, it would be interesting to see whether or not any of these methods could benefit from the sorting pairing method and whether or not it would yield more accurate results.
Finally it would be interesting to retrieve a computational baseline of the ESPRIT algorithm in terms of the number of flops or pure computations that is required, particularly when it comes to the pairing schemes when dealing with eigenvalues. It would be interesting to see how the traditional, sorting, and exhaustive method would compare to one another. As it was mentioned before in Chapters 3 and 4, the computational complexity increases for the exhaustive method as the number of sources increase as well as the miss rate for the traditional pairing method. It would seem as if the sorting pairing method is somewhere in between, achieving similar composite rms standard deviations and composite rms bias as the exhaustive method but being significantly less computationally taxing. This would be of interest to scientists in the field of oceanography and signal processing where there can be a large number of sources in a large field of buoys, and when the battery life and/or processing power is limited by either the size or required longevity of the buoy. end M = double(N); % Single will give us trouble on indexing. WV = 1:K; % Working vector. lim = K; % Sets the limit for working index. inc = 1; % Controls which element of WV BC = prod(M-K+1:M); % Pre-allocation of return arg. BC1 = BC / ( prod(1:K)); % Number of comb blocks. PN = zeros(round(BC),K,class(N)); L = prod(1:K) ; % To get the size of the blocks. cnt = 1+L; P = perms_loop(K); % Only need to use this once. PN(1:(1+L-1),:) = WV(P); % The first row.