Source Localization and Sensing: A Nonparametric Iterative Adaptive Approach Based on Weighted Least Squares

Array processing is widely used in sensing applications for estimating the locations and waveforms of the sources in a given field. In the absence of a large number of snapshots, which is the case in numerous practical applications, such as underwater array processing, it becomes challenging to estimate the source parameters accurately. This paper presents a nonparametric and hyperparameter, free-weighted, least squares-based iterative adaptive approach for amplitude and phase estimation (IAA-APES) in array processing. IAA-APES can work well with few snapshots (even one), uncorrelated, partially correlated, and coherent sources, and arbitrary array geometries. IAA-APES is extended to give sparse results via a model-order selection tool, the Bayesian information criterion (BIC). Moreover, it is shown that further improvements in resolution and accuracy can be achieved by applying the parametric relaxation-based cyclic approach (RELAX) to refine the IAA-APES&BIC estimates if desired. IAA-APES can also be applied to active sensing applications, including single-input single-output (SISO) radar/sonar range-Doppler imaging and multi-input single-output (MISO) channel estimation for communications. Simulation results are presented to evaluate the performance of IAA-APES for all of these applications, and IAA-APES is shown to outperform a number of existing approaches.


I. INTRODUCTION
The goal of array processing is to estimate the locations and waveforms of sources by combining the received data from multiple sensors so that the desired signal is enhanced, while the unwanted signals, such as interference and noise, are suppressed. In active sensing applications such as radar/sonar range-Doppler imaging, the aim is to find targets present in a region of interest. For channel estimation in communications, the aim is to estimate the non-zero channel taps and their Doppler shifts, which are then fed to the subsequent equalizer for symbol detection. Low numbers of snapshots and low signal-to-noise ratios (SNR) are among the many challenges that array processing systems frequently face. Another challenge is the presence of nearby sources, in terms of location or Doppler, since closely spaced sources are harder to discriminate. Herein both passive array processing and active sensing are considered.
The most basic approach to array processing is the classical delay-and-sum (DAS) method, in which the received signal from each sensor is weighted and delayed so as to focus on different points in space. However, this method suffers from low resolution and high sidelobe levels. There is a vast amount of literature on methods that provide superior performance over the DAS approach when certain assumptions are met [1]. The well-known standard Capon beamformer (SCB) [2] and multiple signal classification (MUSIC) [3,4] methods provide superresolution when the sources are uncorrelated and the number of snapshots is high. Many extensions to these methods have been proposed to deal with modelling errors, such as steering vector mismatches. (See, e.g., [5]- [9].) However, none of these methods is able to cope with very low snapshot numbers, coherent or highly correlated sources, or severe noise.
Only a few snapshots are available when the environment being sensed by the array is stationary for a short duration of time. Moreover, to avoid smearing, i.e., losing resolution because of wide main beamwidths in the array response, averaging can only be done over a small bandwidth [10,11]. Therefore the number of available snapshots, which is directly related to the time-bandwidth product, can be very small, sometimes as small as 3, for applications such as underwater array processing. Furthermore, as discussed later, the data models for single-input single-output (SISO) radar/sonar range-Doppler imaging and multi-input single-output (MISO) channel estimation in communication problems are similar to the model used in array processing with a single snapshot and an arbitrary array geometry (dictated by the probing waveforms).
The array processing problem has been carried into the sparse signal representation area by noticing that the number of actual sources is usually much smaller TABLE I Notation used in the Text k ¢ k 0`0 -norm k ¢ k 1`1 -norm k ¢ k 2`2 -norm k ¢ k F Frobenius norm of a matrix the Hadamard (elementwise) matrix product tr (¢) trace of a matrix (¢) T transpose of a vector or matrix (¢) H conjugate transpose of a vector or matrix than the number of potential source points that can be considered [12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29]. Sparsity-based techniques have also been used in spectral estimation [30], image processing, and array design [31][32][33], among many other application areas. Sparse signal representation algorithms can deal with a few snapshots (even one). However, they may require large computation times and the fine tuning of one or more hyperparameters. Application areas, ranging from passive source localization to active radar/sonar range-Doppler imaging and channel estimation for communications are addressed here. Section II considers passive sensing applications and describes a weighted least squares-based iterative adaptive approach for amplitude and phase estimation (IAA-APES). The algorithm is named IAA-APES herein since its derivation resembles that of the amplitude and phase estimation (APES) algorithm [34][35][36]. (See Section IIC.) This name also distinguishes it from the maximum likelihood (ML)-based iterative adaptive approach (IAA-ML) discussed in the Appendix to further motivate IAA-APES. Using the Bayesian information criterion (BIC) [37,38], IAA-APES can be extended to yield point source estimates. (This approach is referred to as IAA-APES&BIC.) Next, a parametric relaxation-based cyclic approach, namely RELAX [39,40], is discussed as a way to further refine the results of IAA-APES&BIC. (This approach is referred to as IAA-APES&RELAX.) Section III presents the data models for the aforementioned active sensing applications and emphasizes the similarities to the passive sensing case. IAA-APES is evaluated via comprehensive simulations in Section IV, and the IAA-APES' performance is compared with that of a number of existing approaches.
Notation: We denote vectors and matrices by boldface lowercase and boldface uppercase letters, respectively. The kth component of a vector x is written as x k . The kth diagonal element of a matrix P is written as P k . See Table I for other symbols and their meanings.

II. PASSIVE SENSING
In passive sensing applications such as aeroacoustic noise measurements [41] and underwater acoustic measurements [42], an array of sensors is used to estimate the desired source characteristics. This section first introduces the data model for such applications and then the IAA-APES algorithm in this model.

A. Data Model
Consider the wavefield generated by K sources located at μ, where μ ¢ =[μ 1 , μ 2 , :::, μ K ] and μ k 1 are the location parameters of the kth signal, k = 1,:::, K. In the narrowband, multi-snapshot case, the M £ 1 array output vector of an M element array in the presence of additive noise can be represented as [5,43] y(n) = A(μ)s(n) + e(n), n = 1,:::, N where N is the number of snapshots, A(μ) is the a(μ 2 ), :::, a(μ K )], and s(n) ¢ =[s 1 (n), s 2 (n), :::, s K (n)] T , n = 1,:::, N, is the source waveform vector at time n. The array steering matrix has different expressions, depending on the array geometry and on whether the source is in the near-field or far-field of the array. For instance, the steering vector corresponding to the kth source for a far-field linear array (where μ k represents the impinging angle of source k in this case) is given by a(μ k ) = [e ¡j(2¼f=c 0 )x 1 cos(μ k ) , :::, e ¡j(2¼f=c 0 )x M cos(μ k ) ] T (2) where f is the center frequency, c 0 is the wave propagation velocity, and x m is the position of the mth sensor, m = 1,:::, M. (See Fig. 1.) Note that for a planar array or near-field sources, only the expressions for the steering vectors have to be modified; the algorithms that are presented can be applied without any modifications, since a(μ) is assumed to be a known function of μ. The number of sources, K, is usually unknown; hence, here, K is considered to be the number of scanning points in the region. In other words every point of a predefined grid that covers the region of interest is considered as a potential source whose power is estimated. Consequently, K will be much larger than the actual number of sources present, and only a few components of fs(n)g will be non-zero. This is the main reason why sparse algorithms can be used in array processing applications.

B. Related Work
We focus our attention on array processing algorithms exploiting sparsity, which have gained noticeable interest recently. Sparse signal representation aims at finding the sparsest s, such that y = As is satisfied, i.e., to minimize ksk 0 , such that y = As, where A is known and y is measured. The problem in its original form is a combinatorial problem and is nondeterministic polynomial-time (NP) hard, making it impractical [12]. Fortunately, when s is sufficiently sparse [12][13][14][15], ksk 0 can be replaced by ksk 1 , which leads to a convex optimization problem that can be solved much more easily by using, for instance, the least absolute shrinkage and selection operator (LASSO) [16] or basis pursuit (BP) [17] algorithms. Alternatively, the focal underdetermined system solution (FOCUSS) algorithm [18], which is derived using Lagrange multipliers, can be used to iteratively solve the sparse problem. A Bayesian approach, such as sparse Bayesian learning (SBL) [19,20] or the approach in [21], can also be used to estimate s. The algorithms mentioned up to this point are for the single-snapshot case only. The extensions of FOCUSS and SBL to the multiple-snapshot case are M-FOCUSS [22] and M-SBL [23], respectively. Another algorithm, the`1-SVD (singular value decomposition) algorithm [24,25] is similar to BP or LASSO, but this algorithm can work with multiple snapshots. (For the single-snapshot case,`1-SVD becomes a LASSO and BP type of method.) M-FOCUSS requires the tuning of two hyperparameters, which might affect the performance of the algorithm significantly.`1-SVD requires the tuning of a hyperparameter and an estimate for the number of sources. Moreover, implementing`1-SVD requires convex optimization software, such as SeDuMi [44]. M-SBL does not require any hyperparameters. However, M-SBL converges quite slowly in its original form [19,20,23].
Besides the above algorithms, which are the focus of our attention in the numerical examples, there are other sparsity based approaches worth mentioning. Reference [26] adds an additional spatial sparsity regularizing term (an`2-norm constraint) to the`1-norm constraint, and it minimizes a cost function similar to that of`1-SVD. However, this method has two hyperparameters, assumes that the source waveforms can be represented by a sparse basis, and has high computational complexity [26]. Fuchs [27,28] uses a sparsity-constrained deconvolution approach that assumes the sources are uncorrelated and that the number of snapshots is large. The sparsity-constrained solution is obtained with a LASSO or BP type of algorithm. Reference [29] introduces two hyperparameter free deconvolution algorithms exploiting sparsity: a sparsity-based extension to the deconvolution approach for the mapping of acoustic sources (DAMAS) [45] (which is similar to [27] and [28] and widely used in practice) and a sparsity based covariance-matrix fitting approach. Extensions to the correlated source case are also provided. However, the methods in [29] are based on the sample covariance matrix, and hence, these methods do not work well with a limited number of snapshots.

