Calibration of Aperture Arrays in Time Domain Using the Simultaneous Perturbation Algorithm

Online calibration is desired in antenna arrays of ultrawide bandwidth. This study proposes a time-domain calibration method based on the simultaneous perturbation (SP) algorithm. Two objective functions were established: power of the received signal at array output or combination of power and correlation coefficient between the signal at array output and a target signal. For both criteria, the convergence settings require only two measurements at each iteration. One advantage of the method is that the entire signal operation for calibration is performed in the time domain. This is achieved by resolving the effects of distortion on time delay of each channel, which accounts for both amplitude and phase distortions at different frequencies. Therefore, the proposed method significantly increased the calibration efficiency for ultrawideband antenna arrays. Since time delay coefficients for calibration associated with array elements were determined independently due to characteristics of the SP, the estimation accuracy of the method is tangential to the number of elements in the array and is mainly dependent on the convergence conditions. This gives the method an additional distinct advantage for calibrating large-scale antenna arrays with ultrawide bandwidth. An estimation accuracy of 99% on time delay adjustments has been achieved and demonstrated.

, [4]. However, calibration of large arrays becomes increasingly more challenging due to an increasing number of parameters to be determined, particularly for very large-scale aperture arrays such as the Square Kilometre Array. This becomes even more challenging at higher operational frequencies for wider bandwidths, due to time delay drift at every array element becoming more sensitive to changes in the environment, and active radio frequency (RF) components. As a result, the number of measurements, in particular at on-site antenna assembly, on the objective function needed for calibration also rises significantly [5], [6], [7], [8].
Stochastic approximation (SA) has been considered to be an effective method for solving problems consisting of a large number of unknown parameters such as the problem faced in the calibration of large-scale antenna arrays [4], [9]. The "simultaneous perturbation" (SP) concept was later introduced to reduce the number of measurements required on the loss function or gradient per iteration for stochastic search algorithms [10], [11]. The SP is also shown to accelerate the convergence of the SA using the adaptive SP (ASP) [12]. It leads to a more efficient adaptive algorithm than traditional finite difference (FD) methods and has a potential on a wide range of practical implementations. However, its actual implementation has not been sufficiently successful because it is extremely costly to estimate the gradients as a large number of loss function measurements involved in the case of gradient-free approach. In the gradient-based scenario, the calculation of gradient usually requires full knowledge of the relationship between the parameters being optimized and the loss function. However, in the optimization process for some applications, gradient is not available or is difficult to compute, such as for array calibration. In [13] and [14], a new simultaneous perturbation SA (SPSA) consensus algorithm for distributed tracking under unknown-but-bounded disturbances is proposed, and it is suitable for distributed problems to estimate time-varying parameters and compensate them.
Determining complex signals (consisting of amplitude and phase) from each array element and compensating the elementto-element variations due to imperfect conditions is the key objective of the array calibration process [15]. Effective phased array calibration methods are based on toggling element phases such that in situ element fields can be measured. However, since these methods must measure not only amplitudes but also phases of the array responses, they are difficult to apply due to practical difficulties, such as large-scale antenna arrays with increased bandwidth where obtaining accurate  phase measurements in real time is still extremely challenging. A novel amplitude-only measurement method was developed in [16]. However, it requires a large number of measurements to extract the calibration coefficients.
Conventionally, calibrating antenna arrays aims to achieve the equalization of amplitude and phase for all antenna channels of the array in the frequency domain [17]. In the calibration process, the characteristics of the array are represented by a manifold vector, which consists of amplitudes as well as phases that occur due to the distance from the source to each element in the array. Once the calibration has been completed, the system produces a synthesis of the coherent signals at the array output. The simplest aperture array calibration technique strives to calculate the optimal phase shifts in all channels by maximizing the received power at the array output [18]. For example, the phase adaptive method was utilized to tune the phase shifters, but this does not guarantee an optimal solution in antenna array calibration problems, even with highresolution phase-shifting circuits [19].
More robust calibration algorithms using phase measurements reduce the implementation scope to laboratory conditions. The rotating element electric-field vector (REV) was developed in 1982, and since then, a few modifications have been made to enhance its capability and it requires measurement at the element level [20]. The most recent calibration algorithms are the ones that do not require measurements at the element level. Hence, no or minimum disruption is caused to the array operation [18]. In these methods, a finite number of power measurements at array output are carried out, and the measured data are then used for estimating unknown phase shifts in order to compensate them for each channel. However, these algorithms are currently limited to narrowband systems as they are essentially implemented in the frequency domain.
There are many techniques developed, mostly in frequency domain, for array calibration. They were tentatively classified into six categories and summarized in Table I. As indicated in Table I, most of the calibration methods operate in the frequency domain. It can be seen in Table I that the number  of measurements required for the calibration is significantly high for some techniques, i.e., a few times more than the total number of elements in the array.
The SPSA has been proven to be an effective stochastic optimization method [10]. The biggest advantage of this method is that it is very effective when the loss function measurements include noise. In addition, it does not need gradient measurements on the loss functions. This method was initially proposed to solve multivariate optimization problems [11]. The method has attracted considerable attention due to its easy implementation. The underlying gradient approximation requires only two measurements of the loss function regardless of the dimension of the optimization problem. Random direction is another gradient search procedure [21], in which the perturbation variables are considered to be uniformly distributed over the surface of a unit sphere. However, it is computationally expensive to obtain the perturbation of random variables. Therefore, in this study, the SPSA approach was adopted as it can handle the optimization of a large number of variables without demanding prior information of gradients, as in the case of array calibration in time domain.
A novel calibration method in time domain is presented in this article. This algorithm does not require to access channel signals or powering down channels for time delay estimation for each array element (signals from array elements are processed through channels in the receiver). The gradients in each iteration were approximated by a highly efficient and easily implemented "SP," and only two measurements are required at each iteration. The total number of measurements is significantly reduced, which is required to resolve the time delay offset for each element of the array in the compensation process. Using the same analogy, in addition to the power measurement approach, the correlation coefficient between the signal based on time delay perturbed element channels and the desired signal is used as a factor for objective function construction. Both forms of objective function in implementing the algorithm for array calibration are investigated in this study.
This is the first attempt to calibrate aperture arrays in the time domain by exploiting the SPSA method, to the best of the authors' knowledge. The proposed method not only avoids seeking calibration coefficients for a large number of frequency points required in wideband array calibrations but also opens a new door for efficient calibration of arrays consisting of a large number of elements, such as the Square Kilometre Array. By adopting the SP technique, estimation errors on time delay distortions of element channels were controlled by a unified convergence criterion and are independent of each other. Hence, the total number of array elements has a negligible effect on accuracy. It is the number of iterations needed to reach convergence that rises with an increasing number of elements in the array and not the estimation error as will be outlined in Section V.
This article is organized as follows. Section II gives the definitions of symbols used in this study. Section III describes the SA theories and convergence conditions. Section IV explains the proposed antenna array calibration method in the time domain. Section V makes a numerical analysis on the algorithm. Section VI introduces the experiments and gives the results. Discussions on the algorithm are made in Section VII. Finally, Section VIII presents the conclusion.

