A New Gain-Phase Error Pre-Calibration Method for Uniform Linear Arrays

In this paper, we consider the gain-phase error calibration problem for uniform linear arrays (ULAs). Based on the adaptive antenna nulling technique, a new gain-phase error pre-calibration method is proposed, requiring only one calibration source with known direction of arrival (DOA). In the proposed method, a ULA with M array elements is divided into M−1 sub-arrays, and the gain-phase error of each sub-array can be uniquely extracted one by one. Furthermore, in order to obtain the accurate gain-phase error in each sub-array, we formulate an errors-in-variables (EIV) model and present a weighted total least-squares (WTLS) algorithm by exploiting the structure of the received data on sub-arrays. In addition, the solution to the proposed WTLS algorithm is exactly analyzed in the statistical sense, and the spatial location of the calibration source is also discussed. Simulation results demonstrate the efficiency and feasibility of our proposed method in both large-scale and small-scale ULAs and the superiority to some state-of-the-art gain-phase error calibration approaches.


Introduction
Array signal processing has been intensively investigated over the last several decades [1,2]. As one of the hotspots in array signal processing, the direction of arrival (DOA) estimation problem of plane waves impinging on an antenna array has drawn widespread attention. Numerous DOA estimation algorithms have been developed to achieve high-resolution performance, e.g., multiple signal classification (MUSIC) [3], estimating signal parameters via rotational invariance technique (ESPRIT) [4], maximum likelihood (ML) approach [5], compressed sensing (CS) method [6] and sparse Bayesian learning (SBL) [7]. These estimators work well under the assumption that prior knowledge of the array manifold is available. However, in practical applications, this assumption is unrealistic since a series of array imperfections including gain-phase errors, mutual coupling and sensor location errors, may affect the ideal array manifold and degrade significantly the DOA estimation accuracy [8][9][10][11][12]. Therefore, it is highly desirable to estimate and calibrate array imperfections before DOA estimation. We focus the paper on estimating the gain-phase errors of uniform linear arrays (ULAs). Typically, the gain-phase errors are caused by the following factors: gradual changes of the behavior of the sensor itself and the internal circuits (due to thermal effects, aging of components, etc.), changes in location of the array elements induced by the vibrating wing of an aircraft or a hydrophone array behind a ship, and changes to the environment around the sensor array [13]. Herein, we consider these three main origins of gain-phase errors and assume that the gain-phase error caused by these changes is in the presence of a constant since all the above-mentioned changes are tiny and slow in a short time.
In recent decades, lots of effort has been made to address the gain-phase errors of antenna arrays [13][14][15][16][17][18][19][20][21][22][23][24][25]. Existing gain-phase error calibration methods can be roughly divided into two categories: self-calibration (without calibration sources) and pre-calibration (with calibration sources). In general, the first type of gain-phase error estimation method, which can jointly estimate the DOA parameters and gain-phase errors, does not require other calibration sources with known directions. Hence, it is more suitable for practice and has evoked much research interest for the past few years [13][14][15][16][17][18]. In [13,14], a self-calibration iterative method based on the eigenstructure is proposed to simultaneously estimate the DOAs of source signals and calibrate gain-phase errors. However, the drawback of this method is that it may suffer from suboptimal convergence. The algorithm proposed in [15] provides a consistent estimate of the gain-phase errors from the ML perspective, but it requires the covariance matrix of the ideal array output to be known. By using the diagonal lines of the covariance matrix, a class of the gain-phase error calibration algorithm is presented for linear equispaced arrays (LEA) in [16], but the perfect covariance matrix is unavailable when the number of snapshots taken to obtain the covariance matrix is small [17]. Z. Liu et al. proposed a novel sparse Bayesian array calibration (SBAC) method for compensating for all the typical array imperfections [18], whereas the online estimation of numerous unknown parameters makes the SBAC more complicated. Many proposals have been conducted on the gain-phase error calibration problem for nonlinear arrays [19][20][21]. In [19], A. Liu et al. presented an eigenstructure DOA algorithm in the presence of gain-phase errors, which is based on the dot product of the array output and its conjugate. Nevertheless, it requires that at least two signals are spatially far from each other. A Hadamard product-based method [20], which performs independently of the phase errors and without the requirement of two signals being spatially far separated from each other, achieves better estimation performance with respect to [19]. However, it suffers from high computational load. In [21], W. Xie et al. extended the amplitude-only measurementsbased technique [19,20] into central symmetric arrays (CSAs) and proposed an algorithm for jointly estimating the DOAs of noncircular sources and gain-phase errors. In addition, for underwater sensing tasks, the random acoustic fluctuations in the medium, such as a highly dynamic ocean, may introduce phase errors in the output of the ULA. To address this issue, a novel L1-norm (absolute-error) maximum projection principal component analysis (PCA) method in [22] is proposed to resist the uncertainty of random acoustic fluctuations. Dubrovinskaya et al. designed acoustic arrays of an arbitrary shape and a practical algorithm to achieve accurate DOA estimates under some array imperfections, such as the spatial ambiguity that must be compensated for [23].
Unlike the self-calibration method, the second type of method, named the pre-calibration method, requires calibration sources with known directions. The eigenstructure-based pre-calibration method with two calibration sources with known directions was developed in [24], which depends on the covariance differencing technique and iterative method. However, since the gain-phase errors are estimated by using the phase of each entry in the steering vector, the possible phase ambiguities need to be checked. Based on the elegant analysis of the Cramer-Rao bounds on calibration and source location accuracies under three different sensor location situations, the authors of [26] indicate that two calibration sources are needed to calibrate the planar array, and the gain-phase estimation errors tend to decrease as the calibrating signal strength or power increases. By utilizing a set of calibration sources in known locations, a ML calibration algorithm is presented in [27] and the array perturbation parameters can be solved by means of the least square method. Nonetheless, the unique determination (identifiability) of the perturbation parameters is not satisfied in this method. Based on the null characteristic of the MUSIC spectrum, the array sensor gain-phase error calibration problem can be formulated as a series of linear equations and solved by a constrained optimization [25]. When the number of the linear equation is larger than that of the gain-phase errors to be estimated, a unique solution can be obtained. However, plenty of calibration sources with known DOAs are needed to solve for the uncertainties of array gain-phase errors. Under the existence of array phase errors, a non-coherent DOA estimation method is proposed in [28,29] based on the modified version of the greedy sparse phase retrieval (GESPR) framework [30]. Nevertheless, multiple calibration sources are required to cope with the phase ambiguity problem, and the global minima cannot be guaranteed [31]. Compared with the self-calibration method, although the pre-calibration method with calibration sources is somewhat limited practically when numerous calibration sources are required, there are two reasons why we take it to be once again in the spotlight. First, the pre-calibration one, in general, is more computationally efficient than the self-calibration since it avoids the alternation between the DOA estimation and array calibration, and the uncertainty calibration procedures are accomplished one at a time. Second, the pre-calibration can provide more satisfactory calibrating accuracies than the self-calibration due to some prior knowledge of the calibration sources [27]. Therefore, to obtain both satisfactory calibrating accuracy and low-cost implementation, we intend to explore an accurate, low-cost, broadly applicable gain-phase error calibrating scheme, which can be available for both small-scale and large-scale ULAs.
The adaptive antenna nulling technique, also known as power-inversion adaptive array [31], was firstly proposed by Appelbaum [32]. Its primitive application was to reduce sidelobe levels in the unknown direction of interferences by weighting the received signal vectors to minimize the interference powers. Alternatively, it is only necessary to have the steering vectors in the directions of interest to be zeros, i.e., null steering. In [33], several constrained null steering algorithms based on adaptive antenna nulling array are introduced, including mainly constrained least-mean-square (CLMS) and QRrecursive least-squares (QR-RLS) algorithms, and some convergence behaviors of the corresponding algorithms are also provided. The infinite impulse response (IIR) adaptive nulling array structure is designed in [34], which consists of a number of shift-invariant subarrays and takes the outputs of previous subarrays as spatial feedback. Moreover, the null steering scheme has been applied in combating spatial acoustic feedback between the hearing aid loudspeaker and microphone(s) [35,36], in which the calculation of null steering beamformer coefficients can be formulated as a least-square (LS) optimization or a min-max optimization so as to be null in the direction of the acoustic feedback. As an extension of [36], the same authors present a soft constraint on designing the null steering beamformer coefficients [37] in order to cope with the incoming acoustic signal preservation and feedback cancellation simultaneously.
In this paper, unlike using nulls of the MUSIC spectrum, we present a new gainphase error pre-calibration method for ULAs by exploiting the adaptive antenna nulling technique [33,34] and null steering algorithm [35][36][37]. Only one known-DOA calibration source is required. We divide a ULA with M array elements into M − 1 sub-arrays, and the gain-phase errors of each sub-array can be uniquely extracted one by one. To be specific, for the mth sub-array, it is able to derive the gain-phase errors as the reciprocal of the cumulative product of a series of complex coefficients, which relate to the sub-array unperturbed null steering vector (SAUNSV). In order to reliably estimate the SAUNSV,w m , we formulate the received data structure of ULAs as an errors-in-variables (EIV) model, based on which a weighted total least-squares (WTLS) algorithm is presented. Furthermore, we show that (1) the optimal location of the unique calibration source is the normal direction of the ULA. (2) When the calibration signal-to-source signal power ratio (CSR) increases, the gain-phase error estimation performance improves. As a consequence, it requires the calibration signal power to be much larger than the source signal. Our main contributions are therefore as follows: (1) By using the adaptive antenna nulling technique, a new gainphase error pre-calibration method with only one calibration source is presented. (2) We formulate an errors-in-variables (EIV) model and propose a WTLS algorithm in order to estimate the SAUNSV, which is a crucial factor for the gain-phase error estimation. (3) Some comparative statistical analyses on the solution to the WTLS problem are derived.
The paper is organized as follows. In Section 2, the signal model for gain-phase error estimation is established and the adaptive antenna nulling technique is reviewed. The proposed method is derived in Section 3. In Section 4, we formulate an EIV model, propose a WTLS algorithm for estimating the SAUNSV and give the statistical analyses on the solution to the WTLS. Section 5 describes the simulation results, and the conclusion is given in Section 6.
Notation: Superscripts (·) H , (·) T and (·) * stand for conjugate transpose, transpose and complex conjugate, respectively. Matrices and vectors are represented by bold uppercase and bold lower-case characters, respectively. The notation diag(a) means forming a diagonal matrix by using the vector a as its main diagonal entries, while diag(Λ) denotes the vector formed by the diagonal of Λ if Λ is a diagonal matrix. The notation E(·) denotes the mathematical expectation, tr(·) stands for the trace operator, || · || and | · | represent the Euclidean norm and absolute value, respectively. ⊗ denotes the Kronecker product. I M denotes a M × M identity matrix, 0 M is a M × 1 column vector with all zero entries and 0 MM is a M × M matrix with all zero entries. vec(·) denotes the vectorization operator, which stacks all columns of a matrix one below the other to form a column vector.