C. IAA-APES
IAA-APES is a data-dependent, nonparametric algorithm based on a weighted least squares (WLS) approach. Let P be a K £ K diagonal matrix, whose diagonal contains the power at each angle on the scanning grid. Then P can be expressed as js k (n)j 2 , k = 1,:::, K: Furthermore, define the interference (signals at angles other than the angle of current interest μ k ) and noise covariance matrix Q(μ k ) to be where R ¢ = A(μ)PA H (μ). Then the WLS cost function is given by (see, e.g., [34]- [36] and [43]) where kxk 2 x and s k (n) represents the signal waveform at angle μ k and at time n. Minimizing (5) with respect to s k (n), n = 1,:::, N, yieldŝ , n = 1,:::, N: (6) This looks like the result that would be obtained by employing APES [34][35][36], but it is actually different than APES since APES obtains Q(μ k ) from the data by forming subapertures, while IAA-APES computes Q(μ k ), as in (4). Moreover, IAA-APES is iterative, but APES is not, and APES cannot be used with arbitrary array geometries. Using (4) and the matrix inversion lemma, (6) can be written aŝ , n = 1,:::, N: This avoids the computation of Q ¡1 (μ k ) for each scanning point, i.e., K times. Moreover, fŝ k (n)g can be computed in parallel for each scanning point, which makes IAA-APES amenable to implementation on parallel hardware. IAA-APES is summarized in Table II. Since IAA-APES requires R, which itself depends on the unknown signal powers, it has to be implemented as an iteration. The initialization is done by a standard DAS beamformer. Our empirical experience is that IAA-APES does not provide significant improvements in performance after about 15 iterations. In IAA-APES, P and, hence R, are obtained from the signal estimates of the previous iteration and not from the snapshots as in conventional adaptive beamforming algorithms, such as SCB, which fails to work properly with coherent or highly correlated sources, or few snapshots [43]. The Appendix provides an alternative derivation of IAA-APES, based on the ML principle. IAA-APES is shown to be an approximation of the IAA-ML algorithm, which is locally convergent due to cyclically maximizing the likelihood function. Hence, the analysis in the Appendix provides an approximate calculation that shows the local convergence of IAA-APES.
Because P is assumed to be a diagonal matrix, one degree of freedom (DOF) is needed to suppress one interfering source. This may lead to a larger than necessary reduction of the number of degrees of freedom when some of the interfering sources are coherent since the cancellation of multiple coherent interfering sources would require only one DOF if the correct structure of P were known. However, we do not assume that the true structure of P is known. Moreover, it is the diagonal structure of P assumed by IAA-APES that makes the algorithm work properly even for low number of snapshot cases and coherent sources.