II. SYMBOLS
This article uses T and its variants as vectors representing time delays of element channels in arrays, and a variable with a hat means that the value is obtained based on approximation, such asĝ k that is a vector consisting of approximated gradients at the kth step. Furthermore, E represents an expectation operator. Rationalized m.k.s units are used, and the main symbols used in this study are defined as follows.
1) s(t) = is the short pulse signal in the time domain (volts). 2) σ 2 p de f = is the variance of the Gaussian pulse signal.
3) w de f = is the constant related to the time width of Gaussian pulse signal used for calibration. 4) P t = is the average power of the transmitted pulse signal (watts). 5) P r = is the average power of the received signal from the array without calibration (watts). 6) s m (t) = is the signal received from the mth element (volts). 7) A m = is the amplitude gain for the mth element. 8) ξ m = is the spatial time delay of the signal received by the mth element due to relative spacing, determined by scan angle and element spacing. 9) T spatial de f = is the vector of spatial time delay related to element spacing and scan angle, which is known to the system. 10) γ m = is the intrinsic time delay generated in random for the mth element to reflect imperfect receiving condition, and it is unknown to the algorithm. 11) T in de f = is the vector of the intrinsic time delays introduced in array element channels, representing impact by environment and imperfect devices, and the values in this vector are the time delays to be calibrated out by the algorithm. 12) T de f = is the vector of time delays for all elements, including spatial time delays and time delays caused by the imperfect receiving condition. 13) τ m = is the approximated time delay of the mth element representing the misalignment of timing between elements, the time delay approximated by the algorithm for calibration. 14)T * de f = is the vector of optimal calibration coefficients of time delay for all element channels approximated by the algorithm. 15) T * = is the scalar loss function determined by the power received at the array output alone or the product of power and the correlation coefficient between the signal from the array output and the ideal reference signal for calibration. 26) a k = is the gain coefficient for the kth iteration. 27) c k = is the perturbation step coefficient for the kth iteration. 28) ki de f = is the ith element of vector k with the total length of M, and k is the current iteration number. 29) k de f = is the perturbation vector with the total length of M, and each element of the vector is generated using a Bernoulli distribution with a probability of 0.5 for ±1.
= is the approximated gradients vector at the kth iteration of the algorithm. 31) r (x, y) = is the correlation coefficient between vectors x and y.