Signal Model
We consider that L narrowband far-field uncorrelated signals from directions θ = [θ 0 , θ 1 , · · · , θ L−1 ] T impinge on a ULA, which consists of M (L < M) isotropic antenna elements. The distance between two adjacent elements is d = λ/2, i.e., a half-wavelength, where λ is the wavelength of the source signal. Without loss of generality, we take the first array element as reference, and hence, the gain-phase errors of the ULA can be modeled as a diagonal matrix, i.e., where Γ 0 = 1 and Γ m = g m e jϕ m , g m and ϕ m are the corresponding gain and phase errors of the (m + 1)th (m = 1, 2, · · · , M − 1) array elements, respectively. In the tth snapshot, the received signal vector is written as where A(θ) = [a(θ 0 ),a(θ 1 ), · · · , a(θ L−1 )] is the array manifold corresponding to the source signals. a(θ l ) = [1,e jπ sin θ l , e j2π sin θ l , · · · , e j(M−1)π sin θ l ] T denotes the steering vector attached to the lth (l = 0, 1, · · · .L − 1) uncorrelated source signal, s l (t). s(t) = [s 0 (t), s 1 (t), · · · , s L−1 (t)] T with equal power σ 2 s which can be defined as E[|s l (t)|] = σ 2 s , and n(t) = [n 0 (t), n 1 (t), · · · , n M−1 (t)] T is the additive noise vector, which is assumed to be spatially and temporally white with zero mean and covariance matrix E[n(t)n(t) H ] = σ 2 n I M , where σ 2 n is the noise variance of each array element. In addition, the noise herein is assumed to be an internal noise model, which means that the noise level in each array element is not affected by the array responses [38,39]. If the external noise model is considered [16], the noise level changes with both the array responses and received signal level.
We herein utilize a unique calibration signal r 0 (t) with known DOA, γ 0 , to estimate and calibrate gain and phase errors. The calibration signal can be turned on and off at our command, and the following assumption associated with the calibration signal holds. Assumption 1. The calibration signal r 0 (t) is uncorrelated with other source signals and measurement noises, γ 0 = θ l .
In this case, the received signal vector of the ULA is expressed as where b(γ 0 ) = [1,e jπ sin γ 0 , e j2π sin γ 0 , · · · , e j(M−1)π sin γ 0 ] T represents the steering vector associated with the calibration source signal. It should be noted that there is a special case when s(t) = 0 L , i.e., no source signals are received, x(t) = Γb(γ 0 )r 0 (t) + n(t), e(t) = n(t), and 0 L is a column vector with L zero entries.

