{'fraction_non_alphanumeric': 0.03253126089772331, 'fraction_numerical': 0.011209086833059333, 'mean_word_length': 4.321712619300106, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 47, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multisource Self-calibration for Sensor Arrays

Calibration of a sensor array is more involved if the antennas have direction dependent gains and multiple calibrator sources are simultaneously present. We study this case for a sensor array with arbitrary geometry but identical elements, i.e. elements with the same direction dependent gain pattern. A weighted alternating least squares (WALS) algorithm is derived that iteratively solves for the direction independent complex gains of the array elements, their noise powers and their gains in the direction of the calibrator sources. An extension of the problem is the case where the apparent calibrator source locations are unknown, e.g., due to refractive propagation paths. For this case, the WALS method is supplemented with weighted subspace fitting (WSF) direction finding techniques. Using Monte Carlo simulations we demonstrate that both methods are asymptotically statistically efficient and converge within two iterations even in cases of low SNR.


I. INTRODUCTION
In this paper we study the calibration of the direction dependent response of sensor array antennas, excited by simultaneously present calibrator sources. The antenna array has arbitrary geometry but identical antennas. The calibration involves the complex gain of the antenna elements towards each source. The source powers are known but we will allow for small deviations in apparent source locations to account for, e.g., refractive propagation paths. This problem is one of the main challenges currently faced in the field of radio astronomy. For low frequency observations (< 300 MHz) this community is building or developing a number of new instruments, for example the low frequency array (LOFAR) [1], the Murchison wide field array (MWA) [2] and the primeval structure telescope (PaST) [3], which are all large irregular phased arrays. These difficulties arise due to the influence of the ionosphere on the propagation of radio waves, which can qualitatively be categorized into the following four regimes depending on the field of view (FOV) of the individual receptors and the baselines between them [4]: 1) All antennas and all lines of sight sample the same ionospheric delays, thus the ionosphere causes no distortion of the array manifold (small FOV, short baselines). 2) Lines of sight from different antennas sample different ionospheric patches, but all sources in the antenna FOV experience the same delay; the ionosphere thus causes an antenna based gain effect (small FOV, long baselines).
Copyright (c) 2008 IEEE. Personal use of this material is permitted. However, permission to use this material for any other purpose must be obtained from the IEEE by sending a request to pubs-permissions@ieee.org.
This work was supported by the Netherlands Institute for Radio Astronomy 3) All lines of sight towards a single source sample the same ionosphere but the delay differs per source. This regime requires source based direction dependent calibration (large FOV, short baselines, see Fig. 1). 4) The ionospheric delays differ per station and per source (large FOV, long baselines).
The first two regimes can be handled under the selfcalibration assumption used in astronomy that all errors are telescope based thereby reducing the estimation problem to a single direction independent gain factor per telescope [5]. The problem of estimating these direction independent gains based on a single calibrator source is treated in [6], while the multiple source case has been discussed in [7]. The problem of finding a complex gain per antenna per source, the fourth regime, is not tractable without further assumptions which allow to parameterize the behavior of the gains over space, time and/or frequency [8]. In this paper we will focus on the third regime, which is sketched in Fig. 1, thus filling the gap in the available literature. It allows for the calibration of individual closely packed groups of antennas such as a LOFAR station or a subarray of the MWA and PaST telescopes and thus forms a valuable step towards the calibration of the whole array.
We will state our problem in general terms, since the problem plays a role in a range of applications varying from underwater acoustics to antenna arrays. This explains the persistent interest in on-line calibration, autocalibration or selfcalibration of sensor arrays [8]- [12]. In most applications the main driver for studies on array calibration is to improve the DOA estimation accuracy. Many studies in this field therefore try to solve for the DOAs and a number of array parameters. In [9], [13], [14] self-calibration schemes are presented which solve for the direction independent gains and the sensor positions, i.e. the directional response of the sensors is assumed known. Other studies assume a controlled environment to calibrate the array by measuring the calibrator sources one at a time [15], [16] or exploit the array geometry, e.g. a uniform linear array (ULA) having a Toeplitz matrix as array covariance matrix [10], [17]. Weiss and Friedlander [18] have presented a technique for almost fully blind signal estimation. Their work, however, focuses on estimation and separation of source signals, not on characterizing the array itself. We thus feel that the problem at hand also forms an interesting addition to the literature available on sensor array calibration in general.
In the next section we introduce the data model and provide a mathematical formulation of the problem. In Sec. III the Cramèr-Rao bound for this estimation problem is discussed. Section IV presents an alternating least squares (ALS) and a weighted alternating least squares (WALS) approach optimizing subsets of parameters iteratively. These algorithms are validated using Monte Carlo simulations in Sec. V. The conclusions are drawn in Sec. VI.
Notation: Overbar (·) denotes conjugation, the transpose operator is denoted by T , the complex conjugate (Hermitian) transpose by H and the pseudo-inverse by † . An estimated value is denoted by E{·}. ⊙ is the element-wise matrix multiplication (Hadamard product), ⊘ is the element-wise matrix division, ⊗ denotes the Kronecker product and • is used to denote the Khatri-Rao or column-wise Kronecker product of two matrices. diag(·) converts a vector to a diagonal matrix with the vector placed on the main diagonal, vecdiag(·) produces a vector from the elements of the main diagonal of its argument and vec(·) converts a matrix to a vector by stacking the columns of the matrix. We exploit many properties of Kronecker products, including the following (for matrices and vectors of compatible dimensions):