D. IAA-APES&BIC
In many applications it is desirable to obtain point estimates rather than a continuous spatial estimate. To achieve this sparsity we incorporate a model-order selection tool, i.e., the BIC [37,38], into IAA-APES. Let P denote a set containing the indices of the peaks selected from the IAA-APES spatial power spectrum estimate. Also let I denote the set of the indices of the peaks selected by the BIC algorithm so far. The IAA-APES&BIC algorithm works as follows: first the peak, from the set P, giving the minimum BIC is selected. Then the second peak, from the set P ¡ I, which together with the first peak gives the minimum BIC, is selected, and so on, until the BIC value does not decrease anymore. 2 The IAA-APES&BIC algorithm is summarized in Table III. BIC i (´) is calculated as follows (see [38]) where´= jIj + 1, jIj denotes the size of the set I, i is the index of the current peak under consideration, and fŝ j (n)g N n=1 is the IAA-APES signal waveform estimate corresponding to angle μ j , j 2 fI [ ig. Note that the second term on the right side of (8) does not matter to peak selection; it matters only when (8) is used to select the number of peaks to retain.

E. IAA-APES&RELAX
RELAX [39,40] is a parametric cyclic algorithm that requires an estimate of the number of sources in the field. The results of IAA-APES&BIC can be used to provide a good initial estimate for the last step of RELAX. This approach is outlined in Table IV. The RELAX iterations can be terminated when the norm of the difference between two consecutive estimates is smaller than a certain threshold (5 £ 10 ¡4 for the examples considered herein). The maximization , n = 1,:::, N end for until (convergence) step of the algorithm can be implemented without much computational effort by means of a fine search only around the peak values of the IAA-APES&BIC result. 3 RELAX can be useful in estimating off-grid sources accurately and for further improving the IAA-APES waveform and angle estimates.

III. ACTIVE SENSING
In active sensing applications, besides the receiver, there are also one or more transmitters. The radar/sonar range-Doppler imaging problem for a SISO system is first investigated in this section. Then the channel estimation problem for MISO communications is discussed.

A. Range-Doppler Imaging
Pulse compression refers to the process of transmitting a modulated pulse and then matched filtering the returned signal, which arrives at the antenna altered by complex coefficients that bear target information [48].
where M is the pulse length. Let  3 The fine search can be implemented efficiently using the fast Fourier transform (FFT) [39,40] or a derivative-free uphill search method, such as the Nelder-Mead algorithm [46,47]. Note that the latter method is available in the MATLAB optimization toolbox with the name of "fminsearch." the Doppler interval of interest is divided into L bins, and the Doppler frequency for the lth Doppler bin is denoted as ! l . Then the M samples of the received signal that is temporally aligned with the return from the range bin of current interest r (see Fig. 2) can be represented by for r = 1,:::, R, where ® r,l (which is proportional to the complex "voltage" radar-cross section (RCS) of the corresponding target) denotes the complex amplitude of the returned signal from the range bin of current interest r and the lth Doppler bin f® m+r,l g denote the complex amplitudes of the returned signals from the adjacent range bins, and e r denotes the noise.
The reflections from nearby range bins are considered to be clutter. The M £ M shift matrix J m takes into account the fact that the clutter returns from adjacent range bins need different propagation times to reach the radar/sonar receiver: for m = 0,:::, M ¡ 1. Equation (12) can be written as Note that (12) is similar to the data model for passive sensing arrays, see (1), but with a single snapshot. The IAA-APES estimate at iteration number i 2 f1, 2, :::g becomes® for l = 1,:::, L, r = 1,:::, R, where Here IAA-APES is applied in a slightly different manner than in the passive array processing case. When ® r is to be estimated, the previous values of f® r+m g m=M¡1 m=¡M+1 are used to estimate R(r). However, only ® r is updated using R(r) and y r . When all f® r g R r=1 are updated in this way, we advance to the next IAA-APES iteration. IAA-APES is initialized by matched filtering, matched to each range and Doppler bin, to obtain the initial f® r,l g.
2) Related Work: Matched filtering has been widely used for pulse compression because it gives optimal signal-power-to-output-power performance in the presence of a single target and white noise. However, in practical radar/sonar systems, matched filter performance is far from desirable since the pulse compression problem is usually clutter limited rather than noise limited (see, e.g., [49] and [50]). Many data independent (see, e.g., [49]- [54]) and data-adaptive (see, e.g., [48] and [55]) approaches have been proposed to achieve improved pulse compression.
Data independent approaches can be designed off-line, and hence, these approaches are convenient for real-time implementations in practical systems. Although receive filters based on data-independent instrumental variables (IV) can be used to achieve excellent pulse compression for the negligible Doppler case, their performance is unsatisfactory in the nonnegligible Doppler case because of high sidelobe level problems [49]. On the other hand data-adaptive approaches result in better performance but at the cost of implementation complexity. Two important data-adaptive methods for pulse compression are the adaptive pulse compression (APC) [48] and the Doppler-compensated APC (DC-APC) [55] algorithms. APC is an iterative minimum mean-squared error (MMSE) based data adaptive approach, and DC-APC is the extension of APC to the nonnegligible Doppler case. However, an IV approach [49] may be preferred to APC since the sidelobe level of the former method can be made arbitrarily low and because the filter coefficients can be computed off-line, whereas APC updates the filter coefficients iteratively and adaptively. In the nonnegligible Doppler case, however, adaptive approaches can perform much better. DC-APC is purported to work with at most one target per range bin, and the method used in [55] to estimate the Doppler values is not appealing from a performance viewpoint. On the other hand, IAA-APES can work with multiple Doppler targets located at the same range bin, and IAA-APES estimates the Doppler values in a robust manner. Moreover, APC and DC-APC require the tuning of hyperparameters, whereas IAA-APES is hyperparameter free. Some other differences are discussed below.
DC-APC assumes that there is at most one target per range bin. Let f® r¡M+1 , :::, ® r+M¡1 g and f! r¡M+1 , :::,! r+M¡1 g denote the complex amplitudes and the Doppler frequencies, respectively, of the targets in the corresponding range bin r. DC-APC estimates the target parameters iteratively as follows for r = 1,:::, R, where Q r is the true noise covariance matrix of e r , which is assumed to be known, and i 2 f1, 2, :::g represents the current iteration number. The initial estimates of f® r+m g M¡1 m=¡M+1 are obtained by using a standard matched filter that neglects the Doppler effect [55]. DC-APC also requires the estimates f! r+m g M¡1 m=¡M+1 of the target Doppler frequencies at each iteration. The Doppler frequency estimation approaches suggested in [55] are ad-hoc and not very accurate, especially for large Doppler shifts. This limitation of DC-APC, though, can be easily corrected by replacing! r in (20) with all possible Doppler frequencies f! l g to deal with multiple targets per range bin and to form range-Doppler images.
The DC-APC iterations defined in (20)-(21) result in the numerical ill-conditioning of R (i¡1) (r). To mitigate this problem, [48] suggests using j® (i¡1) r j ± in (20) instead of j® (i¡1) r j 2 , and ¾ ± instead of ¾ 2 , where 0 · ± · 2 and where the noise e r is assumed to be white with a known variance ¾ 2 . This approach requires the delicate tuning of ± (at each iteration), but [48] and [55] do not provide a clear guideline on how to do this. We remark that if the noise covariance matrix Q r in (21) is set to zero, then APC/DC-APC becomes identical to FOCUSS, with the sparsity parameter p = 0 and the regularization parameteŗ = 0. (See, e.g., (16) in [22]). In addition, with the introduction of ±, APC/DC-APC is still identical to FOCUSS, now with 1 ¡ p=2 = ±=2 and¸= ¾ ± . We also note that both APC/DC-APC and FOCUSS are related to the approaches in [21], [30], and [56]. It is interesting to note that many different ideas lead to the same result.
IAA-APES assumes that multiple targets can exist within the same range bin and that f® r,l g is calculated as in (18), where the term in the denominator can be viewed as the current estimate of j® (i¡1) r j 2 used in (20). This would be best estimate for j® (i¡1) r j 2 , obtained by applying the minimum variance distortionless criterion, had the true covariance matrix been known [5]. One more advantage of IAA-APES is that the parameter estimates are inherently unbiased for the signal of interest, whereas this is not the case in DC-APC (or FOCUSS). Furthermore, in IAA-APES, the statistical properties of e r are assumed unknown and are taken into account implicitly as false targets in the current and adjacent range and Doppler bins; see (19). No matrix inversion problems were encountered with IAA-APES during our numerical simulations because of the unbiasedness property of IAA-APES.

