Critical Review of Basic Methods on DoA Estimation of EM Waves Impinging a Spherical Antenna Array

: Direction-of-arrival (DoA) estimation of electromagnetic (EM) waves impinging on a spherical antenna array in short time windows is examined in this paper. Reﬂected EM signals due to non-line-of-sight propagation measured with a spherical antenna array can be coherent and/or highly correlated in a snapshot. This makes spectral-based methods inefﬁcient. Spectral methods, such as maximum likelihood (ML) methods, multiple signal classiﬁcation (MUSIC), and beamforming methods, are theoretically and systematically investigated in this study. MUSIC is an approach used for frequency estimation and radio direction ﬁnding, ML is a technique used for estimating the parameters of an assumed probability distribution for given observed data, and PWD applies a Fourier transform to the capture response and produces them in the frequency domain. Although they have been previously adapted and used to estimate DoA of EM signals impinging on linear and planar antenna array conﬁgurations, this paper investigates their suitability and effectiveness for a spherical antenna array. Various computer simulations were conducted, and plots of root-mean-square error (RMSE) against the square root of the Cram é r–Rao lower bound (CRLB) were generated and used to evaluate the performance of each method. Numerical experiments and results from measured data show the degree of appropriateness and efﬁciency of each method. For instance, the techniques exhibit identical performance to that in the wideband scenario when the frequency f = 8 GHz, f = 16 GHz, and f = 32 GHz, but f = 16 GHz performs best. This indicates that the difference between the covariance matrix of the signal is coherent and that the steering vectors of signals impinging from that angle are small. MUSIC and PWD share the same problems in the single-frequency scenario as in the wideband scenario when the delay sample d = 0. Consequently, the DoA estimation obtained with ML techniques is more suitable, less biased, and more robust against noise than beamforming and MUSIC techniques. In addition, deterministic ML (DML) and weighted subspace ﬁtting (WSF) techniques show better DoA estimation performance than the stochastic ML (SML) technique. For a large number of snapshots, WSF is a better choice because it is more computationally efﬁcient than DML. Finally, the results obtained indicate that WSF and ML methods perform better than MUSIC and PWD for the coherent or partially correlated signals studied.


Introduction
Direction-of-arrival estimation is a popular topic in various electromagnetic-related fields of study, finding applications in military surveillance, radar, sonar, and mobile communication systems [1][2][3][4][5][6][7][8][9]. Understanding the directions-of-arrival (DoAs) of incoming signals on receiving antenna can be effectively used to localize the positions of the corresponding sources. It facilitates adaptive beamforming of the receiving pattern to improve the sensitivity of the system towards the desired directions of signals and mitigates undesirable interferences. Invariably, it ensures that the array provides the maximum beam towards the desired users and nulls in the directions of interferences, and it improves the performance of mobile and base stations. Therefore, DoA estimation of electromagnetic (EM) waves impinging on an antenna array is very important.
A spherical antenna array (SAA) configuration is useful in obtaining an array with isotropic characteristics. There is also the spherical phased antenna array, which is widely used in spacecraft and satellite communication. Electronic beam scanning by the phased antenna array remains a better choice because it provides a hemispherical scan and almost uniformly distributed gain [3,10]. Another major advantage of the SAA configuration is its three-dimensional symmetry, which is an advantage in the spatial analysis of signals. The geometry and size usually provide the directivity and degree of beam scanning in the elevation and azimuth needed from the array. A linear array, which scans in one axis, and a planar array, which scans in the x-and y-axes, are acceptable for a bounded scan demand of 60-70 • [10]. The projected apertures of these arrays drop by the cosine factor of the scan angle. The array directivity is reduced for a large scan angle. For larger steerability, a multifold planar configuration or, in a general form, an SAA is preferred [10][11][12]. The directivity of an SAA remains almost the same when an active source is moved over the surface of the SAA. Since signals measured with an SAA can be coherent and highly correlated, it becomes important to investigate the effectiveness of the available DoA estimation methods with an SAA.
Various DoA estimation methods have been reported in the literature. The most commonly used methods are: the multiple signal classification (MUSIC) algorithm [13,14], estimation of signal parameters through rotational invariance technique (ESPRIT) [15,16], beamforming, and the maximum likelihood (ML) DoA estimator [17,18]. They have been deployed and applied to antenna arrays with arbitrary geometries and provide an appropriate estimate of DoA. Recently, the use of phased-array-based sub-Nyquist sampling [19] was proposed. The temporal properties of the impinging signals were explored to improve the detection capabilities. To date, linear and planar antenna array configurations have been considered in most DoA estimations. However, recent publications emphasize the importance and rapid development of SAAs, particularly in satellite communication [3,10].
Multipath signals can be coherent or highly correlated in a narrow band of frequency. Coherent signals can cause spectral-based methods to fail in DoA estimation in linear or planar arrays [1,11], and hence, these methods are generally employed with spherical arrays [11,20,21]. The performance of MUSIC and ESPRIT for the estimation of correlated signals can be improved by using spatial smoothing techniques [22][23][24]. In addition, different smoothing techniques have been proposed for general SAAs [21,25]. The smoothing techniques preprocess the array output covariance matrix over time [26], space [21,24], or frequency [25] before the estimation of DoA. In SAA processing, the array is usually divided into subarrays by transforming the SAA into a uniform linear array [21,27]. Another spatial smoothing approach is the formation of a subarray using eigenbeam space rotation [28,29]. The free space response may be different over time (i.e., DoA, amplitude, and phase) due to refractivity and absorption. The frequency resolution is reduced by averaging over the frequency, while averaging over time reduces the temporal resolution of the estimate. The spatial resolution is reduced by averaging over space because the number of elements for each subarray is less than the number of elements in the entire array.
Contrary to MUSIC and other similar methods, ML methods [30,31] with no smoothing can handle coherent signals. The coherent signal model can be U-dimensional (where U is the number of signals due to multipath) when each signal is differently modeled. If a spatial spectrum is computed with a grid size D, where the correct size and sampling are user-defined, then the DoA estimation using the ML function has a dimensionality of U × D, but spectral-based techniques have 1 × D dimensionality. Therefore, the search space quickly becomes very large as the number of signals increases. For ML methods [31], high dimensionality gives rise to the deployment of nonlinear optimization algorithms for reduced computational time. The general limitation of all methods is that the estimable Electronics 2022, 11,208 3 of 25 number of signals must be less than the number of elements I, i.e., I > U [30]. Moreover, in practice, the free space response is always time-and frequency-dependent.
This paper studies the DoA estimation of EM signals in free space trapped with an SAA. Specifically, we consider the case where more than one signal (i.e., in a multipath scenario) is present in a single analysis window. The number of signal points or responses due to multipath at the receiving antenna within a short time window rises with the square of time t. This impact is usually termed the density of reflection [20]. Therefore, with time, it results in a situation where a single analysis window encompasses multiple signals. There is a complete overlap of signals if the distance between the generating source and the receiving antenna is the same for one or more signals. This usually occurs for first-order signals in a symmetric source-receiver. We examined the ML DoA estimation methods within the context of the SAA's spatial response and large sample approximation (LSA), known as weighted subspace fitting (WSF). LSA assumes that the available number of measurements is large. Because of the aforementioned characteristics, the conditions in which one or more signals are present in one analysis window for narrowband and wideband frequencies are considered in this work. The key innovations and contributions of this study are as follow: Five DoA estimation techniques are presented and investigated in this study. MU-SIC, plane-wave decomposition (PWD), weighted subspace fitting (WSF), and both deterministic (D) and stochastic (S) ML techniques were applied in the analysis of EM waves using SAA processing for the first time in this study. WSF is the large sample approximation of SML. II. We adopted ML techniques in the spatial signal analysis with SAA and confirmed their suitability in DoA estimation using an SAA. It is numerically demonstrated how ML techniques outweighed MUSIC and beamforming methods and how WSF produces the best DoA estimation among the ML techniques. III.
The results obtained from ML techniques are compared to beamforming and MUSIC methods throughout the simulation and analysis section. The ML-based DoA estimation is more appropriate, less biased, and more robust against noise than beamforming and MUSIC techniques. In addition, the WSF technique exhibits lower computational complexity than DML and SML. IV.
We explain how DoA estimation performance with a rigid sphere is expressed as a function of frequency. In addition, the radius of the sphere affects the estimation accuracy at a particular frequency. For instance, for a given SAA with a radius of 8.4 cm, the most accurate performance was achieved at 16 GHz. If similar performance were required, for example, at about 8 GHz, then a larger radius of SAA would be required.
The major novelty of our work as compared with Tervo and Politis [11] is that we particularly extended the DoA methods to the antenna array rather than the microphone array, backed with experimental data in electromagnetics vis-à-vis a spherical antenna array, which is, in the end, the ground truth to test and evaluate any method. This is the major contribution of our work.
The rest of this paper is organized as follows. Section 2 present a general description of the SAA configuration. In Section 3, mathematical formulations for EM signals in free space, the space domain, spherical harmonic decomposition, and SAA with spherical harmonic decomposition (SHD) are derived. In Section 4, five DoA estimation methods are analyzed: MUSIC, PWD, WSF, and deterministic (D) ML and stochastic (S) ML methods. The numerical experimental tests, results, and discussion are presented in Section 5, and finally, conclusions based on the findings are given in Section 6.

General Configuration of SAA
The basic configuration and geometry of the SAA are presented in Figure 1, where θ i and φ i represent the elevation and azimuth of the ith element, respectively. θ d and φ d denote the elevation and azimuth of the incident plane wave, respectively, while s d is the Electronics 2022, 11, 208 4 of 25 source signal. The radiating antenna elements are evenly distributed over a sphere of radius r. There is a conventional distribution for an antenna array on a geodesic sphere depending on 1 of the 13 Archimedean solids or the 5 platonic solids [3,11]. Here, the fundamental merit is that all elements encounter the same neighborhood. Conversely, other various distributions could be determined based on the beam-scanning requirement and the number of antenna elements in use. For an icosahedron distribution along the elevation axis, the elements are organized in circular bands. In each circular band, the locations in azimuth are clearly estimated. Furthermore, any far-field observation point can be defined as p(r, θ, φ). The basic configuration and geometry of the SAA are presented in Figure 1, where and represent the elevation and azimuth of the ith element, respectively. and denote the elevation and azimuth of the incident plane wave, respectively, while is the source signal. The radiating antenna elements are evenly distributed over a sphere of radius r. There is a conventional distribution for an antenna array on a geodesic sphere depending on 1 of the 13 Archimedean solids or the 5 platonic solids [3,11]. Here, the fundamental merit is that all elements encounter the same neighborhood. Conversely, other various distributions could be determined based on the beam-scanning requirement and the number of antenna elements in use. For an icosahedron distribution along the elevation axis, the elements are organized in circular bands. In each circular band, the locations in azimuth are clearly estimated. Furthermore, any far-field observation point can be defined as ( , , ). and represent the elevation and azimuth of the ith element, respectively. and denote the elevation and azimuth of the incident plane wave, respectively, while is the source signal.
The SAA is in active mode, where each element is supplied with an amplifier. To receive a beam pattern in a specific direction ( , ), the elements lying within the predetermined cone angle along the beam direction are in active mode. These elements are activated using accurate phases estimated from the equations in [10,32,33]. Furthermore, the SAA configuration (as in Figure 2) is 3-D symmetrical, which is an advantage in the analysis of the spatial signal. Both its size and geometry produce the directivity and beamscanning level in the azimuth and elevation angles needed from the array. For more steerability, a multifold planar configuration or, in a more general form, an SAA is preferred [3][4][5][6][7]. The directivity of an SAA is almost the same when an active source is moved over the surface of the SAA.
(a) Figure 1. Spherical array geometry system for localization. θ i and φ i represent the elevation and azimuth of the ith element, respectively. θ d and φ d denote the elevation and azimuth of the incident plane wave, respectively, while s d is the source signal.
The SAA is in active mode, where each element is supplied with an amplifier. To receive a beam pattern in a specific direction (θ 0 , φ 0 ), the elements lying within the predetermined cone angle along the beam direction are in active mode. These elements are activated using accurate phases estimated from the equations in [10,32,33]. Furthermore, the SAA configuration (as in Figure 2) is 3-D symmetrical, which is an advantage in the analysis of the spatial signal. Both its size and geometry produce the directivity and beam-scanning level in the azimuth and elevation angles needed from the array. For more steerability, a multifold planar configuration or, in a more general form, an SAA is preferred [3][4][5][6][7]. The directivity of an SAA is almost the same when an active source is moved over the surface of the SAA. The basic configuration and geometry of the SAA are presented in Figure 1, where and represent the elevation and azimuth of the ith element, respectively. and denote the elevation and azimuth of the incident plane wave, respectively, while is the source signal. The radiating antenna elements are evenly distributed over a sphere of radius r. There is a conventional distribution for an antenna array on a geodesic sphere depending on 1 of the 13 Archimedean solids or the 5 platonic solids [3,11]. Here, the fundamental merit is that all elements encounter the same neighborhood. Conversely, other various distributions could be determined based on the beam-scanning requirement and the number of antenna elements in use. For an icosahedron distribution along the elevation axis, the elements are organized in circular bands. In each circular band, the locations in azimuth are clearly estimated. Furthermore, any far-field observation point can be defined as ( , , ). and represent the elevation and azimuth of the ith element, respectively. and denote the elevation and azimuth of the incident plane wave, respectively, while is the source signal.
The SAA is in active mode, where each element is supplied with an amplifier. To receive a beam pattern in a specific direction ( , ), the elements lying within the predetermined cone angle along the beam direction are in active mode. These elements are activated using accurate phases estimated from the equations in [10,32,33]. Furthermore, the SAA configuration (as in Figure 2) is 3-D symmetrical, which is an advantage in the analysis of the spatial signal. Both its size and geometry produce the directivity and beamscanning level in the azimuth and elevation angles needed from the array. For more steerability, a multifold planar configuration or, in a more general form, an SAA is preferred [3][4][5][6][7]. The directivity of an SAA is almost the same when an active source is moved over the surface of the SAA.