II. DATA MODEL
We consider an array of p elements located at a single site ("Regime 3"). Denote the complex baseband output signal of the ith array element by x i (t) and define the array signal vector We assume the presence of q mutually independent i.i.d. Gaussian signals s k (t) impinging on the array, which are stacked in a q ×1 vector s(t). Likewise the sensor noise signals n i (t) are assumed to be mutually independent i.i.d. Gaussian signals and are stacked in a p × 1 vector n(t). If the narrow band condition holds [19], we can define the q spatial signature vectors a k (p×1), which describe for the kth source the phase delays at each antenna due only to the geometry.
The q sources are considered calibrator sources, thus we assume that their powers, their nominal positions, and the locations of the antennas (hence a k ) are known. Refractive effects caused by ionospheric phase gradients may shift the apparent locations of the sources, thus we will also consider cases where the a k are only parametrically known.
The sensors are assumed to have the same direction dependent gain behavior towards the q source signals received by the array. This can include the antenna pattern and ionospheric phase effects. They are described by gain factors g 0k and are collected in a matrix G 0 = diag([g 01 , g 02 , · · · , g 0q ]). The direction independent gains and phases (individual receiver gains) can be described as γ = [γ 1 , γ 2 , · · · , γ p ] T and φ = [e jφ1 , e jφ2 , · · · , e jφp ] T respectively with corresponding diagonal matrix forms Γ = diag(γ) and Φ = diag(φ). With these definitions, the array signal vector can be described as where A = [a 1 , · · · , a k ] (size p × q) and G = ΓΦ; for later use we also define g = γ ⊙ φ.
The signal is sampled with period T and N sample vectors are stacked into a data matrix s1 , σ 2 s2 , · · · , σ 2 sq ] T and the noise covariance matrix is In this model, Σ s is known. Since G 0 and Σ s are diagonal matrices, we can introduce Since the direction dependent gains g 0k are not known, the real valued elements of σ = [σ 2 1 , σ 2 2 , · · · , σ 2 q ] T are not known either. This implies that our problem is identical to estimating an unknown diagonal signal covariance matrix without any DOA dependent errors. We may thus restate (5) as and solve for G, Σ and Σ n under the assumption that A is known (or parametrically known). The data model described by (7) is commonly used in papers on sensor array calibration (e.g. [6], [20]). Flanagan and Bell [9], Weiss and Friedlander [13] and See [14] effectively use the same model, but focus on position calibration of the array elements and are therefore more explicit on the form of A. If the source positions and locations of the sensors within the array are known, an explicit formula for a k can be used to compute the nominal spatial signature vectors. Estimation of source locations is required to account for the refractive effects produced by an ionospheric phase gradient.
The ith element of the array is located at r i = [x i , y i , z i ] T . These positions can be stacked in a matrix R = [r 1 , r 2 , · · · , r p ] T (size p × 3). The position of the kth source can be denoted as l k = [l k , m k , n k ] T . The source positions can be stacked in a matrix L = [l 1 , l 2 , · · · , l q ] T (size q × 3). The spatial signature matrix A can thus be described by where the exponential function is applied element-wise to its argument. In the remainder of this paper we will specialize to a planar array having z i = 0 for convenience of presentation but without loss of generality. From (7) we observe that G and Σ share a common scalar factor. Therefore we may impose the boundary condition σ 2 1 = 1. To solve the ambiguity in the phase solution of G we will take the first element as phase reference, i.e. φ 1 = 0 is imposed. 1 When solving for source locations a similar problem occurs: a single rotation of all DOA vectors can be compensated by the direction independent gain phase solution. We will therefore fix the position of the first source.
In this paper we will address four related sensor array calibration problems based on this data model. These are summarized below where the parameter vectors adhering to the aforementioned boundary conditions are stated explicitly.