B. Channel Estimation
The purpose of channel estimation in communications is to provide the subsequent equalizer with an accurate channel estimate so that the transmitted signals can be recovered successfully at the receiver side [57,58]. SISO channel estimation for communications problem is similar to SISO radar/sonar range-Doppler imaging, with the only difference being that the delay for the former is due to single-trip propagation, while the delay for the latter results from round-trip propagation. denote the ith transmitted pulse, and let L denote the total number of Doppler bins. Then similar to (12), the received signal that is temporally aligned with the return from the rth tap, r = 1,:::, R, with R denoting the total number of taps of the channel, can be written as where J m was defined in (13), and d(! l ) was defined in (11). When determining the parameters of the ith channel f® r,l (i)g, the term P L l=1 ® r,l (i)s i (! l ) is the signal term, and all other terms in (23) are considered as clutter and noise.
Consequently, IAA-APES can be applied directly to (23) to estimate f® r,l (i)g. Like in radar/sonar range compression problems, to determine f® r,l (i)g with high accuracy, the transmitted pulses need to have both good auto and cross-correlation properties [59,60].
2) Related Work: Most of the demodulation algorithms used in practice rely heavily on the accurate estimation of the channel impulse response. As in many other applications, sparse signal estimation approaches have also been proposed in this context. The main motivation for this is that underwater communications [61][62][63] and wireless channels are appropriately modelled as sparse channels consisting of only a few non-zero taps [64]. The existing approaches we evaluate for sparse channel estimation, besides IAA-APES, include the matching pursuit (MP), orthogonal matching pursuit (OMP) [65][66][67], and least squares matching pursuit (LSMP) [68] algorithms, which have been used for sparse channel estimation and equalization in many applications [69][70][71][72]. It is difficult to determine the stopping criterion when using matching pursuit algorithms, and user intervention is needed.

