Direction finding using ESPRIT with interpolated arrays

zyxwvuts zyxwvutsrqponmlkjihgfedcbaZY zyxwvutsrq zyxwvutsrqpo zyxwvuts zyxwvuts zyxwvutsr 1473 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 39, NO. 6, JUNE 1991 01 is a complex chi-squared random variable with K complex degrees of freedom. y is also a complex chi-squared random variable. Recall that the columns of U are independent. The columns of B are also independent and distributed identically to the columns of U since B = U H is obtained from U via a unitary transformation (see, e.g., [7, Appendix A, part F]). The columns of Z and q are independent and Gaussian distributed with covariance I; thus, [9, theorem 3.2.121 is applicable to y. This theorem says that y is a complex chi-squared random variable with M - K degrees of freedom and is independent of q. Furthermore, y and a are independent of s so p is independent of s. Now y/a is a complex F distributed random variable. A simple change of variables indicates that p is complex beta distributed P(P) = (M ( M - I)! pK-l(l - l ) ! ( K - l)! -K p)M-K-l . (27) The distribution of 2, given s is thus P ( 4 Is) = M ( M - I)! a2(M - K - l)! (K - l)! zyxwvutsrq The dependence on s is only through U = s Hs so, in principle, the distribution of is obtained by integrating the product of p(2, 1 s) and the distribution of u 2 over a 2 . However, the moments of 2, are, in general, of greater interest than the distribution itself. It is straightforward in this case to obtain the moments of 2, because p and u 2 are independent. For example, the mean of 2, is E{;,} = M - ’ E { a 2 }E { p } = .f- I. S. Reed, J. D. Mallet, and L. E. Brennan, “Rapid convergence rate in adaptive arrays,” IEEE Trans. Aerosp. Electron. Sysr., vol. AES-IO, pp. 853-863, NOV. 1974. J. Capon and N. R. Goodman, “Probability distributions for estimators of the frequency-wavenumber spectrum,” Proc. IEEE, vol. 58, pp. 1785-1786, Oct. 1970. R. A. Monzingo and T. W . Miller, Introduction to Adaptive Arrays. New York: Wiley, 1980, ch. 6. A. B. Baggeroer, “Confidence intervals for regression (MEM) spectral estimates,” IEEE Trans. hform. Theory, vol. IT-22, pp. 534545, Sept. 1976. E. J. Kelly, “An adaptive detection algorithm,” IEEE Trans. Aerosp. Electron. Sysr., vol. AES-22, pp. 115-127, Mar. 1986. E. J. Kelly and K. M. Forsythe, “Adaptive detection and parameter estimation for multidimensional signal model,” M.I.T. Lincoln Lab., Tech. Rep. 848, Apr. 1989. M. W. Ganz, R. L. Moses, and S. L. Wilson, “Convergence of the SMI and the diagonally loaded SMI algorithms with weak interference,” IEEE Trans. Antennas Propagat., vol. 38, pp. 394-399, Mar. 1990. R. J. Muirhead, Aspecrs of Multivariate Statistical Theory. New York: Wiley, 1982, ch. 3. B. D. Van Veen and R. A. Roberts, “Partially adaptive beamformer design via output power minimization,” IEEE Trans. Acoust., Speech, Signal Proc., vol. ASSP-35, pp. 1524-1532, Nov. 1987. K. M. Buckley, “Spatiallspectral filtering with linearly constrained minimum variance beamformers, ” IEEE Trans. Acoust., Speech. Signal Processing, vol. ASSP-35, pp. 249-266, Mar. 1987. L. J. Griffiths and C. W. Jim, “An alternative approach to linearly constrained adaptive beamforming,” IEEE Trans. Antennas Propagat., vol. AP-30, pp. 27-34, Jan. 1982. T. Kailath, Linear Systems. Englewood Cliffs, NJ: Prentice-Hall, 1980, Appendix. J. Capon, “High-resolution frequency-wavenumber spectrum analysis,” Proc. IEEE, vol. 57, pp. 1408-1418, Aug. 1969. A M’ where U: is the variance or power of the signal. Equation (29) indicates that the average MSE associated with the signal presence is directly proportional to the signal power and the number of adaptive degrees of freedom, but inversely proportional to the number of data vectors. Note that the presence of a strong signal results in large MSE. V. DISCUSSION The expected values of the output power and MSE due to the noise are within 3 dB of the optimum values after M = 2 K data vectors, while the expected value of the excess MSE due to the signal presence is down by 3 dB after M = 2 K data vectors. These results clearly indicate the benefits of reducing the number of adaptive degrees of freedom K the beamformer output is defined with fewer data vectors (M > K), and faster convergence to the optimum is obtained. The disadvantage of reducing K is an increase in the asymptotic MSE e and noise output power P,. This increase represents a loss in performance associated with reducing K and indicates that T, should be designed to minimize P , or equivalently e as suggested in [lo]. REFERENCES [I] N. R. Goodman, “Statistical analysis based on a certain multivariate Gaussian distribution,” Ann. Math. Stat., vol. 34, pp. 152-177, Mar. 1963. Direction Finding Using ESPRIT with Interpolated Arrays Anthony J . Weiss and Motti Gavish Abstract-The technique of interpolated arrays is applied to ESPRIT-type direction finding methods. The resulting method uses sensor arrays with an arbitrary configuration, thus eliminating the hasic restrictive requirement of ESPRIT for two (or more) identical arrays. Our approach allows for resolving D narrow-hand signals if the number of sensors is, at least, D 1, while the original ESPRIT method requires at least 2D sensors. Moreover, it is shown that while ESPRIT performs poorly for signals propagating in parallel (or close to parallel) with the array displacement vector, the advocated technique does not exhibit such weakness. Finally, using two subarrays, ESPRIT cannot resolve azimuth and elevation even when the sensors are not collinear. However, the interpolated ESPRIT procedure resolves azimuth and elevation using only a single array. All the above mentioned advantages are obtained with a reasonable increase of computation load, thus preserving the hasic and most outstanding advantage of ESPRIT. We also discuss, and illustrate numerically, the performance of the original ESPRIT when the sensor locations are perturbed. It is shown + Manuscript received May 25, 1990; revised December 16, 1990. A. J. Weiss is with the Department of Electronic Systems, Faculty of Engineering, Tel-Aviv University, Ramat-Aviv 69978, Israel. M. Gavish is with Communication Systems Division, ELTA Electronics Industries Ltd., Ashdod 77102, Israel, and with the Department of Electronic Systems, Faculty of Engineering, Tel-Aviv University, Ramat-Aviv 69978, Israel. IEEE Log Number 9143786. zyxwvutsrq 1053-587X/91/0600-1473$01 .OO 0 1991 IEEE 1474 zyxwv zyxwvutsrqpon zyxwvu IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 39, NO. 6, JUNE 1991 that the application of interpolation, in this case, can significantly improve the original ESPRIT performance. Our approach relies on the full knowledge of the array manifold like most basic array processing techniques and in contrast with the original ESPRIT. Therefore, the advocated technique should not be considered as a competitor of ESPRIT but as an exploitation of the ESPRIT ideas for standard direction finding. I. INTRODUCTION ESPRIT [5]-[lo] is a novel algorithm for the estimation of the direction of amval (DOA) of narrow-band signals, using an array of sensors. ESPRIT has lower computation and storage requirements than other suboptimal efficient eigenstructure techniques, such as the MUSIC algorithm [I I]. This advantage is achieved by using a displacement invariant array, i.e., sensors occurring in matched pairs with identical displacement vectors. Refinements of the original [6]-[8] formulation of the algorithm resulted in a class of ESPRIT-type methods, e.g., TLS-ESPRIT [5], matrix pencil methods [9], and PRO-ESPRIT [IO]. An additional advantage of these techniques is their robustness to array imperfections [5]. ESPRIT is not the only direction finding method that requires special sensor array configuration. For example, the root-MUSIC algorithm solves for the roots of a polynomial instead of performing the conventional MUSIC search and exhibits better performance at reasonable computational load [ 121. The root-MUSIC method can be applied only to uniform linear arrays. Spatial smoothing techniques, introduced to cope with perfectly (or highly) correlated signals, are applicable to a restricted class of array geometries [3]. The recently proposed approach of interpolated arrays [1]-[4] offers a generalization of the above methods to arbitrary array geometries. The basic idea is to estimate the outputs of a virtual array from the real array data, using an interpolation technique. The implementation of root-MUSIC and spatial smoothing with interpolated arrays have been investigated and found promising [ 11-[3]. In this correspondence we combine the idea of virtual interpolated arrays (VIA) with ESPRIT-type methods, obtaining a new class of DOA estimation algorithms, called VIA-ESPRIT. These algorithms preserve the computational advantage of ESPRIT as well as its robustness, and can be used with arbitrary sensor geometry. Furthermore, the use of VIA-ESPRIT eliminates the redundancy in the physical ESPRIT array. VIA-ESPRIT requires D 1 sensors to estimate the directions of D sources, instead of 2 0 sensors generally required by ESPRIT (excepting special cases such as linear uniform arrays). ESPRIT needs a large number of sensors for multiple dimensions estimation, e.g., azimuth and elevation. VIAESPRIT replaces the required additional sensors by virtual interpolated sensors, thus reducing considerably the hardware requirements of ESPRIT. VIA-ESPRIT also allows for the selection of the displacement vector regardless of the actual sensor configuration, thus reducing the performance sensitivity of ESPRIT to signal direction. It is also shown that, when the sensor locations are perturbed, the application of interpolation can significantly improve the original ESPRIT performance. However, VIA-ESPRIT relies on the full knowledge of the array manifold like most basic array processing techniques and in contrast with the original ESPRIT. Therefore, VIA-ESPRIT should not be considered as a replacement of ESPRIT but as an algorithm based on the ESPRIT ideas for standard direction finding. The correspondence is organized as follows. Section I1 discusses the implementation of interpolated array methods to ESPRIT algorithms, resulting in VIA-ESPRIT techniques. Simulation results are shown in Section 111. Finally, Section IV contains our conclusions. 11. VIA-ESPRIT A . Interpolated Arrays The interpolated array technique employed is explained in [ 11[3]. The first step in the design of an interpolated array consists of dividing the field of view of the array into L sectors. The size of the sectors depends on the array geometry and on the desired interpolation accuracy. Typical sector size is 30"-60" [ 11-[3]. Next, a set of angles 8, (the subscript 1 denotes the sector index) are selected for the design of the interpolation matrix for each sector. The optimal selection of the set is still an open problem. For the implementation of the interpolation in this correspondence, as well as in [I 1-[3], a typical number of 50-100 equidistant angles spanning each sector has been used. The interpolated array is defined by placing virtual elements at desired locations. The real array manifold Alo and the interpolated array manifold A,, are constructed from steering vectors associated with the set 8,. (The zero subscript is used to distinguish between the above manifolds and the direction matrix corresponding to the real signals.) Then an interpolation matrix BI is designed to satisfy in the least square sense, the equation zyxwvutsrqp The size of B, is M x M , where M is the number of sensors. The accuracy of the interpolation is examined by comparing the ratio of the Frobenius norms of (A,o - B,A,,) and If this ratio is not sufficiently small, the sector size is reduced and the interpolation coefficients are recalculated. The result of the procedure outlined above is a set of interpolation matrices { B , }, which are computed only once (off-line) for any given array. Given these interpolation matrices, data covariance, noise covariance, and signal subspace of the virtual array can be calculated from the corresponding matrices of the real array. B. Interpolated TLS-ESPRIT An ESPRIT array is composed of matched doublets, and can be grouped into two identical subarrays, X and Y , translationally separated by a constant displacement vector. The basic idea of VIAESPRIT is that any arbitrary array can be augmented to an ESPRIT geometry via interpolation. The original array plays the role of the ESPRIT X subarray, while the Y subarray is a virtual subarray, constructed by an interpolation technique as described above. Then any of the ESPRIT-type methods can be implemented, using only half of the elements. Let D narrow-band plane waves, centered at frequency wo, impinge on an arbitrary planar array composed of M > D elements, from directions { e , , 02, * , e,}. Using complex (analytic) signal representation, and grouping the signals received by the M array elements in the M X 1 vector x ( t ) , we can write zyxwvutsrqpon zyxwv + x(t) = AS(?) + n(t) (2) where ~ ( tis) a D x 1 vector whose kth element represents the kth signal. The matrix A is the M X D direction matrix whose columns (a(0,); k = 1, * * . , D} are the steering vectors of the D wavefronts. Finally, the M X 1 vector n(t) represents additive noise and its covariance, denoted by E, is assumed known. The signal and the noise are assumed to be stationary and ergodic complex-valued random processes with zero mean. In addition, the noise is assumed to be uncorrelated with the signals, and the signals are assumed not to be perfectly correlated. We shall present in the following the implementation of TLSESPRIT with virtual arrays. For each interpolation sector, a virtual array is constructed by translating the real array by a displacement zyxwvutsrqp zyxw zyxwvutsr zyxwvutsrqpon zyxwvutsrq zyxwvu 1475 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 39, NO. 6, JUNE 1991 vector. Note that the displacement vector (length and direction) is not necessarily the same for all sectors and it can be selected to maximize any desired performance criterion. The signal subspace corresponding to the virtual array is obtained by multiplying the real signal subspace Ex by the corresponding interpolatiofl matrix E!. The resulting VIA-ESPRIT algorithm can be summarized as follows: both azimuth and elevation estimates are required, another pair of subarrays (preferably orthogonal to the first pair) is necessary. However, the two angles can be calculated independently (parallel processing can be implemented) and the computational advantage becomes even more pronounced. ESPRIT computational load grows linearly with the dimension (neglecting the effort for the correct pairing of angle estimates), while that of MUSIC grows exponentially. The VIA-ESPRIT technique applied to multiple dimensions uses virtual arrays with different displacements to obtain the necessary sensor pairs. For example, one pair consists of the real array and an interpolated array translated along the x axis, while the second pair is abtained by translating the real array along the y axis. The interpolation sectors are designed in two dimensions, and the acceptable sector size generally decreases. Pairing of the angle estimates is based on [lo], namely the D correct pairs are those for which the steering vector projection onto the noise subspace is minimal. zyxwvu zyxwvutsrqponm Obtain an estimate of the covariance matrix R,, = E{x(t)xH(r)}, denoted Rn, using the measurements. Perform the generalized eigenvalue (GEV) decomposition of {Rxx,E}, i.e., find A and E that satisfy RUE = c En (3) . . . , Ad,is a diagonal matrix of the Xu, and E = [e,, . . . , e M ]is a matrix where A = diag { X I , GEV's, XI 2 . . . 2 consisting of the corresponding eigenvectors. Estimate the number of sources D using one of the available techniques for this purpose. Construct the estimated array signal subspace Ex given by Ex = [e,, . . . , ed]. (4) Compute the signal subspace EY of the interpolated array using Ex and the interpolation matrix E, Ey = (5) B/E,y. Compute the eigendecomposition 5 ... 2 i26, and Calculate the eigenvalues of (- EI2E;') and denote them by {&; k = 1 , . . . , D } . Estimate the DOA's using = Load The interpolated root-MUSIC algorithm [2] is an excellent candidate for direction finding using an arbitrary geometry array. In the following we compare its computational burden to that of VIAESPRIT. Let us look at one dimension (azimuth) estimation within one interpolation sector, and ignore the signal subspace computation, common to both algorithms. The major effort in the root-MUSIC implementation is solving for the roots of a polynomial of degree 2 ( M - I ) , which is equivalent to a 2(M - 1) x 2(M - 1 ) eigenvalue decomposition (EVD), requiring on the order of 120(M operations [ 1 3 ] . TLS-ESPRIT involves a 2 0 X 2 0 EVD, i.e., 120D3 operations as the most computationally expensive step. (Note that for simpler versions of ESPRIT [ 6 ] - [ 9 ] ,a D X D EVD is sufficient). Hence, for D << M, VIA-ESPRIT exhibits significant computational advantages as compared to the interpolated root-MUSIC method. zyxwvutsrq zyxwvutsrq where A = diag {XI, . . . , X2s},XI partition E into D X D submatrices 8, D. Computation arccos {arg ( & + ) / ( T A ) } ; k = 1, ... ,D (7) where A is the translational displacement of the interpolated subarray, expressed in half-wavelength units. The DOA's are defined relative to the direction of the displacement vector. Steps 5-8 should be repeated for each sector. DOA estimates from all sectors are collected, while discarding all out-of-sector estimates. Finally, D bearing estimates are selected. The selected DOA's correspond to eigenvalues &.. that are closest to the unit circle. This criterion is rather ad hoc, however, it was found to perform successfully in the simulations (see the next section). VIA-ESPRIT variations corresponding to other ESPRIT-type algorithms are derived in a similar manner and are not elaborated here. Since the interpolation is not perfect, the advocated algorithm is expected to exhibit errors that cannot be reduced by increasing the number of samples or the SNR. This estimation error can be controlled by selecting the appropriate sector size. The interpolation error increases not only with the sector size, but also with the displacement vector magnitude. This should be taken into account when the displacement is selected. C . Direction Finding in Three Dimensions An ESPRIT array consisting of doublets can be used to estimate DOA's in a plane (e.g., azimuth) containing the doublet axes. If E. Improving the Original ESPRIT in Nonideal Situations In practical implementations of the original ESPRIT method, it is expected that the two subarrays will not be perfectly matched. One of the reasons for mismatching is the difficulty associated with placing several doublets, all having identical displacement vectors. Assuming that the real sensor locations are known, the interpolated array approach can be useful in this case to calculate the data associated with sensors placed at the nominal (desired) locations. At least two approaches may be considered: 1) Use two virtual subarrays, the X subarray and the Y subarray, located at the desired locations; or 2 ) accept the actual X subarray as one array and define the virtual Y subarray by displacing the actual X subarray. In both cases the data of the virtual subarrays is obtained via interpolation of the corresponding actual subarray data. In the first approach one has to interpolate two subarrays while in the second approach only a single subarray has to be interpolated. Once the data of the two subarrays is available, the ESPRIT algorithm can be applied. In this work we used the second approach with excellent results as demonstrated in Section 111. 111. SIMULATION RESULTS Extensive Monte Carlo simulations have been run to evaluate the performance of VIA-ESPRIT. The examples presented in this section use the TLS versions of both ESPRIT and VIA-ESPRIT and the selected performance criterion is the total rms DOA estimation 1476 zyxwvutsrqponmlkjihgf IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 39, NO. 6, JUNE 1991 error over a number of 300 experiments. We define total rms error as the root of the arithmetic mean of the squared errors of all the DOA's. The signals are uncorrelated and corrupted by white Gaussian noise. The number of sources is assumed to be known. The sector size is 45" and 100 equidistant angles covering each sector are employed in the interpdlation matrix design. For example, the interpolation matrix for ther sector [45', 9 0 °] is calculated using the angles (45" + (90" - 45")k/99; k = 0, 1, * . . , 99). We use ESPRIT as a reference for the VIA-ESPRIT performance. However, one should keep in mind that, although employing a smaller number of sensors, VIA-ESPRIT exploits the full knowledge of the array manifold, while ESPRIT does not. The first simulated array (Figs. 1-4) is identical to that of IS]. Five elements are nonuniformly placed on the following x coordinates 0, 1, 3, 5.5, 7 (in half-wavelength units, X/2), defining the ESPRIT X subarray. The Y subarray is obtained by translating the above along the sensor line by A = X/4. As explained in the previous section, VIA-ESPRIT uses a virtual interpolated Y subarray, thus we compare in the following a five-element array to a tenelement true ESPRIT geometry. Fig. 1 shows the DOA rms error of both VIA-ESPRIT and ESPRIT, for two equal SNR sources located at 6 3 " and 7 0 ° , as a function of SNR. Each of the 300 Monte Carlo experiments used N = 100 snapshots for the cbvariance estimation. Fig. 2 illustrates the performance as a function of the snapshots number, at a fixed SNR of 15 dB; the scenario is similar to Pig. 1, except that the sources are now closer, namely, at 6 6 " and 70". The error of ESPRIT improves, approximately, with 1 while the error of VIA-ESPRIT is lower bounded by the interpolation error. However, in many practical situations a small interpolation error can be tolerated. Interesting and somewhat surprising results can be seen for SNR < 30 dB in Fig. 1, and for N < 450 in Fig. 2: VIAESPRIT outperforms ESPRIT, although only half the number of sensors are used. It should be emphasized, however, that such results are not general, but depend on the local interpolation quality and on the combination of DOA's within the sector. These phenomena may be explained by the exploitation of the array manifold information by VIA-ESPRIT and the fact that ESPRIT ignores it. To illustrate the performance of VIA-ESPRIT for various DOA values, the sources have been located at 80" & 6/2, and the source separation 6 has been changed. The results are plotted in Fig. 3. For a better comparison between the two methods which use different numbers of array elements, we present ESPRIT performance not only at the same SNR (15 dB) as VIA-ESPRIT, but also at SNR zyxwvutsrqponm zyxwvutsrqpon zyxwvutsr zyxwvutsrqponml zyxwvutsrqp I 5 10 20 I5 25 30 35 40 45 50 SNR [dB1 Fig. 1. The DOA rms error versus SNR. Linear nonuniform X subarray of 5 elements located on the x axis at 0, 1, 3 , 5.5, 7 ( X / 2 units); A = X/4 along the x axis; DOA = 63" and 70"; number of snapshots = 100. (a) VIA-ESPRIT, sector = [45", 90"J. (b) ESPRIT. zyxwvutsrqpon zyxwvut /a, 0.2 I I 0 200 100 300 zyxwv 400 I 1 I MH) 800 700 NUMBER OF SNAPSliOT9 Fig. 2. The DOA rms error versus snapshots number. Linear nonuniform Xsubarray of 5 elementslocatedon thexaxis at0, 1 , 3 , 5.5,7 (h/2units); A = X/4 along the x axis; DOA = 66" and 7 0 " ; SNR = 15 dB; (a) VIAESPRIT, sector = [45", 90"]. (b) ESPRIT. = 12dB. We investigate now the performance of VIA-ESPRIT when there are signals in different sectors. Two sources of equal SNR = 15 dB are located at 6 6 " and 70" in the [45", 90'1 sector, while a third signal with variable SNR = 15 dB + ASNR is located at 110" in the [90", 135"] sector. The DOA rms error corresponding to the three sources is plotted in Fig. 4, as a function of ASNR, for both ESPRIT and VIA-ESPRIT. No interference between the signals in different VIA-ESPRIT sectors is observed, even for significantly different power levels. ESPRIT-type techniques perform poorly for signals propagating in parallel (or close to parallel) with the array displacement vector. VIA-ESPRIT does not exhibit such weakness, since a preferred displacement direction can be selected for each interpolation sector. This is illustrated in the next example, which simulates the DOA estimation of one signal, using a circular X subarray with A4 = 8 elements. The spacing between sensors is equal to X/2. ESPRIT used a real Y subarray with doublet separation of A = X/4 (along the x axis). VIA-ESPRIT interpolated the Y subarray at the 4 0 8 10 12 14 18 LO 20 DOA SBPARATION [des] Fig. 3 . The DOA rms error versus DOA separation 6. Linear nonuniform Xsubarray of 5 elements located on the x axis at 0, 1, 3 , 5.5, 7 (X/2 units); A = X/4 along the x axis; DOA = 80" f 6/2; (a) VIA-ESPRIT, SNR = 15 dB, sector = [45", 90"]. (b) ESPRIT, SNR = 15 dB. (c) ESPRIT, SNR = 12dB. zyxwv zyxwvutsrqpo zyxwvutsrqponm zyxwvutsrqponmlkjih 1411 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 39, NO. 6, JUNE 1991 -P I I zyx 1 d zyxwvutsrqponm I 0 0 10 20 30 40 50 60 70 80 90 DOA [de&!] Fig. 4. The DOA rms error versus SNR separation ASNR. Linear nonuniform X subarray of 5 elements located on the x axis at 0, 1, 3, 5.5, 7 ( X / 2 units); A = X/4 along the x axis; DOA = 66", 70", and 110"; SNR = 15 dB, 15 dB, and 15 dB ASNR. (a) VIA-ESPRIT, sectors = [45", 90"], [90", 135"]. (b) ESPRIT. + same distance, however, for each sector the displacement vector was chosen to be perpendicular to the midsector direction. The results are shown in Fig. 5, which plots the DOA m s error as a function of the source direction relative to the displacement vector. As expected, ESPRIT exhibits considerable performance degradation at low DOA values, while the VIA-ESPRIT plot is a periodic function with period equal to the sector size. Assume now that an ESPRIT nonuniform linear array with M = 5 doublets is designed as in the first example, but sensor location errors exist. The main requirement for ESPRIT is a constant displacement vector. Thus, for simplicity, only the Y subarray has been perturbed, introducing independent zero-mean Gaussian location errors on both x and y sensor coordinates. The resulting ESPRIT performance degradation as a function of the perturbation standard deviation is depicted in Fig. 6. A virtual array has been placed at the nominal Y subarray location, and its output has been interpolated using the actual (perturbed) array data. The results for the corrected array using interpolation are shown in the figure as well. For the considered error range, the interpolation performed excellently and the effect of location errors was practically eliminated. The reason is that the perturbations are sufficiently small in order to enable almost perfect interpolation, while the original ESPRIT did not use any information about the sensor location errors. The last example of this section illustrates the implementation of VIA-ESPRIT for three-dimensional direction finding using an L-shape array with M = 5 elements. The (x, y) coordinates of the real sensors are (0, 0), (1, 0), (2, 0), (0, l), and (0, 2) in (X/2) units. The first virtual array has been constructed by translating the real array by (h/4) along the x axis; it is exploited to estimate the angles { a k ;k = 1, * * * , D}between the three-dimensional source directions and the x axis. The procedure is identical to azimuth estimation for a planar array and source scenario. The angles Pk relative to the y axis are estimated using a second interpolated array, obtained by translating the real array by the same amount along the y axis. Given a pair of angles (a,/3), the azimuth and elevation can be easily obtained. Angle pairing is performed as explained in the previous section, and a failure at this step results in bias of the estimates. Two sources with equal SNR = 15 dB have been placed at a, = 83", PI = 61" and a2 = 88", p2 = 67", respectively. The interpolation sector size is 10" for each one of the (CY,p) coordinates, e.g., [80", 90°] and [60", 70'1 sectors, respectively. The design of the interpolation matrix has been performed by defining Fig. 5. The DOA rms error versus source direction for one signal. Circular X subarray of 8 elements spaced X/2 apart; A = X/4 along the x axis for ESPRIT, perpendicular to the midsector direction for VIA-ESPRIT; (a) VIA-ESPRIT, SNR = 15 dB, sectors = [k45" - 22.5", k45" + 22.5"], k = 0, 1, . . . , 7 . (b) ESPRIT, SNR = 15 dB. (c) ESPRIT, SNR = 12 dB . zyxwvuts STANDARD DEVIATION OF SENSOR LOCATION [IamMa/Zl Fig. 6. The DOA rms error versus standard deviation U of sensor location. Linear nonuniform X subarray of 5 elements located on the x axis at 0, 1, 3, 5.5, 7 ( X / 2 units); A = X/4 along the x axis; DOA = 63" and 70"; SNR = 15 dB; the Y subarray elements have Gaussian location errors with standard deviation U on each of the x and y coordinates. (a) VIA-ESPRIT (interpolated Y subarray placed at the nominal location using the perturbed array data), sector = [45", 90"]. (b) ESPRIT. zyxwv zyxwvutsrqpo 15 equidistant angles on each coordinate, resulting in 225 interpolation angles. The VIA-ESPRIT estimates (mean and standard deviation), based on 300 experiments of 200 snapshots each, are (&, B , ) = (83.94' f 0.94", 61.61' D2) = (87.09" f 0.95", 66.45" f 0.91') f 0.94") and the total rms error is 1.21 '. Note that in this case the sectors are rather small, and the number of sectors is increased. However, this method requires only 1/3 of the original ESPRIT elements. IV. CONCLUSIONS The technique of interpolated arrays was applied to ESPRITtype direction finding methods. The resulting VIA-ESPRIT techniques eliminate the need for two (or more) identical subarrays, and can be applied to arbitrary geometry arrays. In addition, the displacement vector can be selected regardless of the actual sensor configuration and ESPRIT sensitivity to signal direction is reduced. A version of the proposed algorithm can be used for im- 1478 zyxwvutsrqponmlkjihgf IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 39, NO. 6, JUNE 1991 proving the original ESPRIT method when the sensors are perturbed. All the advantages are obtained with a reasonable increase of computation load. However, the array manifold must be known, in contrast with the original ESPRIT. ACKNOWLEDGMENT The authors wish to thank Dr. R. Roy and the anonymous referees for comments that improved this contribution. REFERENCES [I] B. Friedlander, “Direction finding with an interpolated array,’’ in Proc. IEEE Int. Con5 Acousr., Speech, Signal Processing, Apr. 1990. [2] B. Friedlander, “The interpolated root-MUSIC algorithm,” submitted for publication. [3] B. Friedlander and A. J. Weiss, “Direction finding for correlated signals using spatial smoothing with interpolated arrays,” submitted for publication. [4] T. P. Bronez, “Sector interpolation of nonuniform arrays for efficient high resolution bearing estimation,” in Proc. IEEE Inf. Con$ Acoust., Speech, Signal Processing, Apr. 1988, pp. 2885-2888. [5] R. Roy and T. Kailath, “ESPRIT-Estimation of signal parameters via rotational invariance techniques,” IEbE Trans. Acousr., Speech, Signal Processing, vol. 37, pp. 984-995, July 1989. [6] A. Paulraj, R. Roy, and T. Kailath, “Methods and means for signal reception and parameter estimation,” p n f o r d University, Stanford, CA, 1985, patent application. [7] R. Roy, A. Paulraj, and T . Kailath, “ESPRIT-a subspace rotation approach to estimation of parameters of cisoids in noise,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP-34, pp. 1340-1342, Oct. 1986. [SI A. Paulraj, R. Roy, and T. Kailath, “A subspace rotation approach to signal parameter estimation,” Proc. IEEE, pp. 1044-1045, July 1986. [9] Y. Hua and T. K. Sarkar, “Matrix pencil method and its applications,” in Proc. IEEE Int. Conj Acousr., Speech, Signal Processing, Apr. 1988, pp. 2476-2479. 1101 M. D. Zoltowski and D. Stavrinides, “Sensor array signal processing via a Procrustes rotations based eigenanalysis of the ESPRIT data pencil,” IEEE Trans. Acoust., Speech, Signal Processing, vol. ASSP37, pp. 832-861, June 1989. [ I 11 R. 0. Schmidt, “A signal subspace approach to multiple emitter location and spectral estimation,’’ Ph.D. dissertation, Stanford University, Stanford, CA, 1981. [12] B. D. Rao and K. V. S . Hari, “Performance analysis of root-MUSIC,” IEEE Trans. Acousr., Speech, Signal Processing, vol. ASSP37, pp. 1939-1949, Dec. 1989. [131 G. H. Golub and C. F. Van Loan, Matrix Computations. Baltimore, MD: Johns Hopkins University Press, 1984. solves an open problem of whether or not there exists a solution to the MEIR problem, proves the fast convergence of the MEIR algorithm proposed by Zhuang et al., and, finally, shows that the system of differential equations which is the basis of the MEIR algorithm is of the Lyapunov type. zyx zyxwvutsrqp I. INTRODUCTION This correspondence represents a continuation of the work reported in an earlier published paper by Zhuang ef al. [7].This correspondence discusses the following: 1) it justifies the proposed maximum entropy image reconstruction (MEIR) formulation; 2) proves the existence of the solution to the MEIR problem, an open problem (see [6]); 3) shows the fast convergence of the proposed MEIR algorithm; 4) demonstrates the MEIR algorithm is governed by a Lyapunov system, whose energy function measures the degree of constraint satisfaction. zyxwvutsrq zyxwvutsrqpo zyxwvutsrqpon zyxwvutsrqpon Maximum Entropy Image Reconstruction Xinhua Zhuang, Robert M. Haralick, and Yunxin Zhao This correspondeftFe is organized as follows. Section I1 covers point 1) above. Section I11 covers 2)-4). The final section presents further research directions and a conclusion. 11. JUSTIFICATION OF MEIR FORMULATION Let the required reconstructed image have positive pixel values fi, . . . ,f, which are to be determined, and on which the entropy is defined, where Let observed image data be given by (3) where e,’s are independent, zero mean, U; variance noise terms. We assume U,’ is known and define a constraint satisfaction function as follovs: \ 2 Typical least squares approaches would like to determine those valuesfi, * . . ,fn which minimize Q ( f i , . * ,&). Rather than this, we seek those f,, . . * , f, which maximize the entropy H ( p , , . . . , p,) subject to the constraint Abstract-This correspondence justifies the maximum entropy image reconstruction (MEIR) formulation proposed by Zhuang et al. (1987), Manuscript received September 18, 1989; revised June 6, 1990. This work was supported by the K.C. Wong Education Foundation, Hong Kong. X. Zhuang is with the Department of Electrical and Computer Engineering, University of Missouri, Columbia, MO 6521 1. R. M. Haralick is with the Department of Electrical Engineering, University of Washington, Seattle, WA 98195. Y. Zhao is with Speech Technologies Laboratory, Panasonic Technologies Inc., Santa Barbara, CA 93105. IEEE Log Number 9 14379 1. QUI,. . . ,A,) = f (5) which comes about from the central limit theorem [l], that is, with probability one zyxwvutsrqpon zyxwvutsrq Thus, provided m is large, we would expect the true value off,, . . . ,f, to satisfy (5). The condition (5) determines the set of feasible images each of which satisfies the given statistical test for consistency with the actual image data i d , , * . , d,,,}. In sum- 1053-587X/91/0600-1478$01 .OO 0 1991 IEEE