Next Article in Journal
Multi-Population Enhanced Slime Mould Algorithm and with Application to Postgraduate Employment Stability Prediction
Previous Article in Journal
Design and Analysis of a Voltage-Mode Non-Linear Control of a Non-Minimum-Phase Positive Output Elementary Luo Converter
Previous Article in Special Issue
Direction-of-Arrival Estimation of Electromagnetic Wave Impinging on Spherical Antenna Array in the Presence of Mutual Coupling Using a Multiple Signal Classification Method
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Review

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

by
Oluwole John Famoriji
* and
Thokozani Shongwe
Department of Electrical and Electronic Engineering Technology, University of Johannesburg, P.O. Box 524, Auckland Park, Johannesburg 2006, South Africa
*
Author to whom correspondence should be addressed.
Electronics 2022, 11(2), 208; https://doi.org/10.3390/electronics11020208
Submission received: 7 December 2021 / Revised: 24 December 2021 / Accepted: 28 December 2021 / Published: 10 January 2022
(This article belongs to the Special Issue Advances in Antennas and Wireless Propagation)

Abstract

:
Direction-of-arrival (DoA) estimation of electromagnetic (EM) waves impinging on a spherical antenna array in short time windows is examined in this paper. Reflected 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 inefficient. Spectral methods, such as maximum likelihood (ML) methods, multiple signal classification (MUSIC), and beamforming methods, are theoretically and systematically investigated in this study. MUSIC is an approach used for frequency estimation and radio direction finding, 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 configurations, 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 efficiency 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 fitting (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 efficient 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.

1. 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 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. MUSIC, 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.
  • 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.
  • 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.
  • 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.

2. 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 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 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.

3. Mathematical Formulations

3.1. 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 π c 3 t 2 V
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
t U = t ( ω ) U V ω 8 π 2 c 3
where ω denotes the angular frequency.

3.2. Analysis of the Space Domain

The position of the antenna’s ith element is presented in the 3-D Cartesian coordinate form as
m i = r i [ sin ( θ i ) cos ( ϕ i ) , sin ( θ i ) sin ( ϕ i ) , cos ( θ i ) ] T
where Φ i = ( θ i , ϕ i ) and r i are the general spherical coordinates, azimuth ϕ i [ 180 ,   180 ) ° , radius r i [ 0 , ) , and elevation θ i [ 0 ,   180 ] ° . Each signal event u arriving at the antenna 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
r u = d u [ sin ( θ u ) cos ( ϕ u ) , sin ( θ u ) sin ( ϕ u ) , cos ( θ u ) ] T
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 :
p i ( k ) = u = 1 U h ˙ i u ( k ) s ˙ u ( k ) + n i ( k )    
where h ˙ i u ( 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 τ i u ( k ) , which defines the travel time of the signal to the antenna and the complex amplitude of each event c i u , i.e.,
h ˙ i u ( k ) = c i u ( k ) exp ( i τ i u ( k ) )
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
c i u ( k ) = c ˙ i u ( k ) c u ( k )         i
where c ˙ i u 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
τ i u ( k ) = k r ¯ u T m i
where r ¯ u = [ sin ( θ u ) cos ( ϕ u ) , sin ( θ u ) sin ( ϕ u ) , cos ( θ u ) ] T is the normal to the wavefront. The SAA response and corresponding delay term is presented as
h ˙ i u ( k ) = c ˙ i u ( k ) exp { i k r ¯ u T m i }
and referred to as the steering vector. For instance, c ˙ i u = 1 for an ideal antenna array. The source response and amplitude of the plane wave are given as the product
s u ( k ) = c ¯ u ( k ) s ˙ u ( k )
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]:
p ( k ) = [ p 1 ( k ) p I ( k ) ] = u = 1 U [ h 1 , u ( k ) h I , u ( k ) ] s u ( k ) + [ n 1 ( k ) n I ( k ) ]                                                                                                         = [ h ( Ω 1 ) , , h ( Ω U ) ] [ s 1 ( k ) , , s U ( k ) ] T + n ( k )     = H ( Ω ) s ( k ) + n ( k )
Ω = [ Ω 1 , , Ω U ] is the set of all unknown DoAs.
The signals and noise are assumed to be zero-mean complex Gaussian processes, i.e., [25]
E { s ( k ) s H ( k ) } = R s s ( 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]:
R p p ( k ) = E { p ( k ) p H ( k ) } = H ( Ω ) R s s ( k ) H H ( Ω ) + σ 2 I
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.

3.3. 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]:
p n m ( k , r ) = i = 1 I α i p ( k , r , Φ i ) [ Y n m ( Φ i ) ] *
and
p ( k , r , Ω ) = n = 0 N h m = n n p n m ( k , r ) Y n m ( Ω ) ,   respectively
where Y n m ( · ) is the SH of degree m and order n [34,35]:
Y n m ( Ω ) = 2 n + 1 4 π ( n m ) ! ( n + m ) ! P n m cos ( c o s θ ) e i ω ϕ
p n m 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 n m 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   c m , depending on the antenna array under consideration. In spherical arrays, whenever k r > N h , spatial sampling occurs [25]. The noiseless pressure in the spherical harmonic domain is expressed in matrix form as [26]
P n m ( k r ) = Y H ( Φ ) Γ Y ( Φ ) B ( k r ) Y H ( Ω ) s ( k ) = B ( k r ) Y H ( Ω ) s ( k ) ,
where Γ = d i a g { [ α 1 , α 2 , , α I ] } is the sampling weights, Y H ( Φ ) Γ Y ( Φ ) = I , and Φ = [ Φ 1 , , Φ I ] is the position of the element within the angular coordinates. The left-hand side of Equation (17) is the vectorized form of SH coefficients, i.e.,
P n m ( k r ) = [ p 0 , 0 ,     p 1 , 1 ,   p 1 , 0 ,   p 1 , 1 ,     ,   p N h ,   N h   ] T ,
k and r are not included for brevity’s sake, and SH are presented in matrix form as
Y ( Ω ) = [ Y 0 0 ( Ω 1 ) ,   Y 0 0 ( Ω 2 )   Y 0 0 ( Ω u ) ,   Y 1 1 ( Ω 1 ) ,   Y 1 1 ( Ω 2 )   Y 1 1 ( Ω u ) ,   Y 1 0 ( Ω 1 ) ,   Y 1 0 ( Ω 2 )   Y 1 0 ( Ω u ) ,   Y 1 1 ( Ω 1 ) ,   Y 1 1 ( Ω 2 )   Y 1 1 ( Ω u ) ,   Y N h N h ( Ω 1 ) ,   Y N h N h ( Ω 2 )   Y N h N h ( Ω u ) ,   ] T
Hence, antenna-array-dependent coefficients are described using a diagonal matrix:
B ( k r ) = d i a g { [ b 0 ( k r ) ,   b 1 ( k r ) ,      b 2 ( k r ) , b 3 ( k r ) , , b N ( k r ) ] }
where each array-dependent coefficient in the case of a rigid sphere is [19]
b n ( k r ) = 4 π i n ( j n ( k r ) j n ( k r ) h n ( k r ) h n ( k r ) )
h n , j n , h n , and j n represent the Hankel and Bessel functions and the corresponding derivatives with respect to kr, respectively.