III. CRB ANALYSIS
The Cramèr-Rao bound (CRB) on the error variance for any unbiased estimator states that [22] where C is the covariance matrix of θ, and C is the Fisher Information Matrix (FIM). For Gaussian signals the FIM can be expressed as where F is the Jacobian evaluated at the true values of the parameters, i.e.
1 In [21] it is shown that p i=1 φ i = 0 is the optimal constraint for this problem. This constraint has the disadvantage that the location of the phase reference is not well defined. Furthermore, the choice for the constraint used here simplifies our analysis in combination with the constraints required to uniquely identify the source locations and the apparent source powers.
For the calibration problem defined by the first scenario, the Jacobian can be partitioned in four parts following the structure in θ: where it can be shown that I s is a selection matrix of appropriate size equal to the identity matrix with its first column removed so that the derivatives with respect to φ 1 and σ 2 1 are omitted. In the second scenario the Jacobian can be partitioned in a similar way, the difference being that the expression for F σn given by Eq. (16) should be replaced by In the third and fourth scenarios an additional component F l is added to the Jacobian containing the derivatives of vec(R) with respect to the source position coordinates. This component can be partitioned as these components can be conveniently written as These equations show that the entries of the Jacobian related to derivatives with respect to the land m-coordinates of the sources are proportional to the xand y-coordinates of the array elements respectively. The physical interpretation of this relation is that a plane wave propagating along the coordinate axis of the coordinate to be estimated provides a more useful test signal to estimate the source location than a signal propagating perpendicular to this axis.

A. Generalized least squares formulation
An asymptotically efficient estimate of the model parameters θ can be obtained via the ML formulation. Since all signals are assumed to be i.i.d. Gaussian signals, the derivation is standard and ML parameter estimates for N independent samples are obtained by minimizing the negative log-likelihood function [23] where R(θ) is the model covariance matrix as function of θ and R is the sample covariance matrix N −1 XX H . It does not seem possible to solve this minimization problem in closed form. As discussed in [23] a weighted least squares covariance matching approach is known to lead to estimates that are, for a large number of samples, equivalent to ML estimates and are therefore asymptotically efficient and reach the CRB.
The least squares covariance model fitting problem can be defined as Equivalently, we consider minimization of the cost function . The more general weighted least squares problem is obtained by introducing a weighting matrix W and optimizing The optimal weight is known to be the inverse of the asymptotic covariance of the residuals, E f (θ 0 )f (θ 0 ) H , where θ 0 is the true value of the parameters [23]. The optimal weight for Gaussian sources is thus The Kronecker structure of W opt allows to introduce W c = R − 1 2 and write the weighted least squares (WLS) cost function as As mentioned, this estimator is asymptotically unbiased and asymptotically efficient [23]. We propose to solve the least squares problems for the four cases defined in section II by alternating between least squares solutions to subsets of parameters. The subsets are chosen such that the solution is either available in the literature or can be derived analytically. We will have four subsets of parameters: the complex direction independent gains of the elements g, the apparent source powers σ, the receiver noise powers σ n , and the source locations L. In the following subsections we will develop the unweighted and weighted least squares solutions for these subsets before putting them together to form the alternating least squares (ALS) or weighted alternating least squares (WALS) method respectively.

B. Estimation of direction independent gains
The least squares problem to find the omnidirectional gains g based on the weighted cost function κ W (θ) can be formulated as where we have introduced R 0 = AΣA H ; in this subsection, R 0 is assumed to be known. This problem can be solved using standard techniques by regarding g and g as independent vector parameters and alternatingly solve for them until convergence 2 . In this approach, the solution for g is Since the result for g is simply the complex conjugate of this relation, it is sufficient to apply Eq. (29) repeatedly until convergence. Although ALS ensures that the value of the cost function decreases in each iteration, it does not guarantee convergence to the global minimum, especially if the initial estimate is poor. The number of iterations required for convergence also depends strongly on the vicinity of the initial estimate to the true value. We are thus interested in obtaining a good initial estimate for g. Fuhrmann [20] has proposed to use a suboptimal closed form solution to initialize the Newton iterations used to solve the ML cost function under the assumption that Σ n is known (or actually no noise is present). A similar problem (with unknown Σ n and R 0 a rank 1 matrix, i.e. a single calibrator source) was also studied by us in [6], where a "column ratio method" (COLR) was proposed. Below, we generalize and improve that technique 3 . All these techniques are suboptimal in the case of multiple calibrator sources, but they can be used to provide the starting point for iterative refinement.
The closed form solution by [20] is derived as follows. If the structure in v = g ⊗ g is neglected and Σ n is known or negligible, we can solve for v in the least squares sense by or equivalently Note that the weights cancel since we solve for one parameter for each entry of R. Subsequently, use the structure of v and assume that the estimate of gg H obtained in the first step has rank 1. An eigenvalue decomposition is used to extract g as the dominant eigenvector from v. It is clear however that the noise on specific elements of v may be increased considerably if some entries of R 0 have a small value. Also, the method is not applicable if Σ n is unknown. In [7] another approach is therefore suggested based on the observation that g i g k = R ik /R 0,ik holds for all off-diagonal elements of R, i.e., for i = k. This implies that This relation is similar to the closure amplitude relation in astronomy [5], [24], which states that in an observation on a single point source or, more generally, in the case of a rank 1 model, the amplitude ratios are related as indicated by the second equality sign in (32).
Since the index k can be chosen freely as long as k = i, j, we can introduce c 1,ij being the column vector containing the values R ik R 0,jk and c 2,ij being the column vector containing the values R 0,ik R jk for all possible values of k = i, j. We can now write (32) in its more general form which has the well-known solution All possible gain ratios can be collected in a matrix M with entries M ij = c † 2,ij c 1,ij . Since the model for M is M ij = g i /g j , the matrix M is expected to be of rank 1, and g can be extracted from this matrix using an eigenvalue decomposition: let u 1 be the eigenvector corresponding to the largest eigenvalue of M, then g = αu 1 , where the scaling α needs to be determined separately since the quotient g i /g j is insensitive to modification by a constant scaling factor applied to all gains. This scaling can be found by minimizing for all off-diagonal elements. By introducing the vectors r 0 = vec − (GR 0 G H ) and r = vec − ( R), where vec − (·) operates like the vec(·) operator but leaves out the elements on the main diagonal of this argument, this cost function can rewritten as which can be solved in the least squares sense by α = r † 0 r. Equation (34) extends the column ratio method proposed in [6] to all elements of the matrix. The column ratio method was introduced to reconstruct the main diagonal of a rank 1 source model, which allows to estimate this rank 1 model using an eigenvalue decomposition without distortion of the results due to unknown sensor noise powers. In our case this allows us to neglect the unknown receiver noise powers. The CRB analysis in [25] shows that it does not matter whether one simultaneously estimates the receiver noise powers and the omnidirectional complex gains exploiting all data, i.e. including the autocorrelations, or ignore the autocorrelations and solve for the direction independent complex gains only.
The Monte Carlo simulations presented in [7] suggest that the method outlined above provides a statistically efficient estimate of g; the pseudo-inverse ensures that the noise on the entries of M does not increase dramatically due to small values in either R or R 0 . It also provides a modification of the initial estimate proposed by Fuhrmann [20] which finds a near optimal solution to the least squares gain estimation problem without further optimization using, e.g., the Newton algorithm. This will generally save computational effort. In the simulations in this paper, we will use Eq. (34) to demonstrate that this method gives statistically efficient results without requiring additional iterations using Eq. (29). Moreover, if three or more iterations are required to ensure convergence using Eq. (29), it also requires less computational effort, as discussed in Sec. IV-G.
Although the pseudo-inverse in Eq. (34) ensures robustness against small entries of either R or R 0 , the fact that the method relies on gain ratios may lead to poor performance if there are small entries in g, e.g. due to failing array elements. This risk can be mitigated by rewriting Eq. (33) to By defining we can enumerate i and j and collect all relations in a single matrix equation This suggests that g can be found by searching for the null space of the p(p − 1) 2 × p matrix C 2 − C 1 . By substituting the details of c 1,ij and c 2,ij in C 2 − C 1 , it is easily seen that g lies indeed in the null space and that this null space is one dimensional. However, this is only true for noise free data. In the practical case of noisy data, there will not be a null space and g lies in the noise subspace. Furthermore, finding g will involve a singular value decomposition on a very large matrix which may be computationally prohibitive.
An alternative approach follows from recognizing that g is obtained from the principal right eigenvector of which is easily found by substitution of the definition of c 1,ij and c 2,ij and use of the Moore-Penrose left inverse. Equation (39) shows that g is the principal right eigenvector of C † 2 C 1 , which has eigenvalue 1 in the noise free case. Simulations with completely randomized R 0 and g suggest that the other eigenvalues are considerably lower than 1. Since it is easily demonstrated that the noise on actual data will lower the main eigenvalue, this is an important result; the contrast between the largest and the second-largest eigenvalue determines the susceptibility of this method to noise on the data.
Note that these methods are still insensitive to a constant scaling factor applied to all gains. This ambiguity is solved by Eq. (36). Also note that all methods presented above enumerate over i, j = k to avoid the diagonal entries of the array covariance matrix which are affected by the system noise. By imposing more stringent restrictions on the choice of i and j, this method can also handle cases in which some elements in the array experience correlator noise, i.e. cases in which Σ n is not diagonal but still has sufficient zero entries to allow the above methods to extract the required information for every combination of g i and g j .

C. Source power estimation
Based on the unweighted cost function κ(θ), σ is found by If weighting is applied, σ follows from Using standard relations for Khatri-Rao and Kronecker products and substituting W c = R − 1 2 we obtain This result confirms the observation by Ottersten, Stoica and Roy [23] that although the derivation involves a square root of the array covariance matrix, the final result only depends on the array covariance matrix and its inverse.

D. Estimating receiver noise powers
If Σ n = diag (σ n ) and no weighting is applied to the cost function, then σ n is found by solving This estimation problem is the same as the sensor noise estimation problem treated in [6], so we just state the result, which is found to be If Σ n = σ 2 n I, a similar derivation gives This result is just the average of the sensor noise estimates obtained when they are estimated individually. If the weighted cost function κ W (θ) is used, we can follow a derivation similar to the one that led us to the estimate for σ in the previous section. We first note that Applying the standard Khatri-Rao and Kronecker product relations and inserting W c = R − 1 2 we get A similar derivation for Σ n = σ 2 n I gives The true value of the covariance matrix R is not known in practical situations. Therefore the measured covariance matrix R is generally taken as estimate of R. It can be shown that this conventional approach leads to a noticeable bias in the estimate of σ n for a finite number of samples N . For simplicity of the argument we will prove this statement for the case where R = σ 2 n I, such that GR 0 G H = 0. Inserting this in Eq. (48) gives Since Since E represents the noise on the data, the expected value of tr (E) and tr E H is zero. This implies that σ 2 n ≈ σ 2 n tr (I) + tr EE H −1 tr (I) This shows that the presence of noise systematically lowers the value of the estimate of σ 2 n . Since R −1 asymptotically in N converges to R −1 , E converges to zero. Therefore this result also shows that the σ 2 n estimate asymptotically converges to the true value if the number of samples N approaches infinity.
However, this bias can be avoided by using the best available knowledge of the estimated parameters, i.e., in (48) use In the first iteration of an alternating least squares algorithm, initial estimates of the parameters are used. These values are replaced by increasingly accurate estimated values in consecutive iterations.

