Rui Feng, Yu Liu, Jie Huang, Jian Sun, Cheng-Xiang Wang,,＊, George Goussetis

1 Shandong Provincial Key Lab of Wireless Communication Technologies, Shandong University, Jinan, Shandong 250100, China

2 Institute of Sensors, Signals and Systems, School of Engineering & Physical Sciences, Heriot-Watt University, Edinburgh, EH14 4AS, U.K

I. INTRODUCTION

With the development of wireless communications, several key technologies such as millimeter wave (mmWave) and massive multiple-input multiple-output (MIMO) were proposed to provide large bandwidth and to increase spectrum efficiency, etc [1]–[4]. Wireless channels using these key technologies exhibit some new channel characteristics. For example, mmWave channels show larger attenuation and higher delay resolution than the conventional channel [5]; in massive MIMO channel, spherical wavefront instead of plane wavefront, as well as the non-stationarity in space domain should be considered [6]–[8]. In order to analyze these new channel characteristics, extensive channel measurements should be carried out and channel parameter estimation algorithms are required to extract both temporal and three-dimensional (3D) spatial information from measurement data. As signals emitted from a transmitter (Tx) arrive at the receiver (Rx) through multipath components (MPCs), the extraction of temporal and spatial information corresponds to estimate delay, azimuth angle of arrival (AAoA), and elevation AoA (EAoA) for each MPC. Considering that there are many array signal processing(ASP) methods, this work aims to summarize the commonly used channel parameter estimation methods in estimating AAoA and EAoA(and delay), and compare their performance in both synthetic and real channel environments.

In this paper, the commonly used array signal processing (ASP)methods to estimate channel parameters are summarized.

Originally, the ASP methods were proposed to estimate spatial information, namely AoA only. As shown in figure 1, they can be divided into three categories, i.e., spectral estimation,parametric subspace-based estimation, and deterministic parametric estimation [9]. Spectral estimation techniques estimate AoAs by forming spatial spectrum and searching for the angles that correspond to the spectrum peaks.This method can be further classified into beamforming and subspace-based techniques[10], [11]. Conventional (Bartlett) beamforming and Capon's beamformer [12] belong to the beamforming technique. The Bartlett beamforming forms spatial spectrum by using weighting vector that maximizes the output power at a certain direction, while Capon's beamformer using weighting vector not only maintains the power constant to the desired direction, but also minimizes the output power from undesired direction. Thus, Capon's beamformer can provide better spatial resolution than the Bartlett beamforming. Multiple signal classification (MUSIC) [13] is a representative algorithm of the subspace-based technique. Different from the beamforming technique, MUSIC algorithm separates the spatial covariance matrix into signal and noise subspaces, and employs the orthogonality between noise subspace and signal subspace to form spatial spectrum. Because the calculation of spatial spectrum employs only noise subspace rather than the whole covariance matrix,MUSIC algorithm is less time consuming than the beamforming technique. In the parametric subspace-based estimation kind of method,estimation of signal parameter via rotational invariance techniques (ESPRIT) [14], least square (LS) ESPRIT, total least square (TLS)ESPRIT, and Unitary ESPRIT [15]–[18] are representative algorithms. The ESPRIT-type algorithms separate the whole Rx array into several sub-arrays, and calculate AoAs directly by employing the invariance properties among the sub-arrays. Unitary ESPRIT is an extension of ESPRIT algorithm which can transform the complex-valued calculation into real-valued domain, thus, to provide reduced computation burden [19]. For the aforementioned spectral estimation and parametric sub-space-based estimation methods, it is required that the MPCs are uncorrelated and the total MPC number is smaller than the Rx array number. However, deterministic parametric estimation kind of method is free from these constraints. Maximum likelihood (ML) technique estimates AoAs by searching for angles that maximize the ML function. It can achieve high accuracy but with very large complexity. Iterative quadratic maximum likelihood(IQML) was proposed to estimate AoAs with the cost function being the weighted LS cost functions [20]. However, it requires multiple snapshots of measurement data and suffers from initialization problem.

Fig. 1. Classification of channel parameter estimation methods.

To further improve the accuracy of parameter estimation and distinguish more MPCs,delay estimation was taken into consideration.Some of aforementioned algorithms were extended to estimate both temporal and 3D spatial information. For example, temporal-spatial-temporal MUSIC (TST-MUSIC) algorithm separates the delay and angle searching procedure into three steps, i.e., two temporal-MUSIC to estimate delay and one spatial-MUSIC to estimate AoA [21]. Considering the pairing of MPC parameters, joint angle and delay estimation (JADE) scheme was proposed which can be combined with MUSIC [22] and ESPRIT-type algorithms to estimate delay and AAoA jointly [23], [24]. TST-MUSIC and JADE scheme break the constraints on the first two categories of ASP methods, and can be used when the total MPC number is larger than the Rx array number. However, they can only be used when multiple snapshots are available. For the case of only one snapshot is applicable, e.g., channel measurements carried out in frequency domain, shift-invariance (SI)JADE algorithm was proposed based on the LS-ESPRIT algorithm to estimate delay and AAoA using a uniform linear array (ULA)[25]. It was also extended to be used in two ULAs to estimate delay, AAoA, and EAoA jointly. Expectation-maximization (EM) [26],space-alternating generalized EM (SAGE)[27], [28], and Richter's maximum-likelihood estimation (RiMAX) algorithms are ML based methods. They can estimate parameters by iteratively searching for the delays and angles that can maximize the ML functions with high accuracy. However, they are usually time consuming due to the iterative searching procedure. SAGE algorithm is an improvement of the EM algorithm by separating the estimated parameters into several subsets in the maximization procedure, thus, to reduce the complexity of the EM algorithm [29], [30]. RiMAX algorithm was proposed to extract not only the MPCs, but also the diffuse multipath components (DMCs) [31].

In terms of wireless channel measurements,for example carried out in frequency domain using vector network analyzer (VNA), one key problem to be solved is to extract interested parameters in a timely and accurately manner.Many research groups used to apply SAGE algorithm to process the measurement data [32].However, in real time channel measurement,algorithms which can provide accurate result with low complexity are preferred. Thus, some researchers resort to classical subspace based methods to reduce the complexity compared with SAGE and RiMAX algorithms [33]. In[34], the performance of MUSIC and SAGE algorithms was compared to estimate AoAs.It was indicated that MUSIC algorithm could provide comparable performance with SAGE algorithm, but with highly reduced processing time. In [35], the Unitary ESPRIT algorithm was introduced to process the massive MIMO channel measurement data. It was indicated that the Unitary ESPRIT algorithm works well for the massive MIMO channel. However, these comparisons are confined in the estimation of AoA. Thus, considering that the comparison of joint temporal and 3D spatial information estimation algorithms is rare, this paper aims to fill this gap and give a thorough summarization of ASP methods in estimating wireless channel parameters. The main contributions of this paper are listed as follows:

1) Performance of several ASP methods to estimate AAoA and EAoA is compared, i.e.,MUSIC, LS-ESPRIT, and Unitary ESPRIT algorithms;

2) Performance of several ASP methods to jointly estimate delay, AAoA, and EAoA is compared, i.e., JADE MUSIC, JADE ESPRIT, SI-JADE, and SAGE algorithms;

3) An improved SI-JADE algorithm is proposed by incorporating with the Unitary ESPRIT algorithm. It can provide reduced computation complexity compared with conventional SI-JADE algorithm;

4) The MUSIC, LS-ESPRIT, Unitary ESPRIT, SI-JADE, improved SI-JADE, and SAGE algorithms are used to the 60 GHz real channel measurement data processing;

5) As the development of wireless communication poses some challenges to the parameter estimation schemes, several future directions are given regarding to the measurement data processing of fifth generation(5G) wireless channel.

The remainder of this paper is organized as follows. System model is given in Section II. Basic concepts of some channel parameter algorithms are introduced in Section III. Simulation results and measurement analysis are given in Section IV. Future directions are given in Section V. Finally, conclusions are drawn in Section VI.

II. SIGNAL MODEL

Fig. 2. Configuration of the Rx array.

We assume the Rx is configured with a horizontal uniform rectangular array (URA) lying in the x-y plane, as shown in figure 2. The total antenna element number is M = Mx×My,where Mxand Mydenote antenna numbers in the x and y axes, respectively. In wireless channel propagation, signals arrive at Rx through MPCs formed by reflection, diffraction, and scattering mechanisms. Considering a narrowband far-field wireless communication scenario, the output at Rx array contributed by the l-th MPC can be expressed as

where Θl=[βl, τl, φl(shuí), θl] is the parameter vector of the l-th MPC, with βl, τl, φl(shuí), and θldenote complex amplitude, delay, AAoA,and EAoA, respectively. L is the total MPC number, s( t) is the transmitted signal, and Nl(t) is the complex white Gaussian noise.The a(ul, vl) is a steering vector representing the Rx array response to the l-th MPC and can be written as

whereanddenote the Kronecker product and transpose operator, respectively. Considering plane wavefront, i.e., MPCs impinge at different array elements with the same AAoAs and EAoAs, phase shifts uland vlcaused by antenna displacements in x and y axes are shown in figure 3. They can be expressed as

where λ denotes the wavelength. The equidistance between two adjacent antenna elements is dxalong x-axis and dyalong y-axis. They are usually assumed equal to or less than half of the wavelength [19].

The total output signal at Rx array is the summation of all MPCs, which is the replicas of the transmitted signal with different complex amplitudes, delays, AAoAs, and EAoAs.It can be expressed as

whereis the s t e e r i n g m a t r i x,and N(t) is the composited noise. The K-point sampled output signal is Y (k ) = ABS (k ) + N(k ) (1 ≤k≤K ),where S(k) and N(k) are the sampled transmitted signal and noise, respectively. For the convenience of understanding and following parameter estimation, the matrix structure of Y(k) is shown in figure 4. The output signals contributed by the l-th MPC at two adjacent data blocks along dimension 1, namely at two antenna elements along x-axis of Rx array,have phase shift e?julcaused by dx. Similarly, along dimension 2 (y-axis) have phase shift e?jvlcaused by dy. It can be seen that the output signal along y-axis is stacked into that along x-axis according to (2). Data blocks in dimension 3 represents K sample points.

The subspace-based method relies on the separation of eigenvectors of the spatial covariance matrix, which can be given as follows,

where [·]Hdenotes the Hermitian operator and P = BS (k )SH(k )BHis the source covariance matrix. If the MPCs are assumed to be uncorrelated, the rank of the source covariance matrix r(P)=L when L＜M. Therefore,r( APAH)=L＜M. According to the properties of matrix, if noises are assumed to be uncorrelated, the sum of two matrices R has(M?L) zero eigenvalues when the number of observation data is in finite. However, the number of observation data is finite in practice.Thus, these zero eigenvalues are corresponding to (M?L) minimum eigenvalues of,whereis the sampled spatial covariance matrix calculated based on the finite observation data. Therefore,can be separated into two segments,

Fig. 3. Spherical coordinates and phase shifts.

Fig. 4. Matrix transposition structure of Y(k).