III. SA THEORIES
To identify the appropriate SA algorithms for array calibrations, convergence condition is an important factor to consider. In this section, we give a brief review of SA theories and their corresponding convergence conditions.
The basic SA algorithm is essentially a stochastic difference equation with a small step size, and the prime concern is its qualitative behavior after many iterations, such as convergence and rate of convergence. Proofs of convergence and the derivation of the rate of convergence have been given for a variety of equations and stochastic processes over a period of several decades. In the early development of stochastic search and optimization-SA, the Robbins-Monro algorithm (also known as RM) for root finding is the basic approach. Recursive procedure for finding the root of a real-valued function g(θ ) was developed. Suppose that the values of g(θ) were not known, but "noise corrupted" observations can be taken with respect to the values of θ. By proposing where a k is an approximate gain sequence satisfying and Y n is the noisy estimation of the value g(θ). The decreasing step sizes a k would provide an implicit averaging of the observations and are essential to ensure the convergence for this algorithm. For another type of basic problems, to seek the value of the variable at the minimum of a smooth functionsometimes-the form of the function is unknown such as for an array calibration problem, and instead, "noise corrupted" observations or measurements can be obtained. Then, the Kiefer-Wolfowitz procedure (also known as KW algorithm) can be used, in which gradients for directing the recursive computation are estimated via FDs using the noisy measurements, and the step sizes are also small. Assume that the value of parameter θ is to be found for the function E F(θ, χ) = f (θ ), where f (·) is continuously differentiable and χ is a random vector; the functions F(·) and f (·) are not known and loosely speaking, and f (θ) is an "estimator" of F(θ, χ). If c k → 0, FD can be expressed as [43] and the iteration can be established as where Y k can be estimated by where −β k is the bias in the FD estimation f θ (θ k ). Accordingly, the iteration relation can be updated as where ψ k is the FD between the estimated f (θ) and the actual F(θ, χ) when θ changes from θ + c k to θ − c k . To ensure the convergence, in addition to the condition for the bias, β k → 0, and the noise term ψ k /(2c k ) is also required to go toward zero (through "average locally"). The KW algorithm, when the noisy estimates of the derivatives are available, is generally based on these derivatives to direct the recursions toward a converged result. For the same minimization problem, when the stochastic gradient is not available and only measurements of loss function are available, an FD approximation to the gradient is proposed. The FD approximation relies on a small change on the variable θ, and the value of loss function is then measured. However, the FD-based SA (FDSA) algorithm can be very costly when the dimension p of the variables is high since we have to measure at least once for each variable to be optimized. With two-sided perturbations, 2 p measurements are required in each recursion. The recursive procedure can be expressed in the general SA form [44] θ k+1 =θ k − a kĝk θ k whereθ k andĝ k are two vectors consisting of the variables for optimization and the corresponding gradient approximations, respectively. The FDSA gradient approximation can be expressed aŝ where ξ i represents a vector with 1 in the ith place the 0's for the rest and c k > 0 defines the steps for changes. The pair (a k , c k ) is the gain sequence for the FDSA algorithm. Compared to the convergence theory for the root-finding RM algorithm, the conditions for convergence of the FDSA algorithm (evolved from the KW algorithm) have an extra gain sequence c k , arising from the bias inĝ k (θ) as an estimator for g k (θ). The conditions for the formal convergence [almost sure (a.s.)] of the FDSA algorithm are given as follows [44]: Clearly, a large number of measurements (2 p times the total number of iterations) are required when p is large. Hence, the SPSA algorithm came into attention, where the required measurements can be significantly reduced by approximating the gradients for all variables simultaneously at each iteration. The gradient approximation for SPSA is generated bŷ In contrast to the FDSA algorithm, with SP, two loss measurements are needed in each approximation cycle instead of 2 p, regardless of the dimension p. The convergence conditions of SPSA algorithm are closely related to that of the FDSA algorithm, in addition to a more stringent requirement on smoothness of the loss function, i.e., three-times differential, two more conditions are added. B.1 Measurement Noise, Relationship Between Measurement Noise and k : This condition together with smoothness of the loss function guarantees that the gradient estimateĝ k (θ) is an unbiased estimate of g k (θ) in the order of O(c 2 k ). B.2 Statistical Properties of the Perturbations: The bounded inverse moments' condition for ki is an important part of SPSA, and the symmetric Bernoulli distribution ±1 satisfies the inverse moments, which ensures thatĝ k (θ) is a nearly unbiased estimate of g k (θ). As mentioned in [44], proofs of the convergence theorems for the FDSA and SPSA have been provided by Fabian (1971), Spall (1992), Dippon and Renz (1997), and so on. The convergence conditions mentioned above illustrate an abstract idea. For solving real problems such as array calibration in this study, we gave specific loss function definitions and checked the convergence conditions for arrays with various numbers of elements. The strategy of selecting optimal parameter coefficients in the algorithm for effective convergences has been studied considering implementations in practice.

