Broadband Extended Array Response-Based Subspace Multiparameter Estimation Method for Multipolarized Wireless Channel Measurements

The clustered delay line channel model, in which each cluster consists of many rays, is widely used for link-level evaluations in mobile communications. Multiple parameters of each ray, including the delay, amplitude, cross-polarization ratio (XPR), initial phases of four polarization combinations, and the azimuth and elevation angles of arrival and departure, must be known and are measured using a channel sounder. The number of rays in every cluster is usually greater than the number of elements in the antenna array of the channel sounder, which represents a challenging issue in multipolarized channel measurements. Based on the broadband extended array response of an electromagnetic vector antenna array, a new subspace estimation method is proposed to resolve a large number of rays. The inter-element spacing of the array can be greater than half the carrier wavelength, which reduces inter-element coupling and simplifies the array design, especially for millimeter-wave bands. The delay of each cluster is first estimated using the reference antenna element. Next, two-dimensional angles of every ray are estimated using the classic rank-deficient multiple signal classification algorithm. Finally, the initial phases, XPR, and amplitude of every ray are estimated. Simulation results validate the proposed method.

propagation characteristics. Therefore, a multiparameter estimation method for in-depth analysis of channel propagation characteristics and channel modeling of the mobile communication systems is important [1], [2]. In the study item for channel models in the 3rd Generation Partnership Project (3GPP), delay, two-dimensional (2D) 1 angles of departure (AoD), 2D angles of arrival (AoA), and cross-polarization ratio (XPR) are essential small-scale parameters [3]- [5]. An accurate and efficient multiparameter estimation method is of great significance in channel measurements. Moreover, ray polarization must be considered in channel modeling. Antenna polarization modeling is also advocated by 3GPP.
Multipolarized multiple-input and multiple-output (MIMO) is a MIMO technique with transmit and receive arrays consisting of multipolarized antennas, e.g., tripolarized or hexapolarized antennas. A tripolarized antenna measures the electric field along three orthogonal directions. A hexapolarized antenna measures the electric field as well as the magnetic field along three orthogonal directions. Previous works have proved that multipolarized MIMO outperforms the corresponding unipolarized MIMO in terms of channel capacity in many situations [6]- [8]. Moreover, dual-polarized antenna configurations outperform single polarization in terms of spatial correlation [9]. The main reason for the capacity gain obtained by multipolarized MIMO relative to unipolarized MIMO is that the correlation between the antenna elements is reduced [10]. Due to these advantages, multipolarized MIMO is considered a promising technique for next-generation mobile communication systems. There has been increasing interest in polarized channel modeling and measurement, especially for multipolarized MIMO systems [10]- [13].
As a typical hexapolarized antenna, the electromagnetic vector antenna (EMVA) has been widely investigated in the literature on multipolarized array signal processing and capacity analysis [7], [14], [15]. However, due to the particularity of the EMVA structure, a collocated EMVA (C-EMVA) is complicated to construct [16]. The major difficulties originate from the antenna feeding method and the elimination of mutual coupling for the collocated elements. Such difficulties will be more apparent in the millimeter-wave multipolarized antenna design. To address these problems, a C-EMVA is 1 The two dimensions are the azimuth and elevation. 0090-6778 © 2021 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
usually reconstructed as a non-collocated EMVA (NC-EMVA) or simplified to a dual-or triple-polarized antenna. The six polarization elements in NC-EMVA are no longer collocated but are placed separately in the form of an array [17], [18]. The NC-EMVA not only simplifies the EMVA design but also suppresses the mutual coupling between the antenna elements. The selection of the inter-element spacing (IES) of non-collocated antenna elements shall consider three aspects. The first one is mutual coupling suppression, which can be simulated by electromagnetic full-wave simulation software. The second one is the placement and heat dissipation of the radio frequency (RF) modules behind the array. The third one is the mobility of the channel sounding system for the practical channel measurement campaign.
Multiparameter estimation methods such as the maximum likelihood (ML) algorithm [19], multiple signal classification (MUSIC) [20], estimation of signal parameters via rotational invariance techniques (ESPRIT) and space-alternating generalized expectation-maximization (SAGE) [21], [22] have been investigated intensively in unipolarized array signal processing. In [21], [23], [24], the delay and AoA are estimated via a subspace-based method, and the number of total rays can be greater than the number of antenna elements but not the number of rays in one cluster . 2 In [22], [25], the delay, AoA, and fading coefficient are estimated using the SAGE algorithm. The SAGE and MUSIC algorithms are compared in [25] using measurement data. The multipolarized array signal processing algorithm has also received substantial attention, particularly for the EMVA [14], [15], [26]- [29].
However, a major difficulty in channel measurements is the parameter estimation of rays using multipolarized arrays when the number of rays in the cluster is greater than the number of array elements. The number of rays in one cluster often exceeds the number of antenna elements in the channel sounder [3], [4]. Existing multiparameter estimation methods usually do not work under such conditions or suffer from low robustness. Although the virtual MIMO can address this problem indirectly by shifting the antenna in space mechanically, it is inconvenient and leads to long measurement times. Therefore, an efficient method to overcome this difficulty with lower hardware complexity is needed.
In this study, a broadband extended array response (BEAR)based subspace estimation method is proposed to address the issue of using the EMVA array (EMVAA). The dimension of the steering vector can be increased substantially via the BEAR. Traditionally, the AoA is obtained in the space domain, which requires the IES to be smaller than half the carrier wavelength. This requirement increases both the complexity of the channel sounding system and the mutual coupling between the antenna elements. The proposed method overcomes the limitation of IES, and the antenna array shape and IES are more flexible than in the traditional method. Furthermore, these advantages are more apparent in millimeter-wave MIMO 2 A cluster is defined as the path that can be resolved in the temporal domain. The number of rays in one path can be larger than the number of antenna elements in the array. channel measurement systems. The main contributions of this work are summarized as follows.
• The BEAR-based subspace estimation method is proposed to estimate the delay, 2D AoD, 2D AoA, initial phases of the four polarization combination states, amplitude, and XPR in multipath environments using a multipolarized antenna array (EMVAA). In particular, the XPR of each ray is estimated in detail, rather than providing an approximate result with the power ratio of each polarized antenna element. Most of the parameters related to channel coefficient generation in the 3GPP protocol can be estimated. • Using the proposed method, the RF module of the channel sounding system is simplified, and the mutual coupling in the antenna array is suppressed because the IES can be larger than half the carrier wavelength. Moreover, 2D angles can be estimated via linear EMVAA with low hardware complexity. Considering the complexity of the broadband loops in antenna design, the proposed method still works well with only dual-or triple-polarized antenna arrays. • The number of rays in one cluster is greater than the number of antenna elements, which is a common difficulty in channel parameter estimation. The proposed BEAR-based method utilizes the response in the spacefrequency-polarization (SFP) domain to solve this problem. The maximum number of rays that the proposed method can resolve is also given. Notations: Upper (Lower) bold-face letters are used to denote matrices (vectors). A K × K identity matrix is denoted as I K . 1 M×N represents the M × N matrix with all elements equal to 1. The kth column vector of an identity matrix is a unit vector denoted as e k . Superscript (·) H denotes the Hermitian transpose, and (·) T denotes the transpose. The Hadamard product, Khatri-Rao product, and Kronecker product are, respectively, denoted by , ⊕, and ⊗. A diagonal matrix with all entries of x on the diagonal is denoted by diag{x}. The phase angle of a complex number x is denoted by arg{x}. The pseudoinverse of X is represented as X † . vec{X} represents the vectorization of X.