Adaptive Antenna Nulling Technique
As shown in Figure 1, the adaptive antenna nulling array [33] is composed of M array elements, which are uniformly spaced. The received signals in the adaptive nulling array are weighted by a set of complex weights, which aim to form nulls in certain directions, e.g., direction of strong interferences and unwanted signals. Moreover, the complex weights can be adjusted by various adaptive updating rules [40,41]. We consider the unperturbed null steering vector ω = [1, ω 1 , · · · , ω M−1 ] T corresponding to the unique calibration source, and the array pattern is expressed as p(γ) = ω H Γb(γ). Based on the adaptive antenna nulling technique, the unperturbed null steering vector can be applied to be null in the calibration signal direction, γ 0 , i.e., where γ 0 ∈ (−90 • , 90 • ) and b (γ 0 ) = Γb(γ 0 ) is defined as the perturbed steering vector with gain-phase errors. Note that the first entry of the null steering vector, ω, is constrained to 1 since we take the first array element as reference.

Proposed Method
We divide a M-element ULA into M − 1 sub-arrays as shown in Figure 2. Each subarray consists of two adjacent sensor elements, the first one of which can be reckoned as the reference. For the mth (m = 1, 2, · · · , M − 1) sub-array, the unperturbed null steering vector and unperturbed steering vector associated with the calibration signal can be defined as w m = [ w m−1 , w m ] T and b m (γ 0 ) = [ e j(m−1)π sin γ 0 , e jmπ sin γ 0 ] T , respectively. Consequently, the nulling array pattern for the mth sub-array can be expressed as where The unperturbed null steering vector in (5) can be rewritten as w m = w m−1wm = w m−1 [ 1,w m ] T , wherew m = [ 1,w m ] T is called the sub-array unperturbed null steering vector (SAUNSV) andw m = w m w m+1 , we also have w 0 =w 0 = 1 and w 1 =w 1 , which can be defined as the first reference entry of the SAUNSV. Thus, for the first sub-array (m = 1), b 1 (γ 0 ) = [ 1, e jπ sin γ 0 ] T , Γ 0 = 1 and Γ 1 can be derived from (5) by For the mth (m = 1) sub-array, the gain-phase error of the mth array element can be derived from (5) as Combining (6) with (7), it yields Since (8) implies that the gain-phase errors are determined on a series of complex coefficientsw 0 ,w 1, · · · ,w M−1 and the known DOA of the calibration source, γ 0 , it is a vital task to obtain the SAUNSV,w m = [ 1,w m ] T , by using the received signals. We propose a WTLS algorithm to estimate them in Section 4. Then, we outline the proposed gain-phase error estimation scheme in Algorithm 1, in which the DOA estimation problem is also considered.
Algorithm 1 Proposed gain-phase error estimation method 1: Turn on the calibration signal, and generate one calibration signal r 0 (t) with known direction, γ 0 . 2: Collect T (T > N) snapshots of receiving signals on the ULA, x(0), x(1), · · · , x(T − 1). 3: Set m = 1, initialize with Γ 0 = 1 and w 0 = 1. 4: for m = 1 to M − 1 do 5: Estimate the second entry of the SAUNSV, w m , by using the WTLS algorithm presented in Section 4. 6: Calculate the estimate of the gain-phase error for the mth array element,Γ m , by using w 0 , w 1, · · · , w m according to (8). 7: end for 8: Turn off the calibration signal (r 0 (t) = 0), collect other T 0 signal snapshots and calibrate the gain-phase errors depending on the estimating results in 4-7, and carry out the DOA estimation by using the MUSIC or other DOA estimation methods.