IV. PROPOSED CALIBRATION METHOD
Due to the growing demands of high-speed communication and high-precision sensing applications, ultrawideband antenna arrays operating at millimeter-wave or THz frequencies have become increasingly more important. It is more efficient to carry out calibration for such arrays in the time domain. As mentioned earlier, the method based on SP in conjuncture with a deterministic approach was proven to be effective in the frequency domain. In this article, the SPSA method is implemented in the time domain for the first time. Fig. 1 shows the setting to calibrate an antenna array with M elements in the time domain. The proposed calibration method is put to a test by employing a short pulse signal in the time domain illuminating the aperture array being measured in a known plane wave. It should be noted that although the calibration waveform is known, the delay between the antenna for transmitting the source signal and the antenna elements of array is not exactly known due to the parasitic random time delays that are caused by the imperfect receiving components and propagation environment. Therefore, the proposed calibration algorithm strives to perform fine tracking on the time delay in each channel and compensate them before acquisition of the received signals at the array output.
Let s(t) denote a time-domain signal, which is a Gaussian pulse as given by where σ 2 p denotes the variance of the Gaussian pulse and w controls the time width. This pulse signal propagated from a transmitter in the far field of an aperture array to be calibrated. The average power of the transmitted pulse signal for calibration can be expressed by Subsequently, the signal received by the mth element of the aperture array before calibration and the synthesizer can be represented by where A m denotes the channel gain between the transmitter of reference and the mth antenna element, consisting of antenna gains and path loss. Moreover, ξ m denotes the spatial time delay of the signal received by the mth element associated with physical distance, and γ m denotes the intrinsic time delay (generated randomly) introduced by other factors from devices and environment. Therefore, the vector T spatial consisting of time delays for M elements is defined as For example, when the angle of incidence, θ = 30 • , and the element spacing of d = 0.015 m, T spatial = [0, 0.025, . . . , 0.375] for M = 16. Random time delays are introduced into each channel to represent the imperfect array construction or distortion caused by the receiver devices, as given by where γ m can be generated by γ m = 2x − 1 and is in the interval (−1, 1) and x = rand(1) is a uniformly distributed random number in the interval (0, 1). The interval for γ m can be adjusted to reflect the actual time delay distortions in real systems. Consequently, the vector of the actual time delays for all channels can be expressed by Let the sampling frequency for the time-domain signal be f s , and hence, the time resolution of the signal is given by Assuming that the total length of the time-domain signal in a processing frame for any channels is L and M is the total number of array elements, then the average received power from the array without calibration can be calculated by The respective time of arrival for the signal received in each channel varies with imperfect array construction and/or other factors. Hence, the synthesized signal in the array receiver without calibration is significantly different from that after all element channels are calibrated.
In order to make an approximation on the objective vector of time delays for calibration,T * , which is used to compensate for the random time delays occurring in the antenna array system within each channel, T in ,T * is defined aŝ where τ m (m = 1, 2, 3, . . . , M) denotes the approximated time delay calibration coefficient for the element m.
The SPSA algorithm has been designed to estimate the distorted time delay in each element channel. Since there are many receiving elements in the array, the problem of calibration has been converted into a multivariate stochastic search and optimization problem. The key problem for the algorithm to address is to seek the values of variants to minimize the value of a loss function. The loss function is crucial for the calibration performance and can be defined in different ways. A more detailed description on loss function definition will be presented in Section V. To illustrate the renewed steps of the SPSA algorithm for array calibration, tentatively, the loss function was assumed to be related to power at array output, and it can be defined as where P k is the average power measured at the array output whenT k is used to compensate the time distortions, T in , and P ideal is the desired power at array output when the time distortions in the element channels are fully eliminated. The step-by-step description given next shows how the SP algorithm was implemented for aperture array calibration in the time domain.
Step 1 (Initialization and Setup Algorithm Coefficients for Perturbation): Set time of arrival of impinging signals for each channel due to difference in propagation path length, T spatial . Then, select the coefficient constants before iteration begins and generate the random time delay vector representing the imperfect channel condition, T in .
Step 2 (Generation of SP Vector, an M-Dimensional Random Perturbation Vector k ): The time delays for all channels in each iteration are decided by the random perturbation vector and the coefficients chosen at the beginning, and a simple choice of the perturbation vector is to use a Bernoulli distribution with probability of 0.5 for each ±1 outcome. An example of k is k = [−1, 1, . . . , −1].
Step 3 (Power Measurements at the Array Output and Evaluations): Obtain two measurements from the power calculation function as shown in (17) based on the SP around the time delays for the current iteration, k, l(T k + c k k ) and l(T kc k k ).
Step 4 (Gradient Approximation): Generate the SP approximation to the unknown gradientŝ where ki is the ith component of the k vector with the total length of M and k is given in Step 2. It is noted that the denominator is the same for all components ofĝ k vector within one iteration and reflects the characteristics of the SP in contrast to other approximation methods.
Step 5 (Updating the Time Delay Estimate): Use the SA form to calculate the time delays for the following iteration step:T Step 6 (Repeat From Step 2 for More Iterations or Terminate): Terminate the algorithm if the condition is met such as l ≤ σ , and l is the loss function as defined in (19) or set a limit for the maximum number of iterations.

V. NUMERICAL ANALYSIS
This section presents the detailed steps, including objective function definitions, coefficient selection strategies, and performance analysis of the proposed calibration algorithm. It was examined for calibration of aperture arrays consisting of various numbers of elements aiming for different levels of accuracy. Once the objective function is defined, the strategy of giving values to the optimization coefficients a k and c k is essential for successful convergence.
The study began with the following assumptions. The source signal given by (10) with 1-ns duration impinged on the planar surface of the aperture array with a particular angle of incidence from a point at far field, and the time of arrival for the signal received on each element varies with the relative distance among elements. The time delays related to the relative distance were calculated by using the parameters, including angle of incidence and the element spacing. Furthermore, the time delay on each channel is distorted by a combination of many factors, including transmit/receive (TX/RX) module impairments, feed circuit variations, mutual coupling, and diffraction due to antenna structures. In the simulation, time delay distortion in each channel is collectively represented by a vector of random numbers, T in , subject to Gaussian distribution with zero mean and within a specified time range.
Definition of the loss function is critically important for the proposed algorithm to be successfully implemented as Sum up the signals from each channel for two separate time delays after applying the simultaneous perturbation, the synthesized signal at array out with backward perturbation, S b (t) = M m=1 s m (t + ξ m +T k (m) + c k km ), whereas with the forward perturbation applied in each element channel, the synthesized signal at array output Obtain the two power measurements with the respective time perturbations, Calculate the gradient vector for the current, the kth, iteration based on the two power measurements from the last step, the gradient vector is,ĝ k (T k ) = Calculate the time delay adjustment vector for the next iteration,T k+1 =T k − a kĝk (T k ), let k = k + 1 12: return Time delay offset for each channel,T * ∼ =T k the outcome of approximation heavily depends on it. Two loss function options were explored in this work: power measurement-based objective function, and product of power measurement and correlation coefficient calculation at array output. The following analysis was made to provide general guidance on efficient implementation of the array calibration method in the time domain based on the SPSA algorithm.