II. SYSTEM MODEL
In this section, the clustered delay line (CDL) channel model, the multipolarized antenna system, and the received signal model are described. Subsequently, the SFP-domain channel response of the wireless channel measurement system is obtained. To clarify the system configuration, the parameters to be used in this work and the corresponding descriptions are listed in Table I.

A. Channel Model
In 3GPP, the CDL model is defined for the full frequency range from 0.5 GHz to 100 GHz, with a maximum bandwidth of 2 GHz, which is used mainly for link-level simulations [3]. In this channel model, the parameters that need to be estimated are the delays of clusters, 2D AoD, 2D AoA, initial phases of the four polarization combination state, amplitude, and XPR of rays. These parameters are defined in this subsection. A total of K rays in L clusters are considered for the multipath environment, where K = L l=1 K l , with K l being the number of rays in the lth cluster. As shown in Fig. 1, the rays in the lth cluster arrive at the receiver with different AoDs and AoAs but approximately the same time delay τ l . Denote the kth ray in the lth cluster as the k l th ray. The azimuth AoA and AoD (AAoA and AAoD) of the k l th ray are denoted by ϕ r,lk and ϕ t,lk . The elevation AoA and AoD (EAoA and EAoD) of the k l th ray are denoted by θ r,lk and θ t,lk . The angle definition is given in Fig. 2 (b). Define Θ t,lk {ϕ t,lk , θ t,lk } and Θ r,lk {ϕ r,lk , θ r,lk }. The direction cosine of AoD is denoted by u t,lk (Θ t,lk ) = [sinθ t,lk cos ϕ t,lk , sinθ t,lk sinϕ t,lk , cos θ t,lk ] T . (1) For the AoA, the direction cosine is u r,lk (Θ r,lk ), which is obtained by substituting the subscript 't' in (1) with 'r'. The ray may change polarization status during propagation. Such a process is represented by lk , ω hv lk , ω vh lk , ω vv lk } being the different initial phases of the horizontal and vertical polarization combinations for the k l th ray and κ lk representing the XPR of the k l th ray. The horizontal polarization is represented using h, and v for the vertical polarization. Note that XPR represents the power ratio of h and v polarization combinations for each ray due to scattering from the physical objects. As shown in Fig. 2 h and v are a set of orthogonal bases in the normal plane of the propagation direction [4].

B. Antenna Pattern and Array Configuration
In this work, a multipolarized MIMO system with M t and M r NC-EMVAs on the transmitting side and receiving side is considered. Fig. 2(a) shows the structure of the applied antenna arrays; the dashed lines depict the C-EMVA. The definitions of the azimuth and elevation angles are shown in Fig. 2(b). The coordinate of the mth C-EMVA can be denoted as r m = [r m,x , r m,y , r m,z ]. In addition, the reference coordinates of the six elements are given as r ex , r ey , r ez , r hx , r hy , and r hz for an NC-EMVA. In this study, the IES of the array can be larger than a half wavelength at the frequency of operation. Thus, the mutual coupling effect is effectively reduced, and the complicated array design is simplified.
In practical channel measurements, the gains of the transmit and receive antenna elements usually vary with the departure and incident angle. The antenna gains are taken into consideration in the received signal modeling progress to de-embed the influence of the antenna pattern. As shown in Fig. 2(a), the dipoles and loops are arranged in different directions for an NC-EMVA array (NC-EMVAA), so each element of an NC-EMVA has a unique antenna gain. Here, the gains of the transmitting NC-EMVA from Θ t are denoted as G t,ex (Θ t ), The gains vector of the transmitting antenna array for the k l th ray is written as Similarly for the AoA, the gains are denoted by substituting the subscript 't' with 'r'. Let a t (f, Θ t,lk ) be the M t × 1 space-frequency domain NC-EMVAA steering vector of the k l th ray for the AoD. It can be expressed as where c is the speed of light. For the AoA, it is denoted by a r (f, Θ r,lk ).
Let A pdt (f, Θ t,lk ) be the 6 × 2 SFP-domain steering matrix of the k l th ray for an NC-EMVA in the transmitter. The 6 rows correspond to the 6 antenna elements of an NC-EMVA, and the 2 columns correspond to the polarization components along the v and h directions. It can be expressed as where D t (f, Θ t,lk ) is the diagonal phase shift matrix of a NC-EMVA for the AoD, i.e., The phase shift vector d t,lk (f, Θ t,lk ) is denoted as (7), shown at the bottom of the page. The space-polarization domain steering matrix for the C-EMVA is expressed as For the AoA, the corresponding expressions are denoted by A pdr (f, Θ r,lk ), Ω(Θ r,lk ), D r (f, Θ r,lk ) and d r,lk (f, Θ r,lk ).

C. Transmit-Receive Signal Model
In this subsection, the transmit-receive signal model is analyzed. First, the transmitting sequence structure is designed. Then, the receive signal model is given for the SFP-domain channel response derivation.
As shown in Fig. 3, there are 6M t transmitting antenna elements respectively numbered (m t , i)th antenna , 3 where m t = 1, . . . , M t and i = 1, . . . , 6. The transmitting sequence of each antenna element is composed of L seq subsequences. Note that L seq ≥ 6M t . The length of each subsequence is N s , which is also equal to the number of discrete Fourier transform (DFT) points. In practical implementations, the DFT uses fast Fourier transform (FFT). Then, the n s th symbol of the l seq th subsequence for (m t , i)th transmitting antenna is denoted as s mt,i [(l seq − 1)N s + n s T ] with n s = 1, . . . , N s and l seq = 1, . . . , L seq , which is expressed as s mt,i,lseq (n s T ) with T being the symbol width for compact expression.
The transmitting sequences for all transmitting antennas in matrix form can be denoted as where S is a 6M t ×L seq N s matrix with 6M t rows corresponding to the 6M t transmitting antenna elements and L seq N s columns corresponding to the L seq subsequences with length N s . The l seq th subsequence for all transmitting antennas can be written as where s lseqns is the n s th symbol of the l seq th subsequence for all transmitting antennas, which is defined as After elaborating the transmitting sequence, the receive signal model is then described. Thus, the n s sample of the l seq subsequence in the nth snapshot of the received signal can be denoted by (12), where z n is the complex white Gaussian noise of the nth snapshot and α n,lk represents the amplitude of the k l th ray in the nth snapshot. The channel is considered as time-varying for different snapshots. The received signal is obtained by summarizing the signal carried by all rays. 3 The (mt, i)th antenna is defined as the ith element of the mtth NC-EMVA.
The remaining variables used in (12) and their corresponding descriptions can be found in Table I.

D. Channel Response in the SFP Domain
In this subsection, the channel response in the SFP domain is derived via the received signal. Let Y nlseq be the received signal of the l seq th transmitting signal in the nth snapshot, i.e., where Y nlseq (f, Θ t , Θ r , Ξ) and y nnslseq (f, Θ t , Θ r , Ξ) are replaced by Y nlseq and y nnslseq for simplicity, respectively. Performing the DFT on each row of Y nlseq ∈ C 6Mr×Ns yields where Y nlseqf ∈ C 6Mr×Ns is the SFP-domain received signal of the l seq th subsequence in snapshot n. There are L seq subsequences in total. DFT is performed on each subsequence, and the SFP-domain received signal in frequency bin f i is selected. Then, the received signal of the nth snapshot in frequency bin f i Y nfi ∈ C 6Mr×Lseq can be obtained as where H nfi ∈ C 6Mr×6Mt and Z nfi respectively represent the channel response matrix and additive white Gaussian noise (AWGN) matrix in the nth snapshot in frequency bin f i . The standard deviation of the noise is σ n with mean zero.S fi ∈ C 6Mt×Lseq is the frequency domain transmitting signal in frequency bin f i for all subsequences, which can be expressed asS where F stands for the DFT matrix and J fi is the selection matrix that selects the signal of the frequency bin f i , which is defined as J fi = I Lseq ⊗ e i . It is easy to design a signal matrixS fi with full row rank. H nfi is the basic component of the SFP channel response in frequency bin f i and snapshot n. Therefore, it is vital to analyze the derivation of H nfi . According to (12) and (15), H nfi is expressed as For simplicity of expression, a r (f i ,Θ r,lk ), g r (Θ r,lk ) and A pdr (f i , Θ r,lk ) will be denoted by a rfi,lk , g r,lk and A pdrfi,lk in the following derivations. Similar simplified expressions are also applied in a tfi,lk , g t,lk and A pdtfi,lk . Define A tfi,lk ∈ C 6Mt×2 as the transmitting array steering matrix of the k l th ray in frequency bin f i , which is denoted as The corresponding receiving array steering matrix is denoted by substituting the subscript 't' in (18) with 'r'. Let can be expressed as where b trfi,lk ∈ C 36MtMr×1 is the steering vector of the k l th ray in frequency bin f i given as Here, A trfi,lk ∈ C 36MtMr×4 is the steering matrix of the k l th ray in frequency bin f i , which is defined as In (20), β lk = vec{T lk } is a column vector. For compact expression, let us rewrite (19) as where and α n is the amplitude vector for all rays in snapshot n denoted as α n = [α n,11 , . . . , α n,1K1 , . . . , α n,LKL , . . . , α n,LKL ] T . (24) In (22), F fi ∈ C K×K is the phase shifting diagonal matrix of frequency bin f i caused by the delay written as Then, h n,fi can be estimated aŝ Thus, the estimated channel response in the SFP domain for all frequency bins and snapshotsĤ ∈ C 36MtMrNs×N can be obtained asĤ which can be further expressed aŝ where B ∈ C 36MtMrNs×K , Γ andẐ respectively represent the steering matrix of the array, the amplitude of rays for N snapshots and the noise matrix, which are given as

A. Maximum Likelihood Estimator
The maximum likelihood estimator is considered an accurate estimator. The likelihood function can be written compactly as Thus, the ML estimate of τ , Θ t , Θ r and Ξ is obtained by maximizing F (τ , Θ t , Θ r , Ξ). Due to the large number of impinging rays and multiple parameters, the ML estimator consumes substantial computational resources. Therefore, it is proposed that these multiple parameters are estimated one after the other to reduce computational complexity.

B. Delay and Angle Estimation
In this subsection, the delay and angle estimation method is derived. The propagation delay of a ray is estimated from the frequency domain channel response of the reference antenna element. In this work, the first antenna element of the first NC-EMVA is selected as the reference antenna element of the transmitting and receiving antenna array. Define a selection matrix J τ = I Ns ⊗ e T 1 , where e 1 is the first column of I 36MrMt . The role of the selection matrix J τ is to select the channel response of the reference antenna element from the SFP-domain channel response. Then, the selection progress is expressed asĤ LetR τ = 1

NĤ τĤ
H τ be the covariance matrix ofĤ τ . The rank ofR τ is the number of clusters without considering the noise. The eigenvalue decomposition (EVD) ofR τ leads to the N s × (N s − L) noise subspace U n,τ . The delay MUSIC estimation can be expressed as where a τ is the steering vector in the frequency domain of the reference antenna element, which can be expressed as With the delay estimated by the above MUSIC technique, the angle estimation is then processed. The AoD can be estimated with the channel response related to the transmitting antenna elements and the reference receiving antenna element. Define a selection matrix J t = I 6NsMt ⊗ e T 1 , where e 1 is the first column of I 6Mr . The role of the selection matrix J t is to select the channel response of the transmitting antennas and the reference receiving antenna element from the SFP-domain channel response. Then, the selected channel response is given asĤ whereB t is the steering matrix of the transmitting array channel response, i.e., Here,B fi is obtained as where J t,ref = I 6Mt ⊗ e T 1 with the length of e 1 being 6M r . The role of the selection matrix J t,ref is to select the channel response of the transmitting antennas from the SFP-domain channel response in any frequency bin. LetR t = 1 NĤ tĤ H t be the covariance matrix ofĤ t . The rank ofR t is equal to the number of rays without considering the noise. The EVD ofR t leads to the 6M t N s × (6M t N s − K) noise subspace U tn . The estimated AoD can be obtained bŷ where P t (Θ t,lk ) is the AoD MUSIC spectrum, Similarly, the AoA is estimated with the channel response related to the receiving antenna elements and the reference transmitting antenna element. Define a selection matrix J r = I Ns ⊗ e T 1 ⊗ I 6Mr , where e 1 is the first column of I 6Mt . Then, the receiving array channel response iŝ whereB r is the steering matrix of the receiving array channel response, i.e., Here,B fi is obtained as where J r,ref = e T 1 ⊗ I 6Mr with the length of e 1 being 6M t . Note thatB fi in AoA estimation is different from that in AoD estimation, which is caused by the different selection matrices J t,ref and J r,ref . LetR r = 1 NĤ rĤ H r be the covariance matrix ofĤ r . The rank ofR r is equal to the number of rays without considering the noise. The EVD ofR r leads to the 6M r N s × (6M r N s − K) noise subspace U rn . The estimated AoA can be obtained bŷ where P r (Θ r,lk ) is the AoA MUSIC spectrum,

C. Frequency Domain Smoothing Technique
In the application of the BEAR-based method, the FFT points are usually too large to perform EVD on the covariance matrix of the channel response. To reduce the effect of the AWGN and the computational complexity, a frequency domain smoothing preprocessing technique is proposed. The whole frequency band is divided into P subfrequency bands, and each subfrequency band corresponds to a virtual subarray. The number of frequency bins in each subfrequency band is N v . And N v should meet the condition min{6M t N v , 6M r N v } > K. The difference in frequency bins in the adjacent subarray is the DFT frequency interval. Then, the channel impulse responses of all the virtual subarrays are summarized. Fig. 4 shows a diagram of the virtual subarray division. The closed-form MUSIC spectra of the AoD and AoA can be derived by means of this technique. See Appendix B for details of the derivation process.

D. Pair Matching of AoD and AoA
The estimated AoD and AoA must be paired. The channel response in any frequency bin, for example, the top 36M r M t rows of H, can be used to perform pair matching. Let H pair = J pair H, where J pair = e T 1 ⊗ I Ns with the length of e 1 being 36M r M t . EVD of H pair yields the noise subspace U pair,n .
The steering matrix for the channel response of the selected frequency bin is constructed as The aim of the pair matching of AoD and AoA is to find the AoD-AoA pair such that P pair = arg min Θ r,lk ,Θ t,lk det A H pair U pair,n U H pair,n A pair . (48)

E. Estimation of Initial Phases, XPR and Amplitude
With the estimated delay and paired angles of the rays, the remaining parameters, XPR, amplitude and initial phases can be estimated. Let For the k l th ray, the XPR can be estimated aŝ With the estimated initial phase and XPR, Ψ can be reconstructed asΨ; then the amplitude of the rays for all N snapshots is estimated asΓ

F. Derivation of the Cramer-Rao Lower Bound (CRLB)
The method of CRLB derivation in [30] is used in this work. The parameters that must be estimated are (ϕ t,lk , θ t,lk , ϕ r,lk , θ r,lk , τ l ), with l = 1, . . . , L and k = 1, . . . , K L . The CRLB is given as where ϑ represents one of the parameters to be estimated, Λ n = diag{α n,lk }, and Υ is given in Appendix C.

G. Implementation of the Bear-Based Subspace Method and Comparisons With Existing Methods
The multiparameter estimation procedure using the BEAR-based subspace method is summarized in Table II. Comparisons with existing methods are presented in Table III. The existing methods cannot resolve rays in the same cluster when the number of rays is greater than the number of antenna elements. For the genuine SAGE algorithm, the resolvable rays are not given clearly, but the low robustness caused by overlapping signals in the temporal domain and algorithm complexity limit its practical implementation. Even so, the SAGE algorithm in [22] is still a widely used algorithm for channel parameter estimation. It is suitable for time invariant channel environments and specular rays extraction. For researchers are interested in the AoA and delay, the SAGE algorithm is a nice candidate approach. However, when researchers are interested in not only the delay and AoA parameters but also the other parameters in 3GPP, BEAR is a potential candidate parameter estimation approach.

H. Complexity Analysis
The computational complexity of the proposed method is calculated in terms of complex multiplications. To reduce the computational burden, the N -point DFT uses the FFT algorithm with computational complexity of N 2 log 2 N . The computational complexity for an EVD of an M × M matrix is O(M 3 ). In Table IV, the complexity of estimating the parameters for all K rays is shown in detail where W τ , W AoD and W AoA represent the number of rank-deficient MUSIC search times, W τ is the one-dimensional search times, which is inversely proportional to the search step, and W AoD and W AoA are 2D search times, which are inversely proportional to the square of the search step. The AoD and AoA estimation is based on the delay of clusters. Thus, the search process needs to be performed L times. To give a more intuitive complexity analysis, the computational complexity of the proposed algorithm is compared with the one of the SAGE algorithm in Fig. 5. The complexity of the SAGE algorithm relies on the number of iteration cycles which is assumed to be 2 here. In actual applications, the number of iteration cycles is affected by the initial values and the convergence speed, which is usually larger than 2. It is easy to find that the complexity of the BEAR is directly related to the virtual subarray size N v .

I. Notifications in BEAR Application
In the above subsections, the BEAR-based subspace multiparameter estimation method is analyzed. Two points warrant further discussion.
The first is the selection of the frequency bins in virtual subarray 1. To avoid correlations between the steering vectors of the rays in different clusters, the minimum frequency bins interval Δf in virtual subarray 1 needs to meet the following conditions: where τ max is the maximum delay of all clusters. The second is the selection of virtual subarray numbers P . The frequency domain smoothing procedure is not an in-phase summation progress in (62) and (63). By smoothing the receiving signals, the SNR gain is expressed by [31] SNR sm (P ) = 1 P Then, the optimal P is obtained by maximizing (57).

IV. NUMERICAL SIMULATIONS AND ANALYSIS
In this section, simulations are performed to show the performance of the proposed BEAR-based subspace estimation method. Multipath propagation channels with additive white Gaussian noises are used. First, the performance of different polarized arrays in the 6×6 MIMO system is given to demonstrate the advantages of the NC-EMVAA. Next, the detailed simulation setup for the comprehensive simulations is given. Then, the comprehensive simulation results of the proposed algorithm are shown. Furthermore, The probability of success and the RMSE of angle estimation are investigated in detail. The XPR, initial phases, and amplitude estimation results are finally given.
In the simulation, high performance of the proposed BEAR-based method is achieved under different SNRs and array sizes. To evaluate the estimation accuracy, the root mean square error (RMSE) of the angle estimation is given as where Q is the number of experiments and ξ represents the azimuth angle or the elevation angle. The condition for judging the success of the estimation is that the sum of the RMSE of the azimuth angle and the elevation angle is less than 1 degree. The probability of success is obtained by the success times in 100 simulations.

