A Complexity-Efficient High Resolution Propagation Parameter Estimation Algorithm for Ultra-Wideband Large-Scale Uniform Circular Array

TL;DR

This paper introduces a novel, complexity-efficient high resolution parameter estimation algorithm for mm-wave channels with large-scale uniform circular arrays, improving accuracy and realism in propagation parameter estimation.

Contribution

It presents a new HRPE algorithm that leverages delay trajectories and phase mode excitation, considering spatial non-stationarity, with validated results from real mm-wave measurements.

Findings

Reduces computational complexity significantly.

Achieves high-resolution estimation of propagation parameters.

Validates effectiveness through real-world mm-wave measurements.

Abstract

Millimeter wave (mm-wave) communication with large-scale antenna array configuration is seen as the key enabler of the next generation communication systems. Accurate knowledge of the mm-wave propagation channels is fundamental and essential. In this contribution, a novel complexity-efficient high resolution parameter estimation (HRPE) algorithm is proposed for the mm-wave channel with large-scale uniform circular array (UCA) applied. The proposed algorithm is able to obtain the high-resolution estimation results of the spherical channel propagation parameters. The prior channel information in the delay domain, i.e., the delay trajectories of individual propagation paths observed across the array elements, is exploited, by combining the high-resolution estimation principle and the phase mode excitation technique. Fast initializations, effective interference cancellations and reduced searching spaces achieved by the proposed schemes significantly decrease the algorithm complexity. Furthermore, the channel spatial non-stationarity in path gain across the array elements is considered for the first time in the literature for propagation parameter estimation, which is beneficial to obtain more realistic results as well as to decrease the complexity. A mm-wave measurement campaign at the frequency band of 28-30 GHz using a large-scale UCA is exploited to demonstrate and validate the proposed HRPE algorithm.

Figures

Tables

Table 1. Table I: Measurement specifications applied in the measurement campaign

Measurement specifications
Tx and Rx Heights	0.84 m	Frequency range	28-30 GHz
Radius of bicoinal array	0.5 m	Frequency points	750
Number of virtual Rxs	720	LoS distance	5 m

Table 2. Table II: Specifications of antennas applied in the measurement campaigns

Antenna specifications
Antenna type	Biconical	Reference horn
Operating frequency range	2-30 GHz	26.4-40.1 GHz
Gain	6 dBi	19 dBi
HPBW in azimuth	Omni	20∘
Polarization	Vertical	Vertical

Equations

$$H_{\ell }(f)={\alpha }_{\ell }{e}^{-j2\pi f{\tau }_{\ell }}$$

$$d_{p,\ell }=\sqrt{{d_{\ell }}^2+r^2-2r{d}_{\ell }\sin {\theta }_{\ell }\cos ({\varphi }_{\ell }-{\varphi }_p)}$$

$$H_{p,\ell }(f)$$

$$\Delta d_{p,\ell }$$

$$Y(p,f;\Theta )=\sum _{\ell =1}^{L}H(p,f;{\Theta }_{\ell })+n(p,f)$$

$$H_{p,\ell }(f)$$

$${\alpha }_p^{\ell }$$

$${\tau }_p^{\ell }$$

$$\Delta d_{p,\ell }$$

$$H_{p,\ell }(f)$$

$$=H_{\ell }(f)\sum _{n=-\infty }^{+\infty }{j}^n{J}_n(2\pi f\frac{r}{c})e^{jn({\varphi }_{\ell }-{\varphi }_p)}$$

$$e^{j\beta \cos \gamma }=\sum _{n=-\infty }^{+\infty }{j}^n{J}_n(\beta )e^{jn\gamma }$$

$$H_{m,\ell }(f)$$

$$=H_{\ell }(f)\sum _{n=-\infty }^{+\infty }{J}_n(2\pi f\frac{r}{c})W_m(f)e^{-jn{\varphi }_{\ell }}\sum _{p=0}^{P-1}\frac{e^{j(m-n){\varphi }_p}}{P}$$

$$\approx H_{\ell }(f)J_m(2\pi f\frac{r}{c})W_m(f)e^{-jm{\varphi }_{\ell }}$$

$$=H_{\ell }(f)e^{-jm{\varphi }_{\ell }}$$

$$\frac{1}{P}\sum _{p=0}^{P-1}{e}^{j(m-n){\varphi }_p}=1$$

$$W_m(f)$$

$$Y(m,f;\Theta )$$

$$G_m(f)=\frac{2}{J_m(2\pi f\frac{r}{c})+J_m^{{}^{\prime }}(2\pi f\frac{r}{c})}$$

$$T_{{\theta }_{\ell }}={\hat{\tau }}_{\ell }-\frac{r}{c}\sin {\theta }_{\ell }\cos ({\hat{\varphi }}_{\ell }-{\varphi }_p)$$

$${\hat{\theta }}_{\ell }=\arg \underset{{\theta }_{\ell }}{\max }C(A_{{\theta }_{\ell }})$$

$$({\hat{\varphi }}_{\ell },{\hat{\theta }}_{\ell },{\hat{d}}_{\ell })=\arg \underset{\varphi ,\theta ,d}{\max }\text{vec}\{\hat{H}(p,f_k;{\Theta }_{\ell }){\}}^T\text{vec}\{W^{\ast }(p,f_k)\}$$

$${\hat{\tau }}_{\ell }=\arg \underset{\tau }{\max }\text{vec}\{\hat{H}(p,f;{\Theta }_{\ell }){\}}^T\text{vec}\{W^{\ast }(p,f)\}$$

$${\hat{\alpha }}_{\ell }=\frac{1}{CK}\text{vec}\{\hat{H}(p,f;{\Theta }_{\ell }){\}}^T\text{vec}\{W^{\ast }(p,f)\}$$

Full text

Index terms— Millimeter-wave, ultra-wideband, large-scale uniform circular array, spatial non-stationarity and channel parameter estimation.

I Introduction

While the next generation communication system (5G) is still on its ascendant, millimeter wave (mm-wave) communication has been seen as the key enabling component due to the vast amount of available spectrum [1, 2, 3]. However, the air attenuation, small antenna aperture etc. would result in severe power loss [4] for the mm-wave propagations, compared to the sub-6 GHz frequency bands with rich multipath components (MPCs) [5]. Nevertheless, the mm-wavelength makes it practical to pack massive antennas [6] in a small area to form large scale antenna arrays. Beamforming [7] and beam-selection [2] techniques can be exploited to compensate the power loss and to enhance the spectrum efficiency through spatial multiplexing. Furthermore, applications such as localization, tracking and surrounding environment reconstruction [8] are promising, e.g., to help vision-disabled people.

To enable the advanced 5G mm-wave techniques and applications, accurate and realistic channel models are fundamental and essential. The establishment of effective mm-wave channel models relies on the comprehensive channel measurements and the channel characteristics extracted from the measured data. In the mm-wave propagations with large scale antenna array configurations, the assumptions applied for the previous generation communication systems (e.g. the Long Term Evolution system) are violated. i) The narrowband assumption [9] is invalid due to the ultra-wide system bandwidth up to several GHz; ii) The two dimensional (2D) propagation assumption [10] was usually assumed in pervious communication systems since the base station was high above, and the propagation distance was large. However, in the mm-wave frequency bands, the propagation distance is generally limited due to the high power loss. It is not practical anymore to assume the 2D propagation. That is, the elevation angels of MPCs [11] should be considered; iii) The mm-wavelength and large array aperture can result in the Fraunhofer far-field distance up to tens or hundreds of meters, which is significantly increased compared to the previous systems. This necessitats the spherical wave propagation [12] model; iv) It is possible that the gain (not just phases) of propagation paths evolve [13, 14] across the array elements. The so-called spatial non-stationarity in path gainmust be considered as well. Thus the mm-wave channel models must provide realistic channel characteristics in multiple parameter domains, i.e., in delay, azimuth, elevation, source distance and amplitude domains. Consequently, accurate, comprehensive yet still complexity-efficient channel parameter estimation algorithms are vital to rapidly promote the 5G progress.

Among the many mm-wave channel spatial profile measurements, the most frequently applied measurement approach is to either mechanically rotate a horn antenna or mechanically move an omnidirectional antenna to form virtual arrays [15]. The vector network analyzer (VNA) is usually used due to its flexibility and easy-configuration. Various estimation algorithms have been proposed for the mm-wave channel estimation based on the empirically measured data. Basically, these algorithms can be classified into four categories, i.e., spectra based approaches, subspace based approaches, sparsity recovery approaches and maximum likelihood approaches as follows. i) For the horn antenna measurements, the measured channels at different steering directions constructed the joint delay-angle power spectra, e.g., as done in [16, 17, 18, 19]. The complexity is of the algorithm is low, while the resulted channel characteristics show dependency on the horn antenna pattern. In [20], the fast fourier transform (FFT) was applied for the extremely large antenna array to distinguishing MPCs. A frequency invariant beamformer was proposed in [21] and extended in [22] for a uniform circular array (UCA) to obtain the delays and azimuths of MPCs in the 3D propagation scenarios. ii) In [23], the multiple signal classification (MUSIC) principle was applied to jointly estimating the delay and azimuth in a 2D mm-wave propagation scenario. Unitary estimation of signal parameter via rotational invariance techniques (ESPRIT) was also exploited in [24] for the 3D mm-wave propagation scenarios with plane wave assumption. iii) Based on the sparsity assumption of mm-wave channels, the authors in [25, 26, 27, 28] recovered the channels and obtained the channel estimations by exploiting the compressive sensing or optimization techniques. iv) The Space Alternating Generalized Expectation-maximization (SAGE) algorithm was applied for the mm-wave channel estimation in [29, 30, 14, 31, 12]. Both the 3D plane-wave propagation scenarios [29, 30, 14] and the 3D spherical propagation scenarios [31, 12] were concerned. However, the fundamental assumption of SAGE is that the received signal contains multiple components that follow orthogonal stochastic measures (OSM) [32, 12]. The violation [33] of the OSM assumption caused by 3D spherical propagation could deteriorate the SAGE performance significantly, e.g., leading to the non-convergence. A maximum likelihood estimator (MLE) based on the expectation-maximization (EM) principle was proposed in [34] aimed for the 3D spherical propagation scenario. The complexity is still considerably high with the 4D parameter searching, though a coarse-to-fine search strategy was applied.