A. Power Measurement-Based Objective Function
The objective function can be decided by the received power at array output as defined in (19). The guidelines from the basic SPSA algorithm for coefficients selection were tested for the specific array calibration problem and presented in the following, but it had been warned that they are guidelines only-may not be the best or even work for every application. In this section, new coefficient selection strategies for effective array calibrations were established and compared with that based on the guidelines from the basic SPSA algorithm.
In order to estimate the magnitude of variations of the loss function at the beginning of the search, the power measurements can be obtained through several sets of T * 0 vectors with its components formed by randomly generated time delays within each channel. The corresponding values for the loss function are calculated by using (19) and arranged in a vector.
1) Basic SPSA Guidelines: As shown in the guidelines of the SPSA method by Spall [11], the choice of these two parameters is related to the specific application of the algorithm, particularly, the magnitudes of the parameters that are to be optimized, i.e., the distorted time delays in each channel for the problem considered in this article. Initial values of a k and c k are suggested in [11] by and where a, c, A, α, and γ are the nonnegative coefficients and can be selected before the optimization iterations begin. Initial values of α and γ are chosen at 0.602 and 0.101, respectively, as suggested in [11]. It was recommended that the parameter c can be approximately equal to the standard deviation of measurement noise. It can be estimated by using the variation of the loss function magnitudes that can be obtained by giving several guessed values to the parameters to be optimized. To reach convergence effectively, prior to implementing iterations, initialization coefficients for the algorithm have to be selected with care as they depend on applications. Section V-A2 addresses the unique coefficient selection process concerning array calibration problems.
2) Optimal Initialization Coefficients Selection: Recall that the relationship between a k and c k must ensure the convergence condition, i.e., ∞ k=0 a 2 k /c 2 k < ∞. Let the ratio between them simply be b k = a k c k (25) and the values of b k were closely observed during estimation steps to ensure the convergence. It is delicate to prevent b k from decaying too fast with an increasing k while maintaining a certain rate of change.
In the early iterations, it is important to keep the change of time delays under control. Before iterations begin, several guessed values can be given to T * 0 to evaluate the magnitude of loss function, approximated gradients, together with the desired magnitude of change in time delay at early iterations, and initialization coefficients can be determined. The detailed procedure to calculate the initialization coefficients based on guesses on T * 0 was presented in Tables II and III. Consequently, a k and c k can then be calculated at each iteration according to (23) and (24). It is noted that the optimal value of a k is related to the magnitude of the loss function, and c k must be small enough to ensure that the perturbation is within  the allowed range of the parameters to be optimized, i.e., the maximum possible time distortion in the system. Therefore, the values ofĝ k have to be closely monitored.
The algorithm terminates when the difference in measured power at the current iteration and the desired power at the array output is close to each other. However, the solution to resolve the time delays based on power measurement alone may need a large number of iterations to converge. Therefore, the scenario-loss function based on the product of power and correlation coefficient-was introduced as another alternative to calculate gradients and establish the termination condition of the algorithm. This will be analyzed in Section V-B.

B. Correlation Coefficient-Related Objective Function
As shown in Section V-A, an objective function based on power measurements as the termination condition for the algorithm requires a large number of iterations to reach convergence. Therefore, an objective function was proposed with an extra parameter-correlation coefficient is introduced to indicate and quantify the similarity of two signals in the time domain. The correlation coefficient is defined to analyze the similarity between two vectors and is given by x = [x 1 , x 2 , . . . , x l , . . . , x L ] and y = [y 1 , y 2 , . . . , y l , . . . , y L ] as wherex andȳ denote the means of vector x and y, respectively. During the approximation process, r (x, y) can be calculated by assuming that x is the synthesized signal from the array output with an ideal receiving condition (i.e., the optimal signal desired at the array output) and y is the synthesized signal vector from the array output when the time offset coefficients, T k , is applied in element channels. The new objective function with an extra factor of correlation coefficient as an indicator for convergence testing is defined as (27) and when the objective function l approaches zero, the synthesized signal at the array output with the corrected time delay for each channel becomes nearly identical to the ideal signal desired out of the receiver. Once this is achieved, the time delay coefficients vector for calibration,T * , takes values of the approximated vector of time delays from the output of the algorithm,T k , which is acquired from the final iteration when the convergence condition is met. Based on the coefficients selection procedure described in Section V-A, the quasi-optimal coefficients for the algorithm were given in Table IV when the number of elements in the arrays was M = 4, 16, . . . , 128. The range of time delay distortions (bounded by a hypothetical maximum time delay distortion) was assumed between −1 and 1 ns. In general, the values for a k , b k , and c k are suggested to decrease as the iteration number, k, increases in order to reach the convergence asymptotically. We found that they can also be assigned with predetermined small constant values, such that when M = 128, a k = 0.1, and c k = 0.01, the convergence condition (l ≤ 0.001) can also be reached.
The implementation of the algorithm with the loss function defined by the combination of power and correlation coefficient shows a faster convergence speed. The comparison of convergence rate using two loss function definitions is shown in Fig. 2. In the simulation, M = 128 and the range of time delay distortion in the 128-element channels was between −1 and 1 ns. All the coefficients were kept the same apart from the calculation of the loss functions.