A. Performance Comparison of Different Polarized Arrays
In Fig. 6, with different polarized antennas in the 6 × 6 MIMO system, the performance comparison of the probability of success is shown. The 6 antenna elements are arranged as a 3 × 2 rectangle array with IES = 15λ c . There are 2 clusters with the number of rays in each cluster being 7. The subsequence number is 100. The rest parameters are given in Table V. It is found that the probability of success achieves 100% for both the uni-polarized array and the cross-polarized array at the SNR of about 9 dB. For the NC-EMVA, the required SNR is only about 6 dB. That is to say, the NC-EMVA has a better performance than uni-polarized and cross-polarized arrays. It is worthy to notice that the XPR cannot be estimated when using the uni-polarized array.

B. Simulation Setup
In this subsection, the setup of the simulations is given in detail in Table V. For the system configuration, the setup of the transmitting antenna array and receiving antenna array are the same, and both arrays are equipped with two EMVAs (12 antenna elements), respectively. The 12 antenna elements are arranged as a rectangle array with size 3 × 4. The carrier frequency is set to be 28 GHz, the IES is larger than half of the carrier wavelength, and the bandwidth of the signals is 500 MHz. To investigate the robustness and accuracy of the proposed method, three kinds of situations (IES= 5λ c , 10λ c , and 15λ c ) are considered where λ c is the wavelength of the carrier frequency. The length of the subsequence is set to N s = 1024, which is equal to the number of FFT points. The number of subsequences is set to L seq = 300. In the frequency domain smoothing procedure, the number of virtual subarrays is set to P = 8. The number of snapshots is set to N = 100, which is larger than the number of rays.
Two clusters are generated in the simulation. The delays of the two clusters are set to τ = [100, 110] ns. Each cluster contains 14 rays. The number of rays is larger than the number of antenna array elements, which is 12 in the simulation. The EAoD and EAoA are set in the range of [0 • , 180 • ]. The AAoD and AAoA are also set in the range of [0 • , 180 • ]. The angles of some rays are set close to each other. These rays are usually scattered or reflected by the closely spaced scatters or reflectors in the propagation environment. The XPRs of all rays are set uniformly from 1 to K to cover the possible range for XPRs. The initial phases are set randomly from −π to π. The amplitude of each ray is set to fluctuate around a fixed value for different snapshots. Some of the rays are set with larger amplitudes, and the rest are set with smaller amplitudes to show the estimation ability for the rays with different powers.

C. Parameter Estimation With Different SNRs
The normalized MUSIC delay spectrum under different SNRs is shown in Fig. 7. Delays can be accurately estimated even under very low SNR conditions. The AoA and AoD estimations are based on the delay estimation results, so accurate delay estimation is essential. The AoD and AoA MUSIC spectra under different SNRs are given in Fig. 8. The rays in one cluster usually result from the same scatter or scatters with a nearby position. Therefore, the angles are distributed in a small scope. Furthermore, some of the rays in one cluster may have the same azimuth or elevation angle. In the simulation of Fig. 8, the azimuth and the elevation of some rays are concentrated within a small degree range. The angle range of the rays in one cluster is [0 • , 180 • ]. The proposed method can be used to estimate the AoD and AoA accurately when the number of rays in one cluster is larger than the number of antenna elements, even if some of the rays have the same azimuths or elevations. As the SNR increases, the peaks of the rays become more clear. The angles can be estimated unambiguously when the space between two antenna elements is larger than half the carrier wavelength, which is different from traditional methods. The reason for the unambiguity of the angle estimation is that the virtual array elements are generated via a broadband array response. Fig. 9 shows the probability of success for the AoD and AoA estimation. The IES is 15λ c , 10λ c , or 5λ c and N v is set to 12 or 128. It can be seen from the results that N v will not have a big impact on the probability of success. This is because the frequency bins in the virtual subarray are of high redundancy for the steering matrix construction. Note that the reduction of N v will lead to a decrease in estimation complexity. Therefore, the selection of the frequency bins should take the estimation complexity constraints and (56) into consideration in the practical implementation. For larger IES, the probability of success approaches 100% under a   lower SNR. This is because the array aperture increases with increasing IES under the same bandwidth. Fig. 10 gives the RMSE of the estimated AAoD and EAoD with different bandwidths. Fig. 11 shows the RMSE of the  estimated AAoA and EAoA with different bandwidths. Here, the N v is set as 128. The results show that the RMSE is affected by the bandwidth. At the same SNR, a wider signal bandwidth leads to better RMSE performance. This is because the larger bandwidth will lead to a larger virtual array aperture. The experiments are performed one hundred times to reduce uncertainty. The CRLB of each situation is given for comparison. The simulation results show that the RMSEs of the proposed method are close to those of the CRLB as the SNR increases. The XPR estimation result of each ray in all clusters, given in Fig. 13, is accurately estimated. The exact value stands for the XPR set in the simulation. The estimated value stands for the XPR estimated by the proposed method. In the simulation, there are two clusters with 14 rays in each cluster. Fig. 12 gives the estimation results of the four initial phases under 10 dB SNR. The four initial phases are randomly generated in the range of [−π, π]. The exact value stands for the initial phases set in the simulations. As shown, the estimated value and the exact value of the four initial phases fit well for all of the rays under 10 dB SNR. Fig. 14 gives the estimation results of the amplitude of all the rays in N snapshots. In the simulation, most of the rays are set with weak power fluctuating in a small range, and only a few  rays are set with strong power. This situation setting is usually practical in a real radio propagation environment. The results show that the amplitude for all of the rays can be estimated accurately under a SNR of 10 dB.
The minimum SNR required to extract all the rays successfully under different L seq with different array sizes is plotted in Fig. 15. The required minimum SNR decreases as L seq and the array size increases because the increase in L seq can suppress the influence of noise and the increase in the array size can enlarge the array aperture. The simulation results can be used to guide the design of a multipolarized channel measurement system that uses the proposed BEAR-based subspace estimation method.