EM Waves in Free Space
The behavior of EM waves in free space can be described as the measured response from a particular signal source to antenna i within a particular space. The plane wave propagates through free space and is trapped at the receiver by multipath. While propagating in free space, the plane wave is affected by various environmental phenomena, such as absorption, diffraction, and reflection [20]. These environmental phenomena alter the phase and amplitude of the propagating plane wave, possibly differently per incidence angle and per frequency.
In addition, the plane wave is altered by the travel distance and attenuated by the channel. Specifically, the absorption rate depends on the composition of the air, distance traveled by the plane wave, and operating frequency [19,20]. Therefore, the absorption by air alone suggests that the signal response varies with time and with frequency and is hence the basis for analysis in the frequency domain. Furthermore, if the signal source is considered a point source, the amplitude of the signal experiences attenuation because the plane wave is spherically spread out. The total number of signals or snapshots U due to multipath per time interval ∆ can be expressed asymptotically as [20] where c is the speed of light, and V is the enclosure volume. This is applicable and correct for all geometries with a homogeneous medium. To reduce the total number of signals within a particular analysis window, a shorter time interval ∆ is usually used. Furthermore, when investigating the response of a signal with a single frequency, the time window length is given as ∆ > 2 ⁄ to ensure that one period of the wavelength can be obtained within the analysis window. Hence, the time instant where U or a smaller number exists in an analysis window of a single frequency can be expressed as where denotes the angular frequency.

EM Waves in Free Space
The behavior of EM waves in free space can be described as the measured response from a particular signal source to antenna i within a particular space. The plane wave propagates through free space and is trapped at the receiver by multipath. While propagating in free space, the plane wave is affected by various environmental phenomena, such as absorption, diffraction, and reflection [20]. These environmental phenomena alter the phase and amplitude of the propagating plane wave, possibly differently per incidence angle and per frequency.
In addition, the plane wave is altered by the travel distance and attenuated by the channel. Specifically, the absorption rate depends on the composition of the air, distance traveled by the plane wave, and operating frequency [19,20]. Therefore, the absorption by air alone suggests that the signal response varies with time and with frequency and is hence the basis for analysis in the frequency domain. Furthermore, if the signal source is considered a point source, the amplitude of the signal experiences attenuation because the plane wave is spherically spread out. The total number of signals or snapshots U due to multipath per time interval ∆t can be expressed asymptotically as [20] U ∆t = 4π where c is the speed of light, and V is the enclosure volume. This is applicable and correct for all geometries with a homogeneous medium. To reduce the total number of signals within a particular analysis window, a shorter time interval ∆t is usually used. Furthermore, when investigating the response of a signal with a single frequency, the time window length is given as ∆t > 2π/ω to ensure that one period of the wavelength can be obtained within the analysis window. Hence, the time instant t U where U or a smaller number exists in an analysis window of a single frequency can be expressed as where ω denotes the angular frequency.

Analysis of the Space Domain
The position of the antenna's ith element is presented in the 3-D Cartesian coordinate form as where has propagated a distance d u from the source to the antenna with a DoA Ω u = (θ u , φ u ) w.r.t. the array origin and can be expressed in Cartesian coordinates as Each signal is considered a plane wave affected by the free space phenomena mentioned above.
A wideband pressure signal p i in the ith antenna for a particular signal source s is depicted in the frequency domain as the addition of all snapshots or EM signals u = 1, . . . , U: where . h iu (k) is the frequency response of a snapshot (an event consists of multiple waves), while n i (k) is a noise variable, which is assumed to be identically distributed and independent for each antenna element. In addition, wave number k = ω/c = 2π f /c. The signal source . s u (k) characterizes the frequency response of the source, which is originally transmitted in a few directions, arriving at the receiving antenna array from direction r u .
For a snapshot with multiple signals, the frequency response depends on the time delay τ iu (k), which defines the travel time of the signal to the antenna and the complex amplitude of each event c iu , i.e., For a plane-wave model to be applied, it is assumed that both signals and sources are in the far field w.r.t. the receiving array. Thus, in line with the plane-wave model, the complex amplitude is expressed as where . c iu is the response of receiving antenna i in the Ω u direction and is assumed to be known a priori, and c u denotes the complex amplitude of the plane wave u, which factors all signal phenomena experienced before reaching the receiving array. In addition, due to the assumption of a short time window analysis and a homogeneous medium, the pathway distance d u is therefore neglected in the directional analysis. As a result, the time delay of the plane wave can be represented with respect to the origin (that is, the mid-point of the sphere in the SAA) as where T is the normal to the wavefront. The SAA response and corresponding delay term is presented as and referred to as the steering vector. For instance, . c iu = 1 for an ideal antenna array. The source response and amplitude of the plane wave are given as the product Electronics 2022, 11, 208 7 of 25 and regarded as the actual signal. For brevity, the array measurements are presented in matrix form. The impulse responses of array elements i = 1, . . . I, i.e., their inputs, are a vector [31]: The signals and noise are assumed to be zero-mean complex Gaussian processes, i.e., [25] E s(k)s H (k) = R ss (k) and E n(k)n H (k) = σ 2 I where E{·} is the expectation, and (·) H is a Hermitian matrix transpose. This results in the following covariance matrix of the array [31]: The assumption of a Gaussian signal in free space is not necessarily true. That is, a signal with a single frequency in a short time window can possibly be deterministic in nature instead of random. The deterministic model of the array covariance matrix is presented in Section 4.5.