C. Analysis on Number of Iterations
The total number of iterations needed in the algorithm to reach the convergence depends on a number of factors: the assigned values for the initialization coefficients, the total    Fig. 3. It indicated that the number of iterations required increases significantly as the convergence condition becomes more stringent. The benefit of a smaller value as convergence condition was potentially compromised by a number of iterations needed. Hence, a balance has to be sought. More novel gradient approximation approaches can be explored to accelerate the convergence, such that when the initial value of c and corresponding c k are fixed, if b k follows a polynomial delay, then a k can be derived according to (25).

D. Analysis on Accuracy of the Algorithm
The accuracy of the algorithm is closely related to the convergence condition. The ratio of discrepancy, between the intrinsic time delays (the hypothesized time delay distortions for real systems) T in and the estimated time delaysT * , with the range of time distortions, is used as a measure for evaluating accuracy of the algorithm (mean{∥T * −T in ∥}/range{T in }). With the number of elements varying between 4 and 128, after the specified convergence was met (l ≤ 0.2, l ≤ 0.1, or l ≤ 0.01), the corresponding estimation errors are shown in Fig. 4. It indicated that average estimation errors in arrays with a different number of elements maintain a similar level as soon as the convergence condition is fixed. This can be better explained by observing the two particular cases (the insets) where the array had 4 and 128 elements. When the convergence condition was met for each case (under l ≤ 0.01), the estimation errors for all element channels distributed evenly of a similar magnitude (with 5% in error) despite the number of the elements in the arrays.

E. Receivers With True Time Delay
Wideband antenna arrays can be dealt with by adopting true time delay (TTD). Initially, switched delay lines were applied in the signal paths, and quantized delays were applied on the element or subarray levels. The recent developments of monolithic microwave integrated circuit (MMIC) chips introduce various amounts of delays over a broad range of frequency. The delay steps are controlled digitally with typical 8-12 bits. However, the range of time delays in a single chip is limited. The cascaded two-stage hybrid transmitter/receiver architecture is shown in Fig. 5. At present, the hybrid analog/digital configuration is preferred in practice to reduce the number of analog-to-digital converters (ADCs) in wideband systems. Unlike the conventional quasi-monochromatic arrays Fig. 4. Estimation error as the number of elements in arrays and convergence conditions varies. The estimation errors were averaged in all element channels and the ratio in percentage between the average time delay errors and the range of time distortions, 2 ns, was given. The inset on the left was for an array with four elements and the one on the right for an array with 128 elements. Both were converged when l was less than 0.01. The loss function was based on the combination of power and correlation coefficient.
where a digital beam can be steered by a simple assignment of phase and amplitude for each element at a single frequency, calibration must be carried out in the time domain for wideband arrays in order to balance the phase and amplitude over the entire bandwidth. Quantization on time delays may create errors as not the exact time delay offsets can be applied during iterations. The estimation was implemented with time delay steps controlled by an 8-bits circuit and the estimation error was compared with that using fully analog time delays in Fig. 6. The convergence (l ≤ 0.01) was reached first with an analog time delay setting, and then, the same number of iterations was used in implementation with an 8-bits control circuit for time delay adjustment.

VI. SIMULATION RESULTS
In order to validate the effectiveness of the proposed method, two scenarios have been considered to carry out calibration in the time domain on finite arrays. In the first scenario, it was assumed that the pulse transmitted for calibration was received without any changes in the element channels, while in the second scenario, pulse distortion caused by transmitting and receiving was taken into consideration. Hence, the influence of pulse distortion on performance of the proposed calibration method was included in the study.

