Direction-of-Arrival Estimation of Electromagnetic Wave Impinging on Spherical Antenna Array in the Presence of Mutual Coupling Using a Multiple Signal Classiﬁcation Method

: A spherical antenna array (SAA) is the conﬁguration of choice in obtaining an antenna array with isotropic characteristics. An SAA has the capacity to receive an electromagnetic wave (EM) with equal intensity irrespective of the direction-of-arrival (DoA) and polarization. Therefore, the DoA estimation of electromagnetic (EM) waves impinging on an SAA with unknown mutual coupling needs to be considered. In the spherical domain, the traditional multiple signal classiﬁcation algorithm (SH-MUSIC) is faced with a computational complexity problem. This paper presents a one-dimensional MUSIC method (1D-MUSIC) for the estimation of the azimuth and elevation angles. An intermediate mapping matrix that exists between Fourier series and the spherical harmonic function is designed, and the Fourier series Vandermonde structure is used for the realization of the polynomial rooting technique. This mapping matrix can be computed prior to the DoA estimation, and it is only a function of the array conﬁguration. Based on the mapping matrix, the 2-D angle search is transformed into two 1-D angle ﬁndings. Employing the features of the Fourier series, two root polynomials are designed for the estimation of the elevation and azimuth angles, spontaneously. The developed method avoids the 2-D spectral search, and angles are paired in automation. Both numerical simulation results, and results from experimental measured data (i.e., with mutual coupling effect incorporated), show the validity, potency, and potential practical application of the developed algorithm.


Motivation and Incitation
Direction-of-arrival (DoA) is a topic of interest in electromagnetic (EM)-related fields in both academia and industry. It finds application in mobile communication systems, vehicular technology, radar, sonar, and military surveillance [1,2]. Knowing the DoA of an incoming electromagnetic wave impinging on a receiving antenna is necessary to effectively localize the positions of the corresponding sources. this enhances the adaptive beamforming of the incoming pattern to boost the system sensitivity towards the wanted signal's direction and to mitigate undesirable interferences. Consequently, the antenna array provides maximum beam to the desired users and nulls in the interference directions to improve the mobile and base station's performance. Hence, DoA estimation in antenna array processing becomes important.
A spherical antenna array (SAA) is important whenever an array with isotropic characteristics is required. A spherical phased antenna array is widely used in spacecraft and satellite communications. The electronic beam scanning from a phased antenna array provides a hemispherical scan and an almost uniformly distributed gain [2,3]. An SAA A spherical antenna array (SAA) is important whenever an array with isotropic characteristics is required. A spherical phased antenna array is widely used in spacecraft and satellite communications. The electronic beam scanning from a phased antenna array provides a hemispherical scan and an almost uniformly distributed gain [2,3]. An SAA configuration (as in Figure 1) has three-dimensional symmetry, which is an advantage in spatial signal analysis. Its size and geometry give the directivity and level of beam scanning in the azimuth and elevation angles required from the array.  [3].
A linear array that scans in one axis, and a planar array that scans in the x-and yaxes, are acceptable for a bounded scan demand of 60-70° [3]. The projected apertures of these arrays drop by the cosine factor of the scan angle. Hence, the array directivity is reduced for a big scan angle. For larger steerability, a multifold planar configuration, or in a general form, an SAA is preferred [3][4][5]. The directivity of an SAA remains almost the same when an active source is moved over the surface of the SAA.