Estimation of the SAUNSV
As can be seen in (3), the received calibration source signal vector Γb(γ 0 )r 0 (t) is perturbed by the source signal plus noise vector e(t) = ΓA(θ)s(t) + n(t). Our primary objective in this section is to formulate an EIV model for each sub-array and propose a WTLS algorithm to estimatew m and Γ m .

EIV Model for the SAUNSV Estimation
For the mth (m = 1, 2, · · · , M − 1) sub-array, the received signal sub-vector, x m (t), with size of 2 × 1 is written by where Suppose, without loss of generality, s(t) = 0 L . By collecting N snapshots of x m (t), three corresponding matrices can be formed as In what follows, we omit the snapshot index t for clarity. According to the nulling array pattern (5) and received signal model (9), we obtain where r 0 = r 0 (t), r 0 (t+1), · · · r 0 (t+N − 1) T .
Note that the first entry ofw m is unit, (13) can be rewritten as where −y and x 0 are two N × 1 column vectors which represent the first and second columns of X H 0 , respectively. In the same way, −e y and e x are the first and second columns of E H x , respectively.
It is obvious that (15) is a standard EIV model for the SAUNSV estimation, and stems from a series of studies by G. H. Golub [42] and S. Van Huffel [43]. The calibration source signal vectors x 0 and y are perturbed by the random error vectors −e x and −e y , respectively, which have to be corrected by the corresponding terms in (15). Furthermore, according to the unconditional-model assumption (UMA) that the source signal s(t) is random [44], we give the following two assumptions.

Assumption 2.
The random vector e 0 = [ e T x e T y ] T is assumed to be a stationary one with zero mean, and its covariance matrix is expressed by where Q e is the covariance matrix of e 0 , Q xx and Q yy represent the auto-covariance matrices of e x and e y respectively , and Q xy denotes the cross-covariance matrix of e x and e y .

Assumption 3.
The matrix E x , which denotes the source signal plus noise, is assumed to be uncorrelated along the columns. This means the vectors e m (t+i) and e m (t+j) (i = j; i, j ∈ {0, 1, 2, · · · , N − 1}) for different snapshots are temporally uncorrelated. Only the row-correlation in E x is considered. Thus, the covariance matrix Q e can be approximated by where σ 2 xx denotes the variance of the source signal plus noise in the mth array element, i.e., e x , σ 2 yy is the variance related to the vector e y , and σ 2 xy is the cross covariance between e x and e y .