A. Unchanged Pulses in the Element Channels
In order to illustrate how to carry out calibration in practice and verify its effectiveness in solving array calibration problems, the SPSA algorithm was adopted to carry out calibration on an aperture array with 16 elements in the time domain. The range of time delay distortion in each element channel is assumed between −1 and 1 ns. In the experiment, the loss function was defined by using the product of received power at the array output and the correlation coefficient between the received signal at the array output (with perturbated time delays in element channels) and the desired signal. The detailed flowchart of running the algorithm for the experiment is denoted in Fig. 7, where l was calculated based on (27) and evaluated against two convergence conditions (l reached 0.2 and 0.005) in each iteration.
The total number of elements in the array for calibration is set to M = 16. A Gaussian pulse with the variance of σ 2 p = 0.01 impinged onto the array from a far field, and the duration of time for the pulse is 1 ns. Initially, a time delay vector, T in , consisting of 16 elements was generated, of which each element was a random number distributed in the range between −1 and 1, representing time delays caused by element variation and imperfect channel condition. Thus, the received signal from each element was defined as: (2)) , . . . , S 16 = s(t − T in (16)), and the combined signal from the 16 channels prior to implementing the calibration procedure will be S no_cal = S 1 + S 2 + · · · + S 16 .
Approximated time delays were introduced in each channel for SP, and iterations were run through until the combined signal from the array output, with perturbation on time delay applied, was similar to the ideal signal (no time delays applied in element channels). When the convergence condition was met, the received signal from the array output was: (16))), and the combined signal can be expressed as During the intermediate steps of the optimization process, T k was renewed at each iteration by introducing time delay perturbation in each channel for the following step; the correlation coefficient between the received signal at the array output (with time delay perturbation in each element channel), S calibrated , and the ideal signal at the array output (without time delays applied in each element channel) was used to establish the condition for convergence.
Initially, a time delay vector T in with 16 elements was generated, a random value between −1 and 1 was assigned for each element to reflect the imperfect receiving condition of the array, and these are the time delays to be calibrated out of the system. The initial values of T in , time delay vector, is given in Table V. As shown in Fig. 7, the SPSA algorithm itself has no prior information about these arbitrary values of time delays. It introduced time delays into channels at each iteration and carried out measurements at the array output. Once the synthesized signal (based on the estimated delay corrections from the algorithm) has a close enough correlation with the combined signals (based on the random time delays at the beginning), i.e., the correlation coefficient between them is greater than a set value, the convergence criterion was met and the time delay vector,T * , from the last iteration is the time delays to be used to make time delay adjustments. The synthesized signals based on T in before and after calibration are compared in Fig. 8. The time delays for all receiving channels are to be compensated based onT * , and they are used to cancel T in . The approximated values ofT * from the algorithm were compared to the intrinsic time delays randomly   Table V, and the estimation error is much lower than 0.01 ns when the range of time deviation is between −1 and 1 ns, which verified the effectiveness of the proposed method.