3.4. 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 Y ^ ( Ω ) 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
Y ^ H ( Ω ) Y H ( Ω )
Equation (17) is the response of a theoretically transformed antenna array, and after equalization, it can be incorporated into the components
Y ^ H ( Ω ) = ( B 1 Y H ( Φ ) Γ ϕ ) ( Y ( Φ ) B Y H ( Ω ) ) = W e n c h ( Ω )
where h ( Ω ) represents the array response vector to the Ω direction or the steering vector, and W e n c is the encoding matrix that implements the equalization and SHT. Ideally, B 1 Y H ( Φ ) Γ ϕ Y ( Φ ) B = I , and Equation (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]:
W ^ e n c ( k ) = Y H ( Ψ ) Γ ψ H ^ H ( k ) ( H ^ ( k ) Γ ψ H ^ H ( k ) + β ( k ) I ) 1
where Γ ψ is the correct weight for the measurement grid, and H ^ ( k ) = [ h ^ 1 , , h ^ 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 t r is accurately described by ( N t r + 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 are
h ^ ( k , Ω ) = ( Y H ( Ψ ) Γ ψ H ^ T ( k ) ) T Y H ( Ω )
Y ^ H ( k ,   Ω ) = W ^ e n c ( k ) h ^ ( k , Ω ) .

4. 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.

4.1. 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 ( k r ) , which results in
q n m ( k ) = Y H ( Ω ) s ( k ) + z ( k ) ,
where z ( k ) = B 1 Y H ( Φ ) Γ n ( k ) . The covariance matrix of the array in this scenario is expressed as
R a a ( k ) = Y H ( Ω ) R s s ( k ) Y ( Ω ) + R z z ( k )
while the covariance matrix of the noise is
R z z ( k ) = σ 2 I | | B 1 ( k r ) Y H ( Φ ) Γ | | 2
which depends on k, and | | · | | is the Frobenius norm. The corresponding time-domain representation of q n m ( k ) is obtained through the inverse Fourier transform 1 { · } as q n m ( t ) = 1 { q n m ( k ) } . Nevertheless, multiplying by B 1 results in unwanted noise amplification in some frequencies [26]. Hence, the same equalization is used in time-domain smoothing, as given in [26]:
H e q ( k ) = [ σ 2 γ n = 0 N h m = n n | | B 1 ( k r ) Y ( Φ ) Γ | | 2 ] 1 2
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 e q ( 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
R ˜ a a ( k ) = Y H ( Ω ) R ˜ s s ( k ) Y ( Ω ) + R ˜ z z ( k )
where
R ˜ z z ( k ) = | | H e q ( k ) | | 2 R z z ( k )
R ˜ s s ( k ) = R { s ( k ) s H ( k ) = | | H e q ( k ) | | 2 R s s ( k )
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 time-independent (i.e., R z z ( t ) = R ˜ z z = σ ˜ a 2 I , and σ ˜ a 2 is the equalized variance).
The covariance matrix estimate of the array for N snapshots (time instant t 1 , ,   t N ) is
R ^ a a = 1 / N i = 1 N q ˜ ( t i ) q ˜ H ( t i ) ,
q ˜ n m ( t ) = 1 { H e q ( k ) q ˜ n m ( 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 frequency-domain 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.

4.2. 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:
Y ( Ω ) = [ Y 0 0 ( Ω 1 ) , Y 1 1 ( Ω 1 ) , Y 1 0 ( Ω 1 ) , Y 1 1 ( Ω 1 ) , Y 3 3 ( Ω 1 ) Y 0 0 ( Ω 2 ) , Y 1 1 ( Ω 2 ) , Y 1 0 ( Ω 2 ) , Y 1 1 ( Ω 2 ) , Y 3 3 ( Ω 2 ) ]
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
Y ( Ω ) = [ Y 0 0 ( Ω ) ,   Y 1 1 ( Ω ) , Y 1 0 ( Ω ) , Y 1 1 ( Ω ) , ,   Y 3 3 ( Ω ) ] .
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:
Ω ^ = a r g min Ω { P ( Ω ) } .
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 ( Ω ) .

4.3. 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 as
^ y ( Ω ) = 1 K k = 1 K W ( k , Ω ) p ^ n m ( k ) .
where p ^ n m ( k ) is the noisy form of Equation (18), and the weighting is
W ( k , Ω ) = 4 π ( N h + 1 ) 2 Y ( Ω ) B 1 ( k r ) ,
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 as
^ y ( Ω ) = 1 N i = 1 N Y H ( Ω ) q ˜ n m ( t i ) ,
The output energy of the time-domain beamformer at Ω is
P P W D ( Ω ) =   | | y ^ ( Ω ) | | 2 = Y ( Ω ) R ^ a a Y H ( Ω ) ,
where R ^ a a is the estimated covariance matrix of the array defined in Equation (34).

4.4. MUSIC

Another important DoA estimation method is MUSIC [11,16]. The MUSIC spatial spectrum is computed as
P M U S I C ( Ω ) = 1 | | E n H Y H ( Ω ) | | 2
where E n is the noise matrix of the array obtained from eigenvalue decomposition of covariance matrix estimate R ^ a a as
R ^ a a = E λ E H = E s λ s E s H + E n λ n E n H
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.

4.5. 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.

4.5.1. 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
R ˜ a a = Y H ( Ω ) R ˜ s s Y ( Ω ) + σ ˜ a 2 I
For time instant t i = t 1 , ,   t N , the SML probability density function is given as
p ( A ˜ N / Ω , R ˜ s s , σ ˜ a 2 ) = i = 1 N 1 π ( N h + 1 ) 2 | R ˜ a a | exp ( q ˜ H ( t i ) R ˜ a a 1 q ˜ ( t i )
where R ˜ a a 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
l o g ( p ( A ˜ N / Ω , R ˜ s s , σ ˜ a 2 ) ) = l o g | R ˜ a a | + T r { R ˜ a a 1 R ^ a a }
The negative log-likelihood solution with respect to S, σ 2 ([31] provides the details), results in
R ˜ s s ( Ω ) = Y ( R ^ a a σ ˜ a 2 ( Ω ) I ) Y H ( Ω )
and
σ ˜ a 2 ( Ω ) = 1 ( N h + 1 ) 2 U T r { P Y ¯ ( Ω ) R ^ a a } ,
respectively, where
P Y ¯ ( Ω ) = I Y H ( Ω ) Y ( Ω )
and
Y ( Ω ) = ( Y ( Ω ) Y H ( Ω ) ) 1 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
P S M L ( Ω ) = log { | Y ( Ω ) R s s ( Ω ) Y H ( Ω ) + σ ˜ a 2 ( Ω ) I | }

4.5.2. 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
R = E { ( q ˜ ( t i ) E { q ˜ ( t i ) } ) ( q ˜ ( t j ) E { q ˜ ( t j ) } ) H } = σ ˜ a 2 I
when the average signal waveform is derived from the signals. According to [31], the DML likelihood function for different snapshots is expressed as
p ( A ˜ N Ω , S ˜ N , σ ˜ a 2 ) = i = 1 N 1 π σ ˜ a 2 I × exp ( q ˜ ( t i ) Y H ( Ω ) s ˜ ( t i ) ) H ( q ˜ ( t i ) Y H ( Ω ) s ˜ ( t i ) ) σ ˜ a 2 )
where S ˜ N = [ s ˜ ( t i ) , , 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]:
σ ^ a 2 = 1 ( N h + 1 ) 2 T r { P Y ¯ ( Ω ) R ^ a a } .
The localization function and signal estimates can be deduced using the negative log-likelihood, which gives rise to a nonlinear least-square problem [31]:
P D M L ( Ω ) = T r { P Y ¯ ( Ω ) R ^ a a }
and
S ^ N = Y ( Ω ) A ˜ N
respectively.
Asymptotically, MUSIC is equivalent to DML for large sizes of samples and uncorrelated signals [31] (i.e., in the incoherent case).

4.6. 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].
The WSF technique is a larger sample approximation of the SML technique, whose localization function is presented as [31]
P W S F ( Ω ) = T r { P Y ¯ ( Ω ) E ^ s W E ^ s H } ,   W = λ ^ 2 λ s 1
where λ ^ = ( λ s σ ^ 2 I ) , λ s denotes the signal eigenvalue matrix, and σ ^ 2 = 1 ( N h + 1 ) 2 d i a g { λ n } defines the average of the ( N h + 1 ) 2 U lowest eigenvalues. This weighting provides the least asymptotic error variance, as depicted in [31].

4.7. 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 estimate ^ Ω , i.e., E { Ω ^ } = Ω , is bound by
E { ( Ω ^ Ω ) ( Ω ^ Ω ) T } = [ E { 2 log ( p ( A ˜ N / Ω ) ) Ω Ω T } ] 1
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):
{ C R L B s t o } i j = σ ˜ a 2 2 [ { T r { Y j P Y ¯ Y i H R ˜ s s Y R ˜ a a 1 Y H R ˜ s s } } ] 1
where Y i is the partial derivative of Y with respect to the ith element in Ω and defined as
Y i = Y Ω i
Y i elements are the SH partial derivatives with respect to elevation angle θ and azimuth angle ϕ and expressed as
Y n m ( Ω ) θ = ( 2 n + 1 ) ( n m ) ! 4 π ( n + m ) ! e i m ϕ csc ( θ )                                                                                × [ ( n + 1 ) cos ( θ ) P n m ( cos ( θ ) ) + ( m n 1 ) P n + 1 m ( cos ( θ ) ) ]
and
Y n m ( Ω ) ϕ = i m Y n m ( Ω ) ,
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]:
{ C R L B d e t } i j = σ ˜ a 2 2 [ { T r { Y j P Y ¯ Y i H R ˜ s s } } ] 1

4.8. 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].

5. 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):
R M S E = E { ( Ω Ω ^ ) ( Ω Ω ^ ) T } ,
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 S N R = 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., S N R e = | | s ˜ ( k ) | | 2 / σ ˜ a 2 . In all cases considered in this paper, we simulated the noise as space-domain i.i.d. complex Gaussian random variables with variance σ 2 .

5.1. 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 50 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.

5.2. 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
R ˜ s s = k = 1 N | | H u ( k ) | | 2 E e q
and E e q 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 S N R 10   dB . The elevation angle estimation has a smaller RMSE, just as in the CRLB prediction.

5.3. 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 ˜ s s = E e q [ 1 1 1 1 ] ,   R ˜ s s = E e q [ 1 0.85 0.85 1 ] ,   R ˜ s s = E e q [ 1 0.5 0.5 1 ] ,   and   R ˜ s s = E e q [ 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 S N R = 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 almost-incoherent 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 ° .

5.4. 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 k r > N h , spatial aliasing exists in high frequencies [36]. However, WSF and ML techniques are not affected by spatial aliasing, because B ( k r ) 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
R ˜ s s ( k ) = | | H e q ( k ) | | 2 [ 1 1 1 1 ] ,  
and the covariance matrix estimate of the array becomes
R ^ a a ( k ) = q ˜ ( k ) q ˜ H ( k ) .
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.
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.

6. 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.

Author Contributions

Conceptualization, O.J.F. and T.S.; methodology, O.J.F.; software, O.J.F.; validation, O.J.F. and T.S.; formal analysis, O.J.F.; investigation, T.S.; resources, T.S.; writing—original draft preparation, O.J.F.; writing—review and editing, T.S.; supervision, T.S. All authors have read and agreed to the published version of the manuscript.

Funding

The APC was funded by the University of Johannesburg, South Africa.

Data Availability Statement

Data are available upon request from authors.

Acknowledgments

This research is supported in part by the University Research Committee (URC) of the University of Johannesburg, South Africa.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Carlin, M.; Rocca, P.; Oliveri, G.; Viani, F.; Massa, A. Directions-of-Arrival Estimation Through Bayesian Compressive Sensing Strategies. IEEE Trans. Antennas Propag. 2013, 61, 3828–3838. [Google Scholar] [CrossRef]
  2. Pralon, M.G.; Del Galdo, G.; Landmann, M.; Hein, M.A.; Thoma, R.S. Suitability of Compact Antenna Arrays for Direction-of-Arrival Estimation. IEEE Trans. Antennas Propag. 2017, 65, 7244–7256. [Google Scholar] [CrossRef]
  3. Kumar, P.; Kumar, C.; Kumar, S.; Srinivasan, V. Active spherical phased array design for satellite payload data transmission. IEEE Trans. Antennas Propagat. 2015, 63, 4783–4791. [Google Scholar] [CrossRef]
  4. Alibakhshikenari, M.; Virdee, B.S.; Azpilicueta, L.; Naser-Moghadasi, M.; Akinsolu, M.O.; See, C.H.; Liu, B.; Abd-Alhameed, R.A.; Falcone, F.; Huynen, I.; et al. A comprehensive survey of metamaterial transmission-line based antennas: Design, challenges, and applications. IEEE Access 2020, 8, 144778–144808. [Google Scholar] [CrossRef]
  5. Alibakhshikenari, M.; Babaeian, F.; Virdee, B.S.; Aïssa, S.; Azpilicueta, L.; See, C.H.; Althuwayb, A.A.; Huynen, I.; Abd-Alhameed, R.A.; Falcone, F.; et al. A comprehensive survey on “various decoupling mechanisms with focus on metamaterial and metasurface principles applicable to SAR and MIMO antenna systems”. IEEE Access 2020, 8, 192965–193004. [Google Scholar] [CrossRef]
  6. Famoriji, O.J.; Shongwe, T. Source Localization of EM Waves in the Near-Field of Spherical Antenna Array in the Presence of Unknown Mutual Coupling. Wirel. Commun. Mob. Comput. 2021, 2021, 1–14. [Google Scholar] [CrossRef]
  7. Famoriji, O.J.; Shongwe, T. Direction-of-Arrival Estimation of Electromagnetic Wave Impinging on Spherical Antenna Array in the Presence of Mutual Coupling Using a Multiple Signal Classification Method. Electronics 2021, 10, 2651. [Google Scholar] [CrossRef]
  8. Shirkolaei, M.M.; Ghalibafan, J. Scannable Leaky-Wave Antenna Based on Ferrite-Blade Waveguide Operated Below the Cutoff Frequency. IEEE Trans. Magn. 2021, 57, 1–10. [Google Scholar] [CrossRef]
  9. Nadeem, I.; Alibakhshikenari, M.; Babaeian, F.; Althuwayb, A.; Virdee, B.S.; Azpilicueta, L.; Khan, S.; Huynen, I.; Falcone, F.J.; Denidni, T.A.; et al. Circular polarized antennas for existing and emerging wireless communication technologies. J. Phys. D Appl. Phys. 2021, 55, 1–15. [Google Scholar]
  10. Kumar, C.; Kumar, B.P.; Kumar, V.S.; Srinivasan, V.V. Dual Circularly Polarized Spherical Phased-Array Antenna for Spacecraft Application. IEEE Trans. Antennas Propag. 2012, 61, 598–605. [Google Scholar] [CrossRef]
  11. Tervo, S.; Politis, A. Direction of Arrival Estimation of Reflections from Room Impulse Responses Using a Spherical Microphone Array. IEEE/ACM Trans. Audio Speech Lang. Process. 2015, 23, 1539–1551. [Google Scholar] [CrossRef]
  12. Pastorino, M.; Randazzo, A. A smart antenna system for direction of arrival estimation based on a support vector regression. IEEE Trans. Antennas Propag. 2005, 53, 2161–2168. [Google Scholar] [CrossRef]
  13. Schmidt, R. Multiple emitter location and signal parameter estimation. IEEE Trans. Antennas Propag. 1986, 34, 276–280. [Google Scholar] [CrossRef] [Green Version]
  14. Swindlehurst, A.; Kailath, T. A performance analysis of subspace based methods in the presence of model errors. I. The MUSIC algorithm. IEEE Trans. Signal Process. 1992, 40, 1578–1774. [Google Scholar]
  15. Zoltowski, M.D.; Haardt, M.; Mathews, C.P. Closed-form 2-D angle estimation with rectangular arrays in element space or beam space via unitary ESPRIT. IEEE Trans. Signal Process. 1996, 44, 316–328. [Google Scholar] [CrossRef] [Green Version]
  16. Gao, F.; Gershman, B. A generalized ESPRIT approach to direction-of-arrival estimato. IEEE Trans. Signal Process. 2005, 12, 254–257. [Google Scholar] [CrossRef]
  17. Ziskind, I.; Wax, M. Maximum likelihood localization of multiple sources by alternating projection. IEEE Trans. Acoust. Speech Signal Process. 1988, 36, 1553–1560. [Google Scholar] [CrossRef]
  18. Stoica, P.; Gershman, A. Maximum-likelihood DOA estimation by data-supported grid search. IEEE Signal Process. Lett. 1999, 6, 273–275. [Google Scholar] [CrossRef]
  19. Wang, F.; Fang, J.; Duan, H.; Li, H. Phased-Array-Based Sub-Nyquist Sampling for Joint Wideband Spectrum Sensing and Direction-of-Arrival Estimation. IEEE Trans. Signal Process. 2018, 66, 6110–6123. [Google Scholar] [CrossRef] [Green Version]
  20. Rao, N.N. Fundamentals of Electromagnetics for Electrical and Computer Engineering; Prentice Hall: Hoboken, NJ, USA, 2009. [Google Scholar]
  21. Hoffman, M. Conventions for the analysis of spherical arrays. IRE Trans. Antennas Propag. 1963, 11, 390–393. [Google Scholar] [CrossRef]
  22. Ottersten, B.; Kailath, T. Direction-of-arrival estimation for wide-band signals using the ESPRIT algorithm. IEEE Trans. Acoust. Speech Signal Process. 1990, 38, 317–327. [Google Scholar] [CrossRef]
  23. Wang, H.; Kaveh, M. Coherent signal-subspace processing for the detection and estimation of angles of arrival of multiple wide-band sources. IEEE Trans. Acoust. Speech Signal Process. 1985, 33, 823–831. [Google Scholar] [CrossRef] [Green Version]
  24. Zhang, Y.; Ye, Z. Efficient Method of DOA Estimation for Uncorrelated and Coherent Signals. IEEE Antennas Wirel. Propag. Lett. 2008, 7, 799–802. [Google Scholar] [CrossRef]
  25. Rafaely, B. Spherical array configurations. In Fundamentals of Spherical Array Processing; Springer Topics in Signal Processing; Springer: Cham, Switzerland, 2019; Volume 16, pp. 81–102. [Google Scholar]
  26. Huleihel, N.; Rafaely, B. Spherical array processing for acoustic analysis using room impulse responses and time-domain smoothing. J. Acoust. Soc. Am. 2013, 133, 3995–4007. [Google Scholar] [CrossRef]
  27. Hafizovic, I.; Nilsen, C.-I.C.; Holm, S. Transformation Between Uniform Linear and Spherical Microphone Arrays with Symmetric Responses. IEEE Trans. Audio Speech Lang. Process. 2011, 20, 1189–1195. [Google Scholar] [CrossRef]
  28. Kellermann, W.; Kowalczyk, K.; Mabande, E.; Sun, H. Comparison of Subspace-Based and Steered Beamformer-Based Reflection Localization Methods. In Proceedings of the 2011 19th European Signal Processing Conference, Barcelona, Spain, 29 August–2 September 2011. [Google Scholar] [CrossRef]
  29. Goossens, R.; Bogaert, I.; Rogier, H. Phase-Mode Processing for Spherical Antenna Arrays With a Finite Number of Antenna Elements and Including Mutual Coupling. IEEE Trans. Antennas Propag. 2009, 57, 3783–3790. [Google Scholar] [CrossRef]
  30. Krim, H.; Viberg, M. Two decades of array signal processing research: The parametric approach. IEEE Signal Process. Mag. 1996, 13, 67–94. [Google Scholar] [CrossRef]
  31. Ottersten, B.; Viberg, M.; Stoica, P.; Nehorai, A. Exact and Large Sample Maximum Likelihood Techniques for Parameter Estimation and Detection in Array Processing. In Radar Array Processing; Springer: Berlin/Heidelberg, Germany, 1993; pp. 99–151. [Google Scholar] [CrossRef]
  32. Sengupta, D.L.; Smith, T.; Larson, R. Radiation characteristics of a spherical array of circularly polarized elements. IRE Trans. Antennas Propag. 1968, 16, 2–7. [Google Scholar] [CrossRef]
  33. Bondyopadhyay, P.K. Geodesic sphere phased arrays for LEO satellite communications. In Proceedings of the IEEE Antennas and Propagation Society International Symposium. Transmitting Waves of Progress to the Next Millennium. 2000 Digest. Held in Conjunction with: USNC/URSI National Radio Science Meeting, Salt Lake City, UT, USA, 16–21 July 2000; pp. 206–209. [Google Scholar]
  34. Famoriji, O.J.; Akingbade, K.; Ogunti, E.; Apena, W.; Fadamiro, A.; Lin, F. Analysis of phased array antenna system via spherical harmonics decomposition. IET Commun. 2019, 13, 3097–3104. [Google Scholar] [CrossRef]
  35. Stockton, R.; Hockensmith, R. Application of spherical arrays—A simple approach. Proc. IEEE AP-S Int. Symp. 1977, 15, 202–205. [Google Scholar]
  36. Zhu, J.; Han, L.; Blum, R.S.; Xu, Z. Mutiple-snapshot newtonalized orthogonal matching pursuit for line spectrum estimation with multiple measurement vector. Signal Process. 2019, 65, 155–164. [Google Scholar] [CrossRef]
  37. Kumar, L.; Hegde, R.M. Stochastic Cramér-Rao Bound Analysis for DOA Estimation in Spherical Harmonics Domain. IEEE Signal Process. Lett. 2014, 22, 1030–1034. [Google Scholar] [CrossRef]
  38. Pastorino, M.; Randazzo, A. The SVM-Based Smart Antenna for Estimation of the Directions of Arrival of Electromagnetic Waves. IEEE Trans. Instrum. Meas. 2006, 55, 1918–1925. [Google Scholar] [CrossRef]
  39. El Zooghby, A.; Christodoulou, C.; Georgiopoulos, M. Performance of radial-basis function networks for direction of arrival estimation with antenna arrays. IEEE Trans. Antennas Propag. 1997, 45, 1611–1617. [Google Scholar] [CrossRef]
  40. Famoriji, O.J.; Ogundepo, O.Y.; Qi, X. An Intelligent Deep Learning-Based Direction-of-Arrival Estimation Scheme Using Spherical Antenna Array with Unknown Mutual Coupling. IEEE Access 2020, 8, 179259–179271. [Google Scholar] [CrossRef]
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.
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.
Electronics 11 00208 g001
Figure 2. SAA with 64 elements as modeled with CST: (a) side view and (b) top view [7].
Figure 2. SAA with 64 elements as modeled with CST: (a) side view and (b) top view [7].
Electronics 11 00208 g002aElectronics 11 00208 g002b
Figure 3. Root-mean-square error over CRLB and 100 Monte Carlo samples: (a) elevation and (b) azimuth angle versus noise variance σ 2 for a single signal in ϕ = θ = 60 ° . The azimuth angle has a larger CRLB and RMSE than the elevation angle. (a) θ ^ , elevation (b) ϕ ^ , azimuth.
Figure 3. Root-mean-square error over CRLB and 100 Monte Carlo samples: (a) elevation and (b) azimuth angle versus noise variance σ 2 for a single signal in ϕ = θ = 60 ° . The azimuth angle has a larger CRLB and RMSE than the elevation angle. (a) θ ^ , elevation (b) ϕ ^ , azimuth.
Electronics 11 00208 g003
Figure 4. Monte Carlo samples against ϕ 2 when S N R = 20   dB and the second signal is delayed by d samples. (a)   ϕ ^ 1 ,   d = 0 (b)   ϕ ^ 2 ,   d = 0 (c) ϕ ^ 1   d = 0.5 (d)   ϕ ^ 2 ,   d = 0.5 (e)   ϕ ^ 1 ,   d = 1 (f) ϕ ^ 2 , d = 1 (g)   ϕ ^ 1 ,   d = 8 (h) ϕ ^ 2 ,   d = 8 .
Figure 4. Monte Carlo samples against ϕ 2 when S N R = 20   dB and the second signal is delayed by d samples. (a)   ϕ ^ 1 ,   d = 0 (b)   ϕ ^ 2 ,   d = 0 (c) ϕ ^ 1   d = 0.5 (d)   ϕ ^ 2 ,   d = 0.5 (e)   ϕ ^ 1 ,   d = 1 (f) ϕ ^ 2 , d = 1 (g)   ϕ ^ 1 ,   d = 8 (h) ϕ ^ 2 ,   d = 8 .
Electronics 11 00208 g004
Figure 5. (ah) RMSE over CRLB and 100 Monte Carlo samples when SNR = 20 dB for chosen frequencies: (a)   ϕ ^ 1 ,   f = 900   MHz (b)   ϕ ^ 2 ,   f = 900   MHz (c)   ϕ ^ 1 ,   f = 2   GHz (d)   ϕ ^ 2 ,   f = 2   GHz (e)   ϕ ^ 1 ,   f = 4   GHz (f)   ϕ ^ 2 ,   f = 4   GHz (g)   ϕ ^ 1 ,   f = 8   GHz (h)   ϕ ^ 2 ,   f = 8   GHz . (ip) RMSE over CRLB and 100 Monte Carlo samples when SNR = 20 dB for chosen frequencies: (i)   ϕ ^ 1 ,   f = 16   GHz (j)   ϕ ^ 2 ,   f = 16   GHz (k)   ϕ ^ 1 ,   f = 32   GHz (l)   ϕ ^ 2 ,   f = 32   GHz (m)   ϕ ^ 1 ,   f = 40   GHz (n)   ϕ ^ 2 ,   f = 40   GHz (o)   ϕ ^ 1 ,   f = 75   GHz (p)   ϕ ^ 2 ,   f = 75   GHz .
Figure 5. (ah) RMSE over CRLB and 100 Monte Carlo samples when SNR = 20 dB for chosen frequencies: (a)   ϕ ^ 1 ,   f = 900   MHz (b)   ϕ ^ 2 ,   f = 900   MHz (c)   ϕ ^ 1 ,   f = 2   GHz (d)   ϕ ^ 2 ,   f = 2   GHz (e)   ϕ ^ 1 ,   f = 4   GHz (f)   ϕ ^ 2 ,   f = 4   GHz (g)   ϕ ^ 1 ,   f = 8   GHz (h)   ϕ ^ 2 ,   f = 8   GHz . (ip) RMSE over CRLB and 100 Monte Carlo samples when SNR = 20 dB for chosen frequencies: (i)   ϕ ^ 1 ,   f = 16   GHz (j)   ϕ ^ 2 ,   f = 16   GHz (k)   ϕ ^ 1 ,   f = 32   GHz (l)   ϕ ^ 2 ,   f = 32   GHz (m)   ϕ ^ 1 ,   f = 40   GHz (n)   ϕ ^ 2 ,   f = 40   GHz (o)   ϕ ^ 1 ,   f = 75   GHz (p)   ϕ ^ 2 ,   f = 75   GHz .
Electronics 11 00208 g005aElectronics 11 00208 g005bElectronics 11 00208 g005c
Figure 6. RMSE versus SNR comparison test between MUSIC, PWD, DML, SML, and WSF.
Figure 6. RMSE versus SNR comparison test between MUSIC, PWD, DML, SML, and WSF.
Electronics 11 00208 g006
Figure 7. GE performance evaluation of MUSIC, PWD, DML, SML, and WSF against SNR using the measured data in [40].
Figure 7. GE performance evaluation of MUSIC, PWD, DML, SML, and WSF against SNR using the measured data in [40].
Electronics 11 00208 g007
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Famoriji, O.J.; Shongwe, T. Critical Review of Basic Methods on DoA Estimation of EM Waves Impinging a Spherical Antenna Array. Electronics 2022, 11, 208. https://doi.org/10.3390/electronics11020208

AMA Style

Famoriji OJ, Shongwe T. Critical Review of Basic Methods on DoA Estimation of EM Waves Impinging a Spherical Antenna Array. Electronics. 2022; 11(2):208. https://doi.org/10.3390/electronics11020208

Chicago/Turabian Style

Famoriji, Oluwole John, and Thokozani Shongwe. 2022. "Critical Review of Basic Methods on DoA Estimation of EM Waves Impinging a Spherical Antenna Array" Electronics 11, no. 2: 208. https://doi.org/10.3390/electronics11020208

APA Style

Famoriji, O. J., & Shongwe, T. (2022). Critical Review of Basic Methods on DoA Estimation of EM Waves Impinging a Spherical Antenna Array. Electronics, 11(2), 208. https://doi.org/10.3390/electronics11020208

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop