(navigation image)
Home American Libraries | Canadian Libraries | Universal Library | Community Texts | Project Gutenberg | Children's Library | Biodiversity Heritage Library | Additional Collections
Search: Advanced Search
Anonymous User (login or join us)
Upload
See other formats

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

•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*.