B. Distorted Pulses in the Element Channels
In practice, the pulses used for calibration will inevitably be distorted to some extent due to limitations in the antenna array and frequency-dependent fading during propagation. Hence, it is necessary to take the signal distortion effects into   (10) was used as the reference signal (w = 0.3 and σ 2 p = 0.01) and the received signals from the element channels had a similar form but with different variance values (σ 2 p varying between 0.0217 and 0.0625) representing distortion effect in different element channels. Second, the Richard wavelet, the second derivative of the Gaussian pulse as in (10), , was used as the reference signal (α = 24) for calibration, whereas the pulses received from element channels used wavelets of different α values (α varying between 8 and 23 for the 16 elements). Finally, the short pulse signal was generated from the signal generator AWG5204 (from Tektronix) with the pulsewidth of 400 ps in the setting, and the received pulses were taken from the output ports of the power splitter where the random propagation delay in each channel was introduced later on after the respective pulse waveforms were collected. The transmitted pulse and the received pulses in the time domain and frequency domain are shown in Fig. 9. It can be seen that the pulse waveform was modified during propagation and this was more clearly observed in the frequency domain. The imperfect array antenna condition was further illustrated by the phenomenon that the received pulses were not showing up in a regular interval. The source pulse was used as a reference signal to launch the algorithm with the abovementioned received pulse distributing among the element channels with random time delays. The calibration results for the three examples Fig. 9. Illustration of pulses transmitted, the reference pulse waveform for calibration, and the received pulses from a power splitter emulating the pulses received from the element channels, and the source pulse from AWG5204 has a width of 400 ps. (a) Transmitted and received pulses in the time domain, the repetitive pulses at the transmit side were shown indicating a regular interval where the received pulses were deviated from the regular intervals, and the time deviations are to be calibrated. (b) Power spectrum of the transmitted and received pulses. mentioned above are given in Fig. 10. It indicated that the accuracy of the time delays estimated is better than 99% even though the pulses received from the element channels were different from the reference pulse waveform for calibration. The calibration performance was compared with other typical calibration methods performed in the frequency domain, and the main parameters are summarized in Table VI. VII. DISCUSSION An idealized solution to the general problem of calibration is to determine a number of unknowns (presumably constant over a period of time at least) from a minimum number of measured quantities (independent of each other.) In a broadband aperture array with a large number of elements, a number of measurements required to resolve calibration coefficients are significant and difficult to implement practically; in particular, expensive near-field facility testing is infeasible for arrays with large dimensions. The method based on time-domain operations presented in this article is an attempt to cope with Fig. 10. Calibration performance when the received pulses in the element channels were different from the reference pulse. Two sets of random delays for 16 element channels were used for verification, and they were defined as cases 1 and 2. (a) Reference pulse for calibration is a Gaussian pulse or a Richard wavelet, and the received pulses in the element channels were different to the reference pulses (numerical pulse waveforms). (b) Reference pulse is the pulse generated from the signal generator with the pulsewidth of 400 ps, and the received pulses in the element channels were distorted and therefore different from the reference pulse (real pulse waveforms). these challenges. The accurate calibration of aperture array antennas does not require the full capability of a near-field range (NFR) facility to physically sample the RF field in front of the array in half-wavelength increments or the response of each radiating element as is required in most existing calibration methods.
The characteristic of the short pulse employed for calibration has a connection with the operational frequency band of aperture arrays. It should be decided following the principle that a minimum distortion is caused to the pulse while receiving from the array elements. However, the filtering effect of array antennas on the signal for calibration can be mitigated by using the pulse signal received from an array element as the reference signal. The algorithm can cope with pulse distortions without a significant compromise on estimation performance as shown in Fig. 10 where distorted pulses from array elements were considered during calibration.
The same method can be applied in the frequency domain by changing the signal for calibration from a short pulse to a narrowband signal in order to determine optimal phase offset coefficients in element channels. The loss function defined in (19) can be adopted. As this is not the focus of this article, we do not elaborate further on this. However, we note that a very high level of accuracy can be achieved if the method is implemented in the frequency domain using such method.
Aperture array technology is evolving with advances in solid-state microwave integrated circuits. Element-level digital beamforming (DBF) becomes increasingly realistic. However, digital subarray architecture represents a feasible compromise between the number of beams and the number of digital channels. For such arrays as shown in Fig. 5, the algorithm can be launched in each subarray, and it can also be applied to make time delay compensations between the subarrays before beams from subarrays being digitized.
For most time delay RFICs, the adjustable delay range is between 0 and 200 ps. The state-of-the-art MMICs have a maximum time delay in the order of 1000 ps covering various frequency bands. The maximum time distortion in element channels the algorithm dealt with was 2000 ps (between −1 and 1 ns) in the analysis. However, the time delay range can be adjusted based on practical needs. The number of iterations will reduce when the range of time delay distortions occurring in array element channels is smaller as it is easier to converge.
It is worth mentioning that the method proposed in this work is applicable for the calibration of broadband antenna arrays employing time delay units at the element or subarray level. Before the utilization of fully digitized antenna arrays becomes more accessible, broadband array systems with TTDs are being developed to mitigate or eliminate "squint" phenomenon by adopting phased-steered approaches for scanning, which works well for a system with a narrow RF bandwidth. However, for applications where high resolution is desired such as imaging and tracking, more radical changes in technology, including integrated circuits of TTD with a higher time delay resolution, such as 8-12-bits time delay shifters, are expected. These higher cost systems are needed to provide the performance requirements for high-resolution applications and the calibration method proposed in this work can greatly reduce the complexity of calibrating such systems. Various antenna array architecture solutions with integrated TTD units are presented in [45].
In practice, the quantization step size of TTDs plays an important role on the accuracy of the proposed method. When the number of control bits for the time delay shifter reduces, the calibration error will rise as the accurate time delay coefficients in element channels are more difficult to acquire and this leads to ambiguity. One reference case has been examined in the previous discussion where the errors of time delay coefficients for calibration at several element channels (in an array with 128 elements) are approximately 10% with 8-bits resolution TTD units.
Antenna arrays can have complex and extended timedomain responses-impulse responses that have the potential to significantly alter the transmitted or received waveforms, either because of the environment or the array itself. On the other hand, by using phase shifters, a significant squintinduced broadening on the main beam will occur when the array has a wide bandwidth; hence, TTDs have to be adopted to restore beamwidth performance. The required TTD, T , and phase shift, ϕ, have the following relation: For a broadband waveform with a pulsewidth in time τ p ≪ Mdsin(θ )/c, even if dispersion exists in the Tx/Rx paths, which occurs because of the reflections or other effects, as the time resolution is high, the responses irrelevant to calibration can be filtered out in the time domain. In addition, the settings used for calibration are mainly considered to be under lineof-sight (LOS) propagation. As mentioned in Section VI-B, the pulse shapes with distortion due to transmitting, propagation, and receiving have a minor effect on calibration performance as soon as the reference pulse signal for calibration is narrow enough compared to the maximum time delay distortion occurring in the element channels. This suggests that the suitable pulse waveform can be chosen based on the following principles: 1) the time-domain waveform covers the 3-dB bandwidth of the operational frequency of the array and 2) the duration of the pulse is significantly less than the maximum time delay distortion in the element channels. One approach can be achieved by measuring the impulse response of one element in the array and then using this impulse response as the reference pulse for calibration, where a vector network analyzer with time-domain transforming capability can be used.

VIII. CONCLUSION
A new calibration method for aperture arrays in the time domain has been proposed. It exploits the SA algorithm with SP. The proposed method can be successfully applied for calibrating antenna arrays of ultrawide bandwidth in the time domain. The main contributions of this article can be summarized as follows: 1) it gives the technical details of recurrence relation needed to implement the algorithm for the purpose of array calibration in the time domain and the strategy on selecting the initial values of the coefficients to start the algorithm; 2) it demonstrates that two types of loss functions can be employed to establish the convergence criterion based on the availability in practice; 3) it provides guidance for the optimal values of the coefficients and their relationship with the problem; 4) small arrays or large arrays can be treated in different tiers by assigning different values for coefficients; and 5) it shows that the convergence can be reached at ease to ensure an estimation error within ±1% of the total time duration of the short pulse that was employed for calibration. This validated the efficacy of the method for calibrating aperture arrays and opens new options for online calibration of ultrawideband arrays without the need for resorting to expensive and time-consuming near-field facility testing. It holds great potential for millimeter or THz applications where increased frequency bandwidth and real-time processing are in need.
The scope of the method is for array systems with TTD units at the element or subarray level, where the time resolution is high and it should be significantly smaller than the maximum time distortion occurring in the element channels. The future work includes a waveform design for optimal calibration performance and applies the proposed method for wideband beamforming under more complex situations.