Related Works, Associated Problems, and Contributions
Researchers have developed different methods of estimating DoA in a spherical harmonic domain. The subspace algorithms have been formulated in a spherical domain [6]. The estimation of signal parameter via rotational invariance technique (ESPRIT) in spherical harmonic (SH) domain (SH-ESPRIT) has been developed based on a tangent function via an SH recursion relationship. However, the SH-ESPRIT is ill-conditioned if the tangent function tends to infinity close to the equator of the spherical system [7]. The two-stage SH-ESPRIT (TS-SH-ESPRIT) uses two recurrence relations of SH functions to compute the azimuth and elevation separately [8]. A two-stage decoupled method (TSDM) that uses unitary SH-ESPRIT in the initial stage and unitary SH root-MUSIC in the second stage has been developed in [9]. The DoA vector SH-ESPRIT (DoAV-SH-ESPRIT) was employed in order to get the directional variables from three separate recurrence relations of SH functions. These three variables are corresponding to the unconventional parts (three) of the A linear array that scans in one axis, and a planar array that scans in the x-and y-axes, are acceptable for a bounded scan demand of 60-70 • [3]. The projected apertures of these arrays drop by the cosine factor of the scan angle. Hence, the array directivity is reduced for a big scan angle. For larger steerability, a multifold planar configuration, or in a general form, an SAA is preferred [3][4][5]. The directivity of an SAA remains almost the same when an active source is moved over the surface of the SAA.

Related Works, Associated Problems, and Contributions
Researchers have developed different methods of estimating DoA in a spherical harmonic domain. The subspace algorithms have been formulated in a spherical domain [6]. The estimation of signal parameter via rotational invariance technique (ESPRIT) in spherical harmonic (SH) domain (SH-ESPRIT) has been developed based on a tangent function via an SH recursion relationship. However, the SH-ESPRIT is ill-conditioned if the tangent function tends to infinity close to the equator of the spherical system [7]. The two-stage SH-ESPRIT (TS-SH-ESPRIT) uses two recurrence relations of SH functions to compute the azimuth and elevation separately [8]. A two-stage decoupled method (TSDM) that uses unitary SH-ESPRIT in the initial stage and unitary SH root-MUSIC in the second stage has been developed in [9]. The DoA vector SH-ESPRIT (DoA V -SH-ESPRIT) was employed in order to get the directional variables from three separate recurrence relations of SH functions. These three variables are corresponding to the unconventional parts (three) of the unit vector that are directional and point to the DoA [10,11]. DoA V -SH-ESPRIT does not have angular complexity and the accuracy of estimation is independent of the DoA. In [12], the SH-MUSIC (SH-MUSIC) approach is presented by Li et al., but it requires very high computational load due to two-dimensional angle searching [13]. In addition, Kumar et al. developed SH root-MUSIC (SH-R-MUSIC) [14]. The SH-R-MUSIC approach circumvents peak searches; however, it only estimates azimuth via polynomial rooting.
Furthermore, a few scientists attempted the manifold separation method (MSM) for the expansion of the traditional root multiple signal classification towards the estimation of azimuths for arrays with arbitrary configurations, giving birth to different search-free techniques [15][16][17][18]. Joint estimations of elevation and azimuth for arbitrary arrays are transformed into bivariate polynomial rooting via MSM [19]. However, the higher order of the polynomial and the larger aperture causes complexity in the computation for the rooting operation. The fast Fourier transform of multiple signal classification (FFT-MUSIC) transformed the two-dimensional MUSIC cost function to the customary 2D-DFT by MSM, and later used the 2D-FFT for the efficient computation of the spatial spectrum, which mitigates the complexity in the computation of root-MUSIC or 2D MUSIC [19].
With the current trend in the growth of technology and applications, systems are becoming smaller, leading to smaller inter-element distances in an array, leading to strongly coupled arrays. This causes strong mutual coupling, impedance mismatch, and poor radiation characteristics. Addressing DoA in the presence of mutual coupling remains an important topic of interest.
This paper presents a 1D MUSIC (1D-MUSIC) method for the estimation of azimuth and elevation of EM waves impinging on the SAA. The key innovations and contributions of this study are as follow: A new method is proposed for computing a mapping matrix that is responsible for building a link between Fourier series and SH functions. Furthermore, if we consider the Fourier series structure that is associated with the Vandermonde feature, it becomes effortless linking the polynomial roots to data containing the DoA. The elevation is first estimated via a one-dimensional search considering the elevation estimates as initialization, and then two polynomials are designed for the calculation of both the elevation and azimuth accordingly. Furthermore, it is important to state that Huang and Cheng [20], first looked into and developed 1D-MUSIC. The major novelty of our work as compared with [20] is that we particularly extended the 1D-MUSIC method to the antenna array as opposed to 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 gap filled in our work.