Spherical Harmonic Decomposition
Using an SAA, the array covariance matrix and pressure are analyzed within the context of spherical harmonic decomposition (SHD). In this section, the presented formulation follows those given in [26,[32][33][34][35]. The SHD model of the pressure signal p(k, r, Ω) for an antenna array with harmonic order N h and radius r is expressed using spherical Fourier transform (SFT) approximation and the corresponding inverse [26]: and where Y m n (·) is the SH of degree m and order n [34,35]: p nm is the spherical harmonic domain coefficient, Φ i is the coordinates of the antenna, Ω = (θ, φ) is the direction of steering where the inverse SFT is formulated, and α i is a sampling weight used to correct errors due to orthonormality. Moreover, P m n is associated with Legendre polynomials, while (·) * represents the complex conjugate. The sampling scheme defines the sampling weights, with the order defined by the number of elements (check [25] for details). Uniform sampling is assumed throughout this paper; therefore, the weight is reduced to α i = 4π/I. We also adopt harmonic coefficients (N h + 1) 2 = 16, where the array harmonic order is N h = 3, and radius r = 8.4 cm, depending on the antenna array under consideration. In spherical arrays, whenever kr > N h , spatial sampling occurs [25]. The noiseless pressure in the spherical harmonic domain is expressed in matrix form as [26] where is the position of the element within the angular coordinates. The lefthand side of Equation (17) is the vectorized form of SH coefficients, i.e., k and r are not included for brevity's sake, and SH are presented in matrix form as Hence, antenna-array-dependent coefficients are described using a diagonal matrix: where each array-dependent coefficient in the case of a rigid sphere is [19] b n (kr) = 4πi n j n (kr) − j n (kr)h n (kr) h n (kr) (21) h n , j n , h n , and j n represent the Hankel and Bessel functions and the corresponding derivatives with respect to kr, respectively.

Spherical Antenna Array (SAA) with Spherical Harmonic Decomposition (SHD)
The application of theoretical discrete spherical harmonic transformation (SHT) to the receptions of a real antenna array as in Equation (15) leads to serious loss of accuracy in the transformed coefficients, as compared to the theoretical result of the antenna array of equal specifications to the actual one [11,12,25]. This is because of the differences between the real antenna array and the theoretical antenna array model, such as calibration errors between antenna elements and antenna mutual coupling effects and positioning error. Conversely, the SHT performance approximates the theoretical counterpart using real array calibration measurements in an anechoic chamber.
The SHT of Equation (15) is expressed as an array of encoding filters, and the equalization with inverse modal weights B −1 can be optimally obtained within the context of least-squares. This is the same as the challenge faced when trying to obtain the spherical steering vector that is optimal w.r.t. the measured properties of the array. In this work, such a measurement-dependent steering vector is adapted. Section 5.3 presents the SAA measurement procedure under consideration.
LetŶ(Ω) represent the measurement steering vector. In an ideal scenario, by approximation, the steering vector is equal to the theoretical steering vector for a plane wave in the transform domainŶ Equation (17) is the response of a theoretically transformed antenna array, and after equalization, it can be incorporated into the componentŝ where h(Ω) represents the array response vector to the Ω direction or the steering vector, and W enc is the encoding matrix that implements the equalization and SHT. Ideally, (22) holds correctly. In a practical situation, both the interpolated array response for an arbitrary Ω and encoding filters can be calculated from response measurements in L directions Ψ = [Ψ 1 , . . . , Ψ L ] arranged uniformly in a grid. We can then compute the encoding filters using the solution of regularized weighted least-squares from measurements [34][35][36]: where Γ ψ is the correct weight for the measurement grid, andĤ(k) = ĥ 1 , . . . ,ĥ L is the set of array-measured snapshots at frequency k in the direction Ψ, while β is a set of regularization parameters. In this study, we computed the weights of Γ ψ via the areas of spherical Voronoi cells of the measurement grid.
In order to acquire the steering vector h(Ω) of the antenna array with high accuracy in any direction, spherical interpolation was conducted through the expansion of the steering vector as a function of its measured spherical harmonic coefficients, as the SH matrices Y are computed up to an order driven by the measurement number. It is discovered that an accurate truncation order N tr is accurately described by (N tr + 1) 2 = L/4 for a regular grid. By visual inspection of the energy of spherical harmonic coefficients, it can be deduced that coefficients above this truncation order approach the noise level. Combining Equations (24) and (25), the final interpolated spherical harmonic coefficients arrived at arê

DoA Estimation
Five DoA estimation methods are presented in this section. MUSIC, PWD, WSF, and both deterministic (D) and stochastic (S) ML techniques were adapted in the analysis of signals using an SAA for the first time in this work. WSF is the large sample approximation of the SML technique, as mentioned earlier.
There exist major differences between WSF, MUSIC, PWD, and ML techniques. Apart from the fact that PWD and MUSIC are not ML techniques, an advantage is that in WSF and ML, one can freely choose both the noise and signal model. In MUSIC, PWD, and other related techniques, we have to look for U maxima from a spatial spectrum of D size, while in maximum likelihood techniques, we only look for a single minimum (or equivalent maximum) from an estimation function of U × D size. It then emerges that ML techniques suffer from the dimensionality curse very quickly as U rises. An extensive search in this greater-dimensional space is computationally forbidden, but nonlinear optimization algorithms are often used to evaluate the minimum within a considerable period. The nonlinear optimization algorithms for WSF and ML techniques are presented in Section 5.1.
Two scenarios are of interest: analysis of signals in a single frequency and wideband analysis of signals. In the case of wideband analysis, a time-domain smoothing algorithm is applied, as in [26,35], which is briefly reviewed in this section. In narrowband analysis, the smoothing algorithm does not provide any advantages, because the signals are coherent. Hence, analysis is conducted for spherical harmonic coefficients in the frequency domain.

Smoothing in Time Domain
The formulations in [26] are followed in the processing of SH domain coefficients. The normalized spherical harmonic coefficients are achieved by multiplication of the left-hand side of Equation (17) by B −1 (kr), which results in where z(k) = B −1 Y H (Φ)Γ n (k). The covariance matrix of the array in this scenario is expressed as while the covariance matrix of the noise is which depends on k, and ||·|| is the Frobenius norm. The corresponding time-domain representation of q nm (k) is obtained through the inverse Fourier transform F −1 {·} as q nm (t) = F −1 {q nm (k)}. Nevertheless, multiplying by B −1 results in unwanted noise amplification in some frequencies [26]. Hence, the same equalization is used in timedomain smoothing, as given in [26]: where σ 2 γ is the constant normalization factor that is the same for all frequencies. In the case of the antenna array considered in this work, the H eq (k) magnitude response has the shape of a high-pass filter. As a result of equalization, the frequency-domain covariance matrix of the array is where R ss (k) = R s(k)s H (k) = ||H eq (k) || 2 R ss (k) (33) and it represents the equalized forms of variables possessing ·. In our estimation, it is assumed that the noise matrix in the time-domain form is spatially white and timeindependent (i.e., R zz (t) = R zz = σ 2 a I, and σ 2 a is the equalized variance). The covariance matrix estimate of the array for N snapshots (time instant t 1 , . . . , t N ) iŝ q nm (t) = F −1 H eq (k) q nm (k) denotes the equalized time-domain spherical harmonic coefficients.
The performance in single-frequency bins is also of interest in this work. We introduce spatial whitening to the steering vectors and covariance matrix estimate in the frequencydomain estimate as in [26]. Therefore, other formulations in this section are provided for time-domain smoothing and are equal to the frequency-domain smoothing equivalent if each whitened form is a substitute for the steering vectors and the covariance matrix estimate.

Steering Vector
The Hermitian transpose of Equation (19) is the steering matrix or vector employed in all of the estimation techniques. For instance, in WSF and ML techniques, the matrix form of the three-signal scenario is the following: The number of harmonic components is (N h + 1) 2 = 16, because only two signals Ω 1 and Ω 2 are present, in contrast to PWD and MUSIC, whose steering is consistently 1-D with the form The above matrices show the basic difference between the 2-D and 1-D search space used for WSF, ML techniques, and spectral-based techniques, respectively.
For WSF and ML techniques, the localization function P(Ω) is 2U-dimensional because Ω has U azimuth and elevation values. The localization function of the minimum argument for WSF and ML techniques produces DoA estimates: The localization function P(Ω) for spectral-based techniques is 2-D, because Ω has two dimensions: azimuth and elevation. Spectral-based techniques for DoA estimates are the U highest local minima in P(Ω).

Beamforming
For many decades, beamforming has been a well-known technique for the estimation of direction [30]. PWD is a frequently used beamformer [37]. PWD has a maximum rejection characteristic of isotropic spatial noise and is expressed over frequencies k = 1, . . . , K aŝ wherep nm (k) is the noisy form of Equation (18), and the weighting is and the unity beamformer in the look direction is ensured by the first term. This given PWD does not whiten the noise and hence exhibits poor performance in the wideband scenario. However, noise whitening can be implemented as given in [19]. PWD can be implemented here using time-domain smoothing aŝ The output energy of the time-domain beamformer at Ω is whereR aa is the estimated covariance matrix of the array defined in Equation (34).

MUSIC
Another important DoA estimation method is MUSIC [11,16]. The MUSIC spatial spectrum is computed as where E n is the noise matrix of the array obtained from eigenvalue decomposition of covariance matrix estimateR aa aŝ where λ represents the eigenvalue matrix, E incorporates the right eigenvectors, and subscripts n and s are the noise and signal subspace, respectively. MUSIC needs a covariance matrix with a full-rank signal. Hence, when we employ MUSIC in localizing U signals, the covariance matrix of the estimated signal should have U eigenvalues that differ from the noise. However, if there are few snapshots or coherent signals of the array covariance matrix, this assumption does not hold. The snapshot refers to frequency-and time-domain spherical harmonic coefficients.

ML Techniques
Generally, ML techniques are used in many estimation scenarios. Readers are referred to [11] for a basic review of ML methods. The ML technique requires both noise and signal models. These models depend on the parameters estimated by the ML technique. Finding the parameters of the models that will mostly explain the observed data is the problem in ML estimation. This paper applies two earlier developed ML techniques in the subspace domain to SAA processing.

Stochastic ML Technique
Stochastic ML (SML) is the first ML technique with the assumption that signals are stochastic in nature and are stochastically processed. In the case of time-domain smoothing, the array covariance matrix is For time instant t i = t 1 , . . . , t N , the SML probability density function is given as where R aa denotes the modeled covariance matrix of the antenna array, |·| is the matrix determinant, and A N = [ q(t i ), . . . , q(t N )] are equalized SH coefficients in the time domain. Usually, in ML techniques, we find the solution from the negative log-likelihood for SML as log p A N /Ω, R ss , σ 2 The negative log-likelihood solution with respect to S, σ 2 ([31] provides the details), results in and

Y(Ω)
are the orthogonal projector and the pseudo-inverse of Y(Ω) onto the null space of Y H (Ω), respectively. The SML localization function, as in [31], is

Deterministic ML Technique
There are no assumptions regarding signal waveforms for the deterministic ML (DML) model (i.e., average signal waveforms have this form: E{ q(t i )} = Y H (Ω) s(t i )). The covariance matrix of the array depends on the noise term when the average signal waveform is derived from the signals. According to [31], the DML likelihood function for different snapshots is expressed as . . , s(t N )] are signals due to multipath. Setting Ω and S N constant from the negative log-likelihood, the variance can be computed as in [31]: The localization function and signal estimates can be deduced using the negative log-likelihood, which gives rise to a nonlinear least-square problem [31]: respectively. Asymptotically, MUSIC is equivalent to DML for large sizes of samples and uncorrelated signals [31] (i.e., in the incoherent case).