Although efforts have been made to characterize the mm-wave propagation channels, a high resolution propagation parameter estimation (HRPE) algorithm which is capable to obtain the high-resolution estimation results of all the spherical propagation parameters (i.e. azimuths, elevations, delays, source distances and amplitudes of MPCs) yet with low computation complexity is still missing in the literature. The existing algorithms are either deficient in obtaining the high-resolution estimation results of all propagation parameters or with fatal computation complexity especially when considering the ultra-wide bandwidth and large-scale array aperture in mm-wave propagation. Moreover, as far as we are concerned, the realistic spatial non-stationarity across the array elements have never been considered for propagation parameter estimation, although it has been widely observed in the mm-wave channel modeling works in the literature, leading non-realistic estimation results and increasing complexity. To fill the above gaps, a complexity-efficient HRPE algorithm is proposed in this paper for the ultra-wideband large-scale UCA. The UCA is considered as it is capable to uniformly cover the whole 360∘ of azimuth compared to the linear arrays and support for phase mode operations. The main contributions and novelties of this paper include:

• •

The proposed algorithm exploits the high resolution estimation technique in the delay domain, the UCA phase mode excitation technique and the maximum-likelihood estimation principle. Fast initializations, effective interference cancellations and reduced searching spaces substantially decrease the computation complexity while maintaining the ability to gain the whole high-resolution estimation results of the 3D spherical mm-wave propagation channel parameters.

• •

The UCA phase mode excitation technique applied in the 3D spherical mm-wave propagation scenarios is investigated through simulations, in terms of its ability to obtain the fast initializations of delays and azimuths.

• •

Individual propagation paths are firstly identified in the delay domain, so that the fast initializations and effective interference cancellations can be achieved for the spherical propagation parameter estimation. Furthermore, the spatial non-stationarity observed across the array elements is considered for the first time in the literature for the propagation parameter estimation, resulting in more realistic estimation results and significantly reduced complexity.

• •

An ultra-wideband mm-wave channel measurement at the frequency band of 28-30 GHz with large-scale (radius of 0.5 m and 720 elements) UCA applied was conducted. The comparisons among the reference channel measured by using the horn antenna, the MLE (or EM) results obtained by applying 4D parameter searchings [34] and the results obtained by using the proposed HRPE algorithm validate the algorithm performance.

The rest of the paper is structured as follows. Sect. II elaborates the signal model. Sect. III describes the measurement campaign and the observed channel characteristics. Sect. IV elaborates the proposed HRPE algorithm. The algorithm validation and the remarks are included in Sect. V. Finally, conclusions are given in Sect. VI.

II Signal Model for The Uniform Circular Array

As illustrated in Fig. 1, the UCA has $P$ isotropic receiver antenna elements (Rxs) uniformly arranged on its perimeter with radius of $r$, i.e., the azimuth of the $p$th Rx is $\phi_{p}=\frac{2\pi p}{P},p=[0,\cdots,P-1]$.

With a large UCA aperture, the spherical wave signal model has to be assumed due to the violation of plane wave assumption. In the underlying channel model, a finite number of $L$ spherical waves are assumed to impinge into the UCA. Although an arbitrary reference point can be selected, we consider the UCA center as the reference point for the model conciseness. The frequency response $H_{\ell}(f)$ at the UCA center, which is contributed by the $\ell$th path (or wave) with parameter set $\Theta_{\ell}=[\tau_{\ell},\phi_{\ell},\theta_{\ell},d_{\ell},\alpha_{\ell}]$, reads

[TABLE]

where $f=[f_{1},\cdots,f_{K}]$ is the vector containing the $K$ frequency points considered, $\Theta_{\ell}$ is the vector containing the parameters of the $\ell$th path, where $\tau_{\ell}$ is the propagation delay, $\phi_{\ell}$ and $\theta_{\ell}$ represent the azimuth and elevation angles, respectively, $d_{\ell}$ is the propagation distance between the UCA center and the last source point during the propagation route, and $\alpha_{\ell}$ denotes the complex amplitude. The propagation distance between the $p$th Rx and the last source point is then calculated as

[TABLE]

Due to the propagation distance difference between the $p$th Rx and the UCA center, the frequency response $H_{p,\ell}(f)$ contributed by the $\ell$th path at the $p$th Rx consequently reads

[TABLE]

where

[TABLE]

with $c$ the light speed, and $\frac{d_{\ell}}{d_{p,\ell}}$ denotes the amplitude change factor according to the Friis transmission equation. The array output contributed by the $\ell$th path then reads $H(p,f;\Theta_{\ell})=[H_{0,\ell},\cdots,H_{P-1,\ell}]^{T}$, where $[\cdot]^{T}$ denotes the transpose operation of the argument. $H(p,f;\Theta_{\ell})$ is a complex-valued matrix with dimension of $P\times K$, i.e., $H(p,f;\Theta_{\ell})\in\mathbb{C}^{P\times K}$.

The array output contributed by all the $L$ paths is formatted as

[TABLE]

where $n(p,f)\in\mathbb{C}^{P\times K}$ is the zero-mean white Gaussian noise with variance of $\sigma^{2}$, and $\Theta=[\Theta_{1},\cdots,\Theta_{L}]$ is the vector containing all the channel parameters. The objective of various estimation algorithms is to obtain the estimation results of $\Theta$. In the sequel, a novel HRPE algorithm is proposed for the ultra-wideband large-scale UCA, which is demonstrated based on a measurement campaign where the realistic channel characteristics are well presented. The proposed algorithm is capable to gain the estimation of $\Theta$.

III Measurement campaign and channel observations

III-A Measurement campaign

In this section, a recently conducted measurement campaign for a UCA is introduced. The measurement scenario and the measurement setup including the transmitter (Tx), Rx, array aperture, frequency range, bandwidth etc. are described.

Fig. 2 illustrates the top-view sketch of the basement where the measurement campaign was conducted. The floor space was around $7.7*7.9$ m2 with few objects including one metallic stairs and one metallic heater. Fig. 2 illustrates a photo taken during the measurement. A VNA was used in the channel measurement. An omnidirectional biconical antenna was exploited as the Tx and was put on the trolley with the VNA onboard. The height of the Tx to the ground was 0.84 m. An identical biconical antenna was fixed on the turntable as the Rx. A virtual UCA was formed by rotating the Rx clockwise with radius of 0.5 m with 720 steps, i.e. $r=0.5$ and $P=720$.111The Fraunhofer distance, i.e. $\frac{2\cdot(2r)^{2}}{\lambda}$, is calculated as 200 m. Plane wave propagation is considered valid only when the source distance $d_{\ell}$ is far greater than it. The distance space between two neighboring Rx positions was 4.4 mm, which is less than the half wavelength (5 mm) at 30 GHz and hence avoids the spatial aliasing [35]. The Rx height was kept the same with that of the Tx, and the starting point of the virtual UCA was marked as the star in Fig. 2, which means that the azimuth and elevation angles of the line of sight (LoS) path were 180∘ and 90∘, respectively. Moreover, the distance between the Tx and the UCA center was 5 m.

The frequency range $f$ was set from 28 GHz to 30 GHz with 750 frequency sweeping points, i.e, $K=750$, in the VNA. The frequency step was 2.7 MHz which corresponds to a maximal observable propagation distance of 111 m. It is sufficient to record the propagation paths confined in the basement. The diameter (1 m) of the UCA is around 6.7 times the intrinsic distance resolution (15 cm) determined by the system bandwidth of 2 GHz. Therefore, the narrowband assumption [9] or the so-called small-scale characterization [36] is invalid in the ultra-wideband measurement campaign. Moreover, the Fraunhofer far-filed distance calculated for the UCA is around 200 m. Thus the near-filed (i.e. spherical wave) propagation must be taken into account.

Two different scenarios, i.e., LoS scenario and Obstructed-LoS (OLoS) scenario, were considered during the measurement. In the OLoS scenario, a blackboard with a metallic substrate and dimension of $1.2*1.2$ m2 was placed in between the Tx and the UCA to block the paths in LoS direction. Furthermore, to obtain references for the measured channels in the two scenarios, we also exploited the horn antenna to replace the Rx biconical antenna to repeat the measurements. It is worthy noting that except that the horn antenna was fixed at the UCA center, i.e. $r=0$, the other settings and measurement procedure were kept the same for the reference measurements. The specifications of the measurement campaign and of the antennas applied in the measurement campaign are included in Table I and Table II, respectively. The presented antenna gains and half-power beamwidths in Table II are the values evaluated in the considered frequency range 28-30 GHz. Readers may refer to [37] for the detailed description of the measurement system including the calibration, etc.

III-B Channel observations

In this section, we present the measurement results to shed lights on the physical propagation mechanisms for the ultra-wideband large-scale array system.

By applying the inverse discrete Fourier transform (IDFT) to the array output $Y(p,f)$ with respect to $f$, the channel impulse responses (CIRs) $h(p,\tau)$ across the $P$ antennas can be obtained. Here, we omit $\Theta$ for a compact notation of $Y(p,f;\Theta)$. Fig. 3 and Fig. 3 illustrate the measured concatenated power delay profiles (CPDPs), i.e., $\left|h(p,\tau)\right|^{2}$ for the LoS scenario and OLoS scenario, respectively. The further delay range is not shown as no paths are present with the selected dynamic range. By observing and comparing Fig. 3 and Fig. 3, we have the following observations.

• •

Finite trajectories exist in both Fig. 3 and Fig. 3. The delay shifts of one path across virtual Rxs can be clearly observed. However, limited by the IDFT operation, sidelobes along these trajectories are also obvious, e.g., along the LoS path as illustrated in Fig. 3.

• •

