# Full text of "Comparison of superresolution algorithms with different array geometries for radio direction finding"

## See other formats

•LEY KNOX LIBRARY Al POSTGRADUATE SC! •JTFREY C/ VS101 NAVAL POSTGRADUATE SCHOOL Monterey, California THESIS COMPARISON OF SUPERRESOLUTION ALGORITHMS WITH DIFFERENT ARRAY GEOMETRIES FOR RADIO DIRECTION FINDING by Ku-Ting Lin September 1998 Thesis Advisor: Second Reader: David C. Jenn Phillip E. Pace Approved for public release; distribution is unlimited. 11 Approved for public release; distribution is unlimited. COMPARISON OF SUPERRESOLUTION ALGORITHMS WITH DIFFERENT ARRAY GEOMETRIES FOR RADIO DIRECTION FINDING Ku-Ting Lin Chung-Shan Institute of Science and Technology, Taiwan B.S.EE., Chung-Cheng Institute of Technology, Taiwan, 1990 Submitted in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE IN ELECTRICAL ENGINEERING from the NAVAL POSTGRADUATE SCHOOL September 1998 ABSTRACT The objective of this thesis is to investigate and evaluate the effectiveness of modern estimation methods with different array geometries as they apply to the problem of bearing estimation. The algorithms were selected from those that apply to multidimensional case, including MUSIC, PHD, minimum norm, and Capon' s beam- former. These four techniques were chosen based on their high resolution capability, and their ability to deal with three-dimensional non-uniform arrays and to estimate both azimuth and elevation angle of arrival (AOA) . Computer simulations were run for linear arrays, circular arrays, and combinations of the two. The test conditions included (1) two closely spaced emitters, and (2) various levels of additive white Gaussian noise. v VI TABLE OF CONTENTS I . INTRODUCTION 1 A. OVERVIEW 1 B. THESIS OUTLINE 6 1 1 . ANTENNA ARRAY 7 A. UNI FROM LINEAR ARRAY 7 B. NON-NUI FORM LINEAR ARRAY (MINIMUM-REDUNDANCY) .. 10 C . CIRCULAR ARRAY 11 D. SPIRAL ARRAY 14 E. COMBINATION OF LINEAR AND CIRCULAR ARRAYS 14 III . SPATIAL SPECTRUM ESTIMATION 17 A. MUSIC 17 B. ESPRIT 20 C . BEAMFORMING 22 1 . Conventional Beamf ormer 22 2 . Capon' s Beamf ormer 24 D. PHD 25 E . MINIMUM NORM METHOD 2 6 IV. IMPLEMENTATION AND COMPUTER CODES 29 A. PARAMETERIC DATA MODEL AND FLOWCHART 2 9 1. Input Parameters 29 VII 2. Signal and Noise Models 31 3 . Antenna Array Model 32 4. Correlation Matrix and Spatial Spectrum Methods 34 B. THE IMPLEMENTATION OF ARRAY MATRIX DECOMPOSITION 35 C. GRAPHICAL USER INTERFACE 37 V. SIMULATION RESULTS AND ANALYSIS 39 A. ANGULAR RESOLUTION of ONE-DIMENSIONAL DF 39 B. TWO-DIMENSIONAL DIRECTION FINDING 40 1 . Azimuth Angular Resolution 40 2 Elevation Angular Resolution 47 3. Multiple Signal Sources 47 4 . Comparison of Azimuth and Elevation Angular Resolution 52 5. Improvement of Angular Resolution of Two- dimensional DF 52 VI . CONCLUSIONS 61 APPENDIX A - THESIS MAIN PROGRAM 65 APPENDIX B - SPECTRUM CALCULATION FOR CIRCULAR ARRAY 75 APPENDIX C - SPECTRUM CALCULATION FOR SPIRAL ARRAY 79 APPENDIX D - SPECTRUM CALCULATION FOR CIRCULAR+LINEAR ARRAY 8 3 APPENDIX E - SPECTRUM CALCULATION FOR LINEAR ARRAY 87 Vlll APPENDIX F - SPECTRUM CALCULATION FOR NON-UNIFORM LINEAR ARRAY 8 9 LIST OF REFERENCES 91 INITIAL DISTRIBUTION LIST 93 IX :■: ACKNOWLEDGMENT I would like to acknowledge my thesis advisor, Professor David C. Jenn for his guidance and encouragement. Most of all, I would like to thank my wife, Seli and my daughter, Grace for their love, support and understanding . XI Xll I . INTRODUCTION A. OVERVIEW The angle of arrival (AOA) is a valuable parameter to be used in deinterleaving signals because a target can not change its position rapidly. Even an airbone radar can not significantly change its position in the few milliseconds of the waveform pulse repetition interval (PRI) time. As a result, the AOA measured by an intercept receiver on the radar is relatively stable value. It requires that a large number of antenna elements and receivers must be matched, either in amplitude or in phase in order to measure the AOA. Direction finding (DF) is the area of electronic- warfare support whose objective is to find the AOA [Ref . 1] . DF systems passively examine the spectrum use by hostile emitters and process the signals to obtain enemy AOA information. The two common approaches to measure AOA are based on amplitude and phase comparison [Ref. 1] . The phase comparison system usually can generate AOA resolution of ±1 degree, which satisfies the modern EW requirement. The number of antennas that can be used in an AOA system is very low in comparison with electronically scanning antennas in a radar system. In an airbone system, the maximum number of elements might be ten. In a ship-based system, the number might be larger. For the phase comparison method at least two separate antenna elements with a predetermined space between them are required. Since the distances from the emitter to the two elements are not the same (except in the broadside case) , the incident wave arrives at the two elements after traveling uneven paths, and thus it arrives with a different phase. The phase delay is proportional to the antenna spacing, wavelength and angle of incidence. Since the first two factors are constant, the direction information can be obtained. DF systems often receive signals from several sources in the same frequency range at the same time and from many different angles of arrival. To resolve multicomponent wavefields, modern techniques of spectral estimation are used. Most recent digital AOA studies have concentrated on high-resolution approaches, such as the MUSIC, ESPRIT and minimum norm methods [Refs. 2-4] . In this thesis the performance of several modern estimation methods is evaluated for different array geometries. This research deals with both the one- and two-dimensional problems: azimuth angle of arrival for linear arrays, and azimuth and elevation angles for a planar array. In addition, an underlying assumption is that the receiver is far enough from the transmitter and the antenna interactions with the platform are controlled so that the arriving wavefront is essentially planar. The main function of a direction finder is to determine the direction of arrival of an incident electromagnetic wave relative to the coordinates of the direction finding site [Ref . 1] . Figure 1 shows a representative direction finding spatial coordinate system with the DF antenna located at the origin. The angle of arrival is specified by the azimuth angle measured from the x-axis, and the elevation angle measured from the x-y plane. Under ideal circumstances the electromagnetic field incident on the DF antenna exhibits a locally plane wave structure with linear polarization. The direction of propagation is indicated by Poynting's vector, P P = -Re{E i xH i *} (l) where * denotes complex conjugation, Ej is the incident electric field intensity, Hj is the incident magnetic field intensity, H,-^. (2, where k is unit vector in the direction of propagation and Z is impedance of free space. In practice, the incident electromagnetic field is usually nonplanar with phase-front distortion created by multipath and scattering, and depolarization produced by nonuniform ionospheric propagation effects if such a path exists. A full-capacity generic DF should measure angle of L l/\ Hi e k Ej +yf k, w y Figure 1. Direction- finding spatial coordinate system with the array located at the origin arrival in three-dimensional space. However, many operational situations require only measurement of the direction component in the azimuth plane. A plane wave, incident at an angle other than the normal to the baseline between the two elements, arrives at one element earlier than the other. Figure 2 illustrates the basic phase delay concept. An incident plane wave arrives at an angle G at antenna 1 inducing a voltage which can be expressed in complex notation as V x exp(y'co/) . After propagrating the distance t/sin9 , the incident plane wave induces a voltage V 2 exp(y(or - <t>) in antenna 2, where the phase delay given by O = (2dn/A) sin0 . (3) Therefore the bearing angle is encoded as a function of phase delay. For the phase difference technique, the bearing angle is computed by using Eq. (3) , where phase difference Figure 2. Phase-difference DF Parameters is measured and the baseline distance d and wavelength A. are known . In the case of linear or planar arrays, a phase difference can be measured or computed for each combination of elements in the array. Then an array correlation matrix can be defined and used in one of several spectral estimation methods to estimate the AOA. B. THESIS OUTLINE Chapter II summarizes several common antenna array- geometries for one-dimensional DF, including uniform and non-uniform linear arrays. The combination of circular and linear arrays for 2-D DF is also investigated. Chapter III covers the .mathematical development of the important estimation methods, including multiple signal classification (MUSIC) , the estimation of signal parameter via rotational invariance technique (ESPRIT) , minimum norm, conventional and Capon's beamforming, and Pisarenko harmonic decomposition (PHD) . Also, DF for multiple emitter AOAs is introduced by using the array covariance matrix. Chapter IV describes the development of a computer code to simulate direction finding antenna performance for signals with white noise present . Matlab is used to implement the estimation methods and perform the array matrix calculations. A graphical user interface (GUI) was developed to allow a variety of parameters to be changed with a minimum of effort on the part of the user. Chapter V shows the simulation results, and validates the code by comparing the data with previously published results. Finally, Chapter VI presents a summary, conclusion and recommendations for future work. II. ANTENNA ARRAYS Antennas can be viewed as spatial filters; they enhance signals in desired directions while simultaneously suppressing signals in other directions. In a phased array, the selectivity is achieved by phase shifting and then superimposing the outputs of all elements [Ref . 1] . For a passive array, the improvement in the signal-to-noise ratio (SNR) is potentially equal to the number of elements. In an array of identical elements, there are two related factors that must be considered. The first is the geometrical arrangement of the array (linear, circular, etc.) . Second is the relative displacement between the elements, once the element arrangement is set . A. UNIFORM LINEAR ARRAY In most EW applications, it is desirable to use as few elements as possible. For example, DF arrays on a ship are linear arrays placed in the horizontal plane to measure the azimuth angle. One linear array can theoretically cover up to 180 degrees of azimuth angle, although it is often limited to 120 degrees to avoid operating in the end-fire mode [Ref. 1] . If the elevation angle is also of interest, a linear array in the vertical direction has to be added. As far as the AOA measurement is concerned, the two linear arrays are usually treated separately. Figure 3 shows a linear array with N elements indexed from q-0,1,..., N-l along the x direction. Let us assume that there is only one plane wave of frequency / and that a wavefront arrives as shown. If the input signal is a sinusoidal wave, then the q th element will have the output g(q,t) = Acos[2nf(t-t )] ' (4) where A is the amplitude of the incident signal, and t ( is the time delay at the q th element with respect to the reference antenna element (q=0) . Since the distance between the antenna elements is d , the delay time can be written as -qd sm§ r q = (5) c where 9 is the incident angle of the input signal, and c is the speed of light. (N-2)d sinG V V cr d — ►«- 6 T £ X N-l Figure 3. A linear array with a plane wave incident The minus sign in Eq. (5) is due to the fact that the equiphase plane arrives at the q th element before it reaches the reference element. Substituting Eq. (5) into Eq. (4) gives qdsinQ g(q,t) = Acos[2M* + )] • ( 6 ) c Let k be a unit vector pointing in the direction of the wave propagation and x a unit vector along the direction of array axis. Eq. (4) can be written as g(q, t) = A cos(2/bi - 2nfqk ■ x— ) ( 7 ) c r „ tz where k • x = cos( — i- 0) = -sin0 . 2 If there are M signals, the output from the gth antenna element is a superposition of all the received signals g(q,t) = Z An cos(2nf m t - Mf^s±) ( 8 , m=l C where /„, and k m are the frequency and direction of the m th signal . Often exponential form is used to express this result as g(q,t) =ZA m exp(j(2nf m t - 2 ^qdk m -x )} ( g } m=l C In order to obtain digital data, the output is sampled at discrete times and therefore the time t will be replaced by integers n =0 , I , . . . , b - 1 , where b is the number of time sample points. If we assume the frequencies of all signals are the same, the factor e 1 ' m ' can be suppressed, and a standard phasor quantity results. For an AOA measurement, the elements cannot be spaced too far apart or an ambiguity will result [Ref . 1] . The shortest distance between two antennas must be less than half of the wavelength of the highest frequency in order to fulfill the Nyquist theorem [Ref. 1] . Therefore, one must sample the incident wavef ront twice per cycle . B. NON-UNIFORM LINEAR ARRAY (MINIMUM -REDUNDANCY) By making use of a theorem by Caratheodory [Ref. 5] it is shown that, 'for a given number of elements N , there exits a distribution of element positions which results in superior spatial spectrum estimators than the uniform linear array. In addition to less hardware, certain nonuniform arrays are capable of large dynamic range (spectral peak to background level), lower sidelobes, and the relatively small estimate bias values. A special type of nonuniform array is the so-called minimum redundancy array (MRA) , in which integer spacings of the base spacing occur only once. A MRA is capable of resolving N(N-l)/2 multiple (simultaneous) signals, whereas a conventional array can only resolve N-l . Some minimum- redundancy array configurations are shown in Table 1. A dot in column 3 represents an element and the 10 integers the number of base spacings between them. (The base spacing is typically A./2 . ) A disadvantage of minimum-redundancy array is that its resolution is not easily increased except by adding more antennas and rearranging the array to the optimum configuration for the new number of antennas. If there are M signals, the output from the q th antenna element is M Mm d a k m ' X }) (10) g(q,t)= I A m exp0(2nf m t ' m q m -)) where d q is the distance between first antenna and the qth antenna . Table 1 . Some minimum redundancy array configurations . Number of elements. jV Array length Configuration 3 3 .1.2. 4 6 .1.3.2. 5 9 .1.3.3.2. 6 13 .1.5.3.2.2. 7 17 .1.3.6.2.3.2. 8 23 .1.3.6.6.2.3.2. 9 29 .1.3.6.6.6.2.3.2. 10 36 .1.2.3.7.7.7.4.4.1 C. CIRCULAR ARRAY A planar array can be used to measure both the azimuth and elevation angles. These types of arrays occupy more space than the linear array but they might be necessary for some special aircraft and shipboard application. Planar arrays include circular and rectangular (of which crossed linear arrays can be considered a subset) . A circular 11 antenna array is shown in Figure 4 . The radius of the circle is R and there are N elements from q =0 to N-l. The array is in the x-y plane and the first antenna element is the reference element iq=0) and it is on the x-axis. The location of the elements are at integer multiples of the angle § = 2n / N . Assume that a plane wave with amplitude A } and frequency / ; is incident on this array with azimuth angle <)), and 6; (the angle between the incident ray and the normal of the array) . The output induced at gth element can be written as gfqJJ^AjexpOVnfjft-^))) c 11) Figure 4 . Circular antenna array where the space delay x , referenced to the center of the array is 12 x q i =-# sine, cos (f)^ . (12) In this case, because the antenna locations are equally spaced in azimuth, the angle <|) 7 can be written as b q i=1$o-bi- (13) Therefore R sin Qj cos( — - - <j> 7 ) g(q,t) = A, exrtj2nf,(t + # )) 2.7ZQ R sin 6 ; cos(— — - <j> y ) = A, exp(j2n(f,t + —^ )) . (14) A.; The incident direction of the input signal can be represented by a unit vector k } in Cartesian coordinates as kj = -(sinOj cos^/X + sinO/ s\n§ l y + cosQ ] z) (15) where x, y and £ are unit vectors in the Cartesian system. The minus sign in front means the vector is pointing to the origin. The position vector of the q th antenna element is r where r q =cos(q$ )x + sm(q$ )y (16) so that kj -r q =-smQ 1 cos(q^ -^ l ) . (17) Extending this to cover multiple signals ( m = 1,2,..., M ) gives 13 m R k -r g(q.t) = 14, ex P 02n(f m t - ° m q )) . ( i 8 m = l K r ■ m The linear array can only provide azimuth angle information because the corresponding propagation vector is two-dimensional. The circular array can provide both azimuth and elevation information since the corresponding propagation vector is three-dimensional. D. SPIRAL ARRAY We examine a variation of the circular array, Eq. (18) , where the radial distances of the elements increase in steps of R /(N-1). Therefore R = R q / N is the distance between q th antenna element and the center point. The element locations lie on an Archemedian spiral. In a circular array, all antenna elements have the same R q . The output of qth antenna element is as follows: m R k ■ r„ g(q,t)= ZA m exp02n(fj- q m q )) . (19) E. COMBINATION OF LINEAR AND CIRCULAR ARRAYS A three dimensional array can be constructed by combining planar and linear arrays. Figure 5 shows a circular array in the x-y plane and a linear array along the z axis. We examine the special case where the elements of the circular array have § = 90 , which results in two elements on the x and y axes equally spaced d from the 14 center. Let there be 3 the elements along the z axis, also spaced d . Therefore, the array has seven elements, all equally spaced. y Figure 5. There dimensional seven element array The positions of the 3 elements of linear array are given by the position vector r h =(h-l)dz where h = 1,2,3, so that k m -r h =-cos(Q m (h-l)d). (20) If the array receives M signals, the output of h th vertical antenna element of the linear array is g linear (h,t) = % A m exp0(2nfj ^^)) m=l 21) 15 Another four antenna elements comprise the planar array which has an output of the form m dk m ■ r Z circular (<l> = Z A„ QXp(j2n(f m t - — H -)) (22) m=l k m where q=l,...,4, M is number of signals, and k m -r is the same as Eq. (18) . This array combines linear and circular arrays to provide 2-D direction finding. From Eq. (21), it is evident that the vertical array only provides elevation angle information. We have described the various array geometries. The elements comprise the data acquisition system which samples the incident signal wavefronts. The next step is to carry out the estimation process. In the next chapter, we introduce spectral -based algorithmic solutions to estimate the AOA. 16 III. SPATIAL SPECTRUM ESTIMATION Many different high-resolution approaches can be used to estimate AOAs from digitized input data. In this chapter, five high resolution methods will be discussed: (1) Multiple Signal Classification (MUSIC) [Ref. 2.]. (2) Estimation of Signal Parameter via Rotational invariance (ESPRIT) [Ref. 3] . (3) Capon's beamformer (Maximum Likelihood Method) [Ref. 6] . (4) Pisarenko Harmonic Decomposition (PHD) [Ref. 7] . (5) Minimum norm [Ref.l]. All of these are cabable of (1) high resolution, (2) estimating both azimuth and elevation AOA, and (3) handling three-dimensional non-uniformly spaced array signals. A. MUSIC MUSIC is a technique used to determine the parameters of multiple wavefronts arriving at an antenna array from measurements made on the signals received at the individual elements [Ref. 2 ] . The waveforms received at the N array elements are linear combinations of the incident wavefronts from M narrowband signals and noise. Thus, the multiple signal classification approach begins with the following model for characterizing the receive vector G as 17 G=AF+V (23) where V a is noise vector. Expanding Eq. (23 S: a,{Q,) a,(Q 2 ) a 2 (Qj) a 2 (Q 2 ) >/" >/" <*i®mY F 2 v 2 <*2<Pm) + a N^M)_ -M_ 7n_ 24 SN _ The size of each matrix is as follows: G : vector of element outputs, N by 1 A : matrix of mode vectors, N by M F : vector of signals, M by 1 V : noise vector, A^ by 1 The incident signals are represented in amplitude and phase at some arbitrary reference point, for instance, the origin of the coordinate system, by the complex vector F . The elements of G and A are also complex in general. The a n (§ m ) are functions of the signal arrival angles and the array element locations. That is, a n (Q m ) depends on the nth array element, its position relative to the origin of the coordinate system, and its response to a signal incident from the direction of the mth source. The m th column of A can be considered as a x mode' vector of responses to the direction of arrival Q m of the m th signal. This N hy 1 mode vector will be denoted by a(Q) . 18 The M by M covariance matrix of the G vector is S = E[GG*] = GG* =AFF*A* + VV' 25 where E is expectation operator and the asterisk denotes complex conjugation. Eq. (25) can be rewritten as A S = APA' + AS () , where P = FF' is a diagonal M by M matrix and S is the noise covariance matrix. When the number of incident wavefronts M is less than the number of array elements N , then APA' is singular. It has a rank less than N ; therefore APA = LS-/LS 1 = This equation is only satisfied with X (not to be confused with wavelength) equal to one of the eigenvalues of S in the matrix of S . Therefore A can only be the minimum eigenvalue A. mm . Observe that any vector orthogonal to A is an eigenvector of S with the eigenvalue a' , the noise variance. Hence, we can write S = UAU H =U S A S U? +U N A N U% 26 where A = diag{ Xj,/. 2 ,^-m) ^ s a diagonal matrix of real eigenvalues. The superscript H represents the Hermitian, U s is the signal eigenvectors and U N is the noise eigenvectors. The noise eigenvectors can be used to form an estimator for the spatial spectrum 19 PMu( * )= am»U N UJ!aQ) 27) Although MUSIC is highly robust it requires characterization of the array response and a spatial search through all possible angles of arrival, G. B. ESPRIT ESPRIT [Ref. 3] eliminates the undesirable features of MUSIC, but only under certain constraints. Assume that there are N points of data (N antennas), which can be divided into two groups, each consisting of N-l points. The first group has an output signal x(t) = As(t) + n ; (t) and the second group y{t) = ABs(t) + n 2 (t) , where A is a steering vector (not related to A in Eq. (24)), s(t) denotes the baseband waveforms, and B is a diagonal matrix B = v l v 2 . v M 28 2tz where v ; = exp(-y — <isin0 ( ) and d is antenna spacing. These data A. will be used to make two matrices by the covariance approach. 20 Let R xx be the auto-covariance matrix of x(t) and R xy the cross-covariance matrix of x(J) and y(t) . If there are M signals then R xx =E[x(t)x*(t)]=ARA* +a 2 I (29) R xy =E[x(t)y(t)]=ARB*A*+cr 2 J ] (30) where ^=[s(7)s (/)] is the signal covariance matrix under the assumption of uncorrelated noise, / is the identity matrix and Jj is an N by A^ matrix with ones along the first lower diagonal off the major diagonal and zeros elsewhere. For a pairwise matched array with doublets the cross covariance matrix becomes R xy =ARB'A\ (31) In the absence of coherent signals, eliminating cr~ , we obtain Cxx=Rxx-c 2 / = ARA* C X y=Rxy- cr 2 J } = ARB'A' which gives C X x- 7C xy = AR{\-yB')A'. (32) The singular values of Eq. (32) are given by the roots of \l-yB'\ = 0. The desired singular values are 21 Vk : = v k , k=l,2, . . ,M. To find the angles we note that for a uniform linear array, the angles of the eigenvalues are equal to 2^7rsinG so that 8* =sirT / (^-) (33) where S k is obtained from y k = exp(y9 A .) . Thus, the direction of arrival can be obtained without a search technique. In this respect computation and storage costs are reduced considerably. C . BEAMFORMING 1 . Conventional Beamf ormer The first attempt to automatically localize signal sources using antenna arrays was through beamforming techniques. The idea is to steer the array beam in one direction at a time and measure the output power. The steering locations which result in maximum power yield the DOA estimates. The array response is steered by forming a linear combination of the antenna outputs y(t) = £w* gi (t) = W H G(t) (34) 1=1 where g t {t) is the output of the /th element, and w,-(0 the corresponding complex weight factor. The output power for b time samples is 22 P(W) = -Y\y{t)\ =-tw H G{t)G H {t)W = W H RW (35) b t=o b 1=0 where R is the signal covariance matrix defined previously. The conventional beamformer is a natural extension of classical Fourier-based spectral analysis to antenna array processing. Suppose we wish to maximize the output power from a certain direction 9 , given a signal s(t) emanating from direction . A measurement at the element output is corrupted by additive noise n(t) and written as g.(t) = a ! (9)s l (t) + n j (t) . The problem of maximizing the output power is formulated as, max[E{w H G(t)G H (t)w} ] = max[W H E^}(t)G H ( k t)\v ] = maxlE\s(t)\ 2 \w H a(Q)\ 2 + a 2 \W\ 2 \ where the assumption of spatially white noise is used. The resulting solution is 1/(6)0(6) The above weight vector can be interpreted as a spatial filter, which has been matched to the impinging signal. Using the weighting vector, the classical spatial spectrum estimator is P BF (Q) = a H (Q)Ra(Q) (38) where R is the signal covariance matrix. 23 For a uniform linear array of isotropic sensors, the steering vector is a(Q) = 1 e** e 20 . . . e>^ 1J *f 39) where <t> = -kdsinQ , N is number of elements and T denotes transpose. The standard beamwidth [Ref. 6] for a uniform 2tz linear array is — , so sources whose angles are closer than N this beamwidth will not be resolved by the conventional beamformer, regardless of the available data quality. 2 . Capon' s Beamformer In an attempt to alleviate the limitations of the conventional beamformer, such as its resolving power of two sources spaced closer than a beamwidth, an improved method was proposed by Capon [Ref. 7] . The optimization problem was posed as min{P(W)} , subject to W H a(Q) = 1 . Hence, Capon's beamformer (also known as the minimum variance distorsionless response filter) attempts to minimize the power contributed by noise and any signals coming from directions other than 9, while maintaining a fixed gain in the look direction G . The optimal W can be found using methods such as the technique of Lagrange multipliers [Ref. 6] , resulting in W CAP - /""(?> . (40) 24 Using the above weight vector, we can find the following spatial spectrum p cap(Q) = — r f ; • (41) The power minimization can also be interpreted as sacrificing some noise suppression capability for more focused rejection in the directions where there are other sources present. The spectral leakage from closely spaced sources is therefore reduced, though the resolution capability of the Capon beamformer is still dependent upon the array aperture and on the signal-to-noise ratio though R . D. PHD The PHD is based on the use of a noise subspace eigen- vector to estimate the incident angle and assumes that the process consists of M complex signals in additive complex white noise. It derives the angle of the signals, their powers and the white noise variance from the known autocorrelation sequence. This algorithm accounts for the conjugate symmetry of the eigenvector. This method is based on an estimation of angles by using the orthogonality property between the signal vectors and vectors in the noise subspace. The minimum eigenvector is orthogonal to the signal vectors, thus we have ai H (Q)u N =0, i = l,2,...,M . (42) 25 Once the minimum eigenvector u N is found the from U N in Eq. (26) , the angles of the signals may be determined by factoring the eigenfilter polynomial [Ref. 7]. The roots are guaranteed to be on the unit circle. The a(Q) is a steering vector that will be zero whenever G = G, , one of the angles of the signals. This leads to the angle estimator function 1 PpHD^)~ II a"(Q)u f 7 ' 43) MINIMUM NORM METHOD This method and the MUSIC method are somewhat alike. The basic idea is to find a vector D that is a linear combination of eigenvectors in the noise subspace. This vector D can be written as Z) = [8 5; . . . 8^] T = [l 5 ; . . . d N f (44) where the superscript T signifies the transpose and TV the number of antenna elements. In this equation d is assigned to be unity. The square of the norm of Eq. (44) is 8 = Z 8? i=0 where M is the number of sources. The minimum norm method minimizes this quantity and hence the name. The procedure is as follows: 26 1. Find the eigenvectors U of the covariance matrix R. Large eigenvalues correspond to signals and small eigenvalues correspond to noise. 2. The vector D can be found from the signal subspace U s in Eq. (26) . 3. The vector D can also be obtained from the noise subspace. The signal subspace U can be partitioned as u oo U 0l U 0N-1 u s = U 10 Ujj U 1N-1 = _ U M0 U M1 U MN-1_ 45) where r\ H is the first row and U' s consists of the remaining rows of U s . These two matrices can be written as T] H = [u 00 U 0] U 0N-l\ U s = U lf) Ujj 1 1N-I U M0 U M1 ■ ■ ■ U MN~1 (46: The D vector can be obtained as D = -U s t\, V-T\ H T\) 47 4. Once the vector D is obtained, the angle estimator function can be written as / Pmn(8) = a(Q)DD H a(Q) H 27 48) So far several high resolution techniques have been applied ■ to the problem of resolving the directions of arrival of incoming signals. From a statistical point of view, the signals arriving at the array can be regarded as random, and thus the signals and noise are uncorrelated. In the next chapter, we implement these methods based on finite observations . 28 IV. IMPLEMENTATION AND COMPUTER CODES A. PARAMETRIC DATA MODEL AND FLOWCHART The AOA estimation methods described in Chapter III were implemented in a computer code. Output data were gathered for a variety of arrays and signals, and a performance comparison made. This chapter describes the structure and operation of the computer simulation, which is written in Mat lab. The approach is model -based in the sense that it relies on certain assumptions made on the observation data. Since the natural causes responsible for signals and noise are often unrelated, it is customary to assume that the signals and noise are uncorrelated with each other. Figure 6 shows the flowchart of the computer program, with the main steps described in the following sections. 1 . Input Parameters For an antenna array with N elements, the size of correlation matrix (correlation between /th and j th element) is N by N . In this case, the maximum number of sources that can be detected cannot exceed (iV-1) . Because the effectiveness of these spatial spectrum methods are a function of signal-to-noise (SNR) ratio, we evaluate the performances of each antenna array by changing SNR at the input. Also, we must consider the spatial resolution capability of each estimation method when applied to each array geometry. For example, test signals that have the same azimuth or elevation angle are generally harder to resolve 29 1 . Number of elements 4. Angles of arrival 2. Number of sources 5. Element spacing 3. Sample points 6. SNR IT 1. Signal model: Laplacian distribution random processing 2. Complex white noise model I Antenna Array Model — Steering vector Uniform Linear Array Non- Uniform Linear Circular Array £ Spiral Array Linear and Circular Array Array output matrix iz Covariance Matrix T Singular Value Decomposition JL Spatial spectrum method MUSIC PHD Min-norm Capon's MLM ESPRIT 1Z lz Iz 2Z Plot: Angle of arrival vs. Power ±z. Angle Figure 6 . Computer simulation flowchart than those not located in a common principal plane. We also examine the effects of the number sampling points. The more 30 sampling points, the more accurate the AOA estimate but it occurs at the expense of increased processing time. Generally the antenna spacing must be less than 0.5 A in order to avoid the ambiguity problem. 2 . Signal and Noise Models In practice, the white noise and signals are modeled as a independent and identically distributed random processes. That means we interpret the received signals as the sample of some waveform at certain specified instants of time, and is referred to as a realization of a discrete time random process. In all cases the signals that arrive can be regarded as random, we assume that signals s(n) are Laplacian probability distributed processes. Because all arrays are not diverse in polarization (i.e., they have all elements identically polarized) , we assume that the frequency and polarization of all received signals are the same. Multi-path and re-radiators in the vicinity of receiving site which may have a mutual coherence are not considered. The signal and white noise generation process are described by the following Matlab pseudo-code: u{ri) = rand(rc) - 0.5 s( n ) = -(u(n)) * log(l - 2* | u(n) |) * A N p w(n) = —^v(n) V2 v(n) = randn(«) + j* randn(rc) (49) 31 where n =0, 1 , . . . , 4095, A is the signal amplitude, N the peak amplitude of white noise and the noise v(n) is a complex Gaussian random variable with both the real and imaginary parts having a variance equal to 1. The Matlab functions rand and randn can be used to generate Gaussian distributed random values. The voltage ratio A/N determines the SNR; in dB it is 2 log ]0 ( A/ N ) . For convenience we let A =1 so that N p =\0- sm/2 ° . (50) 3 . Antenna Array Model ■Based on the discussion in Chapter II the steering matrices for five array geometries can be determined. The steering matrix for the uniform linear array takes the form 1 ... i AQ) = exp(j2ii(N - l)d sin (J) 7 ) exp(j2n(N - l)d sin (j) M ) 51) where N is the number of antenna elements, M is the number of sources and § M is the AOA of the Mth source, - ft/ < rh < ft/ The steering matrix for a the non-uniform linear array is 32 AQ) = Qxp(j2nd x _]S'm§i) Qxp(j2nd iW _,sm^ M ) 52 where d l is the spacing between the reference element and the (/+l)th antenna. For a minimum redundancy array, the d t are specified in Table 1. The steering matrix for circular array is exp(J2nR cos(-(f) / )sinO / ) • • • exp(j2nR cos(-(j) / )sinG A/ ) ,4(9, 40 Qxp(j2nR cosC^^sinG/) exp(j2nR cos((|) NM ) sin M ) 53 where <j>„ m - n§ - (j) m , Q m is the elevation angle of mth source, (J) m is the azimuth angle of mth source and R is the radius of the circle. For the present analysis it is assumed that R is equal to 0.5 A, < ty m < 2n, and < m < ^ . The steering matrix for the spiral array is Qxp(j2nRj cos(-$ I )smQ ! ) • • • exp(j2nR; cos(-(j) ; )sin0 A/ ) A(Q,$) = Qxp(j2nR N cos((j) N1 ) sin l ) exp(./27t# v cos(<|) ^ ) sin M ) 54 where R n is the distance between nth antenna element and center point. In our model, the maximum R n is 0.5/1. 33 The steering matrix for the combination linear and circular array can be partitioned into sub-matrices associated with the linear and circular parts of the array ^(0,*) 1 Qxp(j2nds'mQ l ) exp(j4nd sinO ,) exp(j2nd cos(-§ j)smQ j) exp(j2nd cos((j) v/ ) sin 9 , ) 1 exp( j 2nd smQ M ) exp(j4nd s'mQ M ) exp(j2nd cos(-<b M )smQ M ) z*p(j2nd cos(<{> NM ) sin G M ) 55 where d is . 4 A . The first three rows in Eq . (55) are for the linear array elements. For these terms there is no dependence on arrival azimuth angle. 4. Correlation Matrix and Spatial Spectrum Methods The output of the array can be calculated by using the signal, noise, and antenna models in Eq. (24) . First the Matlab function is used to find the singular value decomposition. The result is substituted into the various algorithms (Eq. (27), (41), (43), or (48)) with the appropriate antenna array steering matrix (Eq. (51) , (52) , (53), (54), or (55)). The outputs for MUSIC, PHD, minimum norm and Capon's beamformer are angle of arrival vs. power. The power must be computed at every angle. We must steer the antenna array in all azimuth directions (0 to 360 degrees) and elevation direction (0 to 90 degrees) in order to find the AOAs . The output of ESPRIT is the direction of arrival 34 which can be obtained without a search. This method has only- been applied to a uniform linear array. B. THE IMPLEMENTATION OF ARRAY MATRIX DECOMPOSITION The signal measured at the output of any element in the array differs from the signal actually received by an amount attributed to noise. Thus, we have an observed signal, which for each element of the array, consists of the actual received signal plus narrow-band noise. Calculations used 6=4096 time sample points for each antenna element, so the size of received signal matrix is N by 4096, where N is the number of elements. The size of signal -in-space matrix is M by 4096, where M is the number of sources. The observed signal matrix is equal to the received signal matrix times the propagation delay matrix. The observed signal matrix is the antenna output matrix G , from which we find the covariance matrix of the observed signal matrix defined by R = E [G* G T ] . The covariance matrix R has two important special properties. First, it is Hermitian, that is, it is equal to the conjugate transpose of itself. Second, it is positive semi-definite if element noise is present. Thus, for computer coding purposes R = E [G* G T ]= E [G G H ] . (56) 35 The size of R matrix is N by N . The Matlab function SVD is used for singular value decomposition. It has the form [U, S, U] = SVD(R) where R = USU H (Eq. (26)). The eigenvector U is an N by N orthogonal signal and noise subspace matrix, and S is an TV by N diagonal matrix. The singular value decomposition is accurately computed by this function. The columns of U are the eigenvectors of RR , and the entries along the diagonal of S are the correspondingly ordered nonnegative square roots of the eigenvalues of R R , which are also the eigenvalues of RR . The numbers in S are the singular values of the matrix R , which is determined by the number of sources. After calculation of U and S , they are substituted into the spectral estimators in Eqs . (27), (43), and (48) to determine the AOAs . The MUSIC and PHD methods make use of the noise subspace matrix; the minimum norm method uses the signal subspace matrix. Capon's beamformer does not use the singular value decomposition of the covariance matrix. For 2-D direction finding, it takes a significant amount of CPU time to search space in order to find the AOAs. The default search range in azimuth angle is 90 degrees. 36 c. GRAPHICAL USER INTERFACE To compare the relative performance of the various arrays and algorithms, the simulation must be run for many- cases. A graphical user interface (GUI) was developed to vary the input parameters, as shown in Figure 7. Also included is an on-line help feature. The help button provides brief explanations of the input parameters. The GUI provides default values to demonstrate the angular resolution of MUSIC for a uniform linear array and a minimum redundancy array with different signal-to-noise ratios. antenna number: 7 spacing: f^" source number: Signal to Noise Ratio(dB): CLOCK CLOSE HELP sample points: pioiT 1-D DIRECTION FINDING -90 ...90 **■ 'y ■*■* * * i i i>- ■> Unifiorm and Non-uniform Linear Array 2-D DIRECTION FINDING: EL angle:0...90deg. AZ angle: 0. AZ angle: j 16 | 22 | 28 | 34 MUSIC-resolution- MOVIE Linear Array with ESPRIT EL angle: | 10 j 15 j 20 [ 25 2-D DIRECTION FINDING CIRCULAR ARRAY Antenna Geometry Uniform Linear Array Non-uniform linear array 2-D DF Array 2-D DIRECTION FINDING SPIRAL ARRAY 2-D DF LINEAR+CIRCDLARARRAY LA-spacing: g.4 LA-ant No: ~3~ Antenna Pattern Linear* circular array Circular array Figure 7 . Graphic user interface In the next chapter, we use this GUI to present a performance comparison of superresolution spectral 37 algorithms for elevating and analyzing their effectiveness in several antenna and emitter geometries. 38 V. SIMULATION RESULTS AND ANALYSIS In this chapter the results of the computer simulation are presented. The performances of the various algorithms with different array and emitter geometries are examined. Computer simulations were run for the conditions of: (1) two closely spaced AOAs for testing azimuth and elevation angular resolution, (2) various levels of additive white noise, (3) multiple signal sources, and (4) changes in the antenna element spacing. By comparing the test results, we can determine which combination of array and algorithm operates best under the given set of conditions. The power spectra are plotted as surfaces. For all cases, the z -axis represents magnitude of the power, the x-axis represents as azimuth angle and the ;y-axis the elevation angle. The majority of the cases use 7 antenna elements; the exceptions are Figures 3 5 and 3 7 which are for 10 elements. A. ANGULAR RESOLUTION OF ONE -DIMENSIONAL DF Initial simulations were run to determine the comparative performance of algorithms in resolving closely spaced AOAs. The signals were assumed to be of equal frequency, polarization and mutually incoherent. Simulations were run for two signals with angles of arrival of 15 and 18 degrees; the corresponding spectral estimates are shown in Figure 8 . The rows correspond to spectra estimates for the PHD, MUSIC, minimum norm and Capon's beamformer algorithms, and the columns correspond to 39 the uniform array (left) and non-uniform linear array (right) . From this figure, we can see that the non-uniform array can achieve better angular resolution. All algorithms except the Capon's beamformer provide sharp distinguishable peaks, but MUSIC has a smoother spectrum (the separation of the peaks is not as pronounced) . Figure 9 shows an interference simulation for 6 sources. The non-uniform (minimum redundancy) array provides better dynamic range for all algorithms. B. TWO-DIMENSIONAL DIRECTION FINDING 1. Azimuth Angular Resolution A range of SNRs were simulated. Figures 10, 11, and 12 show the results for a high SNR (=25 dB) , at an elevation angle of 15 degrees and AAZ=5 degrees. The MUSIC, PHD and min-norm methods have similar results for all three 2-D array geometries except the Capon's beamformer. However, when the SNR is lower (15 dB) , the circular array performs better, as indicated by a comparison of Figures 13, 14, and 15. Even though two peaks exist, there is no notch between them and therefore it may be difficult to correctly resolve the two spectral peaks. Figures 16, 17, and 18 show the results for SNR=3 dB, and two sources at the same elevation angle but separated in azimuth ( AAZ = 3 degrees) . All arrays can resolve the two peaks. However, for the lower SNR level case, the circular array is generally superior. 40 Uniform linear array m PHD | T3 (D -50 A 1 Q. .mn — ^^ Non-unifonn linear array MUSIC 1 \ -50 -100 -50 50 AOA(degree) -50 -100 -50 AOA(degree) 100 100 100 100 Figure 8. Comparison of uniform (left) and non-uniform (right) linear arrays , 2 sources , 7 elements , separation angle 3 degrees (15° and 18°) . SNR=25 dB Uniform linear array -100 -50 50 AOA(degree) Non-uniform linear array 100 m ■a 100 CD 100 -50 -100 -100 -50 -50 -100 -100 -50 50 100 MUSIC L." 100 Mini-norm 100 100 -50 50 AOA(degree) 100 Figure 9. Unifrom (left) and non-uniform (right) linear arrays, SNR=25 dB, 7 antennas, 6 sources: 15°, 25°, 35°, 45°, 55°, and 65° 41 MUSIC PHD 50 EL(deg) 50 AZ(deg) Mini-norm EL(deg) AZ(deg) Capon MLM 100 EL(deg) AZ(deg) EL(deg) AZ(deg) Figure 10. Linear+circular array, SNR=25 dB, 2 sources at (16°. 15°) and (21°, 15°) . 7 elements MUSIC EL(deg) AZ(deg) Mini-norm PHD ~ 60- QQ Z. 3 40- 1 20- S. o. -20 J EL(deg) AZ(deg) Capon MLM EL(deg) AZ(deg) 100 50 AZ(deg) Figure 11. Circular array, SNR=25 dB, 2 sources at (16°, 15°) and (21°. 15°) . 7 elements 42 MUSIC PHD EL(deg) AZ(deg) Mini-norm 50 50 EL(deg) AZ(deg) Capon MLM 50 AZ(deg) 50 AZ(deg) Figure 12. Spiral array, SNR=25 dB, 2 sources at (16°, 15°) and (21°, 15°) , 7 elements MUSIC PHD 60 i m 40- -a |f 20- 1 o. -20- ~ 60- CD .„ % 40- 1 20- ° ■ a. u -20 J EL(deg) AZ(deg) Mini- norm EL(deg) AZ(deg) Capon MLM EL(deg) AZ(deg) EL(deg) AZ(deg) Figure 13. Circular array, SNR=15 dB, 2 sources at (16°, 15°) and (21°, 15°), 7 elements 43 MUSIC PHD EL(deg) AZ(deg) Mini-norm EL(deg) AZ(deg) Capon MLM EL(deg) AZ(deg) 100 50 AZ(deg) Figure 14. Linear+circular array, SNR=15 dB, 2 sources at (16°, 15°) and (21°,16°) / 7 elements MUSIC PHD _..••■' 60- m 40 • f 20- 1 0- -20 - • • ^im 50^ -"~"~50 E L(deg) a: ^(deg) norm CD "C3 EL(deg) AZ(deg) Capon MLM EL(deg) AZ(deg) 100 50 AZ(deg) Figure 15. Spiral array, SNR=15 dB, 2 sources at (16°, 15°) and (21°, 15°), 7 elements 44 MUSIC PHD EL(deg) AZ(deg) Mini-norm EL(deg) AZ(deg) Capon MLM EL(deg) AZ(deg) EL(deg) AZ(deg) Figure 16. Circular array, SNR=30 dB, 2 sources at (16°, 15°) (19°, 15°) , 7 elements and MUSIC EL(deg) AZ(deg) Mini-norm EL(deg) AZ(deg) PHD EL(deg) AZ(deg) Capon MLM EL(deg) AZ(deg) Figure 17. Linear+circular array, SNR=30 dB, 2 sources at (16°,15°) and (19°. 15°). 7 elements 45 MUSIC PHD 50 Mini-norm 30 j ■ ■•■"'.'.'/:■■■■'"'.'.[ "■•■!!' ■••.. m 60 j ?■ 40- .;..--i i /^^ | 20- §. o. ^^^^H ^^_ 'Lj^^CncSkVBP^^^^^ -20 J 50 80 h m 60- ?- 40- | 20- g. 0- -20 J 50 50 EL(deg) AZ(deg) Capon MLM EL(deg) AZ(deg) EL(deg) AZ(deg) Figure 18. Spiral array, SNR=30 , 2 sources at (16°, 15°) and (19°, 15°) , 7 elements 60 I 20 r 60 • od 40 I 20 -20 MUSIC EL(deg) Mini-norm 60 80 20 40 60 80 EL(deg) PHD 80 60 r m /in -o 40 ■ A: ■ ■ :.... 50 EL(deg) 100 Figure 19. Circular array, SNR=25, 2 sources at (16°, 15°) and (16°, 18°) , 7 elements 46 2 . Elevation Angular Resolution Simulations were run for the two signal case with angles of arrival at 18 and 15 degrees in elevation and 16 degrees in azimuth. For the SNR equal to 25 dB, the corresponding spectral estimates are shown in Figures 19, 20, and 21. The results for the circular and linear+circular arrays are almost the same, but the spiral array cannot correctly resolves the peaks. For a 15 dB SNR the data, which is shown in Figures 22 and 23, indicate that all algorithms, except Capon's beamformer, correctly resolve the two spectral peaks with a circular array. Only two algorithms can resolve the peaks for the circular array, so the elevation angular resolution of the linear+circular array appears to be better for this scenario. 3. Multiple Signal Sources Multiple signal simulations are shown in Figures 24, 25, and 26. There are four signals arriving at azimuth and elevation angles of (10°, 20°) , (55°, 15°) , (30°, 16°) , and (45°, 25°) with equal amplitudes. The circular antenna array and linear+circular array can resolve the four peaks but the circular array has more sidelobes. Only the MUSIC algorithm correctly resolves the signal sources for the spiral array. In this scenario the linear+circular array is best. 47 MUSIC PHD LLJ in O Q. LI I 41 20 -20 I 60 40 20 -20 1TJ ili ::: 20 40 BO EL.(deg) Mini-norm CD i — Hi o a. 60 40 20 -20 50 40 30 20 10 , : I- Stf&5 MftfcmJ^i^i^j^^Bn i i . i 20 40 60 80 EL(deg) Capon MLM 20 40 60 EL(deg) 80 Figure 20. Linear+circular array, SNR=25 dB, 2 sources at (16° / 15°) and (16°, 18°), 7 elements MUSIC PHD 60 S 40 I 20 -20 20 40 60 ELfdeg) Mini-norm 80 _ 60 CD ^40 1 20 : : : : • Ill 1 : : " » -ZAindaeMaa** ' ■ ■ t ■ 20 40 60 80 EL(deg) Figure 21. Spiral array, SNR=25 dB, 2 sources at (16°, 15°) and (16°, 18°), 7 elements 48 MUSIC PHD 40 LLJ 5 20 o -20 20 40 60 80 EL(deg) Mini-norm 20 40 60 80 EL(deg) apon MLM BO ~ 40 oo — 5 20 -20 20 40 60 EL(deg) 80 Figure 22. Circular array, SNR=15 dB, 2 sources at (16°, 15°) and (16°, 18°), 7 elements MUSIC PHD 60 r £T 40 | 20 b -20 20 40 60 80 EUfdeg) Mini-norm EL(deg) Capon MLM 60 ■ m 40 • -a I 20 o -20 20 40 60 EL(deg) 80 100 Figure 23. Linear+circular array, SNR=15 dB, 2 sources at (16°, 15°) and (16°, 18°), 7 elements 49 Mi i:-";ii PHD 60- S 40- powen r j isj CD CD CD 50 "S* 50 EL(deg) AZ(deg) Mini-norm 50 AZ(deg) Capon MLM 50 CD [ ^M 50 AZ(deg) 100 Figure 24. Linear+circular array, SNR=25 dB, 7 elements, 4 sources at (10°, 20°) , (55°, 15°) , (30°, 60°) , and (45°, 25°) MUSIC 50 AZ(deg) Mini-norm PHD EL(deg) AZ(deg) Capon MLM EL(deg) AZ(deg) 50 AZ(deg) 100 Figure 25. Circular array, SNR=25 dB, 7 elements, 4 sources at (10°, 20°) , (55°, 15°) , (30°, 60°) , and (45°, 25°) 50 MUSIC PHD 80 a EL(deg) AZ(deg) Mini-norm 50 AZ(deg) . •'..■!' • ' 50 50 EL(deg) AZ(deg) Capon MLM 50 AZ(deg) Figure 26. Spiral circular array, SNR=25 dB, 7 elements, 4 sources at (10°, 20°) , (55°, 15°) , (30°, 60°) , and (45°, 25°) MUSIC PHD 50 AZ(deg) Mini-norm 50 AZ(deg) 50 50 EL(deg) AZ(deg) Capon MLM 50 AZ(deg) Figure 27. Circular array, SNR=15 dB, 7 elements, 2 sources at (16°, 15°), and (16°, 20°) 51 4. Comparison of Azimuth and Elevation Angular Resolution The estimation methods and arrays were compared for two emitters with a 5 degree separation. Figures 27, 28, and 29 are for SNR=15 dB and AEL=5 degrees; Figures 13, 14, and 15 for SNR=25 dB,' AAZ=5 degrees. We see that the elevation angular resolution is better than the azimuth angular resolution since the sidelobes in azimuth are higher than the sidelobes in elevation for all three array- configurations. Even at lower SNRs , all three arrays employing the MUSIC method can distinguish the two sources in elevation. 5. Improvement of Angular Resolution of Two-dimensional DF Next we investigate changing the radius of the circular array to 1.5 wavelengths. The results are shown in Figures 30 and 31. Comparing Figures 13 and 14 with Figures 30 and 31 (SNR=15 dB, AAZ = 5 degrees) , the angular resolution is better with the larger radius as expected. Larger spacing may be susceptible to ambiguities. From Figures 32 and 33, where the range of azimuth angles is from to 360 degrees, it is evident that there is no serious ambiguity problem. Thus, we can improve the angular resolution by . changing the antenna spacing. 52 MUSIC PHD m 40- t3 *=• 20- a> 1 o. a. -20 ^ EL(dgg) AZ(deg) Mini-norm 50 AZ(deg) 50 50 EL(deg) AZ(deg) Capon MLM 50 AZ(deg) Figure 28. Linear+circular array, SNR=15 dB, 7 elements, 2 sources at (16°, 15°), and (16°, 20°) MUSIC 50 AZ(deg) Mini-norm 50 AZ(deg) PHD 60 -I 50 50 EL(deg) AZ(deg) Capon MLM Figure 29. Spiral array, SNR=15 dB, 7 elements, 2 sources at (16°, 15°), and (16°, 20°) 53 The circular array has a smoother spectrum and less spikes than linear+circular array, since the spacing of linear array part is larger than o . 5 A, . From these simulation results, we conclude that the larger antenna spacing will enhance angular resolution at the expense of more ripples on the spectrum. We can suppress the sidelobes and enhance resolution by increasing the number of elements in the x - y plane and increasing the radius of circle and then adding a linear part on z -axis with a spacing less than 0.5/1. Two elevation cases are considered. For a large elevation angle (80 degrees) , there are four signals arriving at azimuth and elevation angles of (10°, 80°), (90°, 80°) , (180°, 80°) , (270°, 80°) , with a SNR=15 dB . The simulations are shown in Figure 34 for the circular array (7 elements, R n =1.5A), and in Figure 35 for the circular+linear array (3 elements, d =0.4/1 on z -axis) . In both cases the sidelobes increase in the elevation direction, but Figure 35 shows lower sidelobes. From these results we conclude that the elevation angular resolution is worse when the elevation angle is large. This is a well known result that is associated with the increase in beamwidth that occurs with array scanning [Ref . 9] . For small elevation angles (5 degrees) , there are four sources arriving at azimuth and elevation angles of (10°, 5°) , (90°, 5°) , (180°, 5°) , (270°, 5°) , and a SNR=15 dB . The 54 simulations are shown in Figure 36 for the circular array (7 elements, R n =1.5A), and in Figure 3 7 for the circular+linear array (3 elements, d=0.4A). Figure 37 shows lower sidelobes in the plane of the sources. These sidelobes expand in the azimuth direction, implying that the azimuth angular resolution is worse when the elevation angle is small. This case is close to a broadside condition for the linear array and, consequently, the resolution ability is essentially that of the linear array. Thus the resolving ability has been improved at low elevation by adding a linear array part (3 elements on the z axis) . Based on the results of the computer simulations and performance comparisons, we present some conclusions and recommendations in the next chapter. 55 MUSIC PHD EL(deg) AZ(deg) Mm i- norm U 1 , m 40- » 20- ^ n o • Q. -20 J EL(deg) AZ(deg) Capon MLM EL(deg) AZ(deg) 50 AZ(deg) 100 Figure 30. Circular array, R n =1.5/l r SNR=15 dB, 7 elements, 2 sources at (16°, 15°), and (21°, 15°) MUSIC PHD EL(deg) AZ(deg) Mini-norm m 40 EL(deg) AZ(deg) Capon MLM EL(deg) AZ(deg) 1 00 Figure 31. linear+circular array, d =1.5/1, SNR=15 dB, 7 elements, 2 sources at (16°, 15°), and (21°, 15°) 56 EL(deg) AZ(deg) Figure 32. Circular array, R =1.5A,, SNR=15 dB, 7 elements, 2 sources at (16°, 15°), and (21°, 15°), MUSIC method EL(deg) AZ(deg) Figure 33. Linear+circular array, d =1.5 A, SNR=15 dB, 7 elements, 2 sources at (16°, 15°), and (21°, 15°), MUSIC method 57 EL(deg) AZ(deg) Figure 34. Circular array, SNR=15 dB, 7 elements, 4 sources at (10°, 80°) , (90°, 80°) , (180°, 80°) , and (270°, 80°) , MUSIC method EL(deg) AZ(deg) Figure 35. Linear+circular array, SNR=15 dB, 10 elements , 4 sources at (10°, 80°) , (90°, 80°) , (180°, 80°) , and (270°, 80°), MUSIC method EL(deg) AZ(deg) Figure 36. Circular array, SNR=15 dB, 7 elements, 4 sources at (10°, 5°) , (90°, 5°) , (180°, 5°) , and (270°, 5°), MUSIC method EL(deg) 111 AZ(deg) Figure 37. Linear+circular array, SNR=15 dB, 10 elements, 4 sources at(10°,5°) , (90°, 5°) , (180°, 5°) , (270°, 5°) , MUSIC method 59 60 VI. CONCLUSIONS Over the past few decades superresolution techniques have been developed to resolve multicomponent plane wave fields. The most noteworthy of these is the MUSIC algorithm. Others include PHD, Minimum norm, ESPRIT and conventional and Capon's beamforming. Each of the methods has its advantages and disadvantages with regard to hardware and computational aspects (array configurations and element spacing) and emitter parameters (number of emitters and their spatial distribution) . This thesis has examined the performance of four of the estimation methods listed above for several linear, planar and volumetric arrays. Based on the formulations described in Chapter III, Matlab simulations of the methods were used to collect data for linear, circular, spiral and circular plus linear arrays. The emitter locations and spacings were varied, as was the signal-to-noise ratio. The data are presented as mesh surfaces, and the surface characteristics can be interpreted in terms of emitter resolution. For one -dimensional DF, the MRA provides the best dynamic range and resolution for a given number of elements and baseline. Since ESPRIT is generally applied to the spacing which is a known constant value, it is not directly applicable to MRAs . Using ESPRIT, the AOA can be obtained without a search technique, and in this regard is different 61 from the other four methods. Therefore, only the other four algorithms have been compared. The simulation results are in agreement with related two-dimensional DF performance simulations reported by Johnson and Miner [Ref . 10] for the MUSIC and Capon's MLM algorithms. The data in Chapter V leads to the following conclusions with regard to the circular and circular+linear array: (1) for high SNR levels, all four methods have about the same the performance, (2) elevation angular resolution is better than azimuth angular resolution for these three array configuration considered, (3) the azimuth angular resolution of the circular array is better than the spiral or circular+linear, (4) the elevation angular resolution of circular+linear array is better than spiral and circular+linear, (5) for multiple signal sources, the linear+circular array is better when antenna spacing is less than 0.5 wavelength, (6) larger element spacings result in better resolution, but there are more ripples and spikes; ambiguities will occur if the spacing becomes too large, (7) when the elevation angles of two signal sources are small, the azimuth angular resolution is worse than it is at large elevation angles, (because of slow falloff of the pattern in azimuth) 62 (8) when the elevation angles of two signal sources are large, the elevation angular resolution is worse than when elevation angles are small (because of slow falloff of the pattern in elevation) . From the data it appears that the MUSIC algorithm is more stable, almost always correctly resolves individual emitters, and more importantly, has low sidelobes . However, it has been demonstrated that under some conditions MUSIC fails to provide the correct number of peaks [Ref . 11] . The circular array has best azimuth resolution but the results of linear+circular array are close to that of the circular array. Based on multiple signal sources the linear+cicular array is probably the better choice. When the antenna spacing is larger (say R =1.5 wavelength), the circular array is better. Having more elements on a larger diameter does provide some advantages . The distance between two adjacent antennas will be small to avoid high sidelobes and offers better azimuth resolution, especially for small elevation angle. Since we use subspace methods to figure out the AOA, the number of sources must be less than the number of array elements. Future work should consider how to deal with this limitation as well as improvements to the low signal-to- noise-ratio case. Any techniques that would allow a decrease in the search time for 2-D direction finding would also be valuable. 6 3 64 APPENDIX A - THESIS MAIN PROGRAM clear all close cv=version; if str2num (cv ( 1 ) ) <5 figure ( 1 ) , uicontrol ( ... 'Style', 'text', . . . ' Units ',' normalized* , ... 'Position' , [0.01 0.3 0.9 0.1 ], ... 'BackgroundColor ' , [0. 8 0.8 0.8], ... ' HorizontalAlignment ' , ' left ' , ... 'String' , 'This program cannot be run on MATLAB 4. X, please use MATLAB 5 . X ' ) ; pause(5), close, exit end load nsource msens=7; snr=2 5; %nsource=2; uicontrol (' Units ',' points ' , ... 'BackgroundColor ', [0. 8 0.8 0.8], ... ' FontSize' ,10, ... ' ForegroundColor' , [0 1], ... 'Position' , [23.5 288.75 125 17.25], ... ' String ' , ' antenna number : ' , ... 'Style' , 'text ' , ... 'Tag' , 'StaticTextl ' ) ; uicontrol (' Units ', 'points ' , ... 'BackgroundColor' , [0. 8 0.8 0.8], ... ' FontSize' ,10, ... ' ForegroundColor ', [1 0], ... 'Position' , [22.75 273.75 125.75 15], ... ' String' ,' source number:', ... 'Style' , 'text' , ... 'Tag' , 'StaticTextl ' ) ; uicontrol ('Units', 'points', . . . 'BackgroundColor' , [0. 8 0.8 0.8], ... 'FontSize' ,9, ... 'Position' , [22 254 125 19], ... 'String' , 'Signal to Noise Ratio(dB):', ... 'Style' , 'text' , ... 'Tag' , 'StaticText2' ) ; msensl = uicontrol (' Units ',' points ' , . ! — ! 'BackgroundColor ', [0 10], ... 'Callback' , ' inn3 ' , ... 'FontSize* ,12, ... •Position' , [150.5 290.25 27 15], ' String ' , num2str (msens ) , ... 65 'Style' , 'edit' , ... 'Tag' , 'EditTextl' ) ; nsourcel = uicontrol (' Units ', 'points ' , ... ' BackgroundColor ' , [0 10], ... 'Callback' , 'resl' , ... ' FontSize' , 12, ... •Position' , [150.5 274.5 27 15], ... ' String ', num2str (nsource) , ... 'Style' , 'edit' , ... 'Tag' , 'EditTextl* ) ; snrl = uicontrol (' Units ',' points ' , ... 'BackgroundColor ', [0 10], ... •Callback' , ' inn3 ' , ... 'FontSize' ,12, ... •Position' , [150.5 258 27 15], ... ' String' , num2str (snr) , ... 'Style' , 'edit* ) ;%, ... uicontrol (' Units ',' points ' , ... 'BackgroundColor* , [0.8 0.8 0.8], ... ' FontSize' , 10, . . . ' ForegroundColor ' , [1 0], ... •Position', [40.9655 242.69 100.552 17.3793], ... 1 String ' , ' sample points : * , ... 'Style' , 'text' ) ; bnsamp=uicontrol ( ' Units * , ' points ' , ... 'BackgroundColor ', [0 10], ... 'Callback' , ' inn3 ' , ... 'Position' , [150.207 242.069 26.069 14.2759], ... 'String' , '4096' , ... 'Style' , 'edit' ) ; uicontrol ( ' Units ' , 'points ' , ... 'BackgroundColor ', [0.8 0.8 0.8], ... ' FontSize' ,10, ... ' ForegroundColor' , [0 1], ... 'Position', [178.759 287.379 54 16.7586], ... ' String ',' spacing :' , ... 'Style' , 'text' ) ; bdspace= uicontrol (' Units ', 'points ' , ... 'BackgroundColor' , [0 10], ... 'Callback' , ' inn3 ' , ... 'Position' , [232.138 288.621 26.069 14.2759], ... ' String' ,'0.4', ... 'Style' , 'edit' ) ; aw=clock; tname = ['Direction Finding Analysis year:', num2str (aw ( 1 ) month: ', num2str (aw (2 ) ) ,',-• -day: ', num2str (aw (3) ) ] ; set (gcf , ' name ' , tname) nsource=str2num (get (nsourcel , ' string ' ) ) ; bdspacel= uicontrol (' Units ', 'points ' , ... 'BackgroundColor ', [0 10], ... 'Callback' , * inn33 ' , ... •Position', [252 27 27 12], ... 'String' , '0.4' , ... 'Style' , 'edit' ) ; bdl= uicontrol (' Units ', 'points ' , ... 'BackgroundColor ', [0 10], ... 66 'Callback' , 'inn33' , ... •Position' , [252 14 27 12] , 'String' , '3' , ... 'Style' , 'edit' ) ; uicontrol (' Units ',' points ' , ... 'BackgroundColor' , [0. 8 0.8 0.8], ... ' FontSize' ,10, ... ' ForegroundColor ' , [ 1 ] , ... 'FontWeight* , 'bold' , ... 'Position', [11. 1724 222.828 218.25 15.7172], ... 'String' ,' 1-D DIRECTION FINDING -90 ...90 deg', ... 'Style' , 'text ' , ... 'Tag' , 'StaticText3' ) ; uicontrol (' Units ' , ' points ' , . . . 'BackgroundColor' , [0. 8 0.8 0.8], ... ' FontSize' ,10, ... ' ForegroundColor ' , ' r ' , ... 'FontWeight ', 'bold' , ... 'Position' , [16.25 207.5 41.75 15], ... ' String ',' angle :' , ... 'Style' , 'text ' , ... 'Tag' , 'StaticText4 ' ) ; for n=l:8 bbl=ui control (' Units ' , 'points ' , . . . 'BackgroundColor' , [0.8 0.8 0.8], ... 'Position' , [56. 25+(84-56. 25) * (n-1) 207 26.25 14.25], 'String' , ' ' , ... •Style' , 'text' ) ; end %clear bbl %bbl=[] for n=l:nsource bbl (n) =ui control ('Units','points', ... •BackgroundColor ', [0 10], ... ' Callback' ,' inn3 ' , ... 'Position' , [56. 25+(84-56. 25) * (n-1) 207 26.25 14.25], 'String' , int2str(-25+10*n) , ... 'Style' , 'edit' ) ; end size (bbl ) ; bbl; uicontrol (' Units ', 'points * , ... 'BackgroundColor' , [0.8 0.8 0.8], ... ' FontSize' ,9, ... ' ForegroundColor ' , ' r ' , ... 'FontWeight' , 'bold' , ... 'Position' , [2.75 126.75 58.75 12.75], ... 'String', ' AZ angle:', ... 'Style' , 'text' ) ; for n=l 67 bb2(n) = uicontrol (' Units ', 'points ' , ... 'BackgroundColor' , [0.8 0.8 0.8], ... 'Position', [58. 25+(90-62. 25)* (n-1) 126.75 26.25 14.25], ... 'String' , ' ' , ... 'Style' , 'text' ) ; end clear bb2 bb2= [ ] ; for n=l:nsource bb2(n) = uicontrol (' Units ', 'points ' , ... 'BackgroundColor' , [0 10], ... 'Callback' , 'inn3' , ... •Position* , [58.25+(90-62.25) * (n-1) 126.75 26.25 14.25], ... 'String' , num2str (10+6*n) , ... 'Style', 'edit' ) ; end uicontrol (' Units ', 'points ' , ... 'BackgroundColor' , [0.8 0.8 0.8], ... 'FontSize' ,10, ... ' ForegroundColor ' , [0 1], ... ' FontWeight ' , ' normal ' , ... 'Position' , [12 156 207 32.25], ... 'String', '2-D DIRECTION FINDING: EL angle:0...90 deg. AZ angle 0. . . 360 deg', . . . 'Style' , 'text' ) ; uicontrol ( ' Units ' , 'points ' , ... 'BackgroundColor' , [0.8 0.8 0.8], ... 'FontSize' ,9, ... 'FontWeight ', 'bold' , ... 'Position' , [2 108 58.75 12.75], ... 'String', 'EL angle:', ... 'Style' , 'text' ) ; for n=l: 8 bb3 (n) = uicontrol (' Units ', 'points ' , ... 'BackgroundColor' , [0.8 0.8 0.8], ... •Position' , [58+(90.75-63) * (n-1) 108 26.25 14.25], 'String' , ' ' , ... 'Style' , 'text' ) ; end clear bb3 bb3= [ ] ; for n=l:nsource bb3(n) = uicontrol (' Units ', 'points ' , ... 'BackgroundColor' , [0 10], ... •Callback' , ' inn3 ' , ... 'Position' , [58+(90. 75-63) * (n-1) 108 26.25 14.25], ' String' , num2str (5 + 5*n) , ... 'Style' , 'edit' ) ; end bval= uicontrol (' Units ', 'points ' , ... 'BackgroundColor' , [0 1], ... 'Callback' , * inn3 ' , ... 68 ' FontSize' ,10, ... ' ForegroundColor ' , [1 10], ... 'Position', [190.552 153.31 36.6207 19.2414], ... 'String' ,[' 90 ' ; ' 180 ' ; ' 270 ' ; ' 360 ' ] , ... ' Style ' , ' popupmenu ' , ... 'Tag' , ' PopupMenul ' ) ; clear bearingll clear bearingn clear ang msensll=str2num (get ( msensl , 'string ' ) ) ; snr=str2num (get (snrl , 'string ' ) ) ; dspace=str2num (get (bdspace, ' string ' ) ) ; nsamp=str2num (get (bnsamp, ' string ' ) ) ; thrv=get (bval, 'Value' ) if (thrv==l) thr=90; elseif (thrv==2 ) thr=180; elseif (thrv==3) thr=27 0; elseif (thrv==4) thr=360; end for n=l:nsource bearingll (n) =str2num (get (bbl (n) ,' string ') ) ; bearingn (n) =str2num (get ( bb2 (n) , 'string ' ) ) ; ang (n) =str2num (get (bb3 (n) , ' string ' ) ) ; end bearingll bearingn ang uicontrol (' Units ' , ' points ' , . . . 'BackgroundColor' , [0.752941 0.752941 0.752941], . ' ButtonDownFcn ' , ' ctlpanel SelectMoveResize ' , ... 'Position' , [288 151.5 125.25 69], ... ' Style ' , ' frame ' , ... 'Tag' , ' Frame 1' ) ; uicontrol ( ' Units ' , 'points ' , ... ' ButtonDownFcn' ,' ctlpanel SelectMoveResize', ... 'Callback' , 'doam' , ... ' FontSize' ,9, ... 'FontWeight ' , 'bold' , ... 'Position' , [230.25 224.25 182.25 28.5], ... ' String ',' Unifiorm and Non-uniform Linear Array', 'Tag' , 'Pushbuttonl' ) ; uicontrol (' Units ',' points ' , ... ' ButtonDownFcn ',' ctlpanel SelectMoveResize', ... 'Callback* , ' eprit ' , ... 'FontSize' ,9, ... 'FontWeight ', 'bold' , ... 'Position' , [293.25 156 116.25 27], ... 'String', 1 Linear Array with ESPRIT', ... 'Tag' , 'Pushbutton2' ) ; uicontrol (' Units ',' points ' , ... ' ButtonDownFcn' ,' ctlpanel SelectMoveResize', ... 'Callback' , 'mvml' , ... 69 1 FontSize' ,9, ... * FontWeight ' , 'bold' , ... 'Position' , [292.5 186.75 117.75 28.5], ... 'String' , ' MUSIC-resolution- MOVIE', ... 'Tag' , 'Pushbutton3' ) ; uicontrol ( ' Units ',' points ' , ... ' ButtonDownFcn ' , ' ctlpanel SelectMoveResize ' , ... 'Callback' , 'doag4 ' , ... ' FontAngle' ,' italic' , ... 'FontName' , 'TIMES' , ... ' FontSize' ,10, ... ' FontWeight ', 'bold' , ... •Position' , [10.5 72.75 270.75 24.75], ... 'String' , '2-D DIRECTION FINDING CIRCULAR ARRAY', ' Tag ' , ' Pushbutton4 ' ) ; uicontrol ( ' Units ', 'points ' , ... ' ButtonDownFcn ',' ctlpanel SelectMoveResize 1 , ... 'Callback' , 'spira' , ... ' FontAngle' , 'italic' , ... 'FontName' , 'TIMES' , ... 'FontSize' ,10, ... 'FontWeight', 'bold' , ... •Position' , [11.25 41.25 270 27], ... 'String' , '2-D DIRECTION FINDING SPIRAL ARRAY', . ' Tag ' , ' Pushbutton5 ' ) ; uicontrol ( 'Units', 'point s', ... ' ButtonDownFcn' ,' ctlpanel SelectMoveResize', ... . 'Callback' , 'cl' , ... ' FontAngle' , 'italic' , ... ' FontName' , 'TIMES' , ... ' FontSize' ,10, ... 'FontWeight' , 'bold' , ... 'Position' , [10.5 14.25 190 22.5], ... 'String' , '2-D DF LINEAR+CIRCULAR ARRAY', ... 1 Tag ' , ' Pushbutton6 ' ) ; uicontrol (' Units ', 'points ' , ... 'BackgroundColor' , [0. 8 0.8 0.8], ... ' FontSize' ,9, ... ' ForegroundColor ' , [0 1], ... 'Position' , [200.759 22.25 54 16.7586], ... ' String ',' LA-spacing: ' , ... 'Style' , 'text* ) ; uicontrol ( ' Units ' , 'points ' , ... 'BackgroundColor' , [0.8 0.8 0.8], ... 'FontSize* ,9, ... ' ForegroundColor ', [0 1], ... 'Position' , [200.759 8.25 54 16.7586], ... 'String' , ' LA-ant No:', ... 'Style' , 'text' ) ; uicontrol ( ' Units ' , 'points ' , ... ' ButtonDownFcn' ,' ctlpanel SelectMoveResize', ... 'Callback' , 'close' , ... ' FontSize' ,12, ... 'FontWeight', 'bold', ... 'Position' , [291.75 258 52.5 25.5], ... 'String' , 'CLOSE' , ... 70 ' Tag ' , ' Pushbutton7 ' ) ; uicontrol (' Units ', 'points ' , ... ' ButtonDownFcn' , ' ctlpanel SelectMoveResize ' , 'Callback' , ' ! clock' , ... ' FontSize' ,12, ... ' FontWeight ' , 'bold' , ... 'Position' , [291.75 284.25 102.75 24.75], ... 'String', 'CLOCK' , ... 'Tag' , ' Pushbutton7' ) ; uicontrol ( ' Units ',' points ' , ... 'BackgroundColor ' , [0. 8 0.8 0.8], ... 'ButtonDownFcn', 'ctlpanel SelectMoveResize', ' FontName ' , ' times ' , ... 'FontSize' ,13, ... 'FontWeight' , 'bold* , ... ' ForegroundColor ' , [1 0], ... 'Position' , [287.25 127.5 128.25 23.25], ' String ', 'Antenna Geometry', ... 'Style' , 'text ' , ... 'Tag' , 'StaticText5' ) ; uicontrol ( ' Units ',' points ' , ... 'ButtonDownFcn', 'ctlpanel SelectMoveResize', 'Callback' , ' sp2 ' , ... ' FontSize' ,9, ... ' FontWeight ', 'bold' , ... 'Position' , [291 87.75 122.25 22.5], ... ' String ',' Non-uniform linear array', ... 'Tag' , 'Pushbutton8 ' ) ; uicontrol ( ' Units ', 'points ' , ... 'ButtonDownFcn', 'ctlpanel SelectMoveResize', 'Callback' , 'spl' , ... 'FontSize ' , 9, ... 'FontWeight ', 'bold' , ... 'Position' , [290.25 111.75 122.25 22.5], ' String' ,' Uniform Linear Array', ... ' Tag ' , ' Pushbutton8 ' ) ; uicontrol ( ' Units ', 'points ' , ... 'ButtonDownFcn', 'ctlpanel SelectMoveResize', 'Callback' , ' g5 ' , ... 'FontSize' ,10, ... ' FontWeight *, 'bold' , ... 'Position' , [291 64.5 122.25 22.5], ... 'String' , ' 2-D DF Array' , ... ' Tag ' , ' Pushbutton8 ' ) ; uicontrol ( ' Units ',' points ' , ... 'ButtonDownFcn', 'ctlpanel SelectMoveResize', 'Callback' , ' helpwin (str, ttl) ' , ... ' FontSize' ,15, ... ' FontWeight ', 'bold' , ... 'Position' , [346.5 258.75 47.25 24], ... 'String' , 'HELP' , ... ' Tag ' , ' Pushbutton8 ' ) ; uicontrol ( ' Units ' , 'points ' , ... 'BackgroundColor' , [0.8 0.8 0.8], ... 'ButtonDownFcn', 'ctlpanel SelectMoveResize', ' FontName ' , ' times ' , ... 71 ' FontSize' ,13, ... ' FontWeight ' , 'bold' , ... ' ForegroundColor ' , [1 0], ... 'Position' , [289.5 41.25 128.25 23.25], ... 'String' , 'Antenna Pattern', ... •Style', 'text' ) ; uicontrol ( ' Units ' , 'points ' , ... ' ButtonDownFcn ' , 'ctlpanel SelectMoveResize ' , 'Callback' , 'cclap' , ... ' FontSize' ,9, ... 'FontWeight' , 'bold' , ... 'Position' , [292.5 27.75 122.25 22.5], ... 'String',' Linear+circular array'); uicontrol ( ' Units ', 'points ' , ... 1 ButtonDownFcn ' , ' ctlpanel SelectMoveResize ' , •Callback' , 'clap', . . . 'FontSize' ,9, ... 'FontWeight', 'bold' , ... 'Position' , [292.5 3.75 122.25 22.5], ... ' String ',' Circular array'); str = { This program was created by LIN, KU-TING, TAIWAN ' Date:9/98 Thesis topic : ' Comparsion of superresolution algorithms with different array geometry' for radio direction finding' Thesis adviser : Prof . David C . Jenn ' second reader: Prof. Pace Phillip' We evaluate the effectiveness of various array designs and processing ' method for two-dimensional direction finding and resolving multiple ' emitters. Using five different spatial spectrum methods to predict the' direction of arrival (DOA) ' Discussion of Topics:' 1. spatial spectrum estimation methods' (a) MUSIC (multiple signal classification' (b) ESPRIT' (c) PHD (Pisarenko harmonic decomposition) ' (d) minmum-norm' (e) Capon beamformer ' 2. Array Geometry' (a) uniform linear array(l-D)' (b) non-uniform linear array (1-D) ' (c) circular array(2-D) ' (d) spiral array(2-D)' (e) linear+circular array (2-D) ' 3. Antenna Pattern' (a) circular array' (b) linear+circular array' 4. Angular Resolution' 72 array 1 HOW TO USE?' ' 1. Input parameter' (a) antenna number: less than 10 for non-uniform linear ' (b)source number:less than antenna number' ' (c) SNR (dB) : signal to noise ratio.' (d) sample points: the more points, the more accurate results ! ' length) ' (e) spacing: 1-D DF antenna spacing (normalized to wave- ' 2-D DF the distance from center point to element (X-Y plane) ' ' (f)LA spacing : antenna spacing for linear array part(Z- axis)of linear+circular array' ' (g)LA-ant No:antenna number for linear array part(Z- axis)of linear+circular array' ' (h) 1-D direction finding angle range:0~180 deg ' (i)2-D DF EL angle:0~90 degree' (j)2-D DF AZ angle:0~360 degree" ' note: choose the AZ searching range --> push popup menu moment ! ' up! ' i i moves . ' i Matlab 5.X ' 10 (=7+3) ! ' resolution ' i sidelobe ' only four choices : 90, 180, 270, 360 degree. 2. When you change the "source number", please wait one 3. Check antenna array geometry and antenna pattern' 4. Choose 1-D or 2-D DF --> Push button, See what you get!' *> — Attention--<* ' 1. The 2-D DFs take some time to search the AOAs, ' please be patient to wait, the word "HOLD ON! !" will show 2. ESPRIT method just for uniform linear array: 20 runs' 3. MUSIC-resolution-movie DEMO --> four tagerts fixed, one sample points=200' 4. This program CANNOT be run on Matlab 4.X, at least 5. The default antenna number of linear+circular array is 6. The larger radius of circular array, the better 7. When the elevation angle of AOA is small (5 degree), the ' expand on the azimuth direction. ' ' 8. When the elevation abgle of AOA is large (80 degree), the sidelobe ' ' expand on the elevation direction.' ' 9. When the radius of circular array is larger (1.5 wavelength) , there' ' are more spikes, supressing the spikes by adding linear array on the Z-axis' * ' } ;ttl = 'DOA' ; %helpwin (str, ttl); text (0.65, 1.02,setstr(108),' font name' , 'symbol' , ' font size' ,20) axis ( 'off ) 73 74 APPENDIX B - SPECTRUM CALCULATION FOR CIRCULAR ARRAY format compact bearing=bearingn; msens=msensll ; msource=nsource; def_msens = 8 ; def_msource = 2; smat = zeros (nsamp, msource) ; for -n = l:msource u = rand (nsamp, 1 ) -0.5; smat(:,n)= -sign(u) .* log ( l-2*abs (u) ) ; end % < > bearing = bearing * pi / 180; % degrees to rad phi=2* pi /msens; ql=[0:msens-l] *phi; for nnn=l : msource cira (nnn, : ) =ql; end for nb=l:msens cb ( : , nb) =bearing ' ; end amat2= ( (2*pi*dspace*cos (cira-cb) ) ) ' ; for ns=l:msource a 3 ( : , ns ) =amat2 ( : , ns ) *sin (ang (ns) *pi/180) ; end amat3=exp ( j *a3) ; % the noisefree sensor signal smat2 = smat*amat3 . ' ; [nsamp, msens] = size (smat2 ) ; nsource=msource; %snr=input ( ' snr= ' ) ; namp=l/ (sqrt (2) v 10" (snr/20) ) ; for n=l:msens nmat ( : , n) =randn (nsamp, 1 ) +j *randn (nsamp, 1 ) ; end %nmat = nmat - ones (nsamp, 1 ) * mean (nmat) ; smat2=smat2 + namp'*'nmat ; ymat=smat2'; % estimate the spatial cross-cumulant matrix cmat = zeros (msens, msens) ; ymat = ymat - ones (nsamp, 1 ) * mean(ymat); % remove mean cmat = conj(ymat' * ymat) / nsamp; % correlation matrix % determine number of sources [umat, smat, vmat] = svd(cmat); svec = diag(smat); format compact % estimate bearing spectra 75 prmat = zeros (msens, msens) ; for i=l:msens if (svec(i) > 0) prmat = prmat + vmat(:,i) * vmat(:,i)' / svec(i); end end prvec = prmat (:,1) ; dtheta=l; theta = [0:dtheta:thr] ' * (pi/180); va=[0:l:90] * (pi/180) ; pp= (length ( theta ) ) ; for nnl=l:pp ciral (nnl, : ) =ql; end for nbl=l:msens cbl ( : , nbl)=theta; end omega = 2*pi*dspace*cos (ciral-cbl ) ; %<■ mth = length (theta) ; mlc = zeros (mth, 1) ; gvec = umat ( 1, 1 : nsource) .' ; gmat = umat (2 : msens, 1 : nsource) ; arvec = [ 1 ; - conj(gmat) * gvec / (1 - gvec'*gvec) ] ; tmat = vmat (:, nsource+1 rmsens) ; %wvmat = tmat * diag (ones (msens-nsource, 1 ) ./ svec (nsource+1 :msens) ) * tmat ' ; tmat = tmat * tmat ' ; nv=length (va) ; if thr==360 for bb=l:nv for k =l:mth steer = exp ( j * omega ( k, : ) ' *sin ( va (bb) ) ) ; smusl(k,bb) =1./ (steer' * tmat * steer); end end else for bb=l:nv for k =l:mth steer = exp(j * omega ( k, :)' *sin (va (bb) )) ; beaml(k,bb) = 1 . /abs (steer ' *inv( cmat)* steer); sminl(k,bb) =1./ abs (arvec. 1 * steer) ^2; smusl(k,bb) =1./ (steer' * tmat * steer); spisl(k,bb) = 1./ abs(steer' * vmat ( : , msens ) )~2; end end end 76 val=[-90: 1: 90] * (pi/180) ; mtl=length ( val ) ; display the angular spectra theta = theta*180/pi; va=va*180/pi; val=val*180/pi; top=0 . 95 ; spacing=0 . 02 ; left=0.80; btnWid=0. 1; btnHt=0.04; bottom=0.05; figure ( i ) uicontrol (' Units ',' points ' , ... ' BackgroundColor' , [0.8 0.8 0.8], ... ' FontSize' ,19, ... 'ForegroundColor' , [0.8 0.8 0.8], ... 'Position' , [178 264 100 19], ... 'String' , 'HOLD ON ! ! ' , ... 'Style' , 'text ' , ... 'Tag' , 'StaticText2 ' ) ; f a=get (0, ' current figure ' ) ; figure ( f a+1 ) , as = 10*logl0 (abs ( (smusl + eps) ' ) ) ; asl=10*logl0 (abs ( (smusl+eps) ' ) ) ; %mesh( theta' , va, 10*logl0 (abs ( (smusl+eps)'))); bb=10*logl0 (abs ( (smusl + eps) ' ) ) ; bias=min (min (bb) ) ; as=max (max (bb) ) ; if thr==360 mesh ( theta ', va, bias + 10*logl0 (abs ( (smusl+eps)'))); %title ( [ 'MUSIC, circular array, SNR= ' ,num2str (snr, 3) , ' , radius-C- spacing: ' , num2str (dspace, 3) ] ) ; zlabel ( 'power (dB) ' ) ; xlabeK 'AZ(deg) ' ) ylabel ( 'EL(deg) ' ) axis([0 90*thrv 90 -20 as]); %xlabel ( [ 'ANTENNA NUMBER= ' , int2str (msens ) ] % ylabel ( [ 'source number = ' , int2str (nsource) zlabel ( ' power ' ) break else as2=10*logl0 (abs ( as3=10*logl0 (abs ( as4=10*logl0 (abs ( subplot (221) ,mesh zlabel ( 'power (dB) ); (spisl+eps ) ' ) ) (sminl+eps) ' ) ) (beaml+eps) ' ) ) theta',va,bias+10*logl0 ) ; 77 abs ( (smusl+eps ) ) title ( 'MUSIC ) ;xlabel ( 'AZ (deg) ' ) ylabeK 'EL (deg) ' ) axis([0 90*thrv 90 -20 as]); %xlabel ( [ ' ANTENNA NUMBER= ' , int2str (msens ) ] ) ; %ylabel ( [ ' source number= ' , int2str (nsource) ] ) ; %zlabel ( 'power (dB) ' ) bb=10*logl0 (abs ( (spisl+eps) ' ) ) ; bias=min (min (bb) ) ; as=max (max (bb) ) ; subplot (222) , mesh ( theta' , va, bias+10*logl0 (abs ( (spisl + eps) * ) ) ) ; title ( 'PHD' ) ; zlabel ( 'power (dB) ' ) ; xlabel ( 'AZ (deg) ' ) ylabel ( 'EL (deg) ' ) axis([0 90*thrv 90 -20 as]); %xlabel ( 'Az-angle' ) ; bb=10*logl0 (abs ( (sminl + eps) ' ) ) ; bias=min (min (bb) ) ; as=max (max (bb) ) ; subplot (223) ,mesh( theta' , va,bias+10*logl0 (abs ( (sminl + eps) ' ) ) ) ; title ( 'Mini -norm' ) ; zlabel ( 'power (dB) ' ) ; xlabel ( 'AZ(deg) ' ) ylabel ( 'EL (deg) ' ) axis([0 90*thrv 90 -20 as]); %ylabel ( 'El-angle' ) ; bb=10*logl0 (abs ( (beaml+eps) ' ) ) ; bias=min (min (bb) ) ; %as=max (max (bb) ) ; subplot (224 ) ,mesh ( theta' , va, abs (bias) +10*logl0 (abs ( (beaml + eps) ' ) ) ) ; title ( 'Capon MLM' ) ; zlabel ( 'power (dB) ' ) ; xlabel ( 'AZ(deg) ' ) ylabel ( 'EL (deg) ' ) end AZ=bearing*130/pi; EL=ang; %title ([ 'ANTENNA NUMBER= ', int2str (msens ), ' SOURCES NUMBER=' , int2str (msource) ] ) %xlabel ( [ 'AZ=' , int2str (AZ (l:msource) ) , ' degree' ] ) % ylabel ( [ ' EL= ' , int2str (EL ( 1 : msource) ) , ' degree' ] ) tname = [' CIRCULAR ANTENNA ARRAY, SNR= ', int2str (snr) ] ; set (gcf , ' name ' , tname) 78 APPENDIX C - SPECTRUM CALCULATION FOR SPIRAL ARRAY format compact bearing=bearingn; def_msens = 8; def_msource = 2 ; msource=nsource; msens=msensll; smat = zeros (nsamp, msource) ; for n = lrmsource u = rand (nsamp, 1 ) -0.5; smat(:,n)= -sign(u) .* log ( l-2*abs (u) ) ; end % the propagation matrix bearing = bearing * pi / 180; phi =2* pi /msens ; ql= [0 :msens-l] *phi; dsp= [1 :msens] *dspace* ( 1/msens) ; %dsp=[l 2 5 11 13 16 18]*1/18; for nnn=l : msource ■ cira (nnn, : ) =ql ; dl (nnn, : ) =dsp; end for nb=l:msens cb ( : , nb ) =bearing ' ; end amat2= ( (2*pi*dl . *cos (cira-cb) ) ) ' ; for ns=l:msource a 3 ( : , ns) =amat2 ( : ,ns) *sin(ang(ns) *pi/180) ; end amat3=exp ( j *a3 ) ; smat2=smat*amat3 . ' ; [nsamp, msens] = size(smat2); nsource=msource; namp=l/ (sqrt (2) *10 A (snr/20) ) ; for n=l :msens nmat ( : , n) =randn (nsamp, 1 ) +j *randn (nsamp, 1) ; end smat2=smat2+namp*nmat ; ymat=smat2;% estimate the spatial matrix cmat = zeros (msens, msens ) ; ymat = ymat - ones (nsamp, 1 ) * mean(ymat); % remove mean cmat = conj(ymat' * ymat) / nsamp; [umat, smat, vmat] = svd(cmat); svec = diag(smat); format compact % estimate bearing spectra prmat = zeros (msens, msens) ; for i=l:msens if (svec(i) > 0) prmat = prmat + vmat ( : , i ) * vmat ( : , i ) ' / svec ( i ) ; end end prvec = prmat ( : , 1 ) ; 79 dtheta=l; theta = [0:dtheta:thr] ' * (pi/180); va=[0: 1: 90] * (pi/180) ; pp= (length (theta) ) ; for nnl=l:pp ciral (nnl, : ) =ql ; d (nnl , : ) =dsp; end for nbl=l:msens cbl ( : , nbl) =theta; end omega = 2*pi*d. *cos (ciral-cbl ) ; mth = ' length (theta) ; mlc = zeros (mth, 1 ) ; gvec = umat ( 1 , 1 : nsource) .' ; gmat = umat (2 :msens, 1 : nsource) ; arvec = [ 1 ; - conj(gmat) * gvec / (1 - gvec'*gvec] tmat = vmat (:, nsource+1 :msens) ; wvmat = tmat * diag (ones (msens-nsource, 1 ) ./ svec (nsource+1 :msens ) ) * tmat'; tmat = tmat * tmat ' ; nv=length (va' % > search angle< % for bb=l:nv for k =l:mth steer = exp(j * omega ( k, :)' *sin (va (bb) )) ; beaml(k,bb) =1./ abs (steer' *inv( cmat) * steer); sminl(k,bb) =1./ abs (arvec.' * steer) ^2; smusl(k,bb) = 1./ (steer' * tmat * steer); spisl(k,bb) = 1./ abs(steer' * vmat ( : , msens ) )~2; end end val=[-90 : 1:90]* (pi/180) ; mtl=length ( val ) ; % display the angular spectra theta = theta*180/pi; va=va*18 0/pi; val=val*180/pi; top=0 .95; spacing=0 . 02 ; left=0.80; btnWid=0. 1; btnHt=0. 04; bottom=0. 05; figure ( 1 ) uicontrol (' Units ', 'points ' , ... 'BackgroundColor' , [0.8 0.8 0.8], •FontSize' ,19, ... ' ForegroundColor' , [0.8 0.8 0.8], .. 'Position', [178 264 100 19], ... 'String' , 'HOLD ON ! ! ' , ... 80 'Style' , 'text ' , ... •Tag' , 'StaticText2' ) ; fa=get(0, ' currentf igure ' ) figure (fa+1) , tname = ['SPIRAL CIRCULAR ANTENNA set(gcf, ' name ', tname) asl=10*loglO (abs ( (smusl+eps ) ' ) ) as2=10 + logl0 (abs ( (spisl+eps ) ' ) ) as3=10*logl0 (abs ( (sminl+eps ) ' ) ) as4=10*loglO (abs ( (beaml+eps ) ' ) ) ARRAY, SNR=' , int2str (snr; bb=10 + loglO (abs ( (smusl+eps! bias=min (min (bb) ) ; as=max (max (bb) ) ; if thr==360 mesh ( theta ', va, bias+10*loglO (abs ( (smusl+eps)'))); %title ([ 'MUSIC, spiral circular array, SNR= ' , num2str (snr, 3! spacing: ' , num2str (dspace, 3) ] ) ; axis([0 90*thrv 90 -20 as]); zlabel ( 'power (dB) ' ) ; xlabel ( 'AZ(deg) ' ) ylabel ( 'EL(deg) ' ) , max- else subplot (221) ,mesh( theta' , va, bias+10*logl0 (abs %title ( 'music ' ) ; axis([0 90*thrv 90 -20 as]); smusl+eps) ' ) ) ) ; bb=10*logl0 (abs ( (spisl + eps) ' ) ) ; bias=min (min (bb) ) ; as=max (max (bb) ) ; zlabel ( 'power (dB) ' ) ; title ( 'MUSIC ) ; xlabel ( 'AZ(deg) ' ) ylabel ( 'EL(deg) ' ) subplot (222 ), mesh ( theta ', va, bias + lO^loglO (abs ( (spisl + eps! title ( 'PHD* ) ; zlabel ( 'power (dB) ' ) ;xlabel ( 'AZ (deg) ' ) ;ylabel ( 'EL (deg) ' ) axis([0 90 + thrv 90 -20 as]); bb=10*logl0 (abs ( (sminl + eps) ' ) ) ; bias=min (min (bb) ) ; as=max (max (bb) ) ; %xlabel ( 'Az -angle' ) ; subplot (223 ), mesh ( theta ', va, bias + 10*logl0 (abs ( (sminl+eps! title ( 'Mini-norm' ) ; zlabel ( 'power (dB) ' ) ; title ( 'MUSIC ) ; xlabel ( 'AZ(deg) ' ) ylabel ( 'EL(deg) ' ) axis([0 90*thrv 90 -20 as]); %xlabel ([' uniform circular array, SNR=' , int2str (snr) ]) ; %ylabel ( 'El-angle' ) ; bb=10*logl0 (abs ( (beaml + eps) ' ) ) ; bias=min (min (bb) ) ; 81 subplot (224 ) , mesh ( theta' , va, abs (bias ) +10* log 10 (abs ( (beaml + eps) ' title ( 'Capon MLM' ) ; zlabel ( 'power (dB) ' ) ; xlabel ( 'AZ(deg) ' ) ylabel ( 'EL(deg) ' ) end AZ=bearing*180/pi; EL=ang; °6title ([ 'SPIRAL ANTENNA 'NUMBER= ', int2str (msens ), ' SOURCES NUMBER=' , int2str (msource) ] ) 82 APPENDIX D - SPECTRUM CALCULATION FOR CIRCULAR+LINEAR ARRAY dspacel=str2num (get (bdspacel, 'string' ) ) lm=str2num (get (bell, 'string ' ) ) bearing=bearingn; msource=2; if nsource>6 msource=6; else msource=nsource; end smat = zeros (nsamp, msource) ; for n = l:msource u = rand (nsamp, 1) -0.5; smat(:,n)= -sign(u) .* log ( l-2 + abs (u) ) ; end % the propagation matrix bearing = bearing * pi / 180; % degrees to radian msensm=msensll ; phi=2*pi/msensm; ql= [0 :msensm-l] *phi; cira=ones (msource, 1 ) *ql ; cb=bearing ' *ones ( 1 , msensm) ; amat2= ( (2*pi*dspace*cos (cira-cb) ) ) ' ; amat22=( (2*pi* [0:lm-l] '*dspacel) ) ; for ns=l:msource a 3 ( : , ns) =amat2 ( : , ns) *sin (ang (ns) * pi/ 180) ; a 32 ( : , ns) =amat22*cos (ang (ns) *pi/180) ; end a3=[a3;a32] ; amat3=exp ( j *a3) ; % the noisefree sensor signal smat2=smat*amat3 . ' ; % add white noise : [nsamp, msens] = size(smat2); nsource=msource ; namp=l/ (sqrt (2) *10 A (snr/20) ) ; for n=l:msens nmat ( : , n) =randn (nsamp, 1 ) + j *radn (nsamp, 1 ) ; end nmat = nmat - ones (nsamp, 1 ) * mean(nmat); smat2=smat2+namp*nmat ; ymat=smat2; estimate the spatial cross-cumulant matrix 83 cmat = zeros (msens, msens ) ; ymat = ymat - ones (nsamp, 1 ) * mean (ymat ) ; % remove mean cmat = conj(ymat' * ymat) / nsamp; % correlation matrix % determine number of sources [umat, smat, vmat] = svd(cmat); svec = diag(smat); s=[l 2 3] ; dtheta=l; format compact estimate bearing spectra prmat = zeros (msens, msens) ; for i=l:msens if (svec(i) > 0) prmat = prmat + vmat(:,i) * vmat ( : , i ) ' / svec(i); end end prvec = prmat ( : , 1 ) ; theta = [0:dtheta:thr] ' * (pi/180); va=[0:l: 90] * (pi/180) ; pp= (length (theta) ) ; ciral=ones (pp, 1) *ql; cbl=theta*ones (l,msensm) ; omega = 2*pi*dspace*cos (ciral-cbl ) ; %< @@@@@*^ mth = length (theta) ; mlc = zeros (mth, 1 ) ; gvec = umat ( 1, 1 : nsource) .' ; gmat = umat (2 : msens, 1 : nsource) ; arvec = [1 ; - con j (gmat) * gvec / (1 - gvec'*gvec) ]; tmat = vmat ( : , nsource+1 :msens) ; wvmat = tmat * diag (ones (msens-nsource, 1 ) ./ svec (nsource + 1 :msens ) ) tmat ' ; tmat = tmat * tmat ' ; nv=length ( va) ; if thr==360 for bb=l:nv for k =l:mth steerl = exp(j * omega ( k, :)' *sin (va (bb) )) ; steer2= exp(j *cos (va (bb) ) *2*pi*dspacel* [0:lm-l]'); steer= [steerl; steer2] ; smusl(k,bb) = 1./ (steer' * tmat * steer); end end else for bb=l : nv for k =l:mth steerl = exp(j * omega (k, :)' *sin (va (bb) )) .; steer2= exp(j *cos (va (bb) ) *2*pi*dspacel* [0:lm-l]'); steer= [steerl; steer2] ; 84 beaml(k,bb) = 1 . /abs (steer ' *inv( cmat)* steer); sminl ( k, bb) = 1./ abs (arvec. ' * steer) ^2; smusl(k,bb) = 1./ (steer' * tmat * steer); spisl(k,bb) = 1 . / abs (steer' * vmat ( : , msens ) ) A 2, end end end val=[-90: 1: 90] * (pi/180) ; mtl=length ( val ) ; % display the angular spectra theta = theta*180/pi; va=va*180/pi; val=val*180/pi; top=0 . 95; spacing=0 . 02 ; left=0.80; btnWid=0. 1; btnHt=0.04; bottom=0. 05; figure ( 1 ) uicontrol (' Units ',' points ' , ... •BackgroundColor' , [0.8 0.8 0. ' FontSize' , 19, . . . ' ForegroundColor' , [0.8 0.8 0.8], 'Position' , [178 264 100 19], 'String' , 'HOLD ON ! ! ' , ... 'Style' , 'text ' , ... 'Tag' , 'StaticText2 ' ) ; f a=get ( 0, ' current figure ' figure ( f a + 1 ) , asl = 10*logl0 (abs ( (smusl+eps) * ) ) ; bb=10*logl0 (abs ( (smusl + eps) ' ) ) ; bias=min (min (bb) ) ; as=max (max (bb) ) ; if thr==360 mesh ( theta ', va, bias+10*logl0 (abs ( (smusl+eps)'))); % title ( [ 'MUSIC, linear + circular array, SNR= ' ,num2str (snr, 3) , ' , C- spacing : ' , num2str (dspace, 3) , ' , L- space : ' , num2str (dspacel, 3 ) , ' , L-ant no : ' , num2str (lm, 3) ] ) ; axis([0 90*thrv 90 -20 as]); %xlabel ( [ 'ANTENNA NUMBER=' , int2str (msens ) ] ) ; %ylabel ([' source number = ' , int2str (nsource) ] ) ; %zlabel( 'power (dB) ' ) ; zlabel ( ' power (dB) ' ) ; xlabel ( 'AZ(deg) ' ) 85 ylabel ( 'EL(deg: else as2=10*logl0 (abs ( (spisl+eps) ' ) ) as3=10*logl0 (abs ( (sminl + eps) ' ) ) as4=10*logl0 (abs ( (beaml+eps) ' ) ) subplot (221) ,mesh ( theta' , va, bias+10*logl0 (abs ( (smusl + eps )• ' ) ) ) ; axis([0 90*thrv 90 -20 as]); zlabel ( ' power (dB) ' ) ; title ( 'MUSIC ) ;xlabel ( 'AZ (deg) ' ) ylabel ( 'EL (deg) ' ) bb=10*logl0 (abs ( (spisl+eps) ' ) ) ; bias=min (min (bb) ) ; as=max (max (bb) ) ; subplot (222) , mesh ( theta' , va, bias+10*logl0 (abs ( (spisl+eps) ' ) ) ) ; zlabel ( 'power (dB) * ) ; title ( ' PHD' ) ;xlabel ( 'AZ (deg) ' ) ylabel ( 'EL (deg) ' ) axis([0 90*thrv 90 -20 as]); bb=10*logl0 (abs ( (sminl+eps) ' ) ) ; bias=min (min (bb) ) ; as=max (max (bb) ) ; subplot (223) , mesh ( theta ' , va, bias+lO^loglO (abs ( (sminl + eps) ' ) ) ) ; title ( 'Mini-norm' ) ; zlabel ( 'power (dB) ' ) ; xlabel ( 'AZ(deg) ' ) ylabel ( 'EL (deg) ' ) axis([0 90*thrv 90 -20 as]); tname= ( [ ' linear+circular, SNR=' , int2str (snr) ] ) ; lylabel ( 'El (degree) ' ) ; bb=10*logl0 (abs ( (beaml+eps) ' ) ) ; bias=min (min (bb) ) ; %as=max (max (bb) ) ; subplot (22 4 ) , mesh ( theta' , va, abs (bias)+10*logl0 (abs ( (beaml+eps) ' title ( 'Capon MLM' ) ; zlabel ( 'power (dB) ' ) ; xlabel ( 'AZ(deg) ' ) ylabel ( 'EL(deg) ' ) end AZ=bearing*180/pi; EL=ang; tname= ( [ ' linear+circular, SNR=' , int2str (snr ) ] ) ; set(gcf, ' name ' , tname) 86 APPENDIX E - SPECTRUM CALCULATION FOR LINEAR ARRAY clear nmat [nsamp, msens] = size (smatl ) ; nsource=msource; namp=l/ (sqrt (2) *10 A (snr/20) ) ; for n=l: msens nmat (:,n) = randn (nsamp, 1 ) +j *randn (nsamp, 1 ) ; end smatl=smatl+namp*nmat ; ymat=smatl ; % estimate the spatial cross-cumulant matrix cmat = zeros (msens, msens ) ; ymat = ymat - ones (nsamp, 1 } * mean (ymat); % remove mean % determine number of sources cmat = conj(ymat' * ymat) / nsamp; % correlation matrix determine number of sources [umat, smat, vmat] = svd(cmat); svec = diag(smat); estimate bearing spectra prmat = zeros (msens, msens) ; for i=l rmsens if (svec(i) > 0) prmat = prmat + vmat(:,i) * vmat(:,i)' / svec(i); end end prvec = prmat (:,!); dtheta=0.5; theta = [-90:dtheta:90] ' * (pi/180); omega = 2*pi*dspace*sin (theta) ; mth = length (theta) ; mlc = zeros (mth, 1 ) ; smus = mlc; spis = mlc; seig = mlc; par mlc; beam = mlc; smin = mlc; gvec = umat ( 1, 1 : nsource) . ' ; gmat = umat (2 : msens, 1 : nsource) ; arvec = [1 ; - con j (gmat) * gvec / (1 - gvec'*gvec) ]; tmat = vmat ( : , nsource+1 :msens) ; 87 wvmat = tmat * diag (ones (msens-nsource, 1 ) ./ svec (nsource+1 :msens; tmat ' ; tmat = tmat * tmat'; for k =1 :mth steer = exp(j * omega(k) * [0 :msens-l ] ' ) ; °obeam(k) = abs (steer' * cmat * steer); beam(k) = 1 . /abs (steer ' * inv(cmat) * steer); smin(k) =1./ abs(arvec.' * steer) "2; smus(k) = 1./ (steer' * tmat * steer); spis(k) = 1./ abs(steer' * vmat ( : , msens) ) A 2; end % display the angular spectra theta = theta*180/pi; spec = abs ( [seig, smus, spis, mlc, par, smin, beam] ); % normalize to abs-max of unity for display only spmax = max ( spec ); spmax = ones (1,7) ./ spmax; spec = spec * diag( spmax); figure (2 ) Isubplot ( 421 ) , stem(svec), title (tname) subplot (421) , plot (theta, 10*logl0 (spec ( : , 3) ) , ' r ' ) title (' Uniform linear array') ylabel ( 'power (dB) ' ) gtext ( 'PHD' ) subplot (423) , plot (theta, 10*logl0 (spec ( : , 2 ) ) , ' r ' ) ylabel ( 'power (dB) ' ) , gtext ( 'MUSIC ) subplot (425) , plot (theta, 10*logl0 (spec ( : , 6) ) , ' r ' ) ylabel ( ' power (dB) ' ) , gtext ( ' Mini-norm' ) subplot (427) , plot (theta, 10*logl0 (spec ( : , 7 ) ) , ' r ' ) , hold ylabel ( 'power (dB) ' ) , gtext ( 'Capon MLM' ) , xlabel ( 'AOA(degree) ' ) tname = ['UNIFORM ANTENNA and NON-UNIFORM ANTENNA LINEAR ARRAY'] set(gcf, ' name ', tname) 88 APPENDIX F - SPECTRUM CALCULATION FOR NON-UNIFORM LINEAR ARRAY clear nmat [nsamp, msens] = size (smat2 ) ; nsource=msource ; %snr=input ( ' snr= ' ) ; namp=l/ (sqrt (2) *10~ (snr/20) ) ; for n=l:msens nmat ( : , n) =randn (nsamp, 1 ) + j *randn (nsamp, 1 ) ; end smat2=smat2+namp*nmat ; ymat=smat2; % estimate the spatial cross-cumulant matrix cmat = zeros (msens, msens ) ; ymat = ymat - ones (nsamp, 1 ) * mean(ymat); % remove mean cmat = conj(ymat' * ymat) / nsamp; % correlation matrix determine number of sources [umat, smat, vmat] = svd(cmat); svec = diag(smat); hold off % estimate bearing spectra prmat = zeros (msens, msens) ; for i=l:msens if (svec(i) > 0) prmat = prmat + vmat(:,i) * vmat(:,i)' / svec(i); end end prvec = prmat ( : , 1 ) ; dtheta=0. 5; theta = [-90:dtheta:90] ' * (pi/180); omega = 2*pi*dspace*sin (theta) ; mth = length (theta) ; mlc = zeros (mth, 1 ) ; smus = mlc; spis = mlc;seig = mlc; par = mlc; beam = mlc; smin = mlc; gvec = umat ( 1, 1 : nsource) . ' ; gmat = umat (2 : msens, 1 : nsource) ; ar.vec = [1 ; - conj (gmat) * gvec / (1 - gvec'*gvec) ]; tmat = vmat (:, nsource+1 : msens ) ; wvmat = tmat * diag (ones (msens-nsource, 1) ./ svec (nsource+1 :msens) ) * tmat ' ; tmat = tmat * tmat ' ; for k =l:mth 89 steer = exp (sqrt (-1) * omega (k) * asp'); obeam(k) = abs (steer' * cmat * steer); beam(k) =1 . / abs (steer 1 *inv( cmat) * steer); smin(k) = 1./ abs(arvec. ' * steer) A 2; smus(k) = 1./ (steer' * tmat * steer); spis(k) = 1./ abs(steer' * vmat ( : , msens) ) A 2; end display the angular spectra theta = theta*180/pi; spec = abs ( [smus, spis, smin, beam] ); % normalize to abs-max of unity for display only spmax = max (spec); spmax = ones (1,4) ./ spmax; spec = spec * diag (spmax); figure (2 ) subplot (422) , plot (theta, 10*logl0 (spec ( : , 2 ) ) ) ; %set (gca, 'position' , [0.06 0.75 0.35 0.17]); title ( 'Non-uniform linear array' ) , y label ( 'power (dB) ' ) , gtext ( ' PHD' ) ; subplot (424) , plot (theta, 10*logl0 (spec ( : , 1) ) ) , ylabel ( 'power (dB) * ) , gtext ( 'MUSIC ) ; subplot (426) , plot (theta, 10*logl0 (spec ( : , 3) ) ) ; ylabel ( 'power (dB) ' ) , gtext ( 'Mini-norm' ) %set (gca, 'position' , [0.06 0.28 0.35 0.17]); %title (' Minimum norm method') subplot (428) , plot (theta, 10*logl0 (spec ( : , 4 ) ) ) ; ylabel ( 'power (dB) ' ) , xlabel ( 'AOA(degree) ' ) , gtext ( 'Capon MLM' ) 90 LIST OF REFERENCES [1] James, T. , Digital techniques for wideband receivers, Artech House, 1995. [2] Schmidt, R. 0., "Multiple emitter location and signal parameter," IEEE Trans, on Antennas and Propagation, vol. AP-34, no. 3, March, 1981. [3] Pillai, S. U. , Array signal processing, New York, Springer-Verlag, 1988. [4] Yaakov, B., and Li, X., Estimation and Tracking: Principle, Techniques , Artech House, 1995. [5] Alan, T. M., "Minimum redundancy linear arrays," IEEE Trans. Antennas and Propagation, vol. AP-16, no. 2, March 1968. [6] Hamid, K. , and Mats, V., "Two Decades of Array Signal Processing Research, " IEEE signal processing magazine, July, 1996. [7] Pisarenko, V. F., "The retrieval of harmonics from a covariance function, " Geophysical Journal of Royal Astronomical Society, no. 33, pp. 374-366, 1973. [8] Monson, H. H., Statistical digital signal processing, John Wiley and Sons, 1980. [9] Balanis, C. A., Antenna Theory, John Wiley and Sons, 1977. 91 [10] Johnson, R. L . , and Miner, G. E., "Comparison of Superresolution Algorithms for Radio Direction Finding, " IEEE Trans. Antennas and Propagation, vol. AES-22, no. 4, July 1986. [11] Roy, R., and Kailath, T., "ESPRIT - Estimation of Signal Parameters Via Rotational Invariance Techniques," IEEE Trans, on Acoustics, Speech, and Signal Processing, vol. 37, no. 7, July 1989. 92 INITIAL DISTRIBUTION LIST Defense Technical Information Center 2 8725 John J. Kingman Rd . , STE 0944 Ft. Belvoir, Virginia 22060-6218 Dudley Knox Library 2 Naval Postgraduate School 411 Dyer Rd . Monterey, California 93943-5101 Chairman, Code EC 2 Department of Electrical and Computer Engineering Naval Postgraduate School Monterey, California 93943-5121 Prof. David C. Jenn, Code EC/Jn 2 Department of Electrical and Computer Engineering Naval Postgraduate School Monterey, California 93943-5121 Prof. Phillip A. Pace, Code EC/Pc 1 Department of Electrical and Computer Engineering Naval Postgraduate School Monterey, California 93943-5121 Chung-Shan Institute of Science Technology 2 P.O.Box#90008-6-14 , Lung- Tang, Tao-Yuan, Taiwan, R.O.C. Ku-Ting, Lin 2 Taipei, Taiwan, R.O.C, IF, 15 , Ln2 ,Hsin-Chung St. 93 6 483NPG TH 10/99 22527-200 im '<A*.