Weighted Subspace Fitting Technique
Subspace fitting techniques are suboptimal approximations of the aforementioned ML techniques. They are useful because they have less computational complexity than maximum likelihood techniques. Furthermore, for larger sample sizes, if certain conditions are met, they are asymptotically the same as ML techniques [31].

Estimation Accuracy: Cramér-Rao Lower Bound
Estimation theory provides performance metrics for techniques by comparing their covariance estimation error against a theoretical lower bound. The Cramér-Rao lower bound (CRLB) is a generally employed bound for the estimation of covariance matrix error [16], [31].
For the scenarios considered in this paper, the error covariance of unbiased estimatê For compactness purposes, Ω and k are excluded from the notation in the following. The above provides the stochastic CRLB by substituting the SML probability density function into Equation (44): where Y i is the partial derivative of Y with respect to the ith element in Ω and defined as Y i elements are the SH partial derivatives with respect to elevation angle θ and azimuth angle φ and expressed as and respectively, and csc(·) denotes the cosecant trigonometric function. All of the partial derivatives are presented in the literature (e.g., MATHEMATICA). The CRLB for a single stochastic source in the spherical harmonic domain is given in [35]. The DML probability density function in Equation (52) can be substituted into Equation (55) to provide the deterministic CRLB [25]:

Signal Detection
Knowledge of the number of signals U is a requirement in DoA estimation. Various techniques used for the computation of the number of signals (i.e., detection) have been given; therefore, readers are encouraged to visit [30] for details. Certain signal detection techniques estimate the subspace dimension from the eigenvalue matrix, for instance, by statistically testing the number of eigenvalues belonging to the noise space [31]. For partially correlated source signal detection, the best methods are model-based methods, e.g., WSF detection and the generalized likelihood ratio test [31]. These methods estimate the DoA and detect the number of U simultaneously. Furthermore, it is expected that these approaches work very well for correlated or coherent signal detection. This paper does not examine signal detection, but U is assumed to already exist as prior knowledge, just as in previous works on this topic, such as [2,13,[36][37][38][39].

Numerical Experiments, Results, and Discussion
Simulations and experiments from measured data with an SAA configuration are described in this section. We used a spherical array with 32 elements on the surface of a rigid sphere in all cases considered. A well-detailed technical description, element positions, etc., for the SAA are given in [3]. The results are examined on the basis of the root-mean-square error (RMSE): and consequently compared against the square root of CRLB. Furthermore, in all of the simulations, we averaged the RMSE over 100 Monte Carlo samples.
The signal-to-noise ratio (SNR) is given in this paper as the space-domain SNR, which implies SNR = 1/σ 2 = σ −2 . Nevertheless, a more defined SNR value is the effective SNR, which can be computed as the relationship between the equalized noise variance and equalized signal variance, i.e., SNR e = s(k) || 2 / σ 2 a . In all cases considered in this paper, we simulated the noise as space-domain i.i.d. complex Gaussian random variables with variance σ 2 .

Using Nonlinear Optimization Techniques for Global Minimum/Maximum Search
The localization functions are nonlinear and consequently require the application of nonlinear optimization techniques if simulation time is required. Similar to [31], this study used a Newton-type search algorithm to find the minimum of WSF, DML, and SML localization functions, where the LM (Levenberg-Marquardt) method is employed. Here, instead of LM, we employed a quasi-Newton (QN) technique with the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, which performs a Hessian estimate by finite differences. The QN technique was selected because it is available in MATLAB as the fminunc function. The QN technique iteration is terminated if the change in estimated parameters and the function value does not exceed 1 × 10 −10 between two or more successive iterations. In addition, the allowable maximum iteration number is fixed at 30 for the QN technique.
The QN technique requires intelligent initial prediction for parameter estimation. Thus, we initialized the parameters so that they were similar to true values by adopting the estimate from the alternating projection (AP) algorithm. In practice, true values may not be available; hence, we propose the approach in [31] as a suggestion. The DML, SML, and WSF initial value was searched using the AP algorithm (i.e., a nonlinear optimization algorithm) [17]. AP always converges to a local minimum and not the global minimum. Moreover, here, we assumed the AP-based local minimum to be the global minimum of the employed search grid. In addition, we used U + 2 iterations in the AP algorithm with initialization and a grid of 5 0 resolution with respect to elevation and azimuth.
In all of the simulation scenarios, the QN technique and BFGS were used to estimate the maxima of PWD together with MUSIC localization functions, and parameters were initialized close to the true values. PWD and MUSIC can be initialized in practice from individual spatial spectra by estimating the U highest local maxima. From the same grid employed above in the alternating projection algorithm, these spatial spectra can be computed.

Single Wideband Signal Scenario
A scenario where there is a single signal in a single time window was first simulated. The length of the window N = 2048 and sampling frequency f s = 48 MHz. Furthermore, the signal was simulated as a wideband plane wave s(t) = δ(t − 1024/ f s ) with harmonic order N h = 3 and coming from Ω = (90, 0) • . In addition, spherical harmonic coefficients in the time domain were computed as presented in Section 4.1 using an N-size discrete Fourier transform. The signal covariance matrix for a plane wave is and E eq ≈ 0.88 is the equalized signal energy in total.
The RMSE of an individual technique over 100 Monte Carlo samples for different noise variances σ −2 is shown in Figure 3. Moreover, Figure 3 shows the theoretical CRLB performance for the deterministic and stochastic cases. For one signal scenario, it can be observed that both techniques have almost the same performance. The RMSE of the techniques adheres to the CRLB when SNR ≤ −10 dB. The elevation angle estimation has a smaller RMSE, just as in the CRLB prediction.

Two Wideband Signals
Four cases were simulated in the second example, where two signals arrive during the same time window. Likewise, as shown in [26], two wideband signals were simulated with one signal delayed by d samples; i.e., the first and second delayed signals are ( ) = ( − 1024 ⁄ ) and ( ) = ( − (1024 + ) ⁄ ) , respectively. When the delay approaches zero, the signal covariance matrix is coherent, and when d is increased, the covariance matrix of the signal is incoherent. The four delays considered here are = 0, = 0.5, = 1, and = 8 samples, and the covariance matrices of the signals are respectively.
The two signals were simulated as plane waves with DoAs Ω = (90,0) and Ω = (90, ) and with harmonic order N = 3, and the second plane-wave direction was changed with 2 × . , where = 0, 0.5, 1, … ,7( ). The value of SNR was then set at 20 dB. The SH coefficients in the time domain were computed as in the previous section and windowed by N = 35 samples of a lengthy rectangular window. The window was concentrated around the first signal at 1024 ⁄ . As a result, the covariance matrix estimate of the array was averaged with the time domain N = 35 samples, as given in Equation (34).
The RMSEs for CRLB and all of the techniques over 100 Monte Carlo samples against when = 20 and the second signal is delayed by d samples are shown in Figure 4. Firstly, it can be observed that PWD performs poorly in the estimation of DoA of multiple signals because its increases as its RMSE increases. This can be interpreted as a bias in estimation: the PWD global maximum is found in between the true values, and the second maximum is 180 apart from the second one; the two signals generate a

Two Wideband Signals
Four cases were simulated in the second example, where two signals arrive during the same time window. Likewise, as shown in [26], two wideband signals were simulated with one signal delayed by d samples; i.e., the first and second delayed signals are s 1 (t) = δ(t − 1024/ f s ) and s 2 (t) = δ(t − (1024 + d)/ f s ), respectively. When the delay approaches zero, the signal covariance matrix is coherent, and when d is increased, the covariance matrix of the signal is incoherent. The four delays considered here are d = 0, d = 0.5, d = 1, and d = 8 samples, and the covariance matrices of the signals are R ss = E eq 1 1 1 1 , R ss = E eq 1 0.85 0.85 1 , R ss = E eq 1 0.5 0.5 1 , and R ss = E eq 1 −0.04 −0.04 0.97 respectively.
The two signals were simulated as plane waves with DoAs Ω 1 = (90, 0) • and Ω 2 = (90, φ 2 ) • and with harmonic order N h = 3, and the second plane-wave direction was changed with 2 1+n k ×0.2 , where n k = 0, 0.5, 1, . . . , 7( • ). The value of SNR was then set at 20 dB. The SH coefficients in the time domain were computed as in the previous section and windowed by N w = 35 samples of a lengthy rectangular window. The window was concentrated around the first signal at 1024/ f s . As a result, the covariance matrix estimate of the array was averaged with the time domain N w = 35 samples, as given in Equation (34).
The RMSEs for CRLB and all of the techniques over 100 Monte Carlo samples against φ 2 when SNR = 20 dB and the second signal is delayed by d samples are shown in Figure 4. Firstly, it can be observed that PWD performs poorly in the estimation of DoA of multiple signals because its φ 2 increases as its RMSE increases. This can be interpreted as a bias in estimation: the PWD global maximum is found in between the true values, and the second maximum is 180 • apart from the second one; the two signals generate a maximum in the spatial spectrum of PWD in between the true values of φ 1 and φ 2 . Due to the large energy that this maximum has, RMSE is lower than CRLB, φ 2 ≤ 3 • with d = 0. In addition, the next PWD spectrum maximum is generated by the steering vector's conjugate value that generated that global maximum. MUSIC has an identical profile to PWD, and its estimation is identically biased when d = 0. Similarly, MUSIC and PWD performance improves slightly when φ 2 > 60 • . Consequently, the segregation is longer than the Rayleigh resolution 180 • /N h = 60 • , and the techniques segregate the two signals. Nonetheless, estimates of PWD and MUSIC are also biased when φ 2 > 60 • .
WSF and ML techniques result in smaller values than CRLB in the case of coherence when φ 2 ≤ 3 • , that is, at d = 0. This is because of the bias in estimation. All of the techniques have their global minimum within the true values and showcase compelling proof for a biased estimate, as in PWD. Hence, the information in second-order covariance matrix moments is highly coherent for an estimation that is not biased when φ 2 ≤ 3 • and d = 0. Firstly, the steering vectors are identical when signals come from various directions that are close to each other. Secondly, the signals are coherent. Obtaining better and quality performance when φ 2 ≤ 3 • would require more antenna elements. When d = 0, for φ 2 > 3 • , WSF and ML techniques clearly perform better than MUSIC or PWD and follow the CLRB.
WSF and ML techniques share the same performance in the partially correlated cases, i.e., d = 0.5 and d = 1. WSF and ML techniques have slightly lower RMSE than MUSIC. This means that, in partially correlated cases, WSF and ML techniques have better performance than MUSIC. ML techniques, WSF, and MUSIC share the same performance in the almostincoherent scenario with d = 8. The performance of PWD when d = 0 is identical to that when d = 0.5, d = 1, and d = 8.
Low variance is assumed by CRLB [26]. Based on [31], RMSE and CRLB agree by approximation when the DoA theoretical standard deviation of the signals is less than the half-angle of separation. This also holds for the cases studied. When φ 2 = 3 • and d = 0, the CRLB is 2.1 • for the first signal and 2.3 • for the second signal, which is higher than half of the separation angle 3/2 = 1.5 • .

Two Signals in One Frequency
Here, the techniques' performance based on a single frequency is examined. The noise is whitened, as stated above, because of the applied frequency-domain processing [25]. That is, the steering vectors and whitened frequency-domain covariance matrix are employed against the time-domain forms (Equation (27)). Spatial aliasing should be considered when estimating DoA with an SAA in a single frequency. Typically, when kr > N h , spatial aliasing exists in high frequencies [36]. However, WSF and ML techniques are not affected by spatial aliasing, because B(kr) does not have any zeros (i.e., continuous), and Y(Ω) and Y(Φ) are the same for all of the frequencies.
When we conduct an analysis of signals in a single frequency, the resulting unavoidably coherent covariance matrix is and the covariance matrix estimate of the array becomeŝ In this simulation, the frequencies chosen are the resonant frequencies that fall within a band. It is worth mentioning that covariance matrices of the array are not averaged over the chosen band, but we examined the performance in single frequencies, as given in Equation (62).
Wideband signals were simulated with d = 0, N = 2048, N h = 3, and 20 dB SNR. The SH coefficients in the frequency domain were estimated by Equation (27) through the discrete Fourier transform (DFT), and frequency bins closest to frequencies f were adopted in this analysis. The resolution of the frequency is f s /N, and the bandwidth of each frequency analyzed is about 2 GHz due to windowing "spectral leakage" or aliasing. Figure 5 shows the simulation results for different frequencies. The techniques exhibit identical performance to that in the wideband scenario when = 8 GHz, f = 16 GHz, and f = 32 GHz, but f = 16 GHz performs best. The reason why CRLB > RMSE when f = 16 GHz, f = 32 GHz, f = 40 GHz, and f = 75 GHz is the same as in the aforementioned wideband scenario. This indicates that the difference between the covariance matrix of signals is coherent, and the steering vectors of signals impinging from that angle are small. MUSIC and PWD share the same problems in the single-frequency scenario as in the wideband scenario when d = 0.
Furthermore, when = 900 MHz, f = 2 GHz, and f = 4 GHz, the estimations from all of the techniques are biased, as RMSE has an average of around 2 • for all φ 2 values. At these frequencies, half of the separation angle is lower than CRLB for all separation angles studied. Therefore, we expect RMSE and CRLB not to meet when there is high error variance. The poor effective SNR causes poor performance in these frequencies. The effective SNR in low frequencies is smaller than in higher frequencies because of the array geometry. Figure 6 shows the simulation result for RMSE against SNR for MUSIC, PWD, DML, SML, and WSF. WSF shows better performance in the low-SNR regime. DML and SML show a high level of similarity to WSF. DML and WSF show greater performance than SML. This is because SML converges to a local minimum but does not exhibit convergence to the global minimum. MUSIC and PWD exhibit worse performance compared to other techniques, as in other cases. WSF and DML exhibit identical behavior. Hence, they are good choices in multiple signal analysis, especially in the single-frequency analysis. Furthermore, WSF is more computationally efficient for large U than DML. In the wideband scenario, there will be no wide gap in the performance of any of the techniques. and = 32 GHz, but = 16 GHz performs best. The reason why CRLB > RMSE when = 16 GHz, = 32 GHz, = 40 GHz, and = 75 is the same as in the aforementioned wideband scenario. This indicates that the difference between the covariance matrix of signals is coherent, and the steering vectors of signals impinging from that angle are small. MUSIC and PWD share the same problems in the single-frequency scenario as in the wideband scenario when = 0. Furthermore, when = 900 , = 2 , and = 4 , the estimations from all of the techniques are biased, as RMSE has an average of around 2 0 for all values. At these frequencies, half of the separation angle is lower than CRLB for all separation angles studied. Therefore, we expect RMSE and CRLB not to meet when there is high error variance. The poor effective SNR causes poor performance in these frequencies. The effective SNR in low frequencies is smaller than in higher frequencies because of the array geometry. Figure 6 shows the simulation result for RMSE against SNR for MUSIC, PWD, DML, SML, and WSF. WSF shows better performance in the low-SNR regime. DML and SML show a high level of similarity to WSF. DML and WSF show greater performance than SML. This is because SML converges to a local minimum but does not exhibit convergence to the global minimum. MUSIC and PWD exhibit worse performance compared to other techniques, as in other cases. WSF and DML exhibit identical behavior. Hence, they are good choices in multiple signal analysis, especially in the single-frequency analysis. Furthermore, WSF is more computationally efficient for large U than DML. In the wideband scenario, there will be no wide gap in the performance of any of the techniques. Generally, in this paper, the results obtained reveal that all of the techniques studied perform well for a single signal in the analysis window. However, if there are two highly  Furthermore, when = 900 , = 2 , and = 4 , the estimations from all of the techniques are biased, as RMSE has an average of around 2 0 for all values. At these frequencies, half of the separation angle is lower than CRLB for all separation angles studied. Therefore, we expect RMSE and CRLB not to meet when there is high error variance. The poor effective SNR causes poor performance in these frequencies. The effective SNR in low frequencies is smaller than in higher frequencies because of the array geometry. Figure 6 shows the simulation result for RMSE against SNR for MUSIC, PWD, DML, SML, and WSF. WSF shows better performance in the low-SNR regime. DML and SML show a high level of similarity to WSF. DML and WSF show greater performance than SML. This is because SML converges to a local minimum but does not exhibit convergence to the global minimum. MUSIC and PWD exhibit worse performance compared to other techniques, as in other cases. WSF and DML exhibit identical behavior. Hence, they are good choices in multiple signal analysis, especially in the single-frequency analysis. Furthermore, WSF is more computationally efficient for large U than DML. In the wideband scenario, there will be no wide gap in the performance of any of the techniques. Generally, in this paper, the results obtained reveal that all of the techniques studied perform well for a single signal in the analysis window. However, if there are two highly Generally, in this paper, the results obtained reveal that all of the techniques studied perform well for a single signal in the analysis window. However, if there are two highly correlated signals, both MUSIC and PWD estimations are biased, and ML techniques are better choices by far. All of the techniques improve when smoothing is utilized, which is depicted in Figure 4. In addition, the performance of all techniques would be improved if frequency-domain smoothing were applied in a manner similar to time-domain smoothing. Furthermore, for lower frequencies (f ≤ 900 MHz), the DoA estimation performance is poor. If correct windowing is used at the lowest frequencies, the estimate suffers the same problems as that at low frequencies. Therefore, an SAA with a larger radius should be used for multiple coherent signals in low frequencies.
As shown here, the DoA estimation performance of many electromagnetic waves with a rigid sphere is a function of the frequency analyzed. In addition, the estimation accuracy in a particular frequency is affected by the sphere radius, as depicted by CRLB. For the present SAA radius of 8.4 cm, the most accurate performance was achieved at 16 GHz. If similar performance is required, for example, at about 4 GHz, then a larger radius is needed. Therefore, better performance is achieved in lower frequencies as the radius becomes larger.
Finally, in order to fulfill advancing technology standards, particular elements are positioned in the systems. The distance between the elements becomes smaller, which leads to strong mutual coupling with poor impedance matching and radiation performance. Incorporating the mutual coupling effect existing between elements, experimentally measured data, which are the ground truth to test any procedure, were used to evaluate all of the methods. Hence, experimental measurement data were further used for performance evaluation and analysis. The SAA was situated at the center of the anechoic chamber, and the source was positioned at 74 DoAs, which were obtained from different combinations of 4 different elevations and 18 different azimuths. We chose the azimuths from 5 degrees to 365 degrees with 20 degrees as a step size. For detailed information on the measurement architecture, readers are referred to the previous paper [40], where the measured data were first published. Gross error (GE) performance evaluation metrics were equally employed to evaluate the methods. The comparison of performance between MUSIC, PWD, DML, SML, and WSF is shown in Figure 7. The WSF exhibits better performance. Both SML and DML show a high degree of similarity to WSF, but WSF and DML perform better than SML. This may be because SML converges to a local minimum but does not converge to the global minimum. PWD and MUSIC have worse performance compared to the other methods. Therefore, DML and WSF are good for the analysis of multiple signals, particularly in single-frequency analysis. accuracy in a particular frequency is affected by the sphere radius, as depicted by CRLB. For the present SAA radius of 8.4 cm, the most accurate performance was achieved at 16 GHz. If similar performance is required, for example, at about 4 GHz, then a larger radius is needed. Therefore, better performance is achieved in lower frequencies as the radius becomes larger.
Finally, in order to fulfill advancing technology standards, particular elements are positioned in the systems. The distance between the elements becomes smaller, which leads to strong mutual coupling with poor impedance matching and radiation performance. Incorporating the mutual coupling effect existing between elements, experimentally measured data, which are the ground truth to test any procedure, were used to evaluate all of the methods. Hence, experimental measurement data were further used for performance evaluation and analysis. The SAA was situated at the center of the anechoic chamber, and the source was positioned at 74 DoAs, which were obtained from different combinations of 4 different elevations and 18 different azimuths. We chose the azimuths from 5 degrees to 365 degrees with 20 degrees as a step size. For detailed information on the measurement architecture, readers are referred to the previous paper [40], where the measured data were first published. Gross error (GE) performance evaluation metrics were equally employed to evaluate the methods. The comparison of performance between MUSIC, PWD, DML, SML, and WSF is shown in Figure 7. The WSF exhibits better performance. Both SML and DML show a high degree of similarity to WSF, but WSF and DML perform better than SML. This may be because SML converges to a local minimum but does not converge to the global minimum. PWD and MUSIC have worse performance compared to the other methods. Therefore, DML and WSF are good for the analysis of multiple signals, particularly in single-frequency analysis.

Conclusions
DoA estimation of electromagnetic waves impinging on an SAA with a small number of snapshots was investigated. This paper studies DoA estimation of electromagnetic (EM) waves impinging on an SAA in short time windows, which finds applications in spacecraft and satellite communication. These techniques can avoid smoothing operations and cope with coherent signals. Multipath signal analysis in a short time window is of particular interest in this paper. Previously, ML techniques have been used within linear and planar antenna array configurations, but this paper presents them within the confines of an SAA configuration. Numerical computer simulations were conducted using the RMSE versus the square root of CRLB to evaluate the performance of MUSIC, PWD beamformer, stochastic and deterministic ML, and the subspace technique WSF (a large sample stochastic ML approximation) in direction estimation using an SAA. The results obtained indicate that the WSF and ML methods perform better than MUSIC and PWD in the coherent or partially correlated signals studied. They correctly estimated the DoA of various signals with adequate resolution, in contrast to PWD and MUSIC. Furthermore, deterministic ML or WSF performed better than stochastic ML because of lower convergence probability. For example, when = 900 MHz, f = 2 GHz, and f = 4 GHz, the estimation from all of the techniques was biased, as RMSE had an average of around 2 • for all φ 2 values. At these frequencies, half of the separation angle was lower than CRLB for all the separation angles studied. Therefore, we expect RMSE and CRLB not to meet when there is high error variance. The poor effective SNR causes poor performance in these frequencies. The effective SNR in low frequencies is smaller than in higher frequencies because of the array geometry. Therefore, ML techniques are good candidates for DoA estimation of electromagnetic waves impinging on an SAA, while WSF is a better choice if the computational complexity factor is of great importance.
In spite of the results presented in this study, there are a number of concerns that should be investigated in the future. For instance, we hope to further extend this work to the calibration and measurement of time delay in various applications, compare them to the state-of-the-art in terms of computational complexity, and evaluate their robustness to confounding factors that inevitably lead to alterations in practical SAA applications. Therefore, ML techniques, particularly DML and/or WSF, may spur positive development in SAA processing and its related applications, such as mobile communication systems, aerospace, and radar systems.