The shapes of trajectories vary due to their different angles-of-arrivals (AoAs). For example, the LoS path in Fig. 3 is dominant with a “(”-alike trajectory. This is consistent with the measurement setup that the virtual UCA elements started at the furthest position to the Tx. It can be expected that these trajectories become straight lines in the narrowband system.

• •

The spatial non-stationarity of the trajectories across the array elements can be clearly observed. That is, one propagation path may present “birth-death” behaviour (or break) across the elements when the array aperture is large. This has also been observed in the other works, e.g. in [13, 14]. However, the non-stationarity has not been considered in the literature regarding the channel parameter estimation to our best knowledge.

• •

In the OLoS scenario as illustrated in Fig. 3, the LoS power attenuates significantly due to the blockage of the blackboard with a metallic substrate. Furthermore, two trajectories disappear in Fig. 3. They are marked as “trajectory 1” and “trajectory 2” in Fig. 3, respectively. It can be inferred that trajectory 1 was contributed by the “wall A” as illustrated in Fig. 2, since its trajectory shape is the same with that of LoS path, and the delay difference is consistent with the measurement setup. Similarly, trajectory 2 contributed by the wall B with an opposite trajectory shape also disappears.

The above observations demonstrate the fact that the parameters (e.g. delays and azimuths) of propagation paths actually can be roughly obtained by observing the CPDPs. This motivates us to come up with the idea that it is possible to distinguish the individual paths by identify their trajectories in the CPDPs figure, so that the channel parameter $\Theta$ can be estimated. However, the IDFT sidelobes blur these paths, resulting in the interference among paths and limiting the resolution. Furthermore, the trajectory shapes of different propagation paths change with respect to the AoAs, which hinders the trajectory identification. To cope with these issues, a novel algorithm is proposed, and its high-resolution ability and low computation cost are also demonstrated.

IV The proposed algorithm

The proposed algorithm basically include three parts, i.e. obtaining the high resolution estimation results of element-wise delays and amplitudes for the UCA, trajectory identification using the phase-mode excitation technique and estimating $\Theta_{\ell}$ according to the maximum likelihood principle.

IV-A Element-wise high resolution estimation

Considering the ultra-wideband nature of the measurement, (3) can be rewritten as

[TABLE]

with

[TABLE]

and

[TABLE]

Two parameters, i.e. $\tau_{p}^{\ell}$ and $\alpha_{p}^{\ell}$ are capable to characterize the $\ell$th path when a single (the $p$th) array element is considered. We denote the parameter vector $\Omega_{p}^{\ell}=[\tau_{p}^{\ell},\alpha_{p}^{\ell}]$, then the $L$ propagation paths impinging at the $p$th antenna element can be characterized by the parameter vector $\Omega_{p}=[\Omega_{p}^{1},\cdots,\Omega_{p}^{L}]$.

It can be observed from Fig. 3 that the conventional spectra-based methods such as the IDFT is not capable to obtain the high-resolution estimation of $\Omega_{p}$. The sidelobes of individual propagation paths exist and interfere with each other, and the weak paths can be buried by the sidelobes of strong paths. Therefore, the SAGE principle [10] or the EM principle[38] is applied to estimating $\Omega_{p}$ element-wise. With an ultra-wide system bandwidth, the SAGE algorithm can resolve two paths with very small relative propagation delay. By updating the parameter estimations sequentially and providing the maximum-likelihood estimation results iteratively, it has been proven practically in [10] that the SAGE algorithm has the ability to resolve two paths with $\Delta\tau\gtrapprox\frac{1}{5B}$. $B$ and $\Delta\tau$ represent the bandwidth and the relative delay between two paths, respectively. In our measurement campaign, paths with propagation distance differences no less than 3 cm can be well resolved by using the SAGE algorithm.

Fig. 4 illustrates the element-wise SAGE estimation results, $\hat{\Omega}=[\hat{\Omega}_{1},\cdots,\hat{\Omega}_{P}]$ for the channel measured in the LoS scenario. In practice, the path number $L$ should be set sufficiently large to fully extract the received power in the channel, which is evaluated as 20 in our case. It can be observed from Fig. 4 that the MPCs are well resolved. Different propagation paths with different trajectories are obvious by visual inspection, e.g. the LoS trajectory is smooth and continuous with the highest power. Furthermore, the channel spatial non-stationarity across the array elements can be clearly observed. The path strength along one trajectory may vary, which can result in the break (“birth-death” behavior) of the trajectory in the extreme case. In the literature, channel stationarity over the array elements is typically assumed for propagation parameter estimation, resulting in unrealistic estimation results and increasing complexity.

IV-B Trajectory identification based on the phase mode excitation

In this section, we firstly introduce the phase mode excitation technique applied for the 3D wideband spherical propagation scenario, and the trajectory identification for the LoS path is presented.

IV-B1 Phase mode excitation

As indicated by (8), $[\tau_{\ell},\phi_{\ell},\theta_{\ell},d_{\ell}]$ determines the delay trajectory of the $\ell$th propagation path in $\Omega$. The prior knowledge of the parameters is helpful to identify the trajectory.

The phase mode excitation technique [39, 40] have been proposed to transform the spatial-domain array response of the UCA, i.e. $Y(p,f)$ into the phase-mode domain so that a frequency-invariant beamformer can be achieved for the joint delay and angle estimation. However, the frequency-invariant beamformer is effective only for the two-dimensional (2D) propagation scenario, i.e. with all the elevations of paths strictly limited to $90^{\circ}$.

According to the Taylor expansion, (4) can be rewritten as

[TABLE]

where $\epsilon_{p,\ell}$ is the remainder of the Taylor series introduced by the spherical wave propagation. To begin with, let us consider the 2D propagation scenario with plane wave assumption, i.e. $\theta_{\ell}=90^{\circ}$ and $d_{\ell}=+\infty$ ($\epsilon_{p,\ell}=0$). The frequency response $H_{p,\ell}(f)$ contributed by the $\ell$th path at the $p$th Rx defined in (3) becomes

[TABLE]

The second equality in (10) holds according to the expansion

[TABLE]

where $J_{n}(\cdot)$ is the first kind of Bessel function with order $n$. By involving the so-called phase mode [40], i.e. the basis function $e^{-jm\phi_{p}}$, and a filter $W_{m}(f)$, the $m$th phase-mode response is formatted as

[TABLE]

The third approximation in (12) holds because

[TABLE]

holds only if $n=m+P\cdot z$ where $z$ is an arbitrary integer, and $J_{n}(2\pi f\frac{r}{c})$ decrease to near 0 when $z\neq 0$. The fourth equability in (12) is achieved by making

[TABLE]

The noiseless array response in the phase-mode domain reads

[TABLE]

Therefore, the joint delay-azimuth estimation can be accomplished by applying the 2D fast Fourier transform (FFT) to $Y(m,f;\Theta)$.

However, in the realistic channel, it is not practical to assume that all the elevations of paths are exactly $90^{\circ}$. Moreover, for the large-scale antenna array, spherical wave propagation cannot be neglected. Consequently, the bessel function part in the second step of (12) becomes $J_{n}(2\pi f\frac{r\sin\theta_{\ell}}{c})$, and the summation part with respect to $p$ in the second step of (12) becomes $\sum_{p}e^{j2\pi\frac{\epsilon_{p,\ell}}{c}}e^{j(m-n)\phi_{p}}$. The null misalignments between $J_{m}(2\pi f\frac{r}{c})$ and $J_{m}(2\pi f\frac{r\sin\theta_{\ell}}{c})$ destroys the phase-mode transform. The additive phases caused by $\epsilon_{p,\ell}$ may also make (13) invalid. To cope with the 3D propagation scenario with plane wave propagation, the filter

[TABLE]

was proposed in [21], where $J^{{}^{\prime}}_{n}(\gamma)$ represents the derivative of $J_{n}(\gamma)$ with respect to $\gamma$. The basic idea is to make use of the cosine-alike oscillation of bessel function $J_{n}(\gamma)$ when $\gamma$ is much larger than $n$. Therefore, the denominator of the filter is approximately the envelope of $J_{n}(\gamma)$. In such a way that the nulls of $J_{n}(\gamma)$ is avoided, while the amplitude is about 2 times of $J_{n}(\gamma)$, resulting in the numerator of (16) as 2.

Fig. 5 illustrates the simulated noiseless delay-azimuth power spectrum $p(\tau,\theta)$ of a single path with $\Theta=[15\,\text{ns},180^{\circ},70^{\circ},5\,\text{m},1]$ by applying the filter defined in (16). The other simulation configurations are the same with that as applied in the measurement campaign. The estimated delay and azimuth angle indicated by the spectrum maxima in Fig. 5 are $15\,\text{ns}$ and $180^{\circ}$, which are identical to the set values. Note that the elevation angle and source distance cannot be estimated in this step. In addition, it can be observed that sidelobes exist in both delay and azimuth domain, and the power is lower than the real path power. The reasons include that i) the null locations and the intervals between these nulls of $J_{m}(2\pi f\frac{r\sin\theta_{\ell}}{c})$ shift with respect to that of $J_{m}(2\pi f\frac{r}{c})$; ii) and that the additive spherical-wave phases are introduced. Nevertheless, the peak of the power spectrum is able to indicate the real $\tau$ and $\theta$ of this path. To investigate the effect of elevation and source distance to the estimation results of delay and azimuth, Fig. 5 illustrates the estimation errors of the path delay with different source distances and elevations, and Fig. 5 illustrates that of azimuths. It can be observed that for the spherical wave propagations with different elevations and source distances, phase mode excitation works well to obtain rough estimation results of azimuth and delay for the dominant path, with low computation complexity of 2D FFT.