E. DOA estimation
The problem of estimating L is that of estimating the direction of arrival of signals impinging on a sensor array, and has been studied extensively. MUSIC [26] and weighted subspace fitting (WSF) [27], [28] are well-known statistically efficient DOA estimation methods applicable to arbitrary sensor arrays. In either case the eigenvalue decomposition of R is interpreted in terms of a noise subspace and a signal subspace, i.e.
where λ 1 > · · · > λ q > λ q+1 ≥ · · · ≥ λ p . The noise eigenvalues of the whitened true covariance matrix are all equal and the number of sources can be derived from the distribution of the eigenvalues of R; we will assume that it is known. In practical situations in which only an estimate of the true covariance matrix is available, methods like the exponential fitting test [29], [30] or information theoretic criteria [31] must be used to determine the number of signals. It is quite straightforward to adapt the WSF method to our needs. The data model assumed in [28] can be described as where we have used A to avoid confusion with A introduced in this paper. To map our data model on the data model described by Eq. (54) we should whiten our array covariance matrix using and then set A = Σ −1/2 n GA. With these substitutions it is straightforward to implement the procedures described in [28].

F. Alternating Least Squares (ALS) and Weighted Alternating Least Squares (WALS)
The ingredients of the previous four subsections can be combined to formulate a (Weighted) Alternating Least Squares solution to the stated optimization problems. We start by introducing an algorithm handling the first two scenarios identified in section II. DOA estimation is then added in a straightforward way.
To estimate g, σ and σ n we propose the following (W)ALS algorithm. 2) Estimate g [i] by an eigenvalue decomposition of the matrix M with its entries obtained by averaging (32) (resp. (34)) over all k = i, j or by an eigenvalue decomposition of the matrix C † 2 C 1 as given by Eq. (39) using A and σ [i−1] as prior knowledge. Note that neither gain calibration approach does require knowledge of the sensor noise powers σ n and that the latter approach is advisable if some elements may have a very low gain.

3) Estimate σ [i]
n using either (44) 1| < δ or i > i max , stop, otherwise increase i by 1 and continue with step 2. The proposed criterion for convergence is based on a measure of the average relative error in all parameters, and will work even if the parameter values differ by orders of magnitude.
The extension with step 5 is referred to as the Extended ALS (xALS) resp. xWALS algorithm.
An algorithm that alternatingly optimizes for distinct groups of parameters can be proven to converge if the value of the cost function decreases in each iteration. In [27], [28] it is demonstrated that WSF minimizes the least squares cost function w.r.t. the parameterization of A, thus providing a partial solution to the least squares problem considered here. In step 2, g is estimated using a method that only provides a near optimal solution to the least squares cost function. Its solution could, however, not be discerned from the true solution in Monte Carlo simulations, as demonstrated later in this paper. Optionally, step 2 could be augmented with one or two iterations of Eq. (29) to assure minimization of the least squares cost function. Alternatively, one could use the proposed estimate in the first iteration and use Eq. (29) in consecutive iterations in which a proper initial estimate is available from the previous iteration. The other parameters are estimated using well known standard solutions for least squares estimation problems. Therefore the value of the cost function is reduced in each step, thus ensuring convergence. Table I summarizes the numerical complexity of different stages of the ALS and WALS algorithms per iteration expressed in the number of complex multiplications. N iter is the number of iterations required for convergence when using Eq. (29). In this table, the notation o(· · · ) is used to denote all the neglected lower order terms. The complexity of the omnidirectional gain estimates is dominated by the eigenvalue decomposition on M, which takes approximately 4p 3 multiplications assuming use of the divide and conquer ALS g 4 1 2 p 3 + pq 2 + o(· · · ) σ 2p 2 q + 12pq 2 + 9q 3 + o(· · · ) σn 0 σ 2 n 0 total 4 1 2 p 3 + 2p 2 q + 13pq 2 + 9q 3 + o(· · · ) WALS g 5p 3 + pq 2 + o(· · · ) g (Eq. (29)) N iter 2 2 3 p 3 + 4p 2 q + 2q 2 p + o(· · · ) σ 2 3 p 3 + 2p 2 q + pq 2 + 2 3 q 3 + o(· · · ) σn method [32]. This step dominates the overall complexity as well, so it may be worthwhile to estimate the number of iterations of the power method [33] required to obtain sufficient accuracy, especially if p is large.

