- Email: [email protected]
Minimum bias spatial filters for beamspace direction-of-arrival estimation
J. Franklin Inst. Vol. 335B, No. 1, pp. 35-52, 1998
~ Pergamon
PII:
S0016--0032(96)00114--7
Copyright © 1997The Franklin Institute Published by ElsevierScienceLtd Printed in Great Britain 0016-0032/98 $19.00+0.00
Minimum Bias Spatial Filtersfor Beamspace Direction-of-arrival Estimation by MOENESS
G. AMINt
and N I T I N
BHALLA
Department o f Electrical and Computer Engineering, Villanova University, Villanova, PA 19085, U.S.A.
ABSTRACT: Spatial filtering is used in array signal processing to favor signals coming from a particular zone in space over other signals. In this paper, we propose the use of the minimum bias (MB) windows, originally devised for multitaper spectral analysis, as spatial preproeessing filters for direction-of-arrival (DOA) estimation. These windows (or tapers) are characterized by their orthonormality, simplicity, and most importantly their close behavior to the well known discrete prolate spheroidal sequence (DPSS) windows. A desired spatial bandwidth can be achieved by adding or removing tapers, without the need of reeomputation. The paper introduces the MB tapers for beamspace processing and compares their DOA estimation performance, using noise-subspace eignestructure methods, to that of the DPSS beams. Further, the MB windows can be well approximated by simple trigonometric functions. These functions have expressions identical to that of the Discrete Sinusoidal Transform (DST). As numerous methods f or fast computations of DST have been devised, the outputs of the approximated MB tapers can be computed efficiently. © 1997 The Franklin Institute. Published by Elsevier Science Ltd
L Introduction
In array signal processing, it is often desired that signals coming from a particular zone be removed, or at least attenuated. These signals may represent intentional jammers or multipath arrivals. A common example is a sensor array mounted on an aircraft which receives signals from a sonobuoy. The scatters of the sonobuoy signal from the aircraft propellers cause a severe multipath effect on the array performance, and therefore, their effect should be mitigated. Due to the coherent signal environment, the adaptive spatial nulling techniques (high resolution DOA spectra) not only fail to form nulls (peaks) in the directions of the interferers and their multipaths, but also tend to cancel the desired signal in the array output. This failure occurs even with the decorrelation effects introduced by the motion of the array on the aircraft (1).
t This work is supported by ONR grant N00014-94-1-1052. 35
36
M. G. Amin and N. Bhalla
Preprocessing techniques such as subarray averaging methods (2) can be employed to mitigate this problem and decorrelate the coherent arrivals prior to both hulling and DOA estimation. These methods, however, reduce the array aperture and are appropriate for point sources. They do not perform well for small arrays, few snapshots and/or under models of distributed sources. An alternative is to use spatial filtering to filter out both uncorrelated and coherent unwanted signals in the spatial domain, which is also known as beamspace processing. This is achieved by lowering the array response over specific spatial sectors which include the jammer arrivals. Accordingly, in a multipath environment, the coherent counterparts of direct path are suppressed instead of being decoupled. The response over the desired sector should be uniform and of unit gain. Beamspace methods have been employed for a variety of applications in the literature. The obvious advantage is the reduction in both the interference and noise whose spatial spectrum lies outside the spatial sector of interest, i.e. increasing the signal-tointerference and noise-ratio (SINR). They also lead to reduced computational complexity, amenability to parallel implementation, low sensitivity to sensor perturbations and deviations from the assumed noise model (3-5), and improving the spatial resolution (6, 7). Multiple orthogonal beams can be generated using spatial filters with different characteristics, each receiving the interference over its stopband. To capture the desired signals, the combined beams should uniformly cover the spatial sector of interest. When properly designed, the outputs of the spatial filters are jammer free and can be used as input for eigenstructure-based methods to perform high-resolution DOA estimation of the desired signals. A popular and widely used set of spatial tapers (filters) is the family of discrete prolate spheroidal sequence (DPSS) (8). Their basic properties were given by Slepian (9). These tapers have much lower sidelobes than rectangular windows and, when used as spatial filters, they yield better performance than DFT beamformers. Recently, a new family of orthonormal basis have been introduced for minimum bias estimates of the power density function in the context of time-series analysis (10). In this paper, we apply the minimum bias (MB) tapers for beamspace processing. The MB tapers can be well approximated by trigonometric functions, and as such, have a much simpler form than the DPSS tapers. A desired frequency bandwidth can be achieved by adding or removing tapers, without the need for recomputations, which is the case with DPSS. The approximated MB tapers, called the sinusoidal tapers, have a closed form analytical expression, which is identical to the sinusoidal basis functions of the Discrete Sinusoidal Transform (DST). Since numerous fast algorithms for DST evaluation exist, efficient computation of outputs corresponding to sinusoidal tapering beams is possible. The minimum bias tapers are discussed in detail in Section II. This section also provides the expressions for the approximated versions of the MB tapers, the sinusoidal tapers. In Section III, we explain the use of the MB tapers for spatial signal processing. The properties of sinusoidal tapers are discussed in detail in Section IV. Section V deals with the important aspect of fast computation of the outputs of the different beams utilizing the sinusoidal tapers. The similarity between the expressions for sinusoidal tapers and DST offers the opportunity to achieve fast computations. Absence of such a feature in DPSS makes the sinusoidal tapers a favorable proposition. The fast DST
M i n i m u m Bias Spatial Filters
37
computations of the data in beamspace are compared with the conventional procedure of inner products, i.e. multiplying the data snapshots with the tapers followed by summing. In Section VI, a comparison of DOA estimation efficacy of MB, sinusoidal, DFT, and DPSS tapers, using subspace methods, is provided to elucidate the behavior of the proposed set of tapers. 11. Minimum Bias Tapers In time-series analysis, the spectrum estimate bias is defined as E[5~(/) - S(f)] = ff~o ([ V ( 9 - f ) I2S O ) - S ( f ) ) d9
(1)
where, S ( f ) is the true spectral density whose estimate 5~(f) is given by S(f) =
~ v(n)x(n) e j2~znl ~
(2)
n=]
The observations {x(n), n = 1,2 . . . . . N} form a time series and V ( f ) is the Fourier transform of the data window (taper) v(n). It is shown in (10) that the minimum value of the bias in (1) is obtained when v(n) is such that S ~ ] V ( f ) ] z f 2 dfis minimized. For discrete-time multiple window spectral analysis, the minimum bias tapers are given by the eigenvectors v°), v~2). . . . , v~ , corresponding to the eigenvalues, sorted in the increasing order, of the matrix C = [Cm,]: l/12
c ..... =
f
(-1)
n= m n-m
2rr 2(n -- m) 2
(3)
nCm
These tapers form an orthonormal set. The taper v~i)provides the minimum possible spectral bias subject to v°3 _1_vv'), j = 1, 2, . . . , i, where 3_ denotes orthogonality. We note that in spatial filtering applications, N represents the number of array sensors rather than the number of observations. In practical situations involving beamspace processing, more than one taper should be used to devise the passband of the spatial filters, for the following reasons: (1) although some a priori knowledge is assumed to be available about the direction of arrival of the desired signals, this information is usually not complete. Using too few tapers may not sufficiently cover the spatial zone within which the desired sources are located; (2) eigenstructure methods for DOA estimation require a minimum array aperture size in accordance with the number of signals whose parameters are to be estimated. These methods are applicable only if the array aperture is of size one more than the number of arrivals. As a result of processing the data by spatial filters, the effective array aperture is reduced. The beamspace array aperture becomes equal to the number of tapers used, which is invariably less than the element space aperture (maximum number of tapers that can be used is equal to the number of sensors, which amounts
M. G. Amin and N. Bhalla
38
to an allpass spatial filter for MB tapers). Hence, if the parameters of P signals are to be estimated by beamspace eigenstructure methods, at least ( P + 1) tapers should be used. The MB tapers can be well approximated by the orthonormal sinusoidal functions, u"), u (z). . . . . u (m (10). The elements of u (k), (u~~), u(2k).... , u~)) T are defined as u~ =
/~
.
~k,~
k = 1,2 . . . . . N
4 N~i- sin Nq_-1
n= 1,2,...,N
(4)
The temporal frequency f is being utilized to elucidate the concepts henceforth. Nevertheless, the results obtained can be directly applied to array processing by appropriately substituting the temporal frequency by the spatial frequency. Spatial response is calculated using directional vectors which have been normalized by the factor x/l/N. Below, we discuss the orthonormal properties of sinusoidal tapers. The Fourier transform of the kth taper for a linear array of N sensors is given by (lO), N
V~k~(f) = ~ v(k)e j2~nl n=l
e --jn[(N+ 1)f-- (k/2)]
j x~(( N + 1)
sin [z~ ( f - 2(NL 1))]
s i n I N n ( f + k2(N+
--(--l)k [~(f+ sin
1))]
t
k
(5)
2(N+ 1))1
sin N ~ I V(k](f) = e
/~[(N+J)f-((k J).'2]]
,jSU l)
(6) sin2 (nf) - sin2 (~2(~+ 1)-)
Equation (5) clearly shows that the frequency response of every taper has two components which have a maximum each at k/2(N+ 1) and - k / 2 ( N + 1), respectively, and thus they are symmetrically placed about the origin. The maximum value (at each of these frequencies) of the kth taper is ~ / 2 . The sidelobes diminish in value as we move away from the maxima. The sidelobes of the two terms are of opposite polarity and tend to cancel each other effectively, resulting in much lower sidelobes for the taper as a whole. The frequency response of the kth taper thus has two peaks, equidistant from the origin. The peaks are slightly shifted from the positions Ik/2(N+ 1)[ because the sidelobes of one term interact with the maximum of the other. This shift is more
Minimum Bias Spatial Filters
39
pronounced for small values of k, where the peaks of the two components are close to each other and hence sidelobe contribution is higher. Orthogonality of tapers is a desired property because it eliminates the redundancy and simplifies computations (6). To prove orthonormality of the sinusoidal tapers, consider the inner product bk,,k, of the two tapers corresponding to kl and k2: (2)~=.
7tkln.
rck2n
1 :
N ~
--
E ( ej[~k'n/(N+l)l-e/[r~(k[n)/(N+l)l)(e/[~k2n/IN+l)l-e
/[~k2,;N+I)I) (7)
n=l Using the results of Appendix A, the above expression simplifies to:
(1)~eJn(kl+ke)--eJ[n(kl+k2)/(N+l)]eln(k'k2)--e b~-,'k~ - --
2(N+ 1).
(
.
. . . l - -. e j[n(k)+k2)/(N+l)]
e Jn(kt - k z ) --
1--e
e -jpt(kl -k2)/(N+ 1)]
1 - - e j[.(k, -k2)/(m+ 1)]
jtn(kl-k2)/(N+ l)]
][Tt(kI-k2)/(N+l)]
e - ~ ( k ~ +k 2) _ e-J[Tt(kl +k2)/(N+ 1 )] "~
+
1 --e -j[~(k~+k2)/(N+ 1)]
I
(8)
In Appendix A, it is shown that for the two cases of (k~ q-k2) being odd or even, the above expression sums to zero, implying orthogonality. For k~ = k2, the sum is unity, indicating normality. IlL M i n i m u m Bias Spatial Filters Let X be the data matrix X = AS + ~/
(9)
where A is the steering matrix associated with the arrivals of both signals and jammers, and S is the signal matrix whose columns are the snapshots of the directional array inputs collected at different time instants. The noise is assumed to be spatially white, i.e. the components of the nondirectional noise vector ~/are uncorrelated. The spatial tapers, whether they represent the DPSS, MB, DFT, or any other beamformer, can be used to define the spatial filter matrix T. Accordingly, the beamformers' outputs are given by Y = THX
(10)
Figure 1 illustrates the transition from element space to beamspace. (For the minimum bias case, the number of these tapers is easily determined as a function of the desired spatial bandwidth.) If the spatial sector of interest is not centered at 0 = O, but rather at 00, then one can steer the array to the desired angular position by modifying the spatial matrix as follows: F = DT
(11)
where F is the spatial filter matrix steered to desired angle, and D = [din,,],
M. G. Amin and N. Bhalla
40
Beam #1
Beam #2
Beam #3
/
Element space (N) Beam space (K)
FIG. I. Elementspace to beamspace.
din,, = 6m, exp [j2rc sin(Oo)dm/2) is the diagonal matrix corresponding to the steering angle 0o. The intersensor spacing of the array is d and 2 is the wavelength of operation. The spatial frequency corresponding to the angle 0o is (d/2) sin 00. The beamspace snapshots (data) constituting the columns of Y can then be used for DOA via subspace methods. The DOA are the peak locations of the spectrum:
S(O) =
1
[a~(O)V.nV~ai(O)l
(12)
Minimum Bias Spatial Filters
41
where, as(0) defines the beamspace steering vector (aj(0) = Tna(0), a(0), being the element space steering vector) and V, is the noise subspace projection matrix whose columns are the eigenvectors corresponding to the minimum eigenvalues of the correlation matrix formed by the filtered data. Matrix H is a suitably selected nonnegative square matrix of size ( K - L), where L is the number of arrivals whose parameters are to be estimated (6). In the case of beamspace MUSIC, H = I, while for beamspace MinNorm, H = ssH/(sHs) 2, where s denotes the ( K - L) × 1 vector Vff(1,0 . . . . . O)r.
IV. Properties of the Sinusoidal Tapers 4.1. General properties In general, a beamformer should have the following two properties, for 0 in the passband (6): (a) it does not introduce superfluous zeros in the passband, i.e. aH(O)TTHa(O) > 0
(13)
This property is important because, if the beamformer itself has nulls in the spatial sector of interest, it will lead to spurious peaks in the spectrum at those locations and also, remove any arrivals at those points, resulting in erratic inferences; (b) the first order differences in the element space directional vectors are translated into first order differences in the beamspace directional vectors, implying Tn~0a(0) =# cTHa(O)
(c = scalar)
(14)
This property ensures that any change in the steering vectors in the element space is not smoothed out due to beamspace processing, but is appropriately reflected in the beamspace steering vectors. In Appendix B, it is shown that these two above mentioned properties hold for the sinusoidal tapers. 4.2. Specific properties In Section II, we briefly discussed the sidelobe and orthonormality properties of the sinusoidal tapers. Below, we elaborate on the overall frequency characteristics of these tapers. The analysis of the combined response of say, K tapers (K < N) will be of importance. From (6), it is evident that, even numbered tapers vanish at f = 0 as the numerator evaluates to zero, implying that they do not contribute at the origin. The derivative of squared response of absolute value of kth taper is evaluated by differentiating (6) w.r.t.f. Summing the above derivative for N = 1, 2 . . . . . K gives the corresponding expression for sum of K tapers: ~r
K
(-
~k 1)k sin2 N + ~
~sin [2(N+ 1)Trf] k= ' {sin2 (Trf) - sin2 ( 2 ~ +
~)}2
42
M. G. Amin and N. Bhalla
rc
I
zk]. 2
K sin2 ( N + l ) r c f - ~ - J s a n N+lZCk
,~r± 1', sin 2xf ~ \~'Tx]
k~l
(15) sin2(ltf)-sin2
2(~
For any K, a maxima or minima exists a t f = 0 as the derivative is zero. It is shown in Appendix C that odd values of K will result in a maximum whereas even values of K provide a minimum spatial response at the origin. Because even tapers do not contribute at the origin, as we increment K, going from an odd to an even value, the response a t f = 0 does not change and in fact, maxima of an odd taper and minima of the next taper, which is even, coincide. 4.2.1. Error band of the spatial filter. Filter design problems have an important parameter--the passband ripple. For the spatial filter designed by sinusoidal tapers, we attempt to derive a similar parameter, the error band. This parameter gives the deviation between the maximum and the minimum points in the filter passband. An obvious method of deriving an expression for the error band is: (a) find the derivative of the sum of squares of magnitudes of frequency responses of tapers; (b) determine the points of maxima and minima in the passband; (c) the difference between the lowest minimum and the highest maximum points gives the error band. Equation (15) gives the expression for the desired derivative. This equation can be simplified by assuming N t o be sufficiently large such that rck/(N+ 1) is small enough for reasonable approximations in the numerator. This allows us to make approximations of the kind sin x ~ x for small x, but of course with inherent approximation errors. The magnitude of this error reduces as N increases. The complexity of the problem requires the use of numerical methods to find its roots. In Fig. 2, a zoomed view of the passband corresponding to 32 sensors and 8 tapers is shown. The vertical grid lines are at frequenciesfp = p/2(N+ 1) where p is an integer (these are the frequencies where individual tapers have the peaks of their components). Interestingly, a close observation of the frequency response of the tapers in Fig. 2 and the values of extrema obtained by MATLAB routines, reveals the following: (1) the extrema lie very close to the frequenciesfp; (2) the lowest minimum in the passband occurs near p = K - 2 ; (3) the highest maxima for even K lies nearfp= l, and for K odd a t f = 0. This pattern motivates to study the behavior of frequency response at frequenciesfp for different K. In Appendix D, it is shown that the response of sums of squares of the absolute values of the first K tapers at fp (k = 1,2 . . . . . K), after employing relevant approximations, is:
Minimum Bias Spatial Filters I I I I I I
~= e
0.995
C/~
~
I I I I
I
0.99
f,\
I
~ 0.985 "0
. -
,i
I
:
..... ...........
- - - 4
t..-
43
: k./ :\
,~. . . . . . . ,~. . . . . . . ~," ," ,I __,,_. . . . • . . . . . . . ,___
i i
i
i i
i
0.98
7° 0.975
i
. . . . . . . . . . . . . . . . . . .
r
. . . . . . . . . . .
. . . . . . . . . . . .
i
-
i i i
i i
'
0.97
i i i
0.965
. . . . . . . . . . . . . . . .
. . . . . .
r - - - r
•
t i
i
0.96
.....
:,, V- _ ~i' . . . . . . .
-,,i. . . . . .
] ~___~___,~ i .....
~,. . . . . . . ~,i . . . ~. . / ]ii
0.955 0.95
-0.106
-0.076
-0.046
-0.015 0 0.015 frequency
0.046
0.076
L
0.106
FIG. 2. Frequency response of sinusoidal tapers for K = 8 and N -- 32.
+
7t2
(p2 _k2)2
zero
p-k
= odd
p-k
= even
(16)
Inferences from the above expression are, for odd p, odd numbered tapers will have null contribution, and for even p, even numbered tapers will have null contribution. At f = 0, the response is, t 8 ( N + 1)(K
~S-
2)/2
1
r=0~' ( 2 r + 1) 2
Keven
(17) 8 ( N + 1) (K- 1)/2 1 L ~7t r=0~ ~( 2 r + l )
Kodd
Table I lists the contribution of different tapers at the origin for N -- 32 and K = 8. An obvious deduction is that the farther tapers have lower contributions at the origin. For an expression of the errorband, it is to be noted that contribution of any taper is zero at frequencies evenps away from its own central frequency. Hence, for instance, the tapers contributing at p = K - 2 are k = K-- 2, K - 1, K-- 3, K-- 5 . . . . . while tapers k = K, K - 4, K-- 6 . . . . have zero contributions at that point. Because, the parity of K plays a major role in determining the response of the tapers, the following analysis is divided into two cases:
44
M . G. A m i n and N. Bhalla T A B L E I.
Contribution o f tapers at the origin (N = 32)
Taper
Value
1 3 5 7
26.75 2.97 1.07 0.55
Case l: K is odd
(Except for K = 1) Maximum occurs at f = 0 and the value is given by (17). As tapers contributing at p = K - 2 are K - 2, K - 1, K - 3, K - 5, . . . , the minimum value is given by:
-{-
It 2
[(K--2)2--4m2] 2
m=l
Case 2: K is even
Maximum occurs at p = 1 with the value:
-I-
7~2
m=l
[1--4m2] 2
(19)
Minimum value is: /N~
)
8(N+ ])(K~_2),,2 (2m+ 1)2 + ~2 ~"~--0 [(K--2)2-- (2m+ 1)2]2
(20)
For, K = 2, the minimum value is 8(N+ 1)/~Z 2. The approximate error band expression for K tapers is the difference between the maximum and the minimum for each of the two cases. We wish to mention that for practical purposes, all the terms in the summations need not be computed since the contribution of tapers away from their central frequency decreases with distance.
V. Fast Computation of Minimum Bias Tapers The DST is widely used in digital signal and image processing such as data compression, digital filtering of signals, image reconstruction, and image coding. Several papers have been published in the context of fast algorithms for DST computations and one of the efficient algorithms is proposed in (11), where Discrete Hartley Transform is applied. Four versions of Discrete Sinusoidal Transform (DST) have been introduced in literature. The expression obtained for the sinusoidal tapers in (4) is identical to the one for the DST of the first kind (12):
45
Minimum Bias Spatial Filters 1200
--
or}
I
I
DsT
i
- - Dir~-'t Multiplicdtion
Z
~1000
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
L
.
,
.
.
.
-J
.
.
.
.
L
___/32
I ~"
/
800
. . . . . . . . . . . . . . . . . . .
,
. . . . . . . . . . . . . . . . . . . . . . .
I.I.
/
0
//: /
--1
. . . . . . . . . . . . . . . . . . . I
200 . . . . . . . . . . .
/ /
i
//~
:'- - - - 2 - / - ~ ' -
. . . . . . . . . . . . . .
- ' - / /. . . .
/
,/
i
400
-F-G"
/
,
i
/
/
i
/ / /
/
~'i i
/
/
/
/
i ,/
/
-~"/ - - - - " - -"U /i
/
. . . . . . . . . . . . . . . . . . . . . . .
~ 24
,i-
/"
./
i
:- 2 - 2 0-
'
/
i
~--< . . . . . .i. . . . . . . . . . . . . . . . . . . . . . . . . . . i
~. ~
- ....
--
i
.~.
12 £//<.
. . . . .
f
~
600
~.- ~28 ~ .
i
i
,
/-
5
10
L i
I
I
15
20
25
30
35
K FIG.
3.
Comparison
[Sin,,] =
of direct multiplication
2~[sin(mn~/N)]m,
and DST.
n = 1, 2 . . . . .
N-
1
(21)
The DST provides the output of the tapers by performing: (N/2) (log2 N ) - N + 1 multiplications,
(22)
and (3N/2)
(log2
N) -2N-3N/8
+ log2 N + 16 additions
{23)
If the direct computation method, which amounts to multiplying the data snapshots with the sinusoidal tapers and summing is used, then N K multiplications will be required for K tapers. Figure 3 compares the number of multiplications needed for K tapers for direct method with those needed using DST for N tapers (N ~> K). Each line for direct method represents the number of multiplications for a particular number of sensors. The number of multiplications corresponding to a particular K for a given N can be obtained by moving along the corresponding line for N and using the x-axis for K. For example, the number of multiplications needed when using 10 tapers for 32 sensors is around 325. Figure 3 clearly illustrates that DST involves fewer multiplications than the direct method. Table II shows percentage savings for different values of N and K. Similar results can be obtained for the required number of additions. It is important to note that when DPSS are used as spatial filters, direct method has to be applied. From the above discussion, this implies more computations than sinusoidal
46
M. G. A m & and N. Bhalla TABLE II.
Efficiency offast computations N
K
DST
Direct method
Saving
8 32
2 8
5 49
16 256
68.75 80.86
tapers. Moreover, DPSS generation is done by eigendecompositon of N x N matrix. This itself creates a burden of computations of the order of N 3 (13). VI. Simulations The characteristics of MB, D F T and DPSS spatial tapers are compared in terms of their passband energy concentration and stopband attenuation. The 32-element array response of a spatial filter formed by a bank of tapers is shown in Fig. 4(a). Eight tapers were considered for both the MB and the D F T beamformers. Both cases offer the same bandwidth of 34 °, which is in turn used to determine the number of the DPSS tapers, as the greatest integer less than [2BN] (9). It is evident from Fig. 4(a) that the DPSS and the MB tapers have similar responses. The first sidelobe of the MB tapers is higher than that of the DPSS, but subsequently DPSS sidelobes have higher peaks. On the other hand, the stopband response of D F T is much higher than the other two. This implies that DPSS and MB tapers have similar rejection performance in their stopbands, which is better than the DFT. Figure 4(b) compares the minimum bias tapers and the sinusoidal tapers, and it is clear that MB tapers are very well approximated by the sinusoidal tapers. In Fig. 5, we compare beamspace MUSIC using the traditional D F T beams, the MB tapers of Eq. (3), the sinusoidal tapers given by Eq. (4), and the DPSS tapers. We consider linear uniform array of 32 sensors and 256 snapshots. For each case, eight spatial beams were formed covering a spatial sector of [0-34°]. The MUSIC spectrum should only be searched over this region, since the projection into the noise subspace of the signal directional vectors at(0) = T/~a(0), incident on the stopband, is small, and subsequently leads to high spectral values. The desired signal at 12 ° and an uncorrelated jammer at - 50 ° have SNR of 20 dB each. All the tapers are able to resolve the desired source lying in the passband. A different scenario is considered in Fig. 6, where the desired signal has SNR of 20 dB and its angle of arrival is 15°. There are six narrowband jammers, each of SNR 40 dB with angles of arrivals - 37.7, - 37.8, - 42.6, - 47.8, - 53.2, - 60.1 °. The temporal characteristics of each jammer at the first reference sensor is identical to the desired signal, presenting a highly coherent environment. It is evident from Fig. 6(a) that the MUSIC algorithm has failed to resolve the direction of arrival of the desired signal for D F T tapers. Figures 6(b)-(d) show that the desired signal is resolved and the performance of minimum bias is very close to that of DPSS. These results imply that in the scenario of multiple coherent jammers, the MB and DPSS reject the jammers effectively and the signal in the passband is resolved by a clear peak. Another inference
Minimum Bias Spatial Filters
47
(a)
o
,
m -10 ....
,......
, .......
'-I
b
*
:
i
g'
-15
- --1 ........ :
5
: :~-20--r-. -251-, t 1
i
,
, i
:
i ..... :
: I
I
' Ii'
,,,
: ....
,--
,,
',,
",',
t~ ,~ .'. ' i
F"~:,~t~ ~ I ~,'kL./ ~,' ^I" -35
'
_401 ~
'¢
tli,ll_~ I ,_'A:',_:,_
;
:'lli~
i
Ili~l]-q/r
: ,--,
~"'.l;~,~' "j lvtl*('f~-a'-~,-, t,,~l~.th,,I,,,
''
' Irlth)'T~l' f['
I/r'~'~
ll,~'!~tl"~' I •
| ~
k,~]*;i~
~'
',.~,,.~,'.fl",'Ir.~l',/i:l~'~
_ _,,,__,
II ',1 ,li
', I it
-40
-20
.....
'
'
" ,
,
: '
~__L, ............
I angle
!,
,'-~,
,,,
' 0
r
', I 20
40
~- :
_-k___
17,1 73
I' ',
.'',, ,
/ :
\
, ,-'--;, ,
;~1,
"
r
l
!
......
:
i,
~,l!;~i;
',
' i' j'ijI I1,,',(~ ',1',~ i I ~
-60
i
:
'~,1
•
,
-, j-'~,~---
~-' '-112~I,
',1
-80
~.,,I ~
'-~',~'t",',, ....
I I:,
,
~
i ~,;~',1,I~
', , , ,,,',,, ~',,' ~t.,,~, ,',K~III,~ i-~---,--~-v,~-i-.'-r~-~H~+l~,-H.k,','l~,-[#'ir~ ,, , , ,,, , ,, ~,,t I , , ,/tll,4 :'t ,, ''
:
.......
,
.I ":~, ...... q¼',l t
~,,tl'
h~ ,~ ,'#
i~
:..... :
--I-:-I-i'-'~i',
I/
~I-~-1-
i
, ........... :
;i--:'--'
' /
;1 i
"
'
,
'i
i
Ill ;,'
n-al
\
I. . . . . .
I i, 60
80
(degrees)
(b) 0
m
l
l
i
i
~--
!
I
Sinusoidal
MB
~" "O
-I0
--
,.
. . . . . . . . . . . . . . . . . . . . . . .
I
I I I
-15
....
-20
..........
. . . . . . . . . . . . . . . . . . .
I
--
r-
-25 ....
,........
',...................
:
-- ~ .......
~ ~/
-.........
.
-401
.
.
.
.
.
........
/
:..-
L___
,
f
1
J
i
i
-80
-60
-40
-20
0
20
angle
.~.,_",
40
i
i
60
80
(degrees)
FIG. 4. A r r a y response of (a) different tapers, (b) M B a n d sinusoidal tapers.
,
48
M. G. Amin and N. Bhalla
(a)
(b) 40 i
30
3O
20
iiiii .... i!
2O m
-o 10i 0
-I0
........
-90
-50
", ,---
~. . . . .
012 DFT
50
_,oiiiii!iiiiiii 90
-90
-50
012
50
90
MB
(d)
(c) 4O
40 ......
i
30
_
20
20
10 . . . . . . ' . . . . . .
1
---i- ---i
.....
. . . . . . . . . .
i
i i
0 -10 I
-90
-50
012 Sinusoidal
50
90
-90
-50
0 12
50
90
DPSS
FIG. 5. Beamspace MUSIC for uncorrelated sources: (a) DFT, (b) MB, (c) sinusoidal, (d) DPSS. that can be drawn is that the sinusoidal tapers are good approximation of MB tapers, since they work equally well. VII. Conclusions In this paper, we proposed the use of minimum bias (MB) tapers, introduced for temporal signal analysis, as spatial preprocessing filters. With MB tapers, a desired bandwidth can be achieved by simply adding or removing tapers, obviating the need for recomputations. The MB tapers can be well approximated by sinusoidal tapers, which have a simple closed form analytical expression. The properties and behavior of the sinusoidal tapers were delineated. An approximate expression for the errorband, i.e. the deviation between the maximum and the minimum values in the passband was derived. The application of the MB tapers to high resolution direction finding eigenstructure methods, specifically the MUSIC technique, was explained and elucidated by a set of simulations. For comparison of stopband attenuation provided by various tapers, the DFT, DPSS, MB and the sinusoidal tapers were employed for beamspace MUSIC and applied to a uniform linear array. The results show that the D F T beamspace MUSIC has the worst performance, whereas the other three sets of tapers perform equally well. This indicates that the MB tapers can be closely approximated by the sinusoidal tapers and, both perform similar to the well-known and extensively used DPSS tapers.
Minimum Bias Spatial Filters (a)
49
(b)
20~
v - - ~ - -
20 . . . . .
i
;' i
i
i
i
t
10 . . . .
. . . .
t~ "O
10
1-t
0
.........
"0
i r b
,, -10
-10 -90
-50
0 15
50
90
-90
-50
0 15
50
90
MB
DFT (c)
(d) i i
i i
i i
20
10
10 m 0
° _,o°
-10 -90
-
iiiiiiiiiiii
"O
-50
0 15
Sinusoidal
50
90
-90
-50
0 15
50
90
DPSS
FIG. 6. Beamspace MUSIC in highly coherent scenario: (a) DFT, (b) MB, (c) sinusoidal, (d) DPSS. The crucial aspect of fast computations when utilizing the proposed set of tapers was explored. It was shown that the expression for the sinusoidal tapers resembles the one for the Discrete Sinusoidal Transform (DST). Since substantial endeavors have been made to develop fast computational algorithms for the DST, the sinusoidal tapers based beamspace method is amenable to efficient implementations and has a clear computational advantage over that of the DPSS.
References (1) F. Haber and M. Zoltowski, "Spatial spectrum estimation in a coherent signal environment using an array in motion", IEEE Trans. Antennas Propaq., Vol. AP-34, No. 3, pp. 301310, 1986. (2) T. J. Shan, M. Wax and T. Kailath, "On spatial smoothing for direction-of-arrival estimation of coherent signals", IEEE Trans. Aeoust., Speech, Signal Processing, Vol. ASSP33, No. 4, pp. 806-811, 1985. (3) G. Bienvenu and L. Kopp, "Decreasing high resolution method sensitivity by conventional beam former preprocessing", "Proc. IEEE Int. Conf. Acoust., Speech, Signal Processing", pp. 33.2.1-33.2.4, 1984. (4) C. L. Byrne and A. K. Steele, "Sector-focused stability for high resolution array processing", "Proc. IEEE Int. Conf. Acoust., Speech, Signal Processing", pp. 234~2343, 1987.
M. G. Amin and N. Bhalla
50
(5) P. Forster and G. Vezzosi, "Application of spheroidal sequences to array processing", "Proc. IEEE Int. Conf. Acoust., Speech, Signal Processing", pp. 2267-2271, 1987. (6) H. Lee and M. Wengrovitz, "Resolution threshold of beamspace MUSIC for two closely spaced emitters", IEEE Trans. Acoust., Speech, Signal Processing, Vol. 38, No. 9, pp. 1545-1559, 1990. (7) X. L. Xu and K. M. Buckley, "Statistical performance comparison of MUSIC in element space and beam-space", "Proc. IEEE Int. Conf. Acoust., Speech, Signal Processing", pp. 2124-2127, 1989. (8) C. Zhou, F. Haber and Q. Shi, "On spatial filtering for angle of arrival estimation in a scattering environment by eigendecomposition-based methods", IEEE Trans. Antennas Propag., Vol. 38, No. 6, pp. 931-934, 1990. (9) D. Slepain, "Prolate spheroidal wave functions, Fourier analysis and uncertainty: V: The discrete case", BellSyst. Tech. J., Vol. 5, pp. 1371-1429, 1978. (10) K. Riedel and A. Sidorenko, "Minimum bias multiple tapering spectral estimation", IEEE Trans. Signal Processing, Vol. 43, No. 1, pp. 188-195, 1995. (11) P. S. Kumar ad K. M. Prabhu, "A set of new fast algorithms for DCTs and DSTs", Int. J. Electronics, Vol. 76, No. 4, pp. 553-568, 1994. (12) Z. Wang, "Fast algorithms for the discrete W transform and for the discrete Fourier transform", IEEE Trans. Acoust., Speech, Signal Processing, Vol. 32, No. 8, pp. 803-816, 1984. (13) G. H. Golub and C. F. Van Loan, "Matrix Computations", Second Edn, The Johns Hopkins University Press, 1989. (14) A. P. Prudnikov et al. "Integrals and Series", Vol. 1, Gordon & Breach, 1988. (15) M. Abramowitz and I. Stegun, "Handbook of Mathematical Functions", Dover Publications, NY, 1970.
Appendix A The product of sinusoidal tapers is written in the equivalent exponential form to facilitate summation over k. The summation in (7) will give four terms represented by: e +_J[nn(kl±k2)/(N+
I)]
Le+-Jt"(k'±k2)/(N+l)l_1 J
n=l
L
J (A1)
For even values of (kz +k2), the above expression yields values of - 1. Hence the terms in the summation cancel to produce a zero. For odd values of (k~ _+k2), the expression takes the form: [~ + e-+/~"(k'±k2'"(N+' __e±j[.(k,+_kz)/(N+l)lJ
=
Fz~(kl + k : ) l
- j c ° t [_2 ~
j
The terms in the summation will be of opposing signs, with the same magnitude and will cancel out each other to give a spatial null. To check normality, set kl = k2 in (7):
2(N+I)
l.=,
=
2 ' A ' + l " { - 2 - 2 N } = 1 (A2)
Minimum Bias Spatial Filters
51
Appendix B Equation (5) shows that the frequency response of the kth taper has two terms of the form sin nNx/sin nx. This structure implies that the sidelobes of one term will have a negligible effect on the main lobe of the other, e.g. the value of one term at the center of the main lobe of the other (where the height is N), is ( - 1 ) k+~. For k = 1, the two components are quite close and hence the sidelobes have higher contributions to the other component, but because the two terms are additive, the result is a nonzero response. The main lobe of one component of the kth taper is centered at k/2(N+ 1) while that of the adjacent tapes is a ( k + I ) / 2 ( N + 1) and their midpoint is ( k + 1/2)/2(N+ 1). The first null of a component of the taper is 1IN from its center, which is farther away from the midpoint of the two tapers, which is 1/4(N+ 1) away from the center. Hence the main lobe of each taper component provides nonzero array response before it intersects with the main lobe of the next one. Consequently, a set of contiguous tapers will provide a nonzero passband. To prove the second property, consider N
sin(nxD eaM~"- l)sinO R
Tna(0) ~ T,~ a(0)
.=, C N
(B1)
~ (n - 1) cos 0 sin(nxk) esu~.
Utt
I )sinO
n=l
where C is a constant, xk = 7tk/(N+ 1), k = 1, 2, . . . , K and M = 2~d/2. For the second property to be true, the summations in (B1) for different ks should yield the same R, which does not happen in this case due to the effect of the terms ( n - 1) and cos 0.
Appendix C The expression for the second derivative of the Fourier transform of sum of K tapers, at f = 0 is: x nk 47~2 xk nk ~ cot 2 cosec 2 21t2 (N+ 1) Y ( - 1)k cot 2 - - + ~ k-% 2 ( N + 1) ~ / + I k = odd 2 ( N + 1) 2 ( N + 1) =
2/r2 & 2 nk F L cot - - / c o s e c N + 1 k=, 2(N+l) L
2
rck { 2 rck )7 + (-- 1) k ( N + 1) 2 2 ( N + 1) _ - c o s e c ~ N + 1);J
(C1)
Here we employ approximation sin x ~ x and simplify the expression further to: 4
8(U+ 1)3 Z
= 8 ( N + 1 ) 3 ~=, ~ ~ ( ~( --l ) k
+ ~ 4: k ~ - - ( - - l ) k4~
}
(C2)
The above expression is solved with the help of the following formulae (14): (--1)k k~l
k~
(--1) m - - ( m _ l ) ! [ f l (m °(1)-(-1)xfl(m-*~(K+l)]
1 -1 m k=,/d' -- ( - - ) ! [0('~-'~(1)-0(m-'~(K+l)]
re=l,2 ....
m = 1,2 . . . .
(C3)
(C4)
M. G. Amin and N. Bhalla
52
fl<")(z) = 1/(2 "+ J){~b"[(z+ 1 ) / 2 ] - ~"(z/2)}, where ~,(z) is the psi or digamma function, defined as ~b(z) = (d/dz)[ln F(z)] = F'(z)/F(z), F(z) being the gamma function. For a given value of K, the solution for (C4) is obtained and it is seen that it takes a negative value for odd K and a positive one for even K, signifying a maxima and a minima, respectively. The recurrence formulae for the psi function (15) are useful in computations. For example, when K = 3,
A
' 1
1
'
~j
Similarly Z3= j [ ( - 1)k/k 4] = --0.9498 and Z3= ~(1/k 4) = 1.0748. Putting these values in (C2) will give the value inside the summation as - 0 . 0 4 0 5 , which implies a maximum.
Appendix D The sum of squared values of magnitude of Fourier response of first K tapers is given by: 1 2 ( N + 1)
x sin2 ( N + l)r~f- - -
sin 2 N + ~
=
.k
k~'
(
~2
(DI)
sin2(nJ) --sin2 \ 2 ( N + 1)]J
After approximations of the kind sin x ~ x for small x, the value atfp is:
+ 2(N+l~,~t~)
1
~<,<~-p~)~
(°2)
{[2(N+ i)}2J For even ( k - p ) , the expression vanishes. For the cases when ( k - p ) is odd, it can be simplified to give (16).
~ Pergamon
PII:
S0016--0032(96)00114--7
Copyright © 1997The Franklin Institute Published by ElsevierScienceLtd Printed in Great Britain 0016-0032/98 $19.00+0.00
Minimum Bias Spatial Filtersfor Beamspace Direction-of-arrival Estimation by MOENESS
G. AMINt
and N I T I N
BHALLA
Department o f Electrical and Computer Engineering, Villanova University, Villanova, PA 19085, U.S.A.
ABSTRACT: Spatial filtering is used in array signal processing to favor signals coming from a particular zone in space over other signals. In this paper, we propose the use of the minimum bias (MB) windows, originally devised for multitaper spectral analysis, as spatial preproeessing filters for direction-of-arrival (DOA) estimation. These windows (or tapers) are characterized by their orthonormality, simplicity, and most importantly their close behavior to the well known discrete prolate spheroidal sequence (DPSS) windows. A desired spatial bandwidth can be achieved by adding or removing tapers, without the need of reeomputation. The paper introduces the MB tapers for beamspace processing and compares their DOA estimation performance, using noise-subspace eignestructure methods, to that of the DPSS beams. Further, the MB windows can be well approximated by simple trigonometric functions. These functions have expressions identical to that of the Discrete Sinusoidal Transform (DST). As numerous methods f or fast computations of DST have been devised, the outputs of the approximated MB tapers can be computed efficiently. © 1997 The Franklin Institute. Published by Elsevier Science Ltd
L Introduction
In array signal processing, it is often desired that signals coming from a particular zone be removed, or at least attenuated. These signals may represent intentional jammers or multipath arrivals. A common example is a sensor array mounted on an aircraft which receives signals from a sonobuoy. The scatters of the sonobuoy signal from the aircraft propellers cause a severe multipath effect on the array performance, and therefore, their effect should be mitigated. Due to the coherent signal environment, the adaptive spatial nulling techniques (high resolution DOA spectra) not only fail to form nulls (peaks) in the directions of the interferers and their multipaths, but also tend to cancel the desired signal in the array output. This failure occurs even with the decorrelation effects introduced by the motion of the array on the aircraft (1).
t This work is supported by ONR grant N00014-94-1-1052. 35
36
M. G. Amin and N. Bhalla
Preprocessing techniques such as subarray averaging methods (2) can be employed to mitigate this problem and decorrelate the coherent arrivals prior to both hulling and DOA estimation. These methods, however, reduce the array aperture and are appropriate for point sources. They do not perform well for small arrays, few snapshots and/or under models of distributed sources. An alternative is to use spatial filtering to filter out both uncorrelated and coherent unwanted signals in the spatial domain, which is also known as beamspace processing. This is achieved by lowering the array response over specific spatial sectors which include the jammer arrivals. Accordingly, in a multipath environment, the coherent counterparts of direct path are suppressed instead of being decoupled. The response over the desired sector should be uniform and of unit gain. Beamspace methods have been employed for a variety of applications in the literature. The obvious advantage is the reduction in both the interference and noise whose spatial spectrum lies outside the spatial sector of interest, i.e. increasing the signal-tointerference and noise-ratio (SINR). They also lead to reduced computational complexity, amenability to parallel implementation, low sensitivity to sensor perturbations and deviations from the assumed noise model (3-5), and improving the spatial resolution (6, 7). Multiple orthogonal beams can be generated using spatial filters with different characteristics, each receiving the interference over its stopband. To capture the desired signals, the combined beams should uniformly cover the spatial sector of interest. When properly designed, the outputs of the spatial filters are jammer free and can be used as input for eigenstructure-based methods to perform high-resolution DOA estimation of the desired signals. A popular and widely used set of spatial tapers (filters) is the family of discrete prolate spheroidal sequence (DPSS) (8). Their basic properties were given by Slepian (9). These tapers have much lower sidelobes than rectangular windows and, when used as spatial filters, they yield better performance than DFT beamformers. Recently, a new family of orthonormal basis have been introduced for minimum bias estimates of the power density function in the context of time-series analysis (10). In this paper, we apply the minimum bias (MB) tapers for beamspace processing. The MB tapers can be well approximated by trigonometric functions, and as such, have a much simpler form than the DPSS tapers. A desired frequency bandwidth can be achieved by adding or removing tapers, without the need for recomputations, which is the case with DPSS. The approximated MB tapers, called the sinusoidal tapers, have a closed form analytical expression, which is identical to the sinusoidal basis functions of the Discrete Sinusoidal Transform (DST). Since numerous fast algorithms for DST evaluation exist, efficient computation of outputs corresponding to sinusoidal tapering beams is possible. The minimum bias tapers are discussed in detail in Section II. This section also provides the expressions for the approximated versions of the MB tapers, the sinusoidal tapers. In Section III, we explain the use of the MB tapers for spatial signal processing. The properties of sinusoidal tapers are discussed in detail in Section IV. Section V deals with the important aspect of fast computation of the outputs of the different beams utilizing the sinusoidal tapers. The similarity between the expressions for sinusoidal tapers and DST offers the opportunity to achieve fast computations. Absence of such a feature in DPSS makes the sinusoidal tapers a favorable proposition. The fast DST
M i n i m u m Bias Spatial Filters
37
computations of the data in beamspace are compared with the conventional procedure of inner products, i.e. multiplying the data snapshots with the tapers followed by summing. In Section VI, a comparison of DOA estimation efficacy of MB, sinusoidal, DFT, and DPSS tapers, using subspace methods, is provided to elucidate the behavior of the proposed set of tapers. 11. Minimum Bias Tapers In time-series analysis, the spectrum estimate bias is defined as E[5~(/) - S(f)] = ff~o ([ V ( 9 - f ) I2S O ) - S ( f ) ) d9
(1)
where, S ( f ) is the true spectral density whose estimate 5~(f) is given by S(f) =
~ v(n)x(n) e j2~znl ~
(2)
n=]
The observations {x(n), n = 1,2 . . . . . N} form a time series and V ( f ) is the Fourier transform of the data window (taper) v(n). It is shown in (10) that the minimum value of the bias in (1) is obtained when v(n) is such that S ~ ] V ( f ) ] z f 2 dfis minimized. For discrete-time multiple window spectral analysis, the minimum bias tapers are given by the eigenvectors v°), v~2). . . . , v~ , corresponding to the eigenvalues, sorted in the increasing order, of the matrix C = [Cm,]: l/12
c ..... =
f
(-1)
n= m n-m
2rr 2(n -- m) 2
(3)
nCm
These tapers form an orthonormal set. The taper v~i)provides the minimum possible spectral bias subject to v°3 _1_vv'), j = 1, 2, . . . , i, where 3_ denotes orthogonality. We note that in spatial filtering applications, N represents the number of array sensors rather than the number of observations. In practical situations involving beamspace processing, more than one taper should be used to devise the passband of the spatial filters, for the following reasons: (1) although some a priori knowledge is assumed to be available about the direction of arrival of the desired signals, this information is usually not complete. Using too few tapers may not sufficiently cover the spatial zone within which the desired sources are located; (2) eigenstructure methods for DOA estimation require a minimum array aperture size in accordance with the number of signals whose parameters are to be estimated. These methods are applicable only if the array aperture is of size one more than the number of arrivals. As a result of processing the data by spatial filters, the effective array aperture is reduced. The beamspace array aperture becomes equal to the number of tapers used, which is invariably less than the element space aperture (maximum number of tapers that can be used is equal to the number of sensors, which amounts
M. G. Amin and N. Bhalla
38
to an allpass spatial filter for MB tapers). Hence, if the parameters of P signals are to be estimated by beamspace eigenstructure methods, at least ( P + 1) tapers should be used. The MB tapers can be well approximated by the orthonormal sinusoidal functions, u"), u (z). . . . . u (m (10). The elements of u (k), (u~~), u(2k).... , u~)) T are defined as u~ =
/~
.
~k,~
k = 1,2 . . . . . N
4 N~i- sin Nq_-1
n= 1,2,...,N
(4)
The temporal frequency f is being utilized to elucidate the concepts henceforth. Nevertheless, the results obtained can be directly applied to array processing by appropriately substituting the temporal frequency by the spatial frequency. Spatial response is calculated using directional vectors which have been normalized by the factor x/l/N. Below, we discuss the orthonormal properties of sinusoidal tapers. The Fourier transform of the kth taper for a linear array of N sensors is given by (lO), N
V~k~(f) = ~ v(k)e j2~nl n=l
e --jn[(N+ 1)f-- (k/2)]
j x~(( N + 1)
sin [z~ ( f - 2(NL 1))]
s i n I N n ( f + k2(N+
--(--l)k [~(f+ sin
1))]
t
k
(5)
2(N+ 1))1
sin N ~ I V(k](f) = e
/~[(N+J)f-((k J).'2]]
,jSU l)
(6) sin2 (nf) - sin2 (~2(~+ 1)-)
Equation (5) clearly shows that the frequency response of every taper has two components which have a maximum each at k/2(N+ 1) and - k / 2 ( N + 1), respectively, and thus they are symmetrically placed about the origin. The maximum value (at each of these frequencies) of the kth taper is ~ / 2 . The sidelobes diminish in value as we move away from the maxima. The sidelobes of the two terms are of opposite polarity and tend to cancel each other effectively, resulting in much lower sidelobes for the taper as a whole. The frequency response of the kth taper thus has two peaks, equidistant from the origin. The peaks are slightly shifted from the positions Ik/2(N+ 1)[ because the sidelobes of one term interact with the maximum of the other. This shift is more
Minimum Bias Spatial Filters
39
pronounced for small values of k, where the peaks of the two components are close to each other and hence sidelobe contribution is higher. Orthogonality of tapers is a desired property because it eliminates the redundancy and simplifies computations (6). To prove orthonormality of the sinusoidal tapers, consider the inner product bk,,k, of the two tapers corresponding to kl and k2: (2)~=.
7tkln.
rck2n
1 :
N ~
--
E ( ej[~k'n/(N+l)l-e/[r~(k[n)/(N+l)l)(e/[~k2n/IN+l)l-e
/[~k2,;N+I)I) (7)
n=l Using the results of Appendix A, the above expression simplifies to:
(1)~eJn(kl+ke)--eJ[n(kl+k2)/(N+l)]eln(k'k2)--e b~-,'k~ - --
2(N+ 1).
(
.
. . . l - -. e j[n(k)+k2)/(N+l)]
e Jn(kt - k z ) --
1--e
e -jpt(kl -k2)/(N+ 1)]
1 - - e j[.(k, -k2)/(m+ 1)]
jtn(kl-k2)/(N+ l)]
][Tt(kI-k2)/(N+l)]
e - ~ ( k ~ +k 2) _ e-J[Tt(kl +k2)/(N+ 1 )] "~
+
1 --e -j[~(k~+k2)/(N+ 1)]
I
(8)
In Appendix A, it is shown that for the two cases of (k~ q-k2) being odd or even, the above expression sums to zero, implying orthogonality. For k~ = k2, the sum is unity, indicating normality. IlL M i n i m u m Bias Spatial Filters Let X be the data matrix X = AS + ~/
(9)
where A is the steering matrix associated with the arrivals of both signals and jammers, and S is the signal matrix whose columns are the snapshots of the directional array inputs collected at different time instants. The noise is assumed to be spatially white, i.e. the components of the nondirectional noise vector ~/are uncorrelated. The spatial tapers, whether they represent the DPSS, MB, DFT, or any other beamformer, can be used to define the spatial filter matrix T. Accordingly, the beamformers' outputs are given by Y = THX
(10)
Figure 1 illustrates the transition from element space to beamspace. (For the minimum bias case, the number of these tapers is easily determined as a function of the desired spatial bandwidth.) If the spatial sector of interest is not centered at 0 = O, but rather at 00, then one can steer the array to the desired angular position by modifying the spatial matrix as follows: F = DT
(11)
where F is the spatial filter matrix steered to desired angle, and D = [din,,],
M. G. Amin and N. Bhalla
40
Beam #1
Beam #2
Beam #3
/
Element space (N) Beam space (K)
FIG. I. Elementspace to beamspace.
din,, = 6m, exp [j2rc sin(Oo)dm/2) is the diagonal matrix corresponding to the steering angle 0o. The intersensor spacing of the array is d and 2 is the wavelength of operation. The spatial frequency corresponding to the angle 0o is (d/2) sin 00. The beamspace snapshots (data) constituting the columns of Y can then be used for DOA via subspace methods. The DOA are the peak locations of the spectrum:
S(O) =
1
[a~(O)V.nV~ai(O)l
(12)
Minimum Bias Spatial Filters
41
where, as(0) defines the beamspace steering vector (aj(0) = Tna(0), a(0), being the element space steering vector) and V, is the noise subspace projection matrix whose columns are the eigenvectors corresponding to the minimum eigenvalues of the correlation matrix formed by the filtered data. Matrix H is a suitably selected nonnegative square matrix of size ( K - L), where L is the number of arrivals whose parameters are to be estimated (6). In the case of beamspace MUSIC, H = I, while for beamspace MinNorm, H = ssH/(sHs) 2, where s denotes the ( K - L) × 1 vector Vff(1,0 . . . . . O)r.
IV. Properties of the Sinusoidal Tapers 4.1. General properties In general, a beamformer should have the following two properties, for 0 in the passband (6): (a) it does not introduce superfluous zeros in the passband, i.e. aH(O)TTHa(O) > 0
(13)
This property is important because, if the beamformer itself has nulls in the spatial sector of interest, it will lead to spurious peaks in the spectrum at those locations and also, remove any arrivals at those points, resulting in erratic inferences; (b) the first order differences in the element space directional vectors are translated into first order differences in the beamspace directional vectors, implying Tn~0a(0) =# cTHa(O)
(c = scalar)
(14)
This property ensures that any change in the steering vectors in the element space is not smoothed out due to beamspace processing, but is appropriately reflected in the beamspace steering vectors. In Appendix B, it is shown that these two above mentioned properties hold for the sinusoidal tapers. 4.2. Specific properties In Section II, we briefly discussed the sidelobe and orthonormality properties of the sinusoidal tapers. Below, we elaborate on the overall frequency characteristics of these tapers. The analysis of the combined response of say, K tapers (K < N) will be of importance. From (6), it is evident that, even numbered tapers vanish at f = 0 as the numerator evaluates to zero, implying that they do not contribute at the origin. The derivative of squared response of absolute value of kth taper is evaluated by differentiating (6) w.r.t.f. Summing the above derivative for N = 1, 2 . . . . . K gives the corresponding expression for sum of K tapers: ~r
K
(-
~k 1)k sin2 N + ~
~sin [2(N+ 1)Trf] k= ' {sin2 (Trf) - sin2 ( 2 ~ +
~)}2
42
M. G. Amin and N. Bhalla
rc
I
zk]. 2
K sin2 ( N + l ) r c f - ~ - J s a n N+lZCk
,~r± 1', sin 2xf ~ \~'Tx]
k~l
(15) sin2(ltf)-sin2
2(~
For any K, a maxima or minima exists a t f = 0 as the derivative is zero. It is shown in Appendix C that odd values of K will result in a maximum whereas even values of K provide a minimum spatial response at the origin. Because even tapers do not contribute at the origin, as we increment K, going from an odd to an even value, the response a t f = 0 does not change and in fact, maxima of an odd taper and minima of the next taper, which is even, coincide. 4.2.1. Error band of the spatial filter. Filter design problems have an important parameter--the passband ripple. For the spatial filter designed by sinusoidal tapers, we attempt to derive a similar parameter, the error band. This parameter gives the deviation between the maximum and the minimum points in the filter passband. An obvious method of deriving an expression for the error band is: (a) find the derivative of the sum of squares of magnitudes of frequency responses of tapers; (b) determine the points of maxima and minima in the passband; (c) the difference between the lowest minimum and the highest maximum points gives the error band. Equation (15) gives the expression for the desired derivative. This equation can be simplified by assuming N t o be sufficiently large such that rck/(N+ 1) is small enough for reasonable approximations in the numerator. This allows us to make approximations of the kind sin x ~ x for small x, but of course with inherent approximation errors. The magnitude of this error reduces as N increases. The complexity of the problem requires the use of numerical methods to find its roots. In Fig. 2, a zoomed view of the passband corresponding to 32 sensors and 8 tapers is shown. The vertical grid lines are at frequenciesfp = p/2(N+ 1) where p is an integer (these are the frequencies where individual tapers have the peaks of their components). Interestingly, a close observation of the frequency response of the tapers in Fig. 2 and the values of extrema obtained by MATLAB routines, reveals the following: (1) the extrema lie very close to the frequenciesfp; (2) the lowest minimum in the passband occurs near p = K - 2 ; (3) the highest maxima for even K lies nearfp= l, and for K odd a t f = 0. This pattern motivates to study the behavior of frequency response at frequenciesfp for different K. In Appendix D, it is shown that the response of sums of squares of the absolute values of the first K tapers at fp (k = 1,2 . . . . . K), after employing relevant approximations, is:
Minimum Bias Spatial Filters I I I I I I
~= e
0.995
C/~
~
I I I I
I
0.99
f,\
I
~ 0.985 "0
. -
,i
I
:
..... ...........
- - - 4
t..-
43
: k./ :\
,~. . . . . . . ,~. . . . . . . ~," ," ,I __,,_. . . . • . . . . . . . ,___
i i
i
i i
i
0.98
7° 0.975
i
. . . . . . . . . . . . . . . . . . .
r
. . . . . . . . . . .
. . . . . . . . . . . .
i
-
i i i
i i
'
0.97
i i i
0.965
. . . . . . . . . . . . . . . .
. . . . . .
r - - - r
•
t i
i
0.96
.....
:,, V- _ ~i' . . . . . . .
-,,i. . . . . .
] ~___~___,~ i .....
~,. . . . . . . ~,i . . . ~. . / ]ii
0.955 0.95
-0.106
-0.076
-0.046
-0.015 0 0.015 frequency
0.046
0.076
L
0.106
FIG. 2. Frequency response of sinusoidal tapers for K = 8 and N -- 32.
+
7t2
(p2 _k2)2
zero
p-k
= odd
p-k
= even
(16)
Inferences from the above expression are, for odd p, odd numbered tapers will have null contribution, and for even p, even numbered tapers will have null contribution. At f = 0, the response is, t 8 ( N + 1)(K
~S-
2)/2
1
r=0~' ( 2 r + 1) 2
Keven
(17) 8 ( N + 1) (K- 1)/2 1 L ~7t r=0~ ~( 2 r + l )
Kodd
Table I lists the contribution of different tapers at the origin for N -- 32 and K = 8. An obvious deduction is that the farther tapers have lower contributions at the origin. For an expression of the errorband, it is to be noted that contribution of any taper is zero at frequencies evenps away from its own central frequency. Hence, for instance, the tapers contributing at p = K - 2 are k = K-- 2, K - 1, K-- 3, K-- 5 . . . . . while tapers k = K, K - 4, K-- 6 . . . . have zero contributions at that point. Because, the parity of K plays a major role in determining the response of the tapers, the following analysis is divided into two cases:
44
M . G. A m i n and N. Bhalla T A B L E I.
Contribution o f tapers at the origin (N = 32)
Taper
Value
1 3 5 7
26.75 2.97 1.07 0.55
Case l: K is odd
(Except for K = 1) Maximum occurs at f = 0 and the value is given by (17). As tapers contributing at p = K - 2 are K - 2, K - 1, K - 3, K - 5, . . . , the minimum value is given by:
-{-
It 2
[(K--2)2--4m2] 2
m=l
Case 2: K is even
Maximum occurs at p = 1 with the value:
-I-
7~2
m=l
[1--4m2] 2
(19)
Minimum value is: /N~
)
8(N+ ])(K~_2),,2 (2m+ 1)2 + ~2 ~"~--0 [(K--2)2-- (2m+ 1)2]2
(20)
For, K = 2, the minimum value is 8(N+ 1)/~Z 2. The approximate error band expression for K tapers is the difference between the maximum and the minimum for each of the two cases. We wish to mention that for practical purposes, all the terms in the summations need not be computed since the contribution of tapers away from their central frequency decreases with distance.
V. Fast Computation of Minimum Bias Tapers The DST is widely used in digital signal and image processing such as data compression, digital filtering of signals, image reconstruction, and image coding. Several papers have been published in the context of fast algorithms for DST computations and one of the efficient algorithms is proposed in (11), where Discrete Hartley Transform is applied. Four versions of Discrete Sinusoidal Transform (DST) have been introduced in literature. The expression obtained for the sinusoidal tapers in (4) is identical to the one for the DST of the first kind (12):
45
Minimum Bias Spatial Filters 1200
--
or}
I
I
DsT
i
- - Dir~-'t Multiplicdtion
Z
~1000
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
L
.
,
.
.
.
-J
.
.
.
.
L
___/32
I ~"
/
800
. . . . . . . . . . . . . . . . . . .
,
. . . . . . . . . . . . . . . . . . . . . . .
I.I.
/
0
//: /
--1
. . . . . . . . . . . . . . . . . . . I
200 . . . . . . . . . . .
/ /
i
//~
:'- - - - 2 - / - ~ ' -
. . . . . . . . . . . . . .
- ' - / /. . . .
/
,/
i
400
-F-G"
/
,
i
/
/
i
/ / /
/
~'i i
/
/
/
/
i ,/
/
-~"/ - - - - " - -"U /i
/
. . . . . . . . . . . . . . . . . . . . . . .
~ 24
,i-
/"
./
i
:- 2 - 2 0-
'
/
i
~--< . . . . . .i. . . . . . . . . . . . . . . . . . . . . . . . . . . i
~. ~
- ....
--
i
.~.
12 £//<.
. . . . .
f
~
600
~.- ~28 ~ .
i
i
,
/-
5
10
L i
I
I
15
20
25
30
35
K FIG.
3.
Comparison
[Sin,,] =
of direct multiplication
2~[sin(mn~/N)]m,
and DST.
n = 1, 2 . . . . .
N-
1
(21)
The DST provides the output of the tapers by performing: (N/2) (log2 N ) - N + 1 multiplications,
(22)
and (3N/2)
(log2
N) -2N-3N/8
+ log2 N + 16 additions
{23)
If the direct computation method, which amounts to multiplying the data snapshots with the sinusoidal tapers and summing is used, then N K multiplications will be required for K tapers. Figure 3 compares the number of multiplications needed for K tapers for direct method with those needed using DST for N tapers (N ~> K). Each line for direct method represents the number of multiplications for a particular number of sensors. The number of multiplications corresponding to a particular K for a given N can be obtained by moving along the corresponding line for N and using the x-axis for K. For example, the number of multiplications needed when using 10 tapers for 32 sensors is around 325. Figure 3 clearly illustrates that DST involves fewer multiplications than the direct method. Table II shows percentage savings for different values of N and K. Similar results can be obtained for the required number of additions. It is important to note that when DPSS are used as spatial filters, direct method has to be applied. From the above discussion, this implies more computations than sinusoidal
46
M. G. A m & and N. Bhalla TABLE II.
Efficiency offast computations N
K
DST
Direct method
Saving
8 32
2 8
5 49
16 256
68.75 80.86
tapers. Moreover, DPSS generation is done by eigendecompositon of N x N matrix. This itself creates a burden of computations of the order of N 3 (13). VI. Simulations The characteristics of MB, D F T and DPSS spatial tapers are compared in terms of their passband energy concentration and stopband attenuation. The 32-element array response of a spatial filter formed by a bank of tapers is shown in Fig. 4(a). Eight tapers were considered for both the MB and the D F T beamformers. Both cases offer the same bandwidth of 34 °, which is in turn used to determine the number of the DPSS tapers, as the greatest integer less than [2BN] (9). It is evident from Fig. 4(a) that the DPSS and the MB tapers have similar responses. The first sidelobe of the MB tapers is higher than that of the DPSS, but subsequently DPSS sidelobes have higher peaks. On the other hand, the stopband response of D F T is much higher than the other two. This implies that DPSS and MB tapers have similar rejection performance in their stopbands, which is better than the DFT. Figure 4(b) compares the minimum bias tapers and the sinusoidal tapers, and it is clear that MB tapers are very well approximated by the sinusoidal tapers. In Fig. 5, we compare beamspace MUSIC using the traditional D F T beams, the MB tapers of Eq. (3), the sinusoidal tapers given by Eq. (4), and the DPSS tapers. We consider linear uniform array of 32 sensors and 256 snapshots. For each case, eight spatial beams were formed covering a spatial sector of [0-34°]. The MUSIC spectrum should only be searched over this region, since the projection into the noise subspace of the signal directional vectors at(0) = T/~a(0), incident on the stopband, is small, and subsequently leads to high spectral values. The desired signal at 12 ° and an uncorrelated jammer at - 50 ° have SNR of 20 dB each. All the tapers are able to resolve the desired source lying in the passband. A different scenario is considered in Fig. 6, where the desired signal has SNR of 20 dB and its angle of arrival is 15°. There are six narrowband jammers, each of SNR 40 dB with angles of arrivals - 37.7, - 37.8, - 42.6, - 47.8, - 53.2, - 60.1 °. The temporal characteristics of each jammer at the first reference sensor is identical to the desired signal, presenting a highly coherent environment. It is evident from Fig. 6(a) that the MUSIC algorithm has failed to resolve the direction of arrival of the desired signal for D F T tapers. Figures 6(b)-(d) show that the desired signal is resolved and the performance of minimum bias is very close to that of DPSS. These results imply that in the scenario of multiple coherent jammers, the MB and DPSS reject the jammers effectively and the signal in the passband is resolved by a clear peak. Another inference
Minimum Bias Spatial Filters
47
(a)
o
,
m -10 ....
,......
, .......
'-I
b
*
:
i
g'
-15
- --1 ........ :
5
: :~-20--r-. -251-, t 1
i
,
, i
:
i ..... :
: I
I
' Ii'
,,,
: ....
,--
,,
',,
",',
t~ ,~ .'. ' i
F"~:,~t~ ~ I ~,'kL./ ~,' ^I" -35
'
_401 ~
'¢
tli,ll_~ I ,_'A:',_:,_
;
:'lli~
i
Ili~l]-q/r
: ,--,
~"'.l;~,~' "j lvtl*('f~-a'-~,-, t,,~l~.th,,I,,,
''
' Irlth)'T~l' f['
I/r'~'~
ll,~'!~tl"~' I •
| ~
k,~]*;i~
~'
',.~,,.~,'.fl",'Ir.~l',/i:l~'~
_ _,,,__,
II ',1 ,li
', I it
-40
-20
.....
'
'
" ,
,
: '
~__L, ............
I angle
!,
,'-~,
,,,
' 0
r
', I 20
40
~- :
_-k___
17,1 73
I' ',
.'',, ,
/ :
\
, ,-'--;, ,
;~1,
"
r
l
!
......
:
i,
~,l!;~i;
',
' i' j'ijI I1,,',(~ ',1',~ i I ~
-60
i
:
'~,1
•
,
-, j-'~,~---
~-' '-112~I,
',1
-80
~.,,I ~
'-~',~'t",',, ....
I I:,
,
~
i ~,;~',1,I~
', , , ,,,',,, ~',,' ~t.,,~, ,',K~III,~ i-~---,--~-v,~-i-.'-r~-~H~+l~,-H.k,','l~,-[#'ir~ ,, , , ,,, , ,, ~,,t I , , ,/tll,4 :'t ,, ''
:
.......
,
.I ":~, ...... q¼',l t
~,,tl'
h~ ,~ ,'#
i~
:..... :
--I-:-I-i'-'~i',
I/
~I-~-1-
i
, ........... :
;i--:'--'
' /
;1 i
"
'
,
'i
i
Ill ;,'
n-al
\
I. . . . . .
I i, 60
80
(degrees)
(b) 0
m
l
l
i
i
~--
!
I
Sinusoidal
MB
~" "O
-I0
--
,.
. . . . . . . . . . . . . . . . . . . . . . .
I
I I I
-15
....
-20
..........
. . . . . . . . . . . . . . . . . . .
I
--
r-
-25 ....
,........
',...................
:
-- ~ .......
~ ~/
-.........
.
-401
.
.
.
.
.
........
/
:..-
L___
,
f
1
J
i
i
-80
-60
-40
-20
0
20
angle
.~.,_",
40
i
i
60
80
(degrees)
FIG. 4. A r r a y response of (a) different tapers, (b) M B a n d sinusoidal tapers.
,
48
M. G. Amin and N. Bhalla
(a)
(b) 40 i
30
3O
20
iiiii .... i!
2O m
-o 10i 0
-I0
........
-90
-50
", ,---
~. . . . .
012 DFT
50
_,oiiiii!iiiiiii 90
-90
-50
012
50
90
MB
(d)
(c) 4O
40 ......
i
30
_
20
20
10 . . . . . . ' . . . . . .
1
---i- ---i
.....
. . . . . . . . . .
i
i i
0 -10 I
-90
-50
012 Sinusoidal
50
90
-90
-50
0 12
50
90
DPSS
FIG. 5. Beamspace MUSIC for uncorrelated sources: (a) DFT, (b) MB, (c) sinusoidal, (d) DPSS. that can be drawn is that the sinusoidal tapers are good approximation of MB tapers, since they work equally well. VII. Conclusions In this paper, we proposed the use of minimum bias (MB) tapers, introduced for temporal signal analysis, as spatial preprocessing filters. With MB tapers, a desired bandwidth can be achieved by simply adding or removing tapers, obviating the need for recomputations. The MB tapers can be well approximated by sinusoidal tapers, which have a simple closed form analytical expression. The properties and behavior of the sinusoidal tapers were delineated. An approximate expression for the errorband, i.e. the deviation between the maximum and the minimum values in the passband was derived. The application of the MB tapers to high resolution direction finding eigenstructure methods, specifically the MUSIC technique, was explained and elucidated by a set of simulations. For comparison of stopband attenuation provided by various tapers, the DFT, DPSS, MB and the sinusoidal tapers were employed for beamspace MUSIC and applied to a uniform linear array. The results show that the D F T beamspace MUSIC has the worst performance, whereas the other three sets of tapers perform equally well. This indicates that the MB tapers can be closely approximated by the sinusoidal tapers and, both perform similar to the well-known and extensively used DPSS tapers.
Minimum Bias Spatial Filters (a)
49
(b)
20~
v - - ~ - -
20 . . . . .
i
;' i
i
i
i
t
10 . . . .
. . . .
t~ "O
10
1-t
0
.........
"0
i r b
,, -10
-10 -90
-50
0 15
50
90
-90
-50
0 15
50
90
MB
DFT (c)
(d) i i
i i
i i
20
10
10 m 0
° _,o°
-10 -90
-
iiiiiiiiiiii
"O
-50
0 15
Sinusoidal
50
90
-90
-50
0 15
50
90
DPSS
FIG. 6. Beamspace MUSIC in highly coherent scenario: (a) DFT, (b) MB, (c) sinusoidal, (d) DPSS. The crucial aspect of fast computations when utilizing the proposed set of tapers was explored. It was shown that the expression for the sinusoidal tapers resembles the one for the Discrete Sinusoidal Transform (DST). Since substantial endeavors have been made to develop fast computational algorithms for the DST, the sinusoidal tapers based beamspace method is amenable to efficient implementations and has a clear computational advantage over that of the DPSS.
References (1) F. Haber and M. Zoltowski, "Spatial spectrum estimation in a coherent signal environment using an array in motion", IEEE Trans. Antennas Propaq., Vol. AP-34, No. 3, pp. 301310, 1986. (2) T. J. Shan, M. Wax and T. Kailath, "On spatial smoothing for direction-of-arrival estimation of coherent signals", IEEE Trans. Aeoust., Speech, Signal Processing, Vol. ASSP33, No. 4, pp. 806-811, 1985. (3) G. Bienvenu and L. Kopp, "Decreasing high resolution method sensitivity by conventional beam former preprocessing", "Proc. IEEE Int. Conf. Acoust., Speech, Signal Processing", pp. 33.2.1-33.2.4, 1984. (4) C. L. Byrne and A. K. Steele, "Sector-focused stability for high resolution array processing", "Proc. IEEE Int. Conf. Acoust., Speech, Signal Processing", pp. 234~2343, 1987.
M. G. Amin and N. Bhalla
50
(5) P. Forster and G. Vezzosi, "Application of spheroidal sequences to array processing", "Proc. IEEE Int. Conf. Acoust., Speech, Signal Processing", pp. 2267-2271, 1987. (6) H. Lee and M. Wengrovitz, "Resolution threshold of beamspace MUSIC for two closely spaced emitters", IEEE Trans. Acoust., Speech, Signal Processing, Vol. 38, No. 9, pp. 1545-1559, 1990. (7) X. L. Xu and K. M. Buckley, "Statistical performance comparison of MUSIC in element space and beam-space", "Proc. IEEE Int. Conf. Acoust., Speech, Signal Processing", pp. 2124-2127, 1989. (8) C. Zhou, F. Haber and Q. Shi, "On spatial filtering for angle of arrival estimation in a scattering environment by eigendecomposition-based methods", IEEE Trans. Antennas Propag., Vol. 38, No. 6, pp. 931-934, 1990. (9) D. Slepain, "Prolate spheroidal wave functions, Fourier analysis and uncertainty: V: The discrete case", BellSyst. Tech. J., Vol. 5, pp. 1371-1429, 1978. (10) K. Riedel and A. Sidorenko, "Minimum bias multiple tapering spectral estimation", IEEE Trans. Signal Processing, Vol. 43, No. 1, pp. 188-195, 1995. (11) P. S. Kumar ad K. M. Prabhu, "A set of new fast algorithms for DCTs and DSTs", Int. J. Electronics, Vol. 76, No. 4, pp. 553-568, 1994. (12) Z. Wang, "Fast algorithms for the discrete W transform and for the discrete Fourier transform", IEEE Trans. Acoust., Speech, Signal Processing, Vol. 32, No. 8, pp. 803-816, 1984. (13) G. H. Golub and C. F. Van Loan, "Matrix Computations", Second Edn, The Johns Hopkins University Press, 1989. (14) A. P. Prudnikov et al. "Integrals and Series", Vol. 1, Gordon & Breach, 1988. (15) M. Abramowitz and I. Stegun, "Handbook of Mathematical Functions", Dover Publications, NY, 1970.
Appendix A The product of sinusoidal tapers is written in the equivalent exponential form to facilitate summation over k. The summation in (7) will give four terms represented by: e +_J[nn(kl±k2)/(N+
I)]
Le+-Jt"(k'±k2)/(N+l)l_1 J
n=l
L
J (A1)
For even values of (kz +k2), the above expression yields values of - 1. Hence the terms in the summation cancel to produce a zero. For odd values of (k~ _+k2), the expression takes the form: [~ + e-+/~"(k'±k2'"(N+' __e±j[.(k,+_kz)/(N+l)lJ
=
Fz~(kl + k : ) l
- j c ° t [_2 ~
j
The terms in the summation will be of opposing signs, with the same magnitude and will cancel out each other to give a spatial null. To check normality, set kl = k2 in (7):
2(N+I)
l.=,
=
2 ' A ' + l " { - 2 - 2 N } = 1 (A2)
Minimum Bias Spatial Filters
51
Appendix B Equation (5) shows that the frequency response of the kth taper has two terms of the form sin nNx/sin nx. This structure implies that the sidelobes of one term will have a negligible effect on the main lobe of the other, e.g. the value of one term at the center of the main lobe of the other (where the height is N), is ( - 1 ) k+~. For k = 1, the two components are quite close and hence the sidelobes have higher contributions to the other component, but because the two terms are additive, the result is a nonzero response. The main lobe of one component of the kth taper is centered at k/2(N+ 1) while that of the adjacent tapes is a ( k + I ) / 2 ( N + 1) and their midpoint is ( k + 1/2)/2(N+ 1). The first null of a component of the taper is 1IN from its center, which is farther away from the midpoint of the two tapers, which is 1/4(N+ 1) away from the center. Hence the main lobe of each taper component provides nonzero array response before it intersects with the main lobe of the next one. Consequently, a set of contiguous tapers will provide a nonzero passband. To prove the second property, consider N
sin(nxD eaM~"- l)sinO R
Tna(0) ~ T,~ a(0)
.=, C N
(B1)
~ (n - 1) cos 0 sin(nxk) esu~.
Utt
I )sinO
n=l
where C is a constant, xk = 7tk/(N+ 1), k = 1, 2, . . . , K and M = 2~d/2. For the second property to be true, the summations in (B1) for different ks should yield the same R, which does not happen in this case due to the effect of the terms ( n - 1) and cos 0.
Appendix C The expression for the second derivative of the Fourier transform of sum of K tapers, at f = 0 is: x nk 47~2 xk nk ~ cot 2 cosec 2 21t2 (N+ 1) Y ( - 1)k cot 2 - - + ~ k-% 2 ( N + 1) ~ / + I k = odd 2 ( N + 1) 2 ( N + 1) =
2/r2 & 2 nk F L cot - - / c o s e c N + 1 k=, 2(N+l) L
2
rck { 2 rck )7 + (-- 1) k ( N + 1) 2 2 ( N + 1) _ - c o s e c ~ N + 1);J
(C1)
Here we employ approximation sin x ~ x and simplify the expression further to: 4
8(U+ 1)3 Z
= 8 ( N + 1 ) 3 ~=, ~ ~ ( ~( --l ) k
+ ~ 4: k ~ - - ( - - l ) k4~
}
(C2)
The above expression is solved with the help of the following formulae (14): (--1)k k~l
k~
(--1) m - - ( m _ l ) ! [ f l (m °(1)-(-1)xfl(m-*~(K+l)]
1 -1 m k=,/d' -- ( - - ) ! [0('~-'~(1)-0(m-'~(K+l)]
re=l,2 ....
m = 1,2 . . . .
(C3)
(C4)
M. G. Amin and N. Bhalla
52
fl<")(z) = 1/(2 "+ J){~b"[(z+ 1 ) / 2 ] - ~"(z/2)}, where ~,(z) is the psi or digamma function, defined as ~b(z) = (d/dz)[ln F(z)] = F'(z)/F(z), F(z) being the gamma function. For a given value of K, the solution for (C4) is obtained and it is seen that it takes a negative value for odd K and a positive one for even K, signifying a maxima and a minima, respectively. The recurrence formulae for the psi function (15) are useful in computations. For example, when K = 3,
A
' 1
1
'
~j
Similarly Z3= j [ ( - 1)k/k 4] = --0.9498 and Z3= ~(1/k 4) = 1.0748. Putting these values in (C2) will give the value inside the summation as - 0 . 0 4 0 5 , which implies a maximum.
Appendix D The sum of squared values of magnitude of Fourier response of first K tapers is given by: 1 2 ( N + 1)
x sin2 ( N + l)r~f- - -
sin 2 N + ~
=
.k
k~'
(
~2
(DI)
sin2(nJ) --sin2 \ 2 ( N + 1)]J
After approximations of the kind sin x ~ x for small x, the value atfp is:
+ 2(N+l~,~t~)
1
~<,<~-p~)~
(°2)
{[2(N+ i)}2J For even ( k - p ) , the expression vanishes. For the cases when ( k - p ) is odd, it can be simplified to give (16).