IV-B2 Trajectory identification

The steps to identify the trajectory for the $\ell$th propagation path in $\hat{\Omega}^{-(\ell-1)}$ based on the phase mode excitation are as follows, where the superscript $-(\ell-1)$ indicates the remainder of $\Omega$ after removing the first $\ell-1$ path trajectories already identified, e.g., $\hat{\Omega}^{0}=\hat{\Omega}$.

Reconstruct the current array output $\hat{Y}^{-(\ell-1)}(p,f)$ according to (6) using $\hat{\Omega}^{-(\ell-1)}$. 2. 2.

Apply the phase-mode transform to $\hat{Y}^{-(\ell-1)}(p,f)$ and then 2D FFT to $\hat{Y}^{-(\ell-1)}(m,f)$ to obtain the angle-delay power spectrum $\hat{z}^{-(\ell-1)}(\tau,\phi)$. The azimuth $\hat{\phi}_{\ell}$ and delay $\hat{\tau}_{\ell}$ of the $\ell$th path are obtained by checking the peak of $\hat{z}^{-(\ell-1)}(\tau,\phi)$. 3. 3.

Define the area $A_{\theta_{\ell}}$ with the trajectory pairs ($T_{\theta_{\ell}}-\Delta\tau$, $T_{\theta_{\ell}}+\Delta\tau$), where