G. Computational complexity
The pseudo-inverse in Eq. (40) can be implemented by a singular value decomposition which is computationally more efficient than direct computation of the Moore-Penrose inverse [33]. Some intermediate results, such as GAΣA H G H , are required more than once. In our calculation we assume that these terms are computed only once and stored for future use. Under this assumption the number of complex multiplications required to compute the noise power reduces to zero. The number of array elements p will generally be considerably larger than the number of source signals q. Therefore the q 3 -, pq 2 -and p 2 q-terms may be considered negligible compared to the p 3 -term in most cases. These results suggest that the better statistical performance of the WALS algorithm compared to the ALS algorithm as demonstrated by the simulation results presented in the next section comes with only a minor increase in complexity.
Both algorithm can be augmented with source position estimates using WSF. In [28] it is demonstrated that WSF can be implemented with O(pq 2 ) complex operations per Gauss-Newton-type iteration of the proposed modified variable projection algorithm once the result from the eigenvalue decomposition is available. The cost of WSF is therefore dominated by the eigenvalue decomposition on R requiring about 4p 3 complex multiplications. This makes the computational cost of source location estimation comparable to the cost of estimating the direction independent gains.

V. SIMULATION RESULTS
The methods proposed in the previous section were tested using Monte Carlo simulations and compared with the CRB. For these simulations a five-armed array was defined, each arm being an eight-element one wavelength spaced ULA. The first element of each arm formed an equally spaced circular array with half wavelength spacing between the elements. The source model used in the simulations is presented in Table II. This source model was created using random number generators to demonstrate that the proposed approach indeed works for arbitrary source models.  Variance of the parameter vector θ = [γ T , φ 2 , · · · , φ 40 , σ 2 2 , · · · , σ 2 5 , σ 2 n ] T estimates obtained in Monte Carlo simulations for the weighted and unweighted versions of the alternating least squares algorithm and compared these results with the CRB.
In the first simulation the direction independent gains, the source powers and the receiver noise powers were estimated assuming an array of equal elements, i.e. the parameter vector was defined as θ = [γ T , φ 2 , · · · , φ 40 , σ 2 2 , · · · σ 2 5 , σ 2 n ] T corresponding with the first scenario mentioned in section II. Data was generated assuming σ 2 n = 10 and N = 10 5 . This value for σ 2 n implies an instantaneous SNR on the strongest source of only 0.1 per array element, a typical situation for radio astronomers. Figure 2 shows the variance of the estimated parameters based on 1000 runs and compares these values with the CRB. As expected the WALS method appears to be asymptotically statistically efficient while the ALS approach does not attain the CRB and shows clear outliers. Figure 3 shows the variance found on a number of representative parameters as function of the number of samples and compares this with the corresponding CRBs. This plot confirms the conjecture raised in the previous paragraph that the WALS method is asymptotically statistically efficient. Figure 4 shows the bias found in the Monte Carlo simulations and compares this with the statistical error of 1σ for a single realization based on the CRB. This result indicates that both methods are unbiased for all parameters. Figure 5 shows the difference of the parameter value at the end of each iteration and its final value obtained from one run of the Monte Carlo simulation for a representative case, i.e. similar results were found for other parameters and other runs. The standard deviation based on the CRB is plotted as  Bias on the estimated parameters (θ = [γ T , φ 2 , · · · , φ 40 , σ 2 2 , · · · , σ 2 q , σ 2 n ] T ) for the weighted and unweighted version of the alternating least squares algorithm in the Monte Carlo simulations. The bias is compared with the standard deviation on the estimates derived from the CRB. reference to facilitate the interpretation of the vertical scale.
The results indicate that only one or two iterations are needed to reach the CRB in this scenario. With regard to other scenarios, this result suggests exponential convergence, i.e. each iteration adds about one significant digit to the parameter estimate, and it indicates that the WALS method converges more rapidly than the ALS method. In these simulation the stop criterion was In the second simulation the direction independent gains, the apparent source powers and source locations and the receiver noise powers were estimated assuming an array of ideal elements. This corresponds to the third scenario described in section II. The parameter vector to be estimated is thus θ = [γ T , φ 2 , · · · , φ 40 , σ 2 2 , · · · , σ 2 5 , σ 2 n , l 2 , · · · , l 5 , m 2 , · · · , m 5 ] T . Data were generated assuming σ 2 n = 10 and N = 1000. This again implies an instantaneous SNR of only 0.1 per array element but in this case the SNR per receiver is only 3.2 Variance on the parameter vector θ = [γ T , φ 2 , · · · , φ 40 , σ 2 2 , · · · , σ 2 5 , σ 2 n , l 2 , · · · , l 5 , m 2 · · · , m 5 ] T estimates obtained in Monte Carlo simulations for the xWALS method, compared with the CRB.
for the strongest source after integration. Figure 6 shows the variance on the estimated parameters based on 1000 runs and compares these results with the CRB. In these simulations we only used the xWALS algorithm, since the previous simulations demonstrated that the ALS method gives inferior results compared to the WALS approach. These results show that the variance on the parameter estimates are close to the CRB. Figure 7 shows the bias on the estimated parameters found in these Monte Carlo simulations. These results indicate that the estimates are unbiased.
The convergence of a representative parameter is shown in Fig. 8. Again one but preferably two iterations are sufficient to get accurate results. However, in this scenario the stop criterion |θ Bias on the estimated parameters (θ = [γ T , φ 2 , · · · , φ 40 , σ 2 2 , · · · , σ 2 5 , σ 2 n , l 2 , · · · , l 5 , m 2 · · · , m 5 ] T ) for the xWALS algorithm obtained in Monte Carlo simulations. The results are compared with the standard deviation on the estimated derived from the CRB.  Fig. 8. The difference between the value of γ 5 at the end of each iteration and its final value is plotted versus the iteration number for one of the Monte Carlo simulation runs. The behavior shown in this plot is representative for the results obtained for other parameters and from other runs. The standard deviation based on the CRB is also shown to demonstrate that one but preferably two iterations are sufficient to obtain a sufficiently accurate result. stops because the maximum number of iterations. which was set to 15, was reached. This indicates that the algorithm tries to interpret the noise as real signal. In the first iteration first the omnidirectional gains, apparent source powers and noise power are estimated. These values are then used in the WSF to find the source locations. The source locations are then used to update the sky model, leading to an update of γ, σ and σ n and the cycle is complete. In the first simulation, the SNR after integration defined as √ N σ 2 i /σ 2 n was a factor 10 better than in this simulation in which the SNR per array element per source after integration is only 2.2 to 3.2.