V. CONCLUSION
A BEAR-based subspace estimation method has been proposed to estimate the delay, azimuth and elevation angles of arrival and departure, amplitude, XPR, and initial phase with EMVAA for multipolarized channel measurements. By extending the array response in the SFP domain, the proposed method can estimate a large number of rays in one cluster and is not limited by the number of antenna elements. The proposed method works even if the rays have similar azimuths and elevations. The XPR, which is usually obtained coarsely based on the power ratio of each polarized antenna element, can be estimated accurately with the proposed method. In addition, the four initial phases of different polarization combinations are estimated accurately. The IES can be larger than half the carrier wavelength without angle ambiguity in the angle estimation. Thus, mutual coupling is effectively suppressed. Simulation results show that larger bandwidth leads to a smaller RMSE of the AoA and AoD. Moreover, larger IES leads to smaller required SNR. The proposed method holds promise for multiparameter estimation in future multipolarized channel measurements.

APPENDIX A DERIVATION OF THE RANK-DEFICIENT MUSIC SPECTRUM OF AOD
The steering matrix of the channel response matrix related to the transmitting antenna array and the reference receiving antenna element isB t , which is related toB fi as (38). According to (20) and (39), the column ofB fi can be rewritten asb fi,lk = (A tfi,lk ⊗ ρ)β lk .
ρ is a 1-by-2 vector, which represents the polarization component of the reference receiving antenna element. In this work, it is denoted by the first row of Ω(Θ r,lk ). Let A t,lk be A t,lk = e −j2πf1τ l A T tf1,lk , . . . , e −j2πfN s τ l A T tfN s ,lk Then, the column ofB t can be written asb t,lk = (A t,lk ⊗ ρ)β lk . In order to decrease the MUSIC searching dimension, the rank-deficient MUSIC technique is applied in AoD estimation [32]. Let A AoD,lk = A t,lk ⊗ ρ, the spectrum can be written as P t (Θ t,lk ) = det A H AoD,lk (Θ t,lk )U tn U H tn A AoD,lk (Θ t,lk ) .
Because the columns of A AoD,lk depend linearly on A t,lk ,  . . . , f Ns−P +p ], whereñ v ranges from 1 to N s −P +1 and P is the number of virtual subarrays. Construct a selection matrix J f,p with size N v -by-N s to select the frequency bins of the pth subarray. The n v th row of J f,p corresponds to theñ v th row of I Ns . Let us define a selection matrix for the pth subarray in the smoothing progress of AoD estimation J t,sm,p = J f,p ⊗ I 6Mt .
EVD of R t,P leads to the 6M t N v × (6M t N v − K) noise subspace U P tn . The smoothed AoD MUSIC spectrum can be expressed as P t (Θ t,lk ) = det Ā H t (Θ t,lk )U P tn U H P tnĀt (Θ t,lk ) where the steering matrixĀ t (Θ t,lk ) is given as A t (Θ t,lk ) = J t,sm,1 A t (Θ t,lk ) (Σ t e lk ⊗ 1 1×2 ).
Note that the antenna pattern is discrete in the space domain and cannot be differentiated. Without loss of generality, it is regarded as a unit gain in the calculation progress. The details of the derivation process for the CRLB are similar to those of [30] and are not presented.