[TABLE]

Note that $d_{\ell}$ is ignored to avoid joint 2D searching. $\Delta\tau$ is aimed to tolerate the delay variation caused by the disturbance and the ignorance of $d_{\ell}$. As an example, Fig. 6 illustrates the path areas $A_{90^{\circ}}$ and $A_{30^{\circ}}$ for the LoS path identification, where $\Delta_{\tau}$ is practically set as $\frac{1}{2B}$. It can be observed that the LoS path trajectory is confined within $A_{90^{\circ}}$ perfectly however not in $A_{30^{\circ}}$, indicating the elevation angle is close to 90 degrees. Hence, the area in which the $\ell$th trajectory is confined can be identified as

[TABLE]

where $C(A_{{\theta_{\ell}}})$ denotes the number of array elements with propagation paths confined in $A_{\theta_{\ell}}$. As an example, Fig. 6 illustrates the $C(A_{{\theta_{1}}})$ with respect to different $\theta_{1}$s for the LoS path. It can be observed from Fig. 6 that $C(A_{{\theta_{1}}})$ is 720 (i.e. $P$) with elevation angles from $60^{\circ}$ to $90^{\circ}$. That is, by applying (18), a rough estimation or a rough initialization of $\theta_{\ell}$ can be obtained in this step. Note that due to the spatial non-stationarity, the maximum value of $C(A_{{\theta_{\ell}}})$ can be less then the elements number $P$. Furthermore, it is possible that for one element, more than one path can be contained in $A_{\hat{\theta}_{\ell}}$. To cope with the issue, the mean power of paths in $A_{\hat{\theta}_{\ell}}$ is calculated, and the path with minimal power difference from the mean power is selected element-wise. The selected paths forming the trajectory in $\hat{\Omega}^{-(\ell-1)}$ are denoted as $\Omega^{\ell}$.

As an example, readers can refer to the identification steps for the LoS trajectory ($\hat{\Omega}^{1}$) as illustrated in Figs. 7(a)-7(c).

IV-C *Maximization estimation for $\Theta_{\ell}$ *

The rough estimation results of delay, azimuth and elevation, i.e., $[\hat{\tau_{\ell}},\hat{\phi}_{\ell},\hat{\theta}_{\ell}]$, have been obtained in the former steps. The estimation of $\Theta_{\ell}$ for the $\ell$th path include three steps as follows.

Reconstruct the array output $\hat{H}(p,f;\Theta_{\ell})$ contributed by this single path according to (6) using $\Omega^{\ell}$ obtained in the trajectory identification. 2. 2.

Update the estimation of $[\hat{\phi}_{\ell},\hat{\theta}_{\ell},\hat{d_{\ell}}]$ by solving the following maximization problem

[TABLE]

where $W(p,f_{k})$ is generated by using (3) with parameter set $[0,{\phi},{\theta},d,1]$, $d$ is no larger than $c\hat{\tau}_{\ell}$, $\text{vec}\{\cdot\}$ denotes the vectorization of the argument matrix, and $(\cdot)^{*}$ represents the conjugation of the argument. Although the maximization problem is solved by 3D searching, the searching space is heavily reduced with known prior information of $\tau_{\ell}$, $\phi_{\ell}$ and $\theta_{\ell}$. Furthermore, without delay considered, only one single frequency point $f_{k}$ is selected for the maximization, which also decreases the computation load significantly. 3. 3.

Update $\hat{\tau}_{\ell}$ by solving the following maximization problem

[TABLE]

where $W(p,f)$ is generated by using (3) with parameter set $[\tau,\hat{\phi}_{\ell},\hat{\theta}_{\ell},\hat{d}_{\ell},1]$. The maximization is achieved by 1D searching with prior information of $\tau_{\ell}$. 4. 4.

Calculate the amplitude $\hat{\alpha}_{\ell}$ as

[TABLE]

where $C$ is the number counted in Sect. IV-B2, i.e. $C=C(A_{\hat{\theta}_{\ell}})$, and the $W(p,f)$ here is generated using (3) with parameter set $[\hat{\tau}_{\ell},\hat{\phi}_{\ell},\hat{\theta}_{\ell},\hat{d}_{\ell},1]$.

IV-D Algorithm implementation

The channel parameters of multiple propagation paths are estimated through an “identification-removing” operation. The procedure is elaborated in the pseudocodes listed in Algorithm 1.

{codebox}\Procname

**Algorithm 1: The channel parameter $\Theta$ is estimated according to the procedure described by the following pseudo-codes:

Input:** Measured array output $Y(p,f)$

Output: The parameter $\Theta$ characterizing the 3D spherical propagation channel \liObtain element-wise high-resolution $\hat{\Omega}$ from $Y(p,f)$; \liLet the path number $\ell$ be 1; \liLet $\hat{\Theta}$ be an empty set. \liwhile True \li if $\ell=1$ then \li  $\hat{\Omega}^{-(\ell-1)}=\hat{\Omega};$ \li else \li  $\hat{\Omega}^{-(\ell-1)}=\hat{\Omega}-\sum_{n=1}^{\ell-1}\hat{\Omega}^{n};$ \li end if \li Obtain $\hat{\Omega}^{\ell}$ from $\hat{\Omega}^{-(\ell-1)}$ according to Sect. IV-B2; \li Obtain $\hat{\Theta}_{\ell}$ from $\hat{\Omega}^{\ell}$ according to Sect. IV-C; \li if $C<C_{s}$ then \li  ** break**; \li end \li$\hat{\Theta}=\hat{\Theta}+\hat{\Theta}_{\ell};$ \li$\ell=\ell+1;$ \liend while