System Model
Considering an SAA with radius r having V isotropic and mutually independent elements with the assumption that no mutual coupling exists among the elements. For an SAA configuration, with the v-th element localization expressed as r[sinΨ v cosψ v , sinΨ v sinψ v , sinΨ v ] T in spherical domain, Ψ v and ψ v represent the elevation and azimuth of v-th element, respectively.
[·] T is transpose operation. Assuming Q number of EM waves arrive on SAA from far field. Let Φ d = (ϑ d , ϕ d ) for d = 1, . . . , D, where ϑ and ϕ represent azimuth and elevation angles of arriving signal. Then the SAA output is expressed in form of matrix as where G(Φ) = [g(Φ 1 ), . . . , g(Φ D )] represents the steering matrix, which is the systemfunction from origin to the SAA, and it carries the data of the DoA. t is index of snapshot.
T is the noise that is Gaussian white in vector form with zero mean and σ 2 n variance. After SH decomposition, the output of SAA is expressed as [21] Electronics 2021, 10, 2651 in far field mode and A = (N + 1) 2 . N represents the highest order of SAA and diag {·} is a diagonal matrix and b n (kr) is described as where j n (kr) and h n (kr) represent the nth order spherical Bessel function and spherical Hankel function of first kind, respectively. The j n (kr), h n (kr) are the defined corresponding derivatives. In open sphere, elements are positioned in the open with no capturing impairment, while in rigid sphere, the elements are placed in rigid sphere. Here, the received EM wave gets scattered by the interface with the rigid sphere [22].
[·] H is the Hermitian transpose operator. k = 2π f /c denotes the signal wave number, f is the operation frequency, and c is the speed of EM propagation. and where Y m n (·) denotes the SH function with m, n representing degree and order [22].
where n ∈ [0, N], m ∈ [−n, n] and P m n (·) is the Legendre polynomial. Furthermore, applying the SFT (spherical Fourier transform) on the output with the consideration of orthogonality relation Y H (Ω)αY(Ω) = I A [23] and α = diag{α 1 , . . . , α V }, where α v denotes the sampling weight of the sampling scheme. From the left hand side of Equation (2), multiply it by Y H (Ω)α, then the model can be obtained in SH domain aŝ Multiplying Equation (6) from the left with B −1 (kr) to ensure frequency independent steering matrix, the system model can be expressed as where The covariance matrix is expressed as Nonetheless, the covariance matrix R c is unseen in practical sense, an approach is to employ sample covariance matrixR c in place of the R cR where T denotes the number of snapshots and X = [ On this note, the output of the SAA has been transformed from spatial to SH domain. The steering vector is only made up of the SH function. If we use the properties of SH function, a mapping matrix can easily be constructed to design a connection between Fourier series and SH function. The SH matrix is decomposed into the multiplication of two matrixes, consisting of the azimuth and elevation information, respectively.
Considering the defined SH function in Equation (5), the SH vector y(ϑ, ϕ) can be described as M R A×A gives a diagonal matrix which consists the coefficient of normalization as If L e ∈ C A×(2N+1) and L a ∈ R A×(2N+1) , then L a matrix will satisfy d(ϕ) = L a f(ϕ) (15) f(ϕ) = e jN ϕ , . . . , 1, . . . , e −jN ϕ T and L e is expressed as q(ϑ) = L e f(ϑ) (16) where f(ϑ) = [ e −jNϑ , . . . , 1, . . . , e jNϑ T [24]. We can establish a relationship between Fourier series and SH function as The matrix is expressed as where I = I A 2 :, 1 : A : A 2 ∈ R A 2 ×A , this implies I can be modeled via extraction from columns I A 2 and the interval of extraction is A. is Kronecker product operator. is the mapping matrix that establish the connection between SH vector and Fourier basis function (2D). The reason for calculation is to determine L a and L e because the L a and L e are not known in Equation (18). Hence, the next event is to the best method for L a and L e representation. Based on the features of d(ϑ) and f(ϑ) vectors, it becomes easier to verify that the L a matrix is of the form where Q z = 0 (2z−1)×A 2 I 2z−1 0 (2z−1)×A 2 T , z = 1, . . . , N + 1 and A 2 = N − z + 1. 0 p×p is zero matrix of p × p size. L e is calculated using two independent recurrence Legendre function relationships that is complex and requires initializing the Legendre function [24]. To ensure accuracy of L e , it is important to consider an easier method to calculate L e . L e is described as the result or solution to Least Square problem arg min 20) || · || F is the Frobenius norm [25]. The optimization problem of Equation (20) is solvable if we choose the Q direction ϑ 1 , . . . , ϑ Q distributed uniformly over (0, π]. The A × Q matrix P = [q(ϑ 1 ), . . . , q(ϑ Q )] can be designed for various angles ϑ 1 , . . . , ϑ Q . In a similar way, A 1 × Q matrix can be formed P = [f(ϑ 1 ), . . . , f(ϑ Q )] and A 1 = (2N + 1). Then L e matrix is expressed below The angle size satisfies Q ≥ 2N + 1. Therefore, inversion of Equation (21) exists. Expressing Equation (18) as a function of Equations (19) and (21), the mapping matrix is obtainable. We can obtain it offline before the estimation of DoA.
Following the derivation of , Equation (7) can be described as The information of DoA can be decomposed into matrixes with Vandermonde structure. Therefore, the root polynomial methods can be developed to compute the azimuth and elevation.

Algorithm for Localization
In this section, the description of developed localization, azimuth, and elevation algorithms is presented. Traditional SH-MUSIC works on noise subspace. If the noise and signal subspace orthogonal properties are explored, the SH-MUSIC algorithm will search the array manifold vector (which is continuous) over ϑ and ϕ for the computation of the cost function D minima as below Based on the above, the traditional SH-MUSIC method accepts the 2-D spectral search with high computational complexity. The reason for proposing 1D-MUSIC is to ensure reduced complexity in the computation of SH-MUSIC. f(ϑ) and f(ϕ) share the features of Vandermonde structures, naturally the information of DoA is associated with polynomial roots. However, the traditional root polynomial techniques in the SH domain only estimates azimuth and elevation via the rooting of the polynomial. Using polynomial rooting for the estimation of both elevation and azimuth, the elevation is first estimated and used as initialization, then we design two root polynomials for the computation of both the azimuth and elevation. Furthermore, the estimation of the azimuth firstly, and using the azimuth estimated as an initialization is equally possible in 1D-MUSIC.

Initialization
The function h(ϑ, ϕ) in Equation (23) is defined to represent cost function as Expressing Equation (17) as a function of Equation (24), h(ϑ, ϕ) is described in another form as where H(ϑ) = I A 1 f(ϑ) H H U n U H n I A 1 f(ϑ) . It becomes simple to say that Equation (25) is a problem of optimization, which can be re-expressed as The elevation can be computed as follows Elevation estimation is achievable via the one-dimensional spectral search.

Estimation of Azimuth
Consider the estimated elevation as the initialization for the computation of azimuth.
where c u can be obtained via mathematics. The polynomial carries 4 N roots. Thus, if l d is roots, then, (l * d ) inversion is automatically roots. Therefore, the roots within and outside the circle are 2 N each. Thel d root closest and within the unit circle has the information of the DoA; hence, azimuth is computed bŷ where ς(·) represents the imaginary part of (·). As we estimateφ d using information from H θ d , bothθ d andφ d become paired.

Estimation of Elevation
We can use azimuth estimates as a pilot to estimate the elevation following the same method that the azimuth was estimated. Equation (25) is then described in a new fashion as , H(φ d ) becomes known after calculating the azimuthφ d . In the same way that the azimuth was estimated, the elevation will also be estimated. Hence, elevation estimation can be achieved usinĝ In conclusion, the elevation can be estimated based on the azimuth and the azimuth can be estimated based on the elevation. It is therefore easy to compute the azimuth and elevation iteratively. A step-by-step 1D-MUSIC procedure is presented in Table 1. The results obtained via simulations indicate the accuracy of the proposed 1D-MUSIC and the estimation does not vary following two consecutive iterations. (1) Calculate L a and L e with Equations (19) and (21).
(2) Put L e , L a into Equation (18) to calculate .

Processes of Initialization
(1) Calculate covariance matrix using Equation (9) to obtain noise subspace and signal subspace.

Evaluation
This section presents the computational complexity evaluation and analysis of 1D-MUSIC in comparison with FFT-MUSIC and SH-MUSIC. Multiplication operation is the major computational factor in the estimators that causes complexity. The azimuth is estimated at first and then used as initialization. In a one-dimensional spectrum search, the time taken to initialize with the azimuth is double that when initializing with the elevation. Hence, elevation is considered for initialization.
In 1D-MUSIC, the computational complexity majorly occurs when elevation is used as initialization. Employing the fast subspace decomposition in [25], the complexity when finding D-dimensional noise subspace U n is o A 2 D . Computing H(ϑ) and corresponding inverse where N e is the value of spectral points searched along the path of elevation. 1D-MUSIC computes both the azimuth and elevation using the two root polynomials when initialization ends. The nature of the complexity involved while solving the roots of the polynomial will not be more than o A 3 ; hence, the complexity of where N a is the spectral points search across the azimuth. In FFT-MUSIC, the coefficient matrix computational complexity and two-dimensional fast Fourier transform is o V 2 D + M 2 t V 2 + N a N e + N a N e log(N a N e ) [26], for M t ≈ 2.5 kR, where R denotes the biggest range between the elements of the SAA. An overview of the comparison between SH-MUSIC, FFT-MUSIC, and the new 1D-MUSIC is as given in Table 2. For N a = 4096, N e = 256, 512, 1024, and 2048; search step in path of the elevation is 0.6 0 , 0.34 0 , 0.16 0 , and 0.07 0 . EM signal considered is 2. For SAA, M t = 16. Table 3 shows the operation of multiplication. Therefore, it is right to say that the proposed method exhibits a lower complexity in computation compared with FFT-MUSIC and SH-MUSIC.

Method
Nature of Complexity

Numerical Experiments, Results, and Discussion
This section presents the results of simulations for the purpose of illustration and performance evaluation of the proposed 1D-MUSIC method in comparison with the MUSIC, TSDM, ESPRIT, and SH-MUSIC methods. An SAA of 1.8 cm radius r, and 32 elements operating at 8.4 GHz, which is uniformly distributed on a rigid sphere, was considered for the simulations. The range between consecutive samples within the uniform sampling framework was constant. The uniform sampling framework causes reduced platonic solids. This only happens for a specific number of elements [23].The SAA has the highest order of N = 4. We conducted simulations with 700 independent Monte Carlo events utilizing MATLAB software, 2018b edition. A 0.1 degree search step was used for SH-MUSIC, MUSIC, and the proposed 1D-MUSIC. The number of the employed iterations was two. The kr was to be less than N to avoid an aliasing problem. In addition, narrowband AM signals (amplitude modulation) in the far field were considered in all the simulated cases.
Case 1: Here, we conducted a study on the fitting error (Equation (20)). We calculated the matrix L e using Equation (21) for a given value of Q. Assuming the elevation is from 0 to 180 degree to verify the validity of Equation (20). The equation for the error is as given below where ϑ ∈ (0, π] [27]. Equation (21) is solvable when Q ≥ 2N + 1. We tested the impacts of various Q on the error. Q is divided into two sections using 9 as a benchmark. Hence, Q is put at 7, 9, 10, and 180, the consequent results of the simulation are depicted in Figure 2. At Q above 9, the estimated error falls below 10 −8 ; this implies that L e matrix is seen as independent of the elevation. If Q is below 9, the error in fitting is big, which renders Equation (21) invalid due to the fact that the inverse of Equation (21) does not always exist. The curve slightly varies with elevation due to the Legendre function that takes cos ϑ as the input. For other cases considered in this paper, Q is set at 180. Case 2: In this case, for various initialization, we evaluated the performance of both the elevation and azimuth. The estimation accuracy of DoA was conducted using RMSEs. Two uncorrelated waves arriving on SAA from (45, 34) and (76, 56) degree. The SNR ranges were from −2 dB to 16 dB, with 200 snapshots. Equation (27) and TS-SH-ESPRIT were utilized for the estimation of elevation to end the initialization process. As shown in Figure 3a, based on a comparison with TS-SH-ESPRIT, higher estimation accuracy of elevation is observed according to Equation (27). Based on the computed elevation, the azimuth is obtained using SH-MUSIC. As depicted in Figure 3b, this indicates a higher accuracy of azimuth estimation based on the initialization via Equation (27). Furthermore, Equation (31) is employed to estimate elevation as a function of the estimated azimuth that is known to be an iteration for the elevation.  Figure 3. RMSEs of (a) elevation estimate; (b) azimuth estimate for various initialization with two sources. Figure 3. RMSEs of (a) elevation estimate; (b) azimuth estimate for various initialization with two sources.
In conclusion, the elevation accuracy of the estimation could be enhanced further with an iteration, which shows that the proposed algorithm runs iteratively. Case 3: This example studies the impact of the number of iterations on the accuracy of estimation of both the azimuth and elevation. We considered 400 snapshots and 3 iterations. Two uncorrelated waves arrived on the SAA from (45, 34) and (76, 56) degrees. The SNR ranged from −2 dB to 16 dB. Based on the results obtained in Figure 4b, there is an improvement in azimuth accuracy of estimation after the second iteration. The impact of iteration on the elevation is depicted in Figure 3a. Since an initialization for elevation has been taken with Equation (31) used for the estimation of elevation that is taken as an iteration, then the elevation estimated is used for initialization and for commencement of the other iteration. The RMSEs result after two iterations is as shown in Figure 4a. It reveals the accuracy level of 1D-MUSIC remains constant after the second iteration. Case 5: This scenario presents the performance at 10 dB SNR with 2 uncorrelated EM waves impinging SAA from (45, 34) and (34, 66) degrees sources. Snapshots range between 100 and 600. The RMSEs decay as we increase the snapshots as indicated in Figure 6, together with CRB (Cramer-Rao bound). This demonstrates how good the proposed 1D-MUSIC performance is.
Case 6: This section verifies the performance of the proposed method by comparing the 1D-MUSIC, MUSIC, SH-MUSIC, TSDM, and ESPRIT algorithms to CRB using the RMSEs. We computed the RMSEs of the methods against the SNR with 100 snapshots. Two uncorrelated signals were impinging on the SAA from (45, 34) and (34, 66) degrees. Case 5: This scenario presents the performance at 10 dB SNR with 2 uncorrelated EM waves impinging SAA from (45, 34) and (34, 66) degrees sources. Snapshots range between 100 and 600. The RMSEs decay as we increase the snapshots as indicated in Figure 6, together with CRB (Cramer-Rao bound). This demonstrates how good the proposed 1D-MUSIC performance is.   Furthermore, the 1D-MUSIC method shares close performance with MUSIC and SH-MUSIC, but reduced complexity in computation. This can be attributed to the fact that both MUSIC and SH-MUSIC require a 2D angle search, but 1D-MUSIC requires a 1D angle search for the computation of the elevation and uses the estimated elevation for initialization. The application of root polynomials to compute the elevation and azimuth reduces the complexity in the computation of SH-MUSIC. Zhuang et al. [19] reports that FFT-MUSIC shares close performance with the traditional MUSIC; this result exhibits identical behavior with the result reported in this paper; however, the proposed algorithm exhibits reduced complexity in computation than the FFT-MUSIC.
Case 7: In addition, using two EM signals the performance of the proposed 1D-MUSIC method is evaluated. The signal length employed is 10 s and sampled at 48 MHz. A short-time Fourier transform is performed of SAA output. The parameters considered for the short-time Fourier transform are: discrete Fourier transform of 128 size, 64 samples hopsize, and the Hamming window. The two signals impinging on the SAA are from (45, 54) and (65, 76) degree. The signal is transformed into a short-time frequency domain via short-time Fourier transform. Each signal frame in the short-time frequency domain is transformed into an SH domain utilizing the SFT. The proposed 1D-MUSIC is then used to compute the azimuth and elevation for all signal frames in the SH domain. The average of the results for all of the frames is considered the final azimuth and elevation. The comparison of performance between DoA V -ESPRIT, 1D-MUSIC, and SH-MUSIC is as shown in Figure 8. 1D-MUSIC exhibits greater performance than DoA V -ESPRIT and SH-MUSIC. Therefore, the proposed 1D-MUSIC algorithm has been demonstrated to be a better candidate in DoA estimation. Case 8: Finally, fulfilling the growing application requirements, a specific number of elements are situated on the systems. The inter-element distance becomes shorter causing strong mutual coupling with poor radiation performance and impedance matching. In order to incorporate the mutual coupling effect that exists between elements, experimental measured data, which are the ground truth to systematically evaluate any procedure, are used. Therefore, experimental measurement data are further used for performance evaluation and analysis. The SAA is positioned in the middle of the chamber, and the source is situated at 74 DoAs, which are obtained from different combinations of 4 different elevations and 18 different azimuths. We selected the azimuths from 5 degrees to 365 degrees with 20 degrees as a step size. For detailed information on the measurement architecture/setup using SAA, readers are referred to the previous paper [4] where the measured data were first published. Gross error (GE) performance evaluation metrics is equally employed to evaluate the methods. The comparison of performance between DoA V -ESPRIT, 1D-MUSIC, and SH-MUSIC is as shown in Figure 9. 1D-MUSIC exhibits greater performance than DoA V -ESPRIT and closer performance to SH-MUSIC, even with unknown mutual coupling. Hence, 1D-MUSIC is a good candidate for estimation of the DoA, particularly in practical cases. Figure 9. GE performance evaluation of methods against SNR using the measured data in [4].

Conclusions
In conclusion, a 1D-MUSIC method for estimation of the DoA of EM waves impinging on the SAA with unknown mutual coupling is presented in this paper. The proposed 1D-MUSIC computes the elevation via a 1D spectral search using the estimated elevations for initialization. Exploring the Vandermonde of Fourier series, we were able to design two root polynomials for elevation and azimuth estimation in iteration without a requirement for pair matching. This exhibits a higher accuracy in estimations in comparison with DoA v -ESPRIT and TSADA. When compared with SH-MUSIC, the proposed 1D-MUSIC exhibits close characteristics; however, the proposed technique exhibits a high reduction in the complexity of computation of the SH-MUSIC. The RMSEs are simulated for the verification of 1D-MUSIC. Furthermore, the estimation accuracy is improved when combining the 1D-MUSIC with subspace fitting, even in the presence of mutual coupling.
In spite of the results presented in this paper, there is a sizeable number of concerns that require investigation in the future. For example, the resolution of 1D-MUSIC will be looked into and improved. Finally, we will further extend this work to calibration and measurement of the time-delay in various applications and compare them to the state-of-the-art methods.