IV. NUMERICAL EXAMPLES
We evaluate the performance of IAA-APES and compare it with various alternative methods in this section. We first focus on passive sensing applications and then shift our attention to active sensing applications.

A. Passive Sensing Examples
This subsection investigates the performance of IAA-APES, M-FOCUSS, M-SBL and`1-SVD for various passive sensing scenarios. Unless noted otherwise, M-FOCUSS is implemented by setting the sparsity parameter p = 0:8 and fine tuning the regularization parameter¸to get the best results. 4 1 -SVD is implemented by assuming that the number of sources is known a priori and by fine tuning the hyperparameter. Since the fine tuning of the hyperparameters assume knowledge of the true source parameters, the so obtained results of M-FOCUSS and 1 -SVD are impractical. M-SBL is implemented by using the alternative update method for the parameters as described in [19], [20], and [23]. For all approaches considered, the scanning grid is uniform in the range from 1 ± to 180 ± , with 1 ± increment between adjacent grid points, unless noted otherwise.
We consider a uniform linear array with M = 12 sensors and half-wavelength interelement spacing.
The far-field narrowband signal waveforms and the additive noise signals are assumed to be circularly symmetric, independent identically distributed (i.i.d.) complex Gaussian random processes with zero mean and variance ¾ 2 , which is varied to obtain various SNR values. Furthermore, each signal waveform is normalized such that (1=N) P N n=1 js k (n)j 2 = P k , k = 1,:::, K 0 , for a given P k value, where K 0 denotes the true number of sources. SNR is defined as 10 log 10 (P k =¾ 2 ), k = 1,:::, K 0 , in decibels (dB), where ¾ 2 is the noise variance. According to the comments made in the Introduction, we consider very low snapshot cases, viz. N = 3 and N = 1.
First we consider three uncorrelated sources at 60 ± , 82 ± , and 90 ± , with 5 dB, 10 dB, and 10 dB powers, respectively, and with N = 3. The noise power is 0 dB, which results in a minimum SNR of 5 dB. Fig. 3 shows the power and location estimates of the algorithms. The circles and the vertical dotted lines that align with these circles represent the true source locations and powers, and the results of 10 Monte-Carlo trials are shown in each plot. DAS clearly suffers from smearing and leakage. IAA-APES provides a much better result than DAS, with low sidelobes and peaks at the true source locations. Moreover, IAA-APES&RELAX indicates the number of sources and their locations and powers accurately. As observed in Fig. 3(d)-(f), the sparse algorithms encounter source splitting, location bias, and power underestimation problems. Also it is hard to tell how many sources are present by solely using the spatial estimates of the sparse algorithms.
Next Fig. 4 considers three coherent sources at 60 ± , 80 ± , and 90 ± , with 10 dB power each and with N = 3. The source waveforms are assumed to be identical for all three coherent sources. The noise power is 0 dB, resulting in a 10 dB SNR. Similar to Fig. 3, the circles and the vertical dotted lines represent the true source locations and the powers, and the results of 10 Monte-Carlo trials are shown in each plot. We observe that IAA-APES is able to resolve the sources successfully and that IAA-APES&RELAX provides accurate point estimates. On the other hand DAS fails to resolve the two closely spaced sources. The performances of the sparse algorithms are similar to those of the previous example.
Finally Fig. 5(a) compares the total mean-squared error (MSE) (sum of each individual MSE) of the angle estimates of each algorithm with the Cramer-Rao bound (CRB) [5], and Fig. 5(b) compares the total angle estimation bias (sum of each individual bias' modulus) of each algorithm for varying SNR. Two uncorrelated sources are placed at 77:51 ± and 90:51 ± , and N is set to 1. (The angle values are picked so that they are not on the size 1 ± grid used by IAA-APES.) To calculate the MSE and bias, we consider only the signals with the two largest powers as the estimated signals for FOCUSS and`1-SVD. We run the sparse algorithms, with a fine grid of size 0:01 ± , around the true source locations. The sparse algorithms are run with K = 950 scanning points, whereas IAA-APES is run with the same resolution as before, i.e., K = 180. IAA-APES&RELAX is used to estimate the off-grid sources. M-SBL results are not shown because of the excessive computation time needed. Each point in Fig. 5 is the average of 100 Monte-Carlo trials. We observe that IAA-APES&RELAX has both better variance and bias characteristics than the other methods for lower SNR.
Note that FOCUSS and`1-SVD both have two hyperparameters, and their selection affects the performances of these two algorithms significantly. Moreover, a different parameter setting should be used depending on the SNR, N, and the source structure, i.e., the number of sources, source spacing, source power levels, and correlation levels. We were able to tune the parameters relatively easily in our simulation scenarios, but when the problem dimensions are large and when there is no prior knowledge of the scenario, it becomes difficult to find good hyperparameters. M-SBL, on the other hand, does not require any hyperparameters, but it takes the longest time to converge. These are the main reasons why BP or LASSO, which are the single-snapshot counterpart of 1 -SVD, and SBL and FOCUSS, which are the single snapshot counterparts of M-FOCUSS and M-SBL, are not considered in the active sensing examples below.
Assuming K À M, the complexity of each IAA-APES iteration is O(M 2 K). The complexity of the BIC extension is negligible compared to that of IAA-APES, and the complexity of the RELAX extension depends on how many sources IAA-APES&BIC determines, the termination condition, and the method chosen for the maximization step. In our experiments RELAX usually converged faster than the time needed for IAA-APES to converge. The complexities of M-FOCUSS and M-SBL are also O(M 2 K) per iteration [23]. The complexity of`1-SVD, on the other hand, is O(K 3 K 3 svd ) [25], where K svd is the estimate of the number of sources. In our simulations IAA-APES was always faster than M-FOCUSS, M-SBL, and`1-SVD, especially when the SNR was low. Note that the actual number of iterations required for convergence depends heavily on the specific scenario. For example for the scenario considered in Fig. 5 at SNR = 2 dB, the times required by FOCUSS and`1-SVD are 3.9 and 5.6 s, respectively, whereas IAA-APES&RELAX requires 0.1 s on average. 5

B. Active Sensing Examples
This section considers examples for the SISO radar/sonar and the MISO channel estimation problems.
We consider three moving targets with 5 dB power and with Doppler shifts of ¡30 ± , ¡25 ± , and 15 ± and six moving targets with 25 dB power and Doppler shifts of ¡70 ± , ¡55 ± , ¡10 ± , 10 ± , 20 ± , and 60 ± , as shown in Fig. 6(a). The background noise is assumed to be a circularly symmetric i.i.d. complex Gaussian random process, with mean zero and a variance of 0 dB. The number of range bins is set to 100, and the number of Doppler bins is set to 37. The resulting minimum SNR is 5 dB. Fig. 6 shows that the matched filter smears the targets in both the Doppler and range domains significantly. Figs. 6(c)-(e) show that MP, OMP, and LSMP results are not satisfactory in this example. The MP, OMP, and LSMP algorithms are all terminated manually to give the best performance. 6 Fig. 6(f) shows that IAA-APES provides a much more useful result (the target locations and Doppler frequencies being easily observable) than the matched filter. IAA-APES&BIC can estimate the strong targets accurately and the weak targets reasonably well, as shown in Fig. 6(g). IAA-APES&RELAX helps improve the results of IAA-APES&BIC even further, as shown in Fig. 6(h).  IAA-APES spatial estimate, but BIC does not modify either the range-Doppler values or the power levels estimated by IAA-APES, whereas RELAX improves upon the IAA-APES estimates.
2) Channel Estimation: The performance of MP, OMP, LSMP, and IAA-APES are also investigated for the channel estimation problem encountered in communications. We consider six transmitters, viz. I = 6, and a single receiver. The probing pulses used in the simulations are obtained by using a cyclic algorithm to ensure both good auto and cross-correlation properties; see [59] and [60] for details. The sampling frequency is assumed to be 24 kHz. M is chosen to be 1024 (which yields a 42.7 ms pulse duration). The number of delay taps are set at 80 to yield a 3.3 ms of delay span, and 13 Doppler bins are used to cover the range from ¡30 Hz to 30 Hz with 5 ± resolution. The six channels are assumed to be independent of each other, and they are simulated as circularly symmetric i.i.d. complex Gaussian random variables, with mean zero and a variance decreasing exponentially with increasing delay, specifically ¾ 2 k = e ¡0:03(k¡1) , where ¾ 2 k is the variance of the channel tap with delay index k, k = 1,:::, 80. The noise is simulated as a circularly symmetric i.i.d. complex Gaussian random process, with mean zero and a ¡20 dB variance. The locations of the non-zero taps of the channels are also simulated randomly and by trying to mimic practical channel functions encountered in underwater communications, see, e.g., [69] and [70]. The power levels of the channel taps were in the range from ¡30 dB to 0 dB, which results in a minimum SNR of ¡10 dB. Fig. 7 illustrates the performances of the algorithms for the first channel; IAA-APES&BIC provides the best result. The results for the other channels, which are not shown for conciseness, are similar to that for the first channel.

V. CONCLUSIONS
This paper has presented IAA-APES in array processing applications. IAA-APES is a nonparametric, hyperparameter free algorithm that is designed to work under severe snapshot limitations and for uncorrelated, partially correlated, and coherent sources, as well as for arrays with arbitrary geometries. Because of the similarities between many active sensing applications and passive array processing, IAA-APES can be applied to these cases as well, without any essential modifications. The BIC can be used in conjunction with the IAA-APES algorithm to yield sparse solutions, which are desirable in many applications. Furthermore, the application of the parametric RELAX algorithm to the IAA-APES&BIC results provides further performance improvement. Simulations showed that IAA-APES and its variations outperform the most prominent methods in the literature in the corresponding application areas. IAA-APES is believed to be a viable candidate for practical applications since it does not require any hyperparameters, has a simple formulation, provides superresolution, facilitates parallel processing, and shows good performance.

APPENDIX. AN APPROXIMATE ML INTERPRETATION OF IAA-APES
We derive a locally convergent ML-based iterative adaptive approach called IAA-ML. IAA-ML is similar to IAA-APES in nature, and we show that IAA-APES can be obtained as a simplified version of IAA-ML.
The negative log-likelihood function of fy(n)g N n=1 can be represented as where ¡ = (1=N) P N n=1 y(n)y H (n) is the so-called sample covariance matrix. It is assumed that the received signal is a multivariate, complex Gaussian random vector with zero-mean and covariance matrix R and that the snapshots are i.i.d. [5]. Minimizing (25) is equivalent to minimizing the Stein's loss, ¡ ln jR ¡1 ¡ j + tr(R ¡1 ¡ ) [75,76], which stresses the covariance fitting character of the cost function, but without the i.i.d. Gaussian assumptions needed by the maximum likelihood loss function.
Assume that Q(μ k ) is known and that the signal power at μ k is to be estimated. Using the fact that jI + ABj = jI + BAj and the matrix inversion lemma together with (4), it can be shown that minimizing (25), with respect to P k , is equivalent to minimizing where Q(μ k ) and a(μ k ) have been replaced by Q k and a k , respectively, for notational simplicity. Setting the first derivative of (26), with respect to P k , to zero, i.e., f 0 (P k ) = 0, gives The second derivative of (26), with respect to P k , is a H k Q ¡1 k a k (1 + P k a H k Q ¡1 k a k ) 2 and hence f 00 (P k ) is f 00 (P k ) = (a H k Q ¡1 k a k ) 2 (1 +P k a H k Q ¡1 k a k ) 2 (28) which is strictly positive. This means thatP k is the unique minimizer of f(P k ). In principle,P k may be negative. However, since P k represents power, it should be nonnegative. The minimizer of f(P k ) subject to the constraint P k¸0 isP k = max(0,P k ) (29) sinceP k is the unique minimizer of f(P k ) (wheñ P k < 0, f(0) · f(P k ) for 8P k¸0 ). By using the matrix inversion lemma in (27) to replace Q k by R, we get ComputingP k requires knowledge of P k and R (recall that R = A(μ)PA H (μ), where P is a diagonal matrix with fP k g on its diagonal). Therefore the algorithm must be implemented iteratively; the initialization is done with DAS. IAA-ML will be locally convergent if R is recalculated after each P k is updated because of the cyclical maximization of the likelihood function. Note that (30) can be written as (assumingP k¸0 ) By the properties of SCB, P k ¼ 1=(a H k R ¡1 a k ) [6,43]. Hence an approximate solution of (31) iŝ P k = a H k R ¡1 ¡ R ¡1 a k (a H k R ¡1 a k ) 2 , k = 1,:::, K: Iterating this equation by building R from the latest estimate of fP k g, we get IAA-APES. This approximation has two advantages. First theP k in (32) is guaranteed to be nonnegative, and it alleviates the need for the procedure in (29). Secondly when fP k g is accurate, the difference P k ¡ 1=(a H k R ¡1 a k ) in (31) is small, as discussed above. When fP k g is inaccurate, however, this difference may not be small, and forcing the difference to zero may provide a better estimate of P k . (IAA-ML tends to work well when the snapshot number is large. It is not considered in the numerical examples because of our focus on cases with few snapshots.) Since IAA-APES can be obtained as a close approximation to IAA-ML, which is locally convergent, IAA-APES is expected to also enjoy local convergence. We have never come across an example where IAA-APES did not converge; however, the search for a convergence proof could be a useful direction for future work.

ACKNOWLEDGMENT
We thank Professor Aaron Lanterman for his useful comments on the revisions of a former version of this paper.