Fig. 7 exemplifies the procedure for the first two propagation paths. Note that $C_{s}$ in line 12 of Algorithm 1 is the threshold for checking the reasonability of $\hat{\Omega}^{\ell}$ identified. The estimation procedure stops if there are not enough elements with propagation paths in the current $\hat{\Omega}^{\ell}$, where it is considered that all the reasonable propagation paths have been identified and estimated. Fig. 8 illustrates all the $\hat{\Omega}^{\ell}$s identified by applying Algorithm 1 to the measured array output in the LoS scenario as illustrated in Fig. 3. In this case, $C_{s}$ is practically chosen as 360, the half of the UCA element number $P$. $\Delta\tau$ is practically chosen as $\frac{1}{2B}$.

V Results and remarks

In this section, the estimation results of the proposed algorithm are firstly presented, followed by the remarks on the proposed algorithm.

V-A Estimation results

Fig. 9(a) and Fig. 9(b) illustrate the delay-azimuth power spectra obtained by using the proposed algorithm and measured with the reference horn antenna, respectively, for the LoS scenario. It can be observed that the estimated delay-azimuth power spectrum is in good consistency with the spectrum measured by using the reference horn antenna. Fig. 9(c) also illustrates the estimated delay-azimuth power spectrum by using the MLE algorithm proposed in [34] with 4D parameter searchings. By comparing Figs. 9(a) and 9(c), it can be observed that the two results are very close. However, due to the fact that the channel non-stationarity across the array elements cannot be considered in the MLE algorithm (or the SAGE algorithm), more than one stationary propagation paths (typically with very similar geometrical parameters yet different complex amplitudes) are estimated to mimic the spatial non-stationarity of one propagation path, which results in more artificial paths estimated and increases the estimation complexity.

Fig. 10 illustrates the source distances for different propagation paths estimated. It can be observed from Fig. 10 that the estimated source distances generally become larger when the propagation distances are larger, which demonstrates the fact that the basement walls were probably smooth and with high reflection effects.

V-B Remarks on the proposed algorithm

The proposed algorithm is essentially a HRPE algorithm based on the maximization-likelihood principle. Compared to the other HRPE algorithms, e.g., in [12, 31, 30, 41, 34], the proposed algorithms have the following advantages.

• •

The proposed algorithm is based on the array-element-wise high resolution estimation in the delay domain as elaborated in Sect. IV-A. In this step, the high resolution ability of the algorithm is guaranteed. Moreover, due to the fact that only delay domain is searched, the complexity in this step is low.

• •

In the trajectory identification step as elaborated in Sect. IV-B, the proposed algorithm exploits the phase mode excitation technique, the delay trajectories of propagation paths can be identified with 2D FFT with low computation load. Note that the trajectory identification procedure also has the effect of interference cancellation among different paths, which means that the considerable iterations applied in [12, 31, 30, 41, 34] are no longer required. Furthermore, the rough estimation results of delay, azimuth and elevation can be obtained.

• •

In the maximization estimation step as elaborated in Sect. IV-C, although a 3D and a 1D parameter searchings are applied, the searching spaces are heavily reduced with the prior initializations obtained in the trajectory identification. Therefore, the computation complexity is decreased significantly compared to the MLE algorithm [34] with exhaustive 4D parameter searchings. Moreover, the joint parameter searching avoids the space-alternating operations in the SAGE algorithm [12, 31, 30], where the violation of OSM can lead to the convergence failure of the SAGE algorithm.

• •

The channel non-stationarity across the array elements is considered in (21), which is firstly addressed in the literature for the propagation parameter estimation. The estimation results are more realistic to reproduce the spherical propagation compared to the SAGE [12, 31, 30] and MLE [34] algorithms where stationarity paths are assumed.

VI Conclusions

In this contribution, a complexity-efficient high resolution parameter estimation (HRPE) algorithm was proposed for the ultra-wideband large-scale uniform circular array (UCA). The algorithm was demonstrated based on a recently conducted measurement campaign in a basement. The delay trajectories of different propagation paths across the virtual array elements were well observed, which can roughly indicate the channel information. The proposed HRPE algorithm takes advantage of the prior channel information obtained in the trajectories and is capable to involve the spatial non-stationarity of realistic channels. Investigations show that the algorithm can obtain the high-resolution results, i.e., the delays, azimuths, elevations, source distances and complex amplitudes for individual spherical propagation paths in the underlying mm-wave channel. Fast initializations and effective interference cancellations reduce the computation load heavily. Furthermore, the consideration of the channel spatial non-stationarity is helpful to obtain more realistic results and to decrease the algorithm complexity, by avoiding the artificial stationary paths being estimated. In addition, it is noteworthy that the basement walls have significant reflection effects for the mm-wave channel at 28-30 GHz.