Proposed WTLS Algorithm
As mentioned in the previous subsection, the random vectors e x and e y are unknown in the EIV model (15), so the problem statement is, given the observation vectors x 0 and y, simultaneously estimate the SAUNSV,w m , vectors e x and e y . Then, we formulate it as the following constrained WTLS optimization problem: where the weighted matrix Q −1 e is the inverse matrix of the covariance matrix Q e . Aiming to solve the problem (18), it is necessary to recast (18) as described in the following lemma.
To attain the stationary point of cost function J m , we take its partial derivative with respect tow * m and set it to 0, i.e., The first partial derivative on the right hand of (20) can be expressed by Due to the fact that ∂z ∂z * = 0 for a complex variable z, the third partial derivative on the right hand of (20) can be given by Based on the derivative property of a complex-valued inverse matrix, i.e., [45], where D is a complex matrix, we have Then, Note that and Thus, (24) can be rewritten as where the 2N × N matrix Q T is the partitioned matrix constructed from the first column to Nth column of Q e .
Inserting (27) into (23), we obtain Then, by combining (21), (22) and (28), (20) can be rewritten by where the N × 1 vector, c m , is defined as Actually, by considering the relation between c m and e 0 in (67), we should note that c m = e x and e y = (Q where the 2N × N matrix Q (2) e is the partitioned matrix constructed from the (N + 1)th column to (2N)th column of Q e . Therefore, based on (29),w m can be estimated iteratively by where w m is the estimate of the second entry of the SAUNSV in the mth sub-array. As can be seen in (32), the SAUNSV estimation requires the covariance matrix Q e to be known a priori, which is unrealistic. However, we can use the sample covariance matrix of the random vector e 0 instead of it. Hence, at the first stage of the proposed method, we must estimate the covariance matrix Q e when the calibration source signals are of no use, i.e., r 0 (t) = 0 in this stage. Specifically, for the mth (m = 1, 2, · · · , M − 1) sub-array, we need to collect K 0 N (K 0 > 1) snapshots of receiving signals on the mth sub-array without the calibration source signals, e m (1), e m (2), · · · , e m (K 0 N), and form K 0 random error vectors e 0 (0), e 0 (1), · · · , e 0 (K 0 − 1). Then, the covariance matrix Q e can be estimated bŷ As a result, the steps of the proposed WTLS are tabulated as Algorithm 2.

5:
Form the matrix B 0 (t) = [ w m (t) T ⊗ I N −I N ].

7:
Estimate w m (t + 1) using (32). 8: The WTLS problem solver is derived from (19), which is only the necessary condition on optimality in WTLS problem. As a consequence, our proposed algorithm may only converge to a local minimum point. In order to attain the strict minimum point, we need to consider the sufficient condition on optimality and the Hessian matrix associated with the cost function J m , which is given by As is known, the sufficient condition of the desired minimum point is that the Hessian matrix H is positive definite. In practice, whether the Hessian matrix is positive definite is unknown since the TLS optimization is actually a non-convex problem [46][47][48][49]. However, fortunately, the LS solution can be commonly taken as the initial value of the iterative algorithm, and this is able to enhance the convergence efficiency. Furthermore, the numerical simulation results in this paper demonstrate empirically that the gain-phase error and DOA estimators work well, even when these local minima are introduced into the proposed gain-phase error estimation scheme.

Solution to the WTLS Problem
Let us consider the special case that s(t) = 0 L , and the WTLS problem is degenerated into the conventional total least square (TLS) problem, i.e., the unconstrained objective function for the mth (m = 1, 2, · · · , M − 1) sub-array in (19) can be expressed as where Q e = σ 2 n I 2N .
On the other hand, it is well known that the solution on the minimization of (35), i.e., the estimated SAUNSV, w m = [ 1, w m ] T , consistently converges to the exact solution provided that lim N→∞ X 0 X H 0 N exists and is positive definite [40]. Accordingly, when N → ∞, and it can be shown that for the estimate of (32), w m → w m .
In the case that there exist source signals, i.e., s(t) = 0 L , the weighting covariance matrix Q e is not diagonal. We rewrite the objective function J m in (19) as where B 0 B H 0 = (1 + |w m | 2 )I N , and the matrix P is defined as Substituting (17) and PP H can be expressed by where the weighting factor κ = (1+|w m | 2 ) [σ 2 xx |w m | 2 +σ 2 yy −2Re(w m σ 2 xy )] . When N → ∞, in the same way as (35), the solution on the minimization of (37), w m , converges to a specific solutionw m asymptotically, wherew m = α 1 (2) α 1 (1), α 1 (1) and α 1 (2) denote the first and second entries of the eigenvector α 1 = [α 1 (1) , α 1 (2) ] T with respect to the smallest eigenvalue of the weighted sample covariance matrix, X 0 PP H X H 0 N, which can be approximated as Recalling Assumption 1, we can obtain where x r = Γ m A m (θ)s(t). From (39)-(41), we can observe that (1) the weighting matrix PP H reduces to a diagonal matrix with the weighting factor κ, which enables the data from different columns of the matrix E H x [ 1,w m ] T to be equally sized. When e x and e y are equally sized and uncorrelated, i.e., σ 2 xx = σ 2 yy and σ 2 xy = 0, one obtains κ = σ −2 xx (For a special case that s(t) = 0 L and E H x only depends on noises, κ = σ −2 xx = σ −2 n ); (2) when N → ∞, the estimated SAUNSV obtained by the WTLS problem (37) hinges on the eigen-decomposition In the following, we shall conduct discussion on the specific solution,w m , in two cases of uncorrelated and correlated source signals, respectively.
(1) Uncorrelated source signal case Let two Hermitian matrices be The eigen-decomposition of Y leads to the eigenvector ff 1 corresponding to the exact SAUNSV in (36), and Z can be expressed as where It is immediately obvious that the matrix Y is perturbed by σ 2 s Z 0 , inevitably causing perturbations on α 1 . Based on the matrix perturbation theory (see [50], pp. 69-70), the perturbed eigenvector α 1 corresponding to α 1 can be given by using a convergent power series where the expression in bracket is a convergent power series provided that σ 2 s is sufficiently small, and µ k (k = 1, 2, · · · , ∞) are the coefficients of the perturbation expansions of α 1 along α 0 . The first-and second-order perturbation expansion coefficients are given as (see Appendix B) and Subsequently, by keeping the first-and second-order terms of the source signal power σ 2 s , we obtain the perturbed term in (43) as and the specific solution to the WTLS problem (37),w m , can be written bȳ where CSR = σ 2 c σ 2 s is called the calibration signal-to-source signal ratio, and ∆α 1 (1) and ∆α 1 (2) denote the first and second entries of ∆α 1 , respectively.

(2) Correlated source signal case
In this case, Z can be partitioned as where R s = σ 2 s I L + F, σ 2 s I L and F contain the diagonal and off-diagonal entries of matrix R s , respectively. MatrixZ 0 can be written as where In the same way as the uncorrelated source case, we can derive the first-order perturbation expansion coefficient as From (50), we can note that the first-order perturbation expansion coefficient, compared to the uncorrelated source case in (44), is determined by both source signal autocorrelation and cross-correlation, which correspond to the first term and second term on the right hand of (50), respectively. Moreover, for the second-order perturbation expansion coefficient, we can find similar results by following (45). As a result, both the specific solution,w m , and the estimated gain-phase error in (8) are perturbed by the cross-correlation between the source signals.

Spatial Location of the Calibration Source
Since the SAUNSV derived from (37) is perturbed from its exact value, this, of course, causes a bias on the null position of the null-spectrum in (5). Let γ 0 = γ 0 + ∆γ 0 be the practical null position when usingw m , and ∆γ 0 is the null position bias. Following (5), we have By expanding b 1 (γ 0 + ∆γ 0 ) in the first-order Taylor series around the ideal null position γ 0 , we have Substituting (53) and ∆w m into (51) yields It is of interest to note from (54) that in the case that the calibration source is located in the normal direction of the ULA (γ 0 = 0 • ), the minimum null position bias, |∆γ 0 | = |∆w m πw m | is obtained.

Simulation Results and Discussion
In this section, several simulations are conducted to verify the validity and effectiveness of the proposed gain-phase error estimation method. We first consider the gain-phase error estimation performance of the proposed method for both small-scale and large-scale ULAs in Section 5.1. Then, in Section 5.2, the proposed method is compared with the eigenstructure method [13], diagonal line method [17], and MUSIC nulling method [25]. At last, in Section 5.3, the DOA estimation performance results are provided, which depend on the aforementioned gain-phase error estimation methods. All the results are obtained by averaging over 500 Monte Carlo trials. Two far-field uncorrelated source signals, with equal power σ 2 s , impinge on the ULA from directions 5 • and 30 • . The power of the unique calibration signal, σ 2 c is set to be 1, which is located at γ 0 = 0 • . The signal-to-noise ratio (SNR) is defined as SNR(dB) = 10 log 10 σ 2 s σ 2 n . In addition, since one calibration signal coexists with the source signals in the gain-phase error estimation stage, we define the calibration signal-to-source signal ratio (CSR) as CSR(dB) = 10 log 10 σ 2 c σ 2 s . The random gain-phase errors, which are proved to be feasible in [51], can be assumed to be Gaussian distribution, i.e., g m ∼ N(0, σ 2 g ), where σ 2 g is the variance of gain errors. In addition, for the phase error, we also have γ m ∼ N(0, σ 2 γ ), where σ 2 γ is the variance of phase errors. In simulations, σ 2 g and σ 2 γ are set to be 0.1 and 36, respectively. The deterministic gain-phase error is also considered in Sections 5.2 and 5.3 for a comparison with the cited algorithms. In order to evaluate the estimation precision of the proposed method, the rootmean-square error (RMSE) of the gain-phase error estimation is used, which is defined as where K denotes the number of trials andΓ (k) is the estimated gain-phase error in the kth trial. Similarly, the RMSE for testing the DOA estimation performance is defined as shown below: whereθ l,k is the estimated DOA of the lth source signal in the kth trial.

Gain-Phase Error Estimation Performance
In this subsection, we study the gain-phase error estimation performance results of the proposed method in terms of CSRs, snapshots, number of array elements, spatial location of the calibration signal, and correlation between the source signals. The simulated RMSEs are calculated depending on the proposed method, while the theoretical RMSE values are derived by using (47). Unless otherwise stated, the total number of snapshot T = 80 and N = 40. Accordingly, we can see that T/N = 2 iterations in WTLS are carried out. Figure 3 shows the gain-phase error estimation RMSE for the proposed method versus CSRs. The numbers of array elements are set to M = 6, 8, 12, 16 for Figure 3a-d, respectively. Similar cases for large-scale ULAs, i.e., M = 64, 128, 256, 512 are illustrated in Figure 4. We can observe that the RMSEs decrease as CSRs increase, and more satisfactory gain-phase error estimation performance results can be achieved when the power of the calibration source is far larger than that of the signal source. This is due to statistical perturbation of the SAUNSVs according to (46), which depends on CSRs. Furthermore, it is evident that the simulation results closely agree with theoretical predictions when the CSR is larger than 20 dB.   In Figure 5, the RMSE of the gain-phase estimation using the proposed method is plotted versus numbers of snapshots T for different array elements M = 1, 050, 100. Two iterations are carried out for the proposed method with N = T/2. It can clearly be seen that the RMSE decreases as T increases and increases by increasing M, which coincides with the explicit expression of the gain-phase error estimate (8). The more estimated coefficients w m involve, the larger the estimated error becomes. Consequently, the gain-phase error estimation accuracy and RMSE degrade with the increase in the number of array elements. The gain-phase error estimation RMSE for different spatial locations of the calibration source, γ 0 , is shown in Figure 6a. We observe that the RMSE performance of the proposed method degrades with the increase in γ 0 , which attributes to the null position bias of the practical null spectrum in (51). Furthermore, the corresponding null position bias, ∆γ 0 , is plotted in Figure 6b, from which the simulation results are nearly close in agreement with the theoretical values calculated by (54). The minimum null position bias is achieved in the case of γ 0 = 0 • .  Figure 7 shows the gain-phase error estimation performance for different correlation factors of two source signals. In this case, the source covariance matrix can be expressed as where ξ is the correlation factor. We can note that the RMSE slightly increases by increasing the correlation factor; however, this adverse impact can be reduced by raising the CSR.

Comparison to Other Calibration Methods
In this subsection, both deterministic and random gain-phase errors are considered. The deterministic gain-phase error is modeled as with M = 5, and the random gain-phase error is also subject to the Gaussian distribution with zero mean and variances given in Section 5.1. For a fair comparision, we also apply only one calibration signal located at 0 • for the MUSIC nulling method [25]. Figure 8a,b compare the gain-phase error estimation performance results of the corresponding methods for different SNRs in 5-element ULA with deterministic gain-phase errors and 16-element ULA with random gain-phase errors, respectively. The eigenstructure method [13] and diagonal line method [17] belong to the method of the self-calibration type, whereas the MUSIC nulling method [25] is a pre-calibration method like our proposed one. Generally, the pre-calibration methods are superior to the self-calibration type because of the application of calibration sources with known directions. The eigenstructure method [13] fails to work no matter how high the SNR is due to its suboptimal solutions and lack of unique properties. For the diagonal line method [17], an insufficient number of snapshots limits the estimation of the covariance matrix, resulting in a weak gain-phase error estimation performance. The proposed method behaves better than other competitive methods. The RMSE performance comparison on the corresponding methods with a different number of array element is shown in Figure 9. It is observed that the RMSEs of the MUSIC nulling and proposed methods increase as the number of array element increases because more parameters are needed to be estimated for larger arrays which causes larger estimating errors. Figure 10 presents the performance comparison for different numbers of snapshot T in 5-element ULA with deterministic gain-phase errors. For the eigenstructure method [13] and proposed method, two iterations are carried out with N = T/2. We can see that the proposed one performs better than others.

DOA Estimation Performance
In this subsection, we compare the DOA estimation performance results of the corresponding methods. The CRB on DOA estimation is also given, which is calculated by Equation (4.6) in ref. [52]. All the simulation scenarios are identical to Section 5.2. For estimating DOA of source signals, another T 0 = 40 signal snapshots are needed to be collected. For effective evaluation and comparison of the proposed method, the DOA estimation performance of the MUSIC algorithm in the absence of gain-phase errors is taken as the benchmark. When applying the eigenstructure method [13], the DOA can be jointly estimated with the gain-phase errors, while the MUSIC estimator is used to obtain the DOA results after the gain-phase errors are calibrated by the diagonal line method [17], MUSIC nulling method [25] and proposed method. Once the gain-phase error estimateΓ is available, the MUSIC spatial spectrum can be obtained by whereÛ e is the gain-phase error-free noise subspace with dimension M × (M − L). It can be obtained after the gain-phase error calibration and extracted by using the eigenvalue decomposition (EVD) of the signal covariance matrix estimate, which is calculated bŷ Finally, we can find the DOA estimateθ by using the peak searching of f (θ). In order to reduce the grid mismatch, we herein make use of a recursive grid refinement searching scheme, which contains three steps: (1) Create the maximum and minimum angular grids, whose grid step sizes are θ s,max and θ s,min , respectively, and obtain the rough DOA estimate by using θ s,max . (2) Decrease the grid step size, θ s,max = θ s,max I and I > 1, and obtain a refined DOA estimate by using the updated θ s,max at a local range around the rough DOA estimate in the last step. (3) Return to step (2) until θ s,max ≤ θ s,min . In our simulations, θ s,max = 1 • , θ s,min = 0.01 • and I = 10. Figure 11a,b show the DOA estimation performance results of the proposed method and MUSIC algorithm versus CSR in 5-element ULA and 16-element ULA, respectively. It can be observed that the proposed method can achieve satisfactory DOA performance, like MUSIC, without gain-phase errors at high CSR.  Figures 12 and 13 show the DOA estimation performance comparison among the corresponding methods in 5-element ULA and 16-element ULA, respectively. It is found that our proposed method offers similar DOA estimation performance results to MUSIC in the absence of gain-phase errors, owing to the fact that the gain-phase errors are correctly estimated and calibrated. Figure 14 provides the DOA estimation performance comparison for different numbers of the array element; we also can see the superiority of the proposed one. Figure 15 offers the DOA estimation performance comparison of different numbers of snapshots in the 5-element array. It can be noted that the RMSE decreases when increasing the number of snapshot T.

Discussion
The simulation results presented in Sections 5.2 and 5.3 show that our method attains more satisfactory gain-phase error and DOA estimation performance with respect to the eigenstructure method [13], diagonal line method [17] and MUSIC nulling method [25]. The simulation was carried out in a MATLAB 2020a environment using an Intel 3.1-GHz processor with 8 GB of RAM and under a Windows 10 operating system. All the simulation results were obtained by averaging over 500 independent Monte Carlo trials. The eigenstructure method requires a two-step iteration, the first one involving the gain-phase error estimation and the second, calculating the DOA estimates. As such, a convergence to a global minimum and a unique solution may not been guaranteed, which typically leads to the performance degradation of the parameter estimation. The diagonal line method makes use of the different diagonal lines of the data covariance matrix to estimate the gain-phase errors. Nevertheless, it is known that the diagonal elements of the covariance matrix are mostly contaminated by the environment noises. As shown in Figure 8, in contrast to the Diagonal line method, the proposed method is not sensitive to the environment noises. It obtains a satisfactory estimation performance of the gain-phase error, even when the SNR is 0 dB, due to one high-power calibration source. The MUSIC nulling method requires multiple calibration sources and needs to solve a set of linear equations. Due to the adjustable CSR and the optimal spatial location of the calibration source, our proposed method can offer relatively better estimates than the MUSIC nulling method, which does not provide results about the adjustable CSR and optimal spatial location. However, the proposed method requires a high CSR (larger than 20 dB), and the estimation of covariance matrix Q e in (33) leads to the increase in the number of signal snapshots as compared with these three competing methods.
Our future work is to validate the effectiveness and efficiency of the proposed method by using NI software-defined radio radio testbed. The testbed will work with real Wi-Fi data under IEEE 802.11az for indoor localization and positioning. In this prototype, the gain and phase errors are mainly caused by imperfections in antennas, feedlines, and RF chains. Thus, integrating the proposed calibration procedure in the testbed is helpful in the application for the DOA estimation problem.

Conclusions
Based on the adaptive antenna nulling technique, a new gain-phase error pre-calibration method is proposed, requiring only one calibration signal with known direction. We divide a ULA with M array elements into M − 1 sub-arrays, and the explicit expression of gainphase errors can be derived exactly by applying the SAUNSVs of each sub-array. In order to estimate the SAUNSVs, we develop a WTLS algorithm, and the solution to the WTLS problem in the statistical sense is also studied. These analyses suggest that (1) the optimal location of a unique calibration source is the normal direction of the ULA, and (2) the proposed one requires the calibration signal power to be much larger than the source signal. Simulation results demonstrate the good agreement between the theoretical analyses and experimental predictions, the efficiencies and superiority of the proposed method in terms of the estimation precision of gain-phase errors and DOAs.
It can be noted that (A10) is an equivalent representation of the Lagrange function (A2) with consideration of the equality constraint (A1).