VI. SUMMARY AND CONCLUSIONS
In this paper we have developed a weighted alternating least squares (WALS) algorithm to solve for direction independent complex gains, apparent source powers and receiver noise powers simultaneously in a snapshot observation by a sensor array. Our solution for finding the direction independent gains extends and improves the available methods in the literature in the case of multiple calibrator sources and provides robustness to small entries in the array covariance matrix and small gain values (due to e.g. failing elements). Although we have assumed independent sensor noise powers, our gain estimation method is easily applied to cases in which some pairs of elements experience correlated noise. We also found that unbiased estimation of the receiver noise powers requires weighting with the best available array covariance matrix model instead of weighting with the measured array covariance matrix, which is common practice in signal processing.
The statistical performance was compared with an unweighted alternating least squares (ALS) algorithm and the CRB and found to provide a statistically efficient estimate thereby outperforming the ALS method. The computational complexity of ALS and WALS scale with 4 1 2 p 3 and about 7p 3 respectively assuming that the number of calibrator sources is much smaller than the number of array elements, i.e. the WALS method comes with only a minor increase in computational burden. Simulations indicate that only one or two iterations are sufficient to reach the CRB and that every iteration adds about one significant digit, the WALS method converging slightly faster than the ALS method.
We extended these methods with source location estimation using weighted subspace fitting. Simulations indicate that this extension to the WALS algorithm provides a statistically efficient simultaneous estimate of omnidirectional gains, apparent source powers, source locations and receiver noise powers. Again, one to two iterations proved to be sufficient to reach the CRB. However, the convergence of the parameter is blocked at some level lower than the CRB by the noise in the measurement, at which point the algorithm keeps jumping from one local optimum to another.