with eigenvectors corresponding to the L larger eigenvalues Λs=diag{γ1,...,γL} form the signal subspace E1,sand eigenvectors corresponding to the (M?L) smaller eigenvaluesform the noise subspace E1,n. Therefore, subspace-based method works only when the total MPC number is smaller than the Rx array number and the MPCs are uncorrelated.

III. 3D CHANNEL PARAMETER ESTIMATION ALGORITHMS

In this section, channel parameter estimation algorithms are introduced in two parts.The first part is the summary of joint AAoA and EAoA estimation methods, i.e., MUSIC,LS-ESPRIT, and Unitary ESPRIT algorithms.The second part is the basic concepts of joint delay, AAoA, and EAoA algorithms, i.e.,JADE MUSIC, JADE ESPRIT, SI-JADE, and SAGE algorithms. All algorithms are introduced to be used for a URA. Furthermore, an improved SI-JADE algorithm is proposed by incorporating with the Unitary ESPRIT algorithm in the second part.

3.1 Joint AAoA and EAoA estimation methods

In table 1, the preliminary comparison of joint AAoA and EAoA estimation algorithms are listed. The detailed procedures of them can be found in the following.

1) MUSIC algorithm: Considering that the noise subspace is orthogonal to the signal subspace theoretically, which spans the same space with A, the following relationship holds

The angle searching problem can be transferred to calculate the Euclidian distance betweenand A in practical, and find the AAoAs and EAoAs that can minimize the distance. Thus, the normalized MUSIC spatial spectrum is defined as the inverse of the Euclidian distance [13], i.e.,

Here, u and v are the phase shifts in x and y directions which corresponding to one certain direction with azimuth angle φ and elevation angle θ. Finally, the AAoAs and EAoAs can be located by finding the local maxima of the spatial spectrum.

2) LS-ESPRIT algorithm: As shown in figure 2, the whole Rx array can be divided into two sub-arrays by selecting the first and last(Mx? 1) columns of antenna elements, i.e.,Sub-arrays 1 and 2. The two sub-arrays are exactly the same except for a known constantdisplacement in x-axis. Thus, the steering matrices of two sub-arrays have the following relationship,

Table I. Properties of joint AAoA and EAoA estimation algorithms.

wherewith I is the unit matrix and 0 is the zero mat r i x,, a n d. As the steering matrix A spans the same space with the signal subspace, the following relationship holds,

where T1is a L×L full rank matrix. Substituting (11) into (10), we can get

where E1,s1= J1E1,sand E1,s2= J2E1,s.Through observation it can be seen that, the eigenvalues ofare au. Another way to get auis the LS method,

Similarly, Rx array can be divided into two sub-arrays with a known constant displacement in y-axis, namely Sub-arrays 3 and 4. The corresponding selecting matrices areandThus, the following relationship holds,

whereand

With auand av, phase shifts for each MPC can be calculated. Then, AAoAs and EAoAs can be calculated according to (3) and (4),

Note that, as a horizontal URA is used,MPCs from above and below cannot be distinguished. But it is able to resolve MPCs with AAoA from 0°to 360°. However, AAoAs calculated through (16) are in the range of 0°to 180°. This problem can be solved by observing the positive/negative signs of uland vl. Then, the calculated AAoAs can be mapped to different quadrants of spherical coordinates. This problem has not been mentioned in other literatures before.

3) Unitary ESPRIT algorithm: Different from the LS-ESPRIT algorithm, Unitary ESPRIT algorithm employs the centrosymmetric property of Rx array and operates in real-valued domain to further reduce the complexity.When Mxand Myare odd, a(ul) and a(vl)can be re-written as

and

respectively, by defining the center of Rx array as the reference antenna. They can be transformed into real-valued by using the following simplest matrix [19],

where 2m+1=M if M is odd and Πmis the exchange matrix. The real-valued steering vectors can be expressed asand. Similar to (2), the steering matrix is stacked into a vector form,

Further considering L MPCs, the steering matrix can be expressed as

By selecting Sub-arrays 1 and 2, the invariance relationship satisfied by A′is [19]

whereandare the real and imaginary parts ofrespectively,and

Similarly, by selecting Sub-arrays 3 and 4,we can get

where Kv1and Kv2are the real and imaginary parts of, respectively, and

Accordingly, the received signal is re-written asLetbeing the L “l(fā)argest” left singular vectors of, whereandare the real and imaginary parts ofis related to A′with an unknown L×L real-valued linear transformation matrix, i.e.,. Then, by substitutingto (23) and (24) we can get

and

where. Therefore, the eigenvaluesof Ψu+jΨvcan be computed by using the joint diagonalization method. Then, uland vlare given as

Finally, the AAoA and EAoA can be calculated as in (15) and (16).

3.2 Joint delay, AAoA, and EAoA estimation methods

In table 2, the preliminary comparison of joint delay, AAoA, and EAoA estimation algorithms are listed. The detailed procedures can be found in the following.

1) JADE MUSIC algorithm: In order to construct a temporal phase shift similar to the spatial phase shift, JADE-based algorithms transform the time domain measurement datainto frequency domain. With Fourier transformation (FT), the temporal phase shift of the l-th MPC is

Table II. Properties of joint AOA and delay estimation algorithms.

where, dfis the frequency interval determined by the maximum and minimum frequency points, i.e.,

By stacking temporal domain samples into the spatial domain, namely stretching the M×K matrix Y(k) into a MK×1 vector,JADE-based algorithms can jointly estimate delay, AAoA, and EAoA. Accordingly, similar to (2), the spatial-temporal manifold can be expressed as

Thus, the frequency domain vectorized Y(k) can be given as

where S and N are the frequency domain vectorized signal and noise.

Considering that JADE MUSIC and JADE ESPRIT can only be used when multiple uncorrelated snapshots are available,Y1can be extended to a MK×Lsmatrixwith Lsdenotes multiple snapshots. It has the matrix structure as shown in figure 5. The received signals contributed by the l-th path at two adjacent data blocks along dimension 1 have phase shift e?jul. Along dimension 2 have phase shift e?jvl, and along dimension 3 have phase shift e?jwlcaused by df.

Fig. 5. Matrix transposition structure of

The sampled spatial-temporal covariance matrix is calculated as

Similar to that of the MUSIC algorithm,haszero eigenvalues theoretically. In practice, these zero eigenvalues are corresponding to (M K?L) minimum eigenvalues of. Eigenvectors corresponding to the L larger eigenvalues form the signal subspace E2,s, and eigenvectors corresponding to the (M K?L) smaller eigenvalues form the noise subspace E2,n. Therefore, JADE-based algorithms overcome the constraint inherent in the subspace-based method, e.g., MUSIC and ESPRIT, and can be used when MPC number is larger than the antenna number as long as L＜MK.

The normalized JADE MUSIC spatial-temporal spectrum is defined as

Here, w is the temporal phase shift caused by delay. Finally, the delays, AAoAs, and EAoAs can be located by finding the local maxima of the spatial-temporal spectrum.

2) JADE ESPRIT algorithm: Similar to the LS-ESPRIT algorithm, JADE ESPRIT algorithm can estimate delay by using the invariance properties of two sub-arrays. The divided two sub-arrays are exactly the same except a constant temporal phase shift.

Based on two frequency do-main sample selection matricesand, the following relationship holds

whereand E2,s6= J6E2,s. Thus, the delay of l-th MPC can be calculated based on wl.

By choosing the selection matrices J1-J4and signal subspaces E2,s1-E2,s4accordingly,auand avcan be calculated similar to (13)and (14), respectively. Regarding the joint estimation of delay, AAoA, and EAoA, joint diagonalization methods proposed in [25] can be referred.

3) SI-JADE algorithm: Instead of assuming multiple snapshots, SI-JADE algorithm employs the spatial and temporal domain smoothing technique to form multiple uncorrelated samples [25]. According to the method proposed in [15], a Hankel matrix can be constructed by left shifting and stacking copies of Y2, which is the FT of Y(k). Assume the selected blocks in spatial and temporal domains are m1and m2, respectively, the reconstructed received data matrix can be expressed as

whereis given as (34) shown in the bottom at this page.

Then, the selection matrices J1-J6can be changed according to the structure of (38)and the corresponding signal subspaces E3,s1-E3,s6can be determined. Finally, Ωu, Ωv,and Ωwcan be calculated similar to (13),(15), and (37).

4) Improved SI-JADE algorithm: Employing the center frequency point as reference, we can re-write (34) in a centrosymmetric form as

The real-valued “temporal manifold” can be acquired as. Similar to(35), the real-valued spatial-temporal manifoldcan be written as

To estimate delays, we employ the shift invariance properties as in (23) and (24), i.e.,

whereKw1and Kw2are the real and imaginary parts ofrespectively. The selection matrix is re-written as

The received signal can be expressed a sB y calculating the signal subspace E3,sof

Then, delays can be calculated as introduced in the Unitary ESPRIT algorithm. Similarly, AAoAs and EAoAs can also be calculated.

5) SAGE algorithm: SAGE algorithm is widely used for channel parameter estimation.It can provide accurate estimation results even in highly correlated scenario. However, the complexity is much higher than other kinds of methods due to the iteration procedure.Complete data (unobservable) Xl(k;Θl) and incomplete data (observable) Y(k) are two key notions in the SAGE algorithm, and the choice of complete data should guarantee the complexity and convergence rate of parameter estimation procedure. As specified in [30],the natural choice of complete data given in(1) leads to a surprisingly simple scheme. Ifis known, the parameters of each path can be estimated by searching parameters that maximize the ML functions.However, as it is unknown in practical, it can be obtained by subtracting the contributions of all the paths except the l-th path from Y(k),

Here, parameters with notation “^” denote the initial assumption or previous estimation of the l′-th path. Equation (43) represents the“E”-expectation step of the SAGE algorithm using parallel interference cancellation method.

Given the complete data acquired above,parameters of MPCs can be derived by searching the values that can maximize the cost functions,

Equations (44)-(47) are referred to the“M”-maximization step. Iteratively carrying out the “E” and “M” steps, a sequence of estimates can be generated until the cost function converges to a stationary point.

where

IV. SIMULATION RESULTS AND MEASUREMENT ANALYSIS

Fig. 6. Angle searching results of MUSIC algorithm.

In this section, we compare the performance of the MUSIC, LS-ESPRIT, and Unitary ESPRIT algorithms to estimate AAoA and EAoA, as well as the performance of JADE ESPRIT, JADE MUSIC, SI-JADE, the improved SI-JADE, and SAGE algorithms to estimate both delay and AoAs in synthetic environments. Estimation results of MUSIC,LS-ESPRIT, Unitary ESPRIT, SI-JADE, the improved SI-JADE, and SAGE algorithms are also compared in 60 GHz real measurement environment.

4.1 Simulation results

1) Performance of joint AAoA and EAoA estimation: As MUSIC algorithm estimates AAoA and EAoA by searching for the spatial spectrum peaks, we assume 4 MPCs to observe the resolution property of MUSIC algorithm. The AAoAs and EAoAs for MPCs 1-4 are set as: {andThe number of sample points is 201, signal-to-noise ratio(SNR) is 20 dB, and Rx array size is set to 4× 4. As shown in figure 6, four MPCs can be isolated with sharp peaks, even though there are slight ambiguities between MPCs 3 and 4.With the increase of SNR and array size, these ambiguities can be suppressed.

With the same simulation settings used in figure 6, the estimated parameters using MUSIC, LS-ESPRIT, Unitary ESPRIT, and SAGE algorithms are shown in figure 7. The MPCs parameters estimated by MUSIC algorithm are determined through finding four local maxima. We can see that all the four algorithms can provide accurate estimation results, and the paring for each MPC is also correct.

In order to compare performance of these algorithms to estimate single MPC with different SNRs and array sizes, we select MPCand Monte Carlo simulation time Mcis set to 500. The selected array sizes are 2× 2 and 4× 4, respectively.Figure 8 shows the root mean square errors(RMSEs) of AoA versus SNR. The RMSE of AoA is calculated as follows

whereandare the estimated AAoA and EAoA for the l-th path. It can be seen that with the increase of SNR and array size, all algorithms show reduced RMSEs of AoA and the RMSEs for SAGE algorithm stay about 0.Under the same array size, MUSIC algorithm shows similar performance at low SNR with LS-ESPRIT and Unitary ESPRIT algorithms,but when SNR is large, MUSIC algorithm shows smaller RMSEs than them. Compare LS-ESPRIT with Unitary ESPRIT algorithm,it can be seen that they have very similar performance at both 2× 2 and 4× 4 scenarios.But Unitary ESPRIT algorithm has less computation complexity.

We assume two MPCs with parameters of MPC 1 given as above and the AAoA of MPC 2 changes from 60°to 75°. Namely, the AAoA difference Δφ of MPC 1 and 2 ranges from 0°to 15°. The differences of EAoAs Δθ are set to 5°and 10°, respectively. Thus,the resolutions of algorithms to separate two MPCs can be compared. As shown in figure 9, with the increase of both AAoA and EAoA differences, all algorithms show improved RMSEs. SAGE algorithm has the lowest RMSEs even with small angle difference. MUSIC has smaller RMSEs than LS-ESPRIT and Unitary ESPRIT algorithms. It can also be observed that the Unitary ESPRIT algorithm shows better resolution and more robust than the LS-ESPRIT algorithm.

Fig. 7. Estimation results of joint AAoA and EAoA estimation algorithms.

Fig. 8. Impacts of SNR and array size on RMSEs of AoA.

2) Performance of joint delay, AAoA,and EAoA estimation: For the joint estimation of delay, AAoA, and EAoA algorithms, one MPC with parameters set asis assumed.The snapshots of JADE MUSIC and JADE ESPRIT algorithms are 20. The selected data blocks for spatial and temporal spatial smoothing of (improved) SI-JADE algorithm are m1=2 and m2=30. The array size is set to 3× 3 and 4× 4, and SNR is ranged from 0 to 35 dB to study the impacts of array size and SNR on the estimation performance of one single path, as shown in figure 10. Figure 10(a) shows the RMSEs of AoA and figure 10(b) shows the RMSEs of delay. The RMSE of delay is calculated as follows whereis the estimated delay for the l-th path.

As can be seen from figure 10(a), all algorithms show decreased RMSEs of AoA with the increase of array size and SNR. Comparing JADE MUSIC with JADE ESPRIT under the same array size condition, RMSEs of JADE MUSIC are larger than JADE ESPRIT at small SNR, however, they decrease rapidly to 0 with the increase of SNR. SI-JADE and the improved SI-JADE algorithms show similar performance, both of them show less satisfactory estimation performance than JADE ESPRIT, which assumes 20 snapshots. SAGE algorithm has the lowest RMSEs of AoA even with small array size.

Fig. 9. Impacts of AAoA and EAoA differences on RMSEs of AoA.

In figure 10(b), RMSEs of delay for JADE MUSIC are less than JADE ESPRIT especially with 4× 4 array. SI-JADE and the improved algorithms show similar results in delay estimation, and they have smaller RMSEs of delay than the JADE ESPRIT algorithm.Still, SAGE algorithm can provide accurate delay estimation results.

The spatial and temporal resolution of algorithms are compared in figure 11 using a 4× 4 array. Here, we assume two MPCs, MPC 1 has the same parameters as used in figure 10.The EAoA of MPC 2 is kept constant to 35°.The AAoA of MPC 2 is changed from 60°to 75°, i.e., AAoA difference between two MPCs ranges from 0°to 15°. The delay differences are set to 5 ns and 10 ns successively. From figure 11(a), most algorithms are not very sensitive to the delay and AAoA differences.SAGE and JADE-MUSIC algorithms can provide very accurate estimation, followed by JADE-ESPRIT algorithm. SI-JADE and the improved algorithm have similar RMSEs,however, the improved version shows more robust performance. From figure 11(b), still JADE-MUSIC and SAGE have the lowest RMSEs. The improved SI-JADE algorithm has lower RMSEs than the original SI-JADE,and it is less sensitive to the delay and AAoA differences. In general, comparing figure 11 with figure 9, the estimation performance is greatly improved with the introduction of delay estimation.

Fig. 10. Impacts of SNR and array size on RMSEs of (a) AoA and (b) delay estimations.

Fig. 11. Impacts of delay and AAoA differences on RMSEs of (a) AoA and (b) delay estimations.

4.2 Measurement analysis

To verify the performance of aforementioned algorithms in real channel measurement data processing, we carry out channel measurements in a typical indoor environment at 60 GHz. The relative location between the Tx and Rx is shown in figure 12. The Tx is configured with a broad beam antenna, which is placed on a positioner (height 2 m) to form a 15× 15 horizontal planar array. At the Rx side, a 25 dBi horn antenna is located on a tripod with 1.6 m height. Detailed system setups can be found in table 3.

In figure 13, the estimated spatial spectrums using MUSIC algorithm are depicted. The assumed MPC number is 200 and the angle steps used to search AAoA and EAoA are both 1°.It shows several distinct spectrum peaks with AAoA centring at certain angles while EAoA distributed from 80°to 90°. However, the spectrum peaks are not sharp enough to locate each closely spaced MPCs. Thus, we can only observe the distribution of several clusters.In order to get estimated AoAs of MPCs, the normalized spectrum values are sorted into descending order. Then, the estimated AoAs of MPCs can be roughly located by selecting angles that correspond to the 200 largest values,as shown in figure 14.

Fig. 12. Relative location of Tx array and Rx.

Fig. 13. Estimated spectrum using MUSIC algorithm.

Table III. System setups.

Fig. 14. Comparison of joint AAoA and EAoA estimation results.

Fig. 15. Comparison of joint delay, AAoA, and EAoA estimation results.

In figure 14, the joint AAoA and EAoA estimation results using MUSIC, LS-ESPRIT,Unitary ESPRIT, and SAGE algorithms are shown. For clarity, the MPC numbers assumed for each algorithm are 200, 20, 20, and 80,respectively. We can see that even though the estimated results are not exactly the same,they have most MPCs with EAoAs ranging from 60°to 90°. Especially, results of MUSIC and SAGE algorithms are consistent at six clusters. According to the real measurement setups as in figure 12, they might be formed by the line-of-sight (LoS) and single bounce components. However, only several MPCs estimated by LS-ESPRIT and Unitary ESPRIT are overlapped with the other two algorithms.

Figure 15 shows the joint delay, AAoA,and EAoA estimation results using SI-JADE,improved SI-JADE, and SAGE algorithms.The estimated MPC numbers for three algorithms are 20, 20, and 80, respectively. The estimated MPCs exhibit as clusters with similar parameters. One cluster with delay around 20 ns is contributed by the LoS components.Refer to the distance between Tx and Rx, the correctness of delay estimation results can be verified. It can be seen that SI-JADE and the improved SI-JADE both can locate several MPCs in each cluster. Within each cluster,more MPCs can be extracted than that shown in figure 14. This may be contributed by the temporal-spatial smoothing technique and the inclusion of delay estimation. However, there are also some MPCs calculated by SI-JADE and improved SI-JADE algorithms that have not been extracted by SAGE algorithm, they might be false MPCs induced by the highly coherent signal, or weak MPCs that have not been detected by SAGE algorithm.

V. FUTURE DIRECTIONS

5.1 Improvements based on the traditional methods

Algorithms mentioned in Section III have their own constraints. For example, subspace-based methods very much rely on the full rank(non-singular) of spatial covariance matrix and ML-based methods are intrinsically time consuming determined by iteration procedure.Considering their pros and cons, hybrid methods can be proposed to exploit the advantage of these algorithms. In addition, the system model is usually assumed to be used in a narrowband scenario. However, with the deployment of new technology such as mmWave, ultra-wideband should be considered. Thus, the phase shift caused by the time delay among antenna elements should also be taken into consideration and the existing algorithms can be further improved.

5.2 Application of machine learning/big data algorithms

Considering the large amount of wireless channel measurement data, algorithms that can estimate channel parameters in an efficient way and be used to the real-time data processing are required. It is well known that machine learning or big data related algorithms are capable of analyzing nonlinear properties of wireless channels and neural network has excellent learning and prediction ability[36]. Therefore, we can incorporate suitable algorithms into channel measurement data processing. In channel characteristics analysis and channel modeling, clustering method such as Kernel-power-density (KPD) algorithm was studied to group MPCs [37]. In [38], Gaussian mixture model (GMM) was used to determine channel parameter distribution and convolutional neural network was used to identify different channels automatically. To estimate time-variant parameters, techniques that can predict the tendency of parameter changes should also be introduced.

5.3 New channel characteristics

With the employment of new technologies,algorithms that can cope with the new channel characteristics are necessary. In massive MIMO channel, spherical wavefront instead of plane wavefront should be considered [39],as well as the non-stationary phenomenon[40]. MmWave frequency bands can provide very large bandwidth and thus higher time resolution than conventional frequency bands.Non-stationarity in frequency domain should also be taken into consideration. In high speed train (HST) and vehicle to vehicle (V2V) communications, the channel show non-stationarity in time domain. Therefore, algorithms that can provide higher time resolution, be used in near- field channel measurement data processing, able to estimate time-variant parameters,and cope with the temporal-spatial-frequency non-stationarity are essential.

VI. CONCLUSIONS

In this paper, we have introduced some commonly used wireless channel parameter estimation methods. MUSIC, LS-ESPRIT, and Unitary ESPRIT algorithms can be used to extract 3D spatial information. JADE MUSIC,JADE ESPRIT, SI-JADE, and SAGE algorithms can be employed to jointly estimate both temporal and 3D spatial information. By incorporating with Unitary ESPRIT algorithm,we have further improved the SI-JADE algorithm to reduce computational complexity.Simulation results have shown that for AAoA and EAoA estimation, MUSIC algorithm has better resolution but larger computation burden than ESPRIT-type algorithms, while Unitary ESPRIT algorithm is more robust than LS-ESPRIT algorithm. For joint extraction of delay, AAoA, and EAoA, SAGE algorithm can provide the most accurate results, followed by JADE MUSIC algorithm, while the conventional SI-JADE algorithm is the worst.The improved SI-JADE algorithm can provide comparable accuracy, but is more robust and less complex than the conventional SI-JADE algorithm. For 60 GHz channel measurement data processing, results have shown that MUSIC algorithm can locate several main clusters of MPCs in AAoA and EAoA estimation. The(improved) SI-JADE algorithm can distinguish more MPCs within each cluster than the LS-ESPRIT and Unitary ESPRIT algorithms in joint delay, AAoA, and EAoA estimation.They have similar estimation results to SAGE algorithm. Also, the feasibilities of complexity reducion methods, such as MUSIC, SI-JADE,and improved SI-JADE algorithms, have been verified in 60 GHz real channel measurement data processing. Finally, we have further discussed a few future research directions.

ACKNOWLEDGMENT

The authors would like to acknowledge the support from the Natural Science Foundation of China (Grant No. 61210002, 61371110),EU H2020 ITN 5G Wireless project (No.641985), EU H2020 RISE TESTBED project(No. 734325), and EPSRC TOUCAN project(Grant No. EP/L020009/1).

[1] C.-X. Wang, F. Haider, X. Gao, et al, “Cellular Architecture and Key Technologies for 5G Wireless Communication Networks,” IEEE Commun. Mag.,vol. 52, no. 2, 2014, pp. 122-130.

[2] G. Q. Mao, “Book Review-5G Green Mobile Communication Networks,” China Commun., vol.13, no. 12, 2016, pp. 297-298.

[3] Z. N. Sun, X. Xu, X. M. Chai, et al, “The New Architecture with Time-Spatial Consistency for 5G Networkss,” China Commun., vol. 13, no. 1, 2016,pp. 68-79.

[4] A. L. Swindlehurst, E. Ayanoglu, P. Heydari, et al, “Millimeter-Wave Massive MIMO: The Next Wireless Revolution?” IEEE Commun. Mag., vol.52, no. 9, 2014, pp. 56-62.

[5] T. S. Rappaport, J. N. Murdock, and F. Gutierrez,“State of the Art in 60-GHz Integrated Circuits and Systems for Wireless Communications,”P(pán)roc. IEEE, vol. 99, no. 8, 2011, pp. 1390-1436.

[6] S. B. Wu, C.-X. Wang, M. Alwakeel, et al, “A Non-stationary 3-D Wideband Twin-cluster Model for 5G Massive MIMO Channels,” IEEE J. Sel. Areas Commun., vol. 32, no. 6, 2014, pp.1207-1218.

[7] S. B. Wu, C.-X. Wang, H. Haas, et al, “A Non-stationary Wideband Channel Model for Massive MIMO Communication Systems,” IEEE Trans.Wireless Commun., vol. 14, no. 3, 2015, pp. 1434-1446.

[8] C.-X. Wang, S. B. Wu, L. Bai, et al, “Recent Advances and Future Challenges for Massive MIMO Channel Measurements and Models,” Sci.China Inf. Sci., vol. 59, no. 2, 2016, doi: 10.1007/s11432-015-5517-1.

[9] B. H. Fleury, M. Tschudin, R. Heddergott, et al,“Channel Parameter Estimation in Mobile Radio Environments using the SAGE Algorithm,” IEEE J. Sel. Areas Commun., vol. 17, no. 3, 1999, pp.434-450.

[10] H. Krim and M. Viberg, “Two Decades of Array Signal Processing Research: The Parametric Approach,” IEEE Sig. Proc. Mag., vol. 13, no. 4,1996, pp. 67-94.

[11] B. D. Van Veen and K. M. Buckley, “Beamforming:A Versatile Approach to Spatial Filtering,” IEEE ASSP Mag., vol. 5, no. 2, 1988, pp. 4-24.

[12] J. Capon, “High-Resolution Frequency-Wavenumber Spectrum Analysis,” Proc. IEEE, vol. 57,no. 8, 1969, pp. 1408-1418.

[13] R. O. Schmidt, “Multiple Emitter Location and Signal Parameter Estimation,” IEEE Trans. Antennas Propag., vol. 34, no. 3, 1986, pp. 276-280.

[14] R. Roy, A. Paulraj, and T. Kailath, “ESPRIT-A Subspace Rotation Approach to Estimation of Parameters of Cisoids in Noise,” IEEE Trans.Acoust., Speech, Sig. Proc., vol. 34, no. 5, 1986,pp. 1340-1342.

[15] M. Haardt and J. A. Nossek, “Unitary ESPRIT:How to Obtain Increased Estimation Accuracy with a Reduced Computational Burden,” IEEE Trans. Sig. Proc., vol. 43, no. 5, 1995, pp. 1232-1242.

[16] P. Strobach, “Total Least Squares Phased Averaging and 3-D ESPRIT for Joint Azimuth-Elevation-Carrier Estimation,” IEEE Trans. Sig. Proc.,vol. 49, no. 1, 2001, pp. 54-62.

[17] R. Roy and T. Kailath, “ESPRIT-Estimation of Signal Parameters via Rotational Invariance Techniques,” IEEE Trans. Acoust., Speech, Sig. Proc.,vol. 37, no. 7, 1989, pp. 984-995.

[18] A. L. Swindlehurst and T. Kailath, “Azimuth/Elevation Direction Finding Using Regular Array Geometries,” IEEE Trans. Aerosp. Electron. Syst.,vol. 29, no. 1, 1993, pp. 145-156.

[19] M. D. Zoltowski, M. Haardt, and C. P. Mathews,“Closed-Form 2-D Angle Estimation with Rectangular Arrays in Element Space or Beamspace via Unitary ESPRIT,” IEEE Trans. Sig. Proc., vol. 4,no. 2, 1996, pp. 316-328.

[20] Y. Bresler and A. Macovski, “Exact Maximum Likelihood Parameter Estimation of Superimposed Exponential Signals in Noise,” IEEE Trans.Acoust., Speech, Signal Processing, vol. 34, no. 1,1986, pp. 1081-1089.

[21] Y. Y. Wang, J. T. Chen, and W. H. Fang, “TST-MUSIC for Joint DOA-Delay Estimation,” IEEE Trans.Sig. Proc., vol. 49, no. 4, 2001, pp. 721-729.

[22] I. Jaafar, R. Amara, H. Boujemaa, et al, “Joint Angle and Delay Estimation of Point Sources,”P(pán)roc. IEEE ICECS'05, Gammarth, Tunisia, Dec.2005, pp. 1-5.

[23] A. J. Van Der Veen, M. C. Vanderveen, and A.Paulraj, “Joint Angle and Delay Estimation(JADE) Using Shift-Invariance Techniques,” IEEE Sig. Proc. Lett., vol. 4, no. 5, 1997, pp. 142-145.

[24] M. C. Vanderveen, A. J. Van Der Veen, and A.Paulraj, “Estimation of Multipath Parameters in Wireless Communications,” IEEE Trans. Sig.Proc., vol. 46, no. 3, 1998, pp. 682-690.

[25] A. J. Van Der Veen, M. C. Vanderveen, and A.Paulraj, “Joint Angle and Delay Estimation Using Shift-Invariance Techniques,” IEEE Trans. Sig.Proc., vol. 46, no. 2, 1998, pp. 405-418.

[26] M. Feder and E. Weinstein, “Parameter Estimation of Superimposed Signals Using the EM Algorithm,” IEEE Trans. Acoust., Speech, Sig. Proc.,vol. 36, no. 4, 1988, pp. 477-489.

[27] B. H. Fleury, R. Heddergott, and M. Tschudin,“Wideband Angle of Arrival Estimation using the SAGE Algorithm,” Proc. IEEE ISSSTA'96,Mainz, Germany, Sept. 1996, pp. 79-85.

[28] R. Feng, J. Huang, J. Sun, et al, “A Novel 3D Frequency Domain SAGE Algorithm with Applications to Parameter Estimation in MmWave Massive MIMO Indoor Channels,” Sci. China Inf.Sci., vol. 60, no. 8, 2017, doi: 10.1007/s11432-017-9139-4.

[29] X. F. Yin, Y. Y. He, Z. N. Song, et al, “A Sliding-Correlator-Based SAGE Algorithm for Mmwave Wideband Channel Parameter Estimation,”P(pán)roc. IEEE EuCAP'14, Hague, Netherlands, Apr.2014, pp. 625-629.

[30] K. I. Pedersen, B. H. Fleury, and P. E. Mogensen,“High Resolution of Electromagnetic Waves in Time-Varying Radio Channels,” Proc. IEEE PIMRC'97, Helsinki, Finland, Sept. 1997, pp.650-654.

[31] A. Richter, M. Landmann, and R. S. Thoma,“RIMAX-A Flexible Algorithm for Channel Parameter Estimation from Channel Sounding Measurements,” COST273 TD(04)045, Athens,Greece, Jan. 2004, http://www.lx.it.pt/cost273/.

[32] J. Huang, C.-X. Wang, and R. Feng, et al,“Multi-frequency MmWave Massive MIMO Channel Measurements and Characterization for 5G Wireless Communication Systems,” IEEE J. Sel. Areas Commun., vol. 35, no. 7, 2017, pp.1591-1605.

[33] R. Feng, Y. Liu, J. Huang, et al, “Comparison of MUSIC, Unitary ESPRIT, and SAGE Algorithms for Estimating 3D Angles in Wireless Channels,”P(pán)roc. IEEE ICCC'17, Qingdao, China, Oct. 2017.

[34] X. F. Yin, L. X. OuYang, and H. W. Wang, “Performance Comparison of SAGE and MUSIC for Channel Estimation in Direction-Scan Measurements,” IEEE Access, vol. 4, no. 1, 2016, pp.1163-1174.

[35] T. Wang, B. Ai, R. S. He, et al, “Two-Dimension Direction-of-Arrival Estimation for Massive MIMO Systems,” IEEE Access, vol. 3, no. 1, 2015,pp. 2122-2128.

[36] J. Huang, C.-X. Wang, L. Bai, et al, “A Big Data Enabled Channel Model for 5G Wireless Communication Systems,” IEEE Trans. Big Data,2017, submitted for publication.

[37] R. S. He, Q. Li, and B. Ai, “An Automatic Clustering Algorithm for Multipath Components Based on Kernel-power-density,” Proc. IEEE WCNC'17,San Francisco, USA, Mar. 2017, pp. 1-6.

[38] H. H. Li, Y. Z. Li, S. D. Zhou, et al, “Wireless Channel Feature Extraction via GMM and CNN in the Tomographic Channel Model,” J. Commun. Inf.Netw., vol. 2, no. 1, 2017, pp. 41-51.

[39] X. F. Yin, S. Wang, N. Zhang, et al, “Scatterer Localization using Large-Scale Antenna Arrays Based on a Spherical Wave-Front Parametric Model,” IEEE Trans. Wireless Commun., vol. 16,no. 10, 2017, pp. 6543-6556.

[40] R. Feng, J. Huang, J. Sun, et al, “Spatial Cross-Correlation Properties of MmWave Massive MIMO Channels,” Proc. IEEE ICCC'17, Qingdao, China, Oct. 2017.