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

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$ m² 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$.¹¹1The 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$ m² 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.

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.

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

with

and

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.

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.

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.

1. 1.

   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

   $\displaystyle(\hat{\phi}_{\ell},\hat{\theta}_{\ell},\hat{d}_{\ell})=\arg\max_{{\phi},{\theta},d}\text{vec}\{\hat{H}(p,f_{k};\Theta_{\ell})\}^{T}\text{vec}\{W^{*}(p,f_{k})\}$

   (19)

   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

   $$\hat{\tau}_{\ell}=\arg\max_{\tau}\text{vec}\{\hat{H}(p,f;\Theta_{\ell})\}^{T}\text{vec}\{W^{*}(p,f)\}$$

   (20)

   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

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

   (21)

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

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.